variable coefficient

Question

Rose 2011-4-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/6097-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!!!!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt Tearle 2011-4-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/6097-variable-coefficient#answer_8481

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 个评论
显示 32更早的评论隐藏 32更早的评论

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

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

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

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-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

%%%%%%%%%%%%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-5-6

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

Matt Tearle 2011-5-8

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

请先登录，再进行评论。

variable coefficient

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

34 个评论
显示 32更早的评论隐藏 32更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

variable coefficient

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

34 个评论 显示 32更早的评论隐藏 32更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

34 个评论
显示 32更早的评论隐藏 32更早的评论