You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
variable coefficient
2 views (last 30 days)
Show older comments
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!!!!
Accepted Answer
Matt Tearle
on 23 Apr 2011
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?
34 Comments
Rose
on 23 Apr 2011
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
on 24 Apr 2011
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
on 25 Apr 2011
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.
Matt Tearle
on 25 Apr 2011
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
on 25 Apr 2011
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....
Oleg Komarov
on 25 Apr 2011
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
on 25 Apr 2011
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
on 25 Apr 2011
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
on 25 Apr 2011
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
on 25 Apr 2011
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
on 26 Apr 2011
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
on 26 Apr 2011
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
on 26 Apr 2011
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
on 26 Apr 2011
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
Matt Tearle
on 26 Apr 2011
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
on 27 Apr 2011
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
on 27 Apr 2011
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]
Walter Roberson
on 30 Apr 2011
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
on 1 May 2011
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
on 2 May 2011
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.
Matt Tearle
on 5 May 2011
[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).
Matt Tearle
on 8 May 2011
What do you mean? Error message? Nothing? Unexpected value? Can you do a whos and report the result?
More Answers (0)
See Also
Categories
Find more on Function Creation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)