How to use ode45 to plot this code?
Mostrar comentarios más antiguos
Hello, I have to use ode45 to plot I on the same figure as the other way I plotted I (using a for loop). I tried doing it this way but the error that keeps coming up is that "Vectors must be the same length". Do you have any suggestions in how to fix that? Thank you.
clc
clear
N=1000;
del_T=0.01; %delta T
B=2; %B=Beta
g=0.5; %g=gamma
R_o=4; %R not
S(1)=999;
I(1)=1;
R(1)=0;
t=0:del_T:30;
for k= 1:length(t)-1
S(1+k)=S(k)-(del_T*B*I(k).*(S(k)/N));
I(1+k)=I(k) + ((B*I(k)*S(k)*del_T)./N)-(g*I(k)*del_T);
R(1+k)=R(k) + g*I(k)*del_T;
tspan = [0 30];
y0 = [S(1) I(1) R(1)];
[t,I_ode]= ode45(@(t,I_ode)odefun(t,I_ode,k,B,g,N,del_T,S,I,R), tspan, y0);
end
plot(t, I)
hold on;
plot (t,I_ode)
function I_ode= odefun(t,I_ode,k,B,g,N,del_T,S, I,R)
I_ode=zeros(3,1);
dSdt= S(k)-(del_T*B*I(k).*(S(k)/N));
dIdt=I(k) + ((B*I(k)*S(k)*del_T)./N)-(g*I(k)*del_T);
dRdt=R(k) + g*I(k)*del_T;
end
Respuestas (1)
Alan Stevens
el 13 de Nov. de 2020
You have mixed up doing your own iterations with allowing ode45 to do its iterations! For ode45, just use the following
N=1000;
del_T=1; %delta T
B=2; %B=Beta
g=0.5; %g=gamma
R_o=4; %R not
S=999;
I=1;
R=0;
tspan = [0 30];
y0 = [S I R];
[t,y]= ode45(@(t,y)odefun(t,y,B,g,N), tspan, y0);
S = y(:,1); I = y(:,2); R = y(:,3);
plot(t, y)
xlabel('time'), ylabel('S, I and R')
legend('S','I','R')
function I_ode= odefun(~,y,B,g,N)
S = y(1); I = y(2); R = y(3);
dSdt= -B*I.*S/N;
dIdt= B*I*S./N-g*I;
dRdt= g*I;
I_ode = [dSdt; dIdt; dRdt];
end
1 comentario
Ethan Hoang
el 13 de Nov. de 2020
Editada: Ethan Hoang
el 13 de Nov. de 2020
Categorías
Más información sobre Ordinary Differential Equations en Centro de ayuda y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!