plot system of 5 odes over time using ode45

Question

Sam 2021-12-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1605800-plot-system-of-5-odes-over-time-using-ode45

编辑： Walter Roberson 2021-12-11

Here are my equations and initial conditions. I'm trying to plot these odes as a function of time with the following initial conditions:

xinit = [10000, 0, 0.01, 0, 0];
f = @(t,y) [(10000 - (0.01*y(1)) - (0.000000024*y(3)*y(1)*20001)); ((0.05*0.000000024*y(3)*y(1)*20001) - y(2)); ((25000*y(2)) - 23); ((0.95*0.000000024*y(3)*y(1)*20001) - (0.001*y(4))); ((15*0.001*y(4)) - (6.6*y(5)))];

I'm looking for oscillatory behavior. Here they are written slightly differently:

      T = 10000;
      Q = 0;
      V = 0.01
      P = 0;
      I = 0;
      dydt(1) = 10000 - 0.01*T - 0.000000024*V*T*20001;
      dydt(2) = 0.05*0.000000024*V*T*20001 - Q;
      dydt(3) = 25000*Q - 23;
      dydt(4) = 0.95*0.000000024*V*T*20001 - 0.001*P;
      dydt(5) = 15*0.001*P - 6.6*I;

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

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2021-12-8

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1605800-plot-system-of-5-odes-over-time-using-ode45#answer_850135

在 MATLAB Online 中打开

tspan = [0 0.184];

xinit = [10000, 0, 0.01, 0, 0];

[t, y] = ode45(@odefun, tspan, xinit);

plot(t, y);

legend({'T', 'Q', 'V', 'P', 'I'}, 'location', 'best');

function dydt = odefun(t, y)

T = y(1); Q = y(2); V = y(3); P = y(4); I = y(5);

dydt = zeros(5,1);

dydt(1) = 10000 - 0.01*T - 0.000000024*V*T*20001;

dydt(2) = 0.05*0.000000024*V*T*20001 - Q;

dydt(3) = 25000*Q - 23;

dydt(4) = 0.95*0.000000024*V*T*20001 - 0.001*P;

dydt(5) = 15*0.001*P - 6.6*I;

end

You cannot go much further. By roughly 1.88 the ode functions refuse to continue as the slopes have become to steep.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Walter Roberson 2021-12-8

编辑：Walter Roberson 2021-12-8

在 MATLAB Online 中打开

Look at the first equation,

and compare to

dydt(1) = 10000 - 0.01*T - 0.000000024*V*T*20001;

The equation starts with the constant s, which matches the 10000 in your code; so far so good.

Then the equation has -

and we can speculate that the

maps to your constant named dt (not a great name considering that you are taking derivatives with respect to t and so have dt in your equations with quite a different meaning... but anyhow.) That plausibly matches the - 0.01*T term, so that's okay.

Then the equation has

. k is a constant, V is the function we have mapped to y(3), T is the function we have mapped to y(1), so the part up to 0.000000024*V*T matches. But after that in the equation is the

part. γ is a constant. But I is the function mapped to y(5) . So the code should not have 20001 at that point, and should instead have (1 + 2*10^4 * I) where I = y(5) is in effect.

All of your 20001 in your code should be (1 + 2*10^4 * I)

You should re-check the code.

And there is no shame at all in putting all of those constants such as s as variables in your function code and using them to make the code clear, instead of hard-coding values like 0.000000024 or 20000 in the equations. Any tiny loss you might take in execution speed for using more variables will be more than made up for in time saved debugging the code.

Walter Roberson 2021-12-8

编辑：Walter Roberson 2021-12-8

You can use tiledlayout() to get larger graphs.

Or just put them into different figures.

Plots 2, 4, 5 look like they might fit well on the same plot, but plots 1 and 3 have a significantly different range and do not fit on the same plot.

Good thing you caught my mistake in the second equation.

Walter Roberson 2021-12-11

编辑：Walter Roberson 2021-12-11

在 MATLAB Online 中打开

tspan = [0:.2:50, 100:100:2000];

xinit = [1000, 0, 0.001, 0, 0]; %OR ;xinit = [1000, 0, 0.001, 0, 0]; OR ;xinit = [1, 0, 0.001, 0, 0]; OR ;xinit = [10000, 0, 0.001, 0, 0];

[t, y] = ode45(@odefun, tspan, xinit);

symbols = {'T', 'T^*', 'V', 'P', 'I'};

figure()

whichplot = [1 3];

N = length(whichplot);

for K = 1: N

subplot(N,1,K);

plot(t, y(:,whichplot(K)));

legend(symbols(K), 'location', 'best');

end

function dydt = odefun(t, y)

T = y(1); T__star = y(2); V = y(3); P = y(4); I = y(5);

s = 10000;

k = 2.4*10^-8;

gamma = 2*10^-4;

f = 0.95;

r = 2.5*10^4;

B = 15;

d_T = 0.01;

d_T__star = 1;

d_V = 0.001;

d_P = 23;

d_I = 6.6;

dydt = zeros(5,1);

dydt(1) = s - d_T*T - k*V*T*(1+gamma*I);

dydt(2) = (1-f)*k*V*T*(1+gamma*I) - d_T__star*T__star;

dydt(3) = r*T__star - d_V;

dydt(4) = f*k*V*T*(1+gamma*I) - d_P*P;

dydt(5) = B*d_P*P - d_I*I;

end

请先登录，再进行评论。

plot system of 5 odes over time using ode45

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

回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

plot system of 5 odes over time using ode45

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

回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

WeChat

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

4 个评论
显示 2更早的评论隐藏 2更早的评论