plot system of 5 odes over time using ode45
5 views (last 30 days)
Show older comments
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 Comments
Answers (1)
Walter Roberson
on 8 Dec 2021
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 Comments
Walter Roberson
on 8 Dec 2021
Edited: Walter Roberson
on 8 Dec 2021
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
on 11 Dec 2021
Edited: Walter Roberson
on 11 Dec 2021
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
See Also
Categories
Find more on Symbolic Math Toolbox in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!