variable coefficient

Hi, I have a system of 4 non linear ODEs. I solved them using ode15s and it worked perfectly. But i used constant coefficients at that time. Now I want to do the same thing when the coefficients are variable. Is there any built-in solver for variable coefficients? If yes, then how do i input the coefficients as functions of time? Any guidance/help will be appreciated. thanks!!!!

Answer 1

[Matt Tearle]

All the solvers can handle variable coefficients - they are completely ignorant of the details of the equations. All they require is a function that returns some f(t,y) (with the dimension of f being the same as the dimension of y).

But perhaps you mean that you want to pass the time-dependent coefficients as parameters?

Can you maybe clarify? Give an example?

[Rose]

thanks once again for responding to my problem, Matt !! I think you developed my interest in MATLAB!!

Yes, the coefficients are time-dependent now.

this is my system of equations:

function dy = funl12(t,y)

global k1 k2 k3 mu d1 d2 r K k alpha gamma n;

dy(1,1) = (k1(t).*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2(t).*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;

dy(2,1) = (mu(t).*(y(2).^n)/(K^n+y(2).^n)).*exp((-y(1)).*k)-(k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d1(t)+gamma.*y(4)).*y(2);

dy(3,1) = (k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d2(t)+gamma.*y(4)).*y(3);

dy(4,1) = r(t).*y(4).*(1-(y(4)./(alpha(t).*(y(2)+y(3)-1e-12))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

all these parameters k1,k2,k3,d1,d2 etc are piecewise constants. So, what I am trying is to make an m-file for each of the parameters. Like, for k1(t), I wrote a sub-routine:

function y = k1(t)

if 0<t<=91.25;

y=0.04226;

elseif 91.25<t<=182.5;

y=0.1593;

elseif 182.5<t<=273.75;

y=0.1460;

else 273.75<t<=365;

y = 0.1489;

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

But this doesn't seem to be working. May be, you can tell me the fault?

[Rose]

let me write the exact code i am using:

%%%%%funl12.m%%%%%%%%%

function dy = funl12(t,y)

global k1 k2 k3 mu d1 d2 r K k alpha gamma n;

dy(1,1) = (k1(t).*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2(t).*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;

dy(2,1) = (mu(t).*(y(2).^n)/(K(t)^n+y(2).^n)).*exp((-y(1)).*k)-(k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d1(t)+gamma(t).*y(4)).*y(2);

dy(3,1) = (k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d2(t)+gamma(t).*y(4)).*y(3);

dy(4,1) = r(t).*y(4).*(1-(y(4)./(alpha(t).*(y(2)+y(3)-1e-12))));

%%%%%%%%%%%%%%%%%%%k1.m%%%%%%%%%%%%

function y = k1(t)

if 0<t<=91.25;

y=0.04226;

elseif 91.25<t<=182.5;

y=0.1593;

elseif 182.5<t<=273.75;

y=0.1460;

else 273.75<t<=365;

y = 0.1489;

end;

%%%%%%%%k2.m%%%%%%%%%%%%%%

function y = k2(t)

if 0<t<=91.25;

y=0.008460;

elseif 91.25<t<=182.5;

y=0.04959;

elseif 182.5<t<=273.75;

y=0.03721;

else 273.75<t<=365;

y=0.04750;

end;

%%%%%%%%%%%%%

etc for all other parameters too.

%%%%%%%main file%%%%%%%%%%%

n=2;

k=0.00003125;

finaltime = 365;

tspan = [0 finaltime];

y0=[0;10000;0;0];%winter

options = odeset('RelTol',1e-12,'AbsTol',1e-12);

sol = ode15s(@funl12,tspan,y0,options);

%[t,y] = ode15s('funpyo',tspan,y0);

X = sol.x;

Y = (sol.y)';

plot(X,Y(:,1),'r-', X,Y(:,2),'g-', X,Y(:,3), X, Y(:,4))

xlabel('time')

ylabel('Y')

legend('m', 'X', 'y', 'M')

figure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

the error i am getting is:

??? Attempted to access k1(0); index must be a positive integer or logical.

Error in ==> funl12 at 4

dy(1,1) = (k1(t).*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2(t).*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;

Error in ==> odearguments at 110

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ==> ode15s at 228

[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> simulate_lis2 at 66

sol = ode15s(@funl12,tspan,y0,options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

i don't know what should i do, Matt!!

[Matt Tearle]

Remove the "global" line from funl12.m

Another approach would be to integrate in pieces, from t=0 to t=91.25, then from 91.25 to 182.5, etc. Use the solution at the end of the previous run as the initial condition to the next run. Then set all the parameters to be constants.

[Rose]

it worked! But can I do it for the several years. I want the code to run for hundred of years. Is it possible?

[Rose]

i also want to convert these piecewise constant functions k1,k2 etc into periodic functions. How can I do that? Can you suggest me something. Thanks a lot!!

[Matt Tearle]

I just looked more closely at your code, which is good, because... it may be running but it's not working as you'd expect. Try it -- enter >> k1(200) Not what you wanted, is it?

This is because of how MATLAB is evaluating the statement "0<t<=91.25". Skipping the details, the short answer is that this is NOT the same as you'd read mathematically. To make a compound inequality, use two conditions, joined with an AND: eg

if (t>0) && (t<=91.25)

However, else/elseif constructs are mutually exclusive, so you don't even need the multiple statements. That is, if t = 200, the first two conditions (t<=91.25 and t<=182.5) will have already been tested and failed, so you *know* t>182.5 -- no need to test it (and doing so just wastes time).

So, bottom line: here's your code rewritten:

if t<=91.25

y=0.04226;

elseif t<=182.5

y=0.1593;

elseif t<=273.75

y=0.1460;

else

y = 0.1489;

end

Now, to the next point: how to do this for multiple years (ie t>365). A simple trick is to do t = mod(t,365) at the beginning of the function. This wraps t to the interval [0,365).

To me it seems like there might be some slight confusion in the way you've set up your inequalities. If t==0, k1=0.1489, but at the very next step (ie any t>0), k1=0.04226. Does that really make sense, given that you're starting at t=0? Up to you to decide, obviously. But if you actually want the intervals [0,91.25), [91.25,182.5), etc, then you can replace all your code with

function y = k1(t)

yv = [0.04226,0.1593,0.1460,0.1489];

y = yv(1+floor(mod(t,365)/91.25));

[Rose]

If you look at my other post, I ended up using

if (t>0) && (t<=91.25)

This is the final code I am using now. Please let me know if this is okay !!

function dy = funl12(t,y)

dy(1,1) = (k1(t).*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2(t).*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;

dy(2,1) = (mu(t).*(y(2).^2)/(K(t)^2+y(2).^2)).*exp((-y(1)).*k(t))-(k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d1(t)+gamma(t).*y(4)).*y(2);

dy(3,1) = (k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d2(t)+gamma(t).*y(4)).*y(3);

dy(4,1) = r(t).*y(4).*(1-(y(4)./(alpha(t).*(y(2)+y(3)-1e-12))));

% Nested k1

function y = k1(t)

y = [0.04226,0.1593,0.1460,0.1489];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

% Nested k2

function y = k2(t)

y = [0.008460,0.04959,0.03721,0.04750];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested k3

function y = k3(t)

y = [0.03384,0.1984,0.1460,0.1900];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested mu

function y = mu(t)

y = [0,500,1500,500];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested d1

function y = d1(t)

y = [0.005263,0.02272,0.04,0.02272];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested d2

function y = d2(t)

y = [0.005300,0.2,0.2,0.2];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested r

function y = r(t)

y = [0.0045,0.0165,0.0165,0.0045];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested K

function y = K(t)

y = [10000,20000,30000,20000];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested k

function y = k(t)

y = [0.00003125,0.00003125,0.00003125,0.00003125];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested alpha

function y = alpha(t)

y = [0.4784,0.4784,0.1321,0.1321];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested gamma

function y = gamma(t)

y = [10^(-6),0.0002,0.0000002,0.0002];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

the main file is:

finaltime = 365;

tspan = [0 finaltime];

y0=[0;10000;0;0];%winter

options = odeset('RelTol',1e-12,'AbsTol',1e-12);

sol = ode15s(@funl12,tspan,y0,options);

%[t,y] = ode15s('funpyo',tspan,y0);

X = sol.x;

Y = (sol.y)';

plot(X,Y(:,1),'r-', X,Y(:,2),'g-', X,Y(:,3), X, Y(:,4))

xlabel('time')

ylabel('Y')

legend('m', 'X', 'y', 'M')

figure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

What I want is:

when t lies between 0 and 91.25, k1=0.04226

when t lies between 91.25 and 182.5, k1=0.1593(but this time, 91.25 is included in the previous interval, so no need to include it again.)

and so on....

[Rose]

this doesn't make any sense that if If t==0, k1=0.1489, but at the very next step (ie any t>0), k1=0.04226.

[Rose]

I entered k1(200), I got an error message :( ======>>>

??? Undefined function or method 'k1' for input arguments of type 'double'.

[Oleg Komarov]

You cannot call k1 which is a nested function from the cmd line.

Also consider, as I already explained in the previous post, the behaviour of histc, an example:

histc(10, [7 8 9 10 11])

will return [0 0 0 1 0], meaning that 10 is in the bin [10 11).

If you want to it to be in bin (9 10], then redefine the call:

histc(10, [7.00001 8.00001 9.00001 10.00001 11.00001])

but I really suggest to adopt the [..., ...) (i.e. 10 <= x < 11)

[Rose]

I am not sure if I got your point or not. Sorry about that. I am really good in programming. But I think the way i set up the code finally with your and Matt's help(above in 7th comment) is okay. Isn't it?? All I want the code to do for is- cover all the days in the year i.e from t=0 to t=365. Don't consider a day two times. If 91.25th day is considered in one interval, it shouldn't be considered again.

One more thing, How can I get to know k1(200) or other values?

[Oleg Komarov]

t = [200, 300,29,0,365];

y = [0.04226,0.1593,0.1460,0.1489];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

[Rose]

I am NOT really good in programming **

[Rose]

Okay !! thanks Oleg!! I will try that. If I am using this nested thing, can I use the same command t = mod(t,365) to run the code for years and years? (above, Matt suggested me that )

[Rose]

i tried the command:

t=mod(t,365)

it didn't work. I am getting error:

Error in ==> odearguments at 110

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ==> ode15s at 228

[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> simulate_lis2 at 43

sol = ode15s(@funl12,tspan,y0,options);

??? Undefined function or variable 't'.

Error in ==> simulate_lis2 at 27

t=mod(t,365);

[Oleg Komarov]

you can use Matt's solution (very clever) or you can use histc.

But at this point I don't really understand what's the problem.

Please edit your first message and paste the whole code you're actually using and format it with code {} button.

Then we can look at it.

[Matt Tearle]

Be very careful about redefining variables if you're using nested functions, because nested functions share workspaces with their parents. So changing t inside a nested function can change t in the parent function (depending on how things are called). That said, you seem to be getting a different error... it's hard to work out why you're getting "undefined function or variable 't'" without seeing the whole code. The odd thing is that it seems to be pointing to "simulate_lis2", which appears to be your main function, not the ode function.

[Matt Tearle]

The issue about the intervals isn't whether a value is included twice, just which interval each value is in. Do you want your intervals to be [0,91.25), [91.25,182.5), etc, or (0,91.25], (91.25,182.5], etc.? That's up to you to determine. But to me, it seems more natural to go with the former, because you're starting at t=0. If your intervals are (0,91.25], ..., (273.75,365], then t=0 (which is equivalent, mod 365, to t = 365) is actually the last point of the last interval from the previous year. In calendar terms, I'd see t=0 as midnight on Jan 1, so t=0 would be part of the first quarter (and the parameters should be the same throughout this quarter). Dec 31 is day 364, right up to 364.999..., which is part of the fourth quarter, then 365 is equivalent to t=0 of the next year.

[Rose]

Okay Matt, Thanks for the excellent explanation of the intervals. I got your point. I also agree with the former one [0,91.25), [91.25,182.5), etc, and i think 'histc' works exactly like that. Isn't it. I will post my code now.wait..

[Rose]

t=mod(t,365)

function dy = funl12(t,y)

dy(1,1) = (k1(t).*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2(t).*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;

dy(2,1) = (mu(t).*(y(2).^2)/(K(t)^2+y(2).^2)).*exp((-y(1)).*k(t))-(k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d1(t)+gamma(t).*y(4)).*y(2);

dy(3,1) = (k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d2(t)+gamma(t).*y(4)).*y(3);

dy(4,1) = r(t).*y(4).*(1-(y(4)./(alpha(t).*(y(2)+y(3)-1e-12))));

% Nested k1

function y = k1(t)

y = [0.04226,0.1593,0.1460,0.1489];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

% Nested k2

function y = k2(t)

y = [0.008460,0.04959,0.03721,0.04750];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested k3

function y = k3(t)

y = [0.03384,0.1984,0.1460,0.1900];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested mu

function y = mu(t)

y = [0,500,1500,500];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested d1

function y = d1(t)

y = [0.005263,0.02272,0.04,0.02272];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested d2

function y = d2(t)

y = [0.005300,0.2,0.2,0.2];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested r

function y = r(t)

y = [0.0045,0.0165,0.0165,0.0045];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested K

function y = K(t)

y = [10000,20000,30000,20000];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested k

function y = k(t)

y = [0.00003125,0.00003125,0.00003125,0.00003125];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested alpha

function y = alpha(t)

y = [0.4784,0.4784,0.1321,0.1321];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

%Nested gamma

function y = gamma(t)

y = [10^(-6),0.0002,0.0000002,0.0002];

idx = logical(histc(t,[0,91.25,182.5,273.75,366]));

y = y(idx);

end

end

%%%%%%%%%%%%main.m %%%%%%%%%%%%%%%%%%%%%

t=mod(t,365);

finaltime = 365;

tspan = [0 finaltime];

y0=[0;20000;0;2000];%winter

options = odeset('RelTol',1e-12,'AbsTol',1e-12);

sol = ode15s(@funl12,tspan,y0,options);

X = sol.x;

Y = (sol.y)';

plot(X,Y(:,1),'r-', X,Y(:,2),'g-', X,Y(:,3), X, Y(:,4))

xlabel('time')

ylabel('Y')

legend('m', 'X', 'y', 'M')

figure

[Rose]

did you want to say that i should prefer separate sub-routines for k1, k2 etc instead of nested things?

[Matt Tearle]

Not necessarily - you just need to make sure you understand the implications. From a quick look at your rate equations, it looks like t shows up only in the coefficients. I assume also that you have one main script and one function for the rate equations, with the coefficients defined in nested functions (children of the main rate equations function). In that case, I'd just put the mod command as the first statement of funl12. That should be all you need.

[Rose]

I tried that, but i got an error:

??? Error using ==> feval

Error: File: funl12.m Line: 23 Column: 1

Function definitions are not permitted in this context.

Error in ==> odearguments at 110

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ==> ode15s at 228

[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> simulate_lis2 at 44

sol = ode15s(@funl12,tspan,y0,options);

[Matt Tearle]

The first line *after* the function declaration. The declaration must always be the first line of code.

function dy = funl12(t,y)

t = mod(t,365);

dy(1,1) = [etc]

[Rose]

it works. But when i ran it for 4000 days, it doesn't give me the same results for every year. It should because the data is the same for every year. I don't know what is the problem actually!!

[Walter Roberson]

4000 days is nearly 11 years -- enough that a minimum of one leap year would be encountered. 2 leap years minimum if the data doesn't happen to cross one of the century numbers that were not multiples of 400 (e.g., 1900 was not a leap year.) Using t = mod(t,365) ignores leap years.

Is it the first and the last years that are different from the rest? If so that might be an issue with running with incomplete years.

[Rose]

I think you mistook my question. I don't have the data for 11 years. I have the data for one year and this data is same for every year(I have the seasonal averages for all the parameters for all the four seasons). You can see the code pasted in the above comments. Suppose, I want to run the code for 4000 years, it should give me the same graph for every year. but it is not!

[Matt Tearle]

Why should it give the same graph for each year? You have a dynamical system where the *rate equations* are periodic, but that doesn't mean the solution has to be. When I run the above code, I get a solution that tends to an equilibrium at [0,0,0,0]. Consequently, each year starts with a different set of initial conditions.

[Rose]

yes, you are right Matt. Thanks a lot for your help. I am so grateful. Thanks to Walter too :-)

[Rose]

f i want to have the value of y(2) after 10 years, how can I do that? What I did was- I ran the code for 10 years, and then I typed y(2) in the command window, but it didn't give me value of y(2) after 10 years.

[Matt Tearle]

[t,y] = odeXX(...) returns a vector t and a matrix y where each column is a component of the solution, and each row corresponds to a time. So y_2 is the second column of y. The last value, then, is y(end,2).

[Rose]

I understand your point but how can I see that value?

[Rose]

I tried y(end,2) but it didn't work!!

[Matt Tearle]

What do you mean? Error message? Nothing? Unexpected value? Can you do a whos and report the result?