plot system of 5 odes over time using ode45

5 views (last 30 days)
Sam
Sam on 8 Dec 2021
Edited: Walter Roberson on 11 Dec 2021
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;

Answers (1)

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

Sign in to comment.

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!