variable coefficient
2 次查看(过去 30 天)
显示 更早的评论
Rose
2011-4-23
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!!!!
采纳的回答
Matt Tearle
2011-4-23
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 个评论
Rose
2011-4-23
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
2011-4-24
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
2011-4-25
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
2011-4-25
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
2011-4-25
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
2011-4-25
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
2011-4-25
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
2011-4-25
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
2011-4-25
I entered k1(200), I got an error message :( ======>>>
??? Undefined function or method 'k1' for input arguments of type 'double'.
Oleg Komarov
2011-4-25
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
2011-4-25
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
2011-4-25
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
2011-4-25
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
2011-4-25
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
2011-4-25
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
2011-4-26
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
2011-4-26
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
2011-4-26
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
2011-4-26
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
2011-4-26
did you want to say that i should prefer separate sub-routines for k1, k2 etc instead of nested things?
Matt Tearle
2011-4-26
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
2011-4-27
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
2011-4-27
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
2011-4-30
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
2011-4-30
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
2011-5-1
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
2011-5-2
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
2011-5-5
yes, you are right Matt. Thanks a lot for your help. I am so grateful. Thanks to Walter too :-)
Rose
2011-5-5
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
2011-5-5
[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
2011-5-8
What do you mean? Error message? Nothing? Unexpected value? Can you do a whos and report the result?
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)