ODEs – Euler Method for SIR Model in MATLAB

dynamical systemsMATLABnumerical methodsordinary differential equations

I know how to solve SIR model with ode45, but I just wanted to try it out using Euler Method.
The SIR model is:

$\dot S=-\beta IS \\ \dot I = \beta IS – \gamma I \\ \dot R = \gamma I$
The code that I wrote is

function Euler
deltaT=1;
s=1;
i=1e-6;
r=1-(s+i);
beta=1.4247;
gamma=0.14286-1e-6;
maxTime=70;
t=[0];
for n=1:maxTime
    t(n+1)=t(n)+deltaT;
    s(n+1)=s(n)-beta*s(n)*i(n)*deltaT;
    i(n+1)=i(n)+(beta*s(n)*i(n)-gamma*i(n))*deltaT;
    r(n+1)=r(n)+gamma*i(n)*deltaT;
end
figure
plot(0:maxTime,s)
hold on 
plot(0:maxTime,i)
hold on 
plot(0:maxTime,r)

when I make the step size to be 1, I get the following figure which looks like the expected behaviour.

enter image description here
But when I make the step size small= 0.001, I get a very different figure.
enter image description here

Why does the behaviour change when I change the step sizes. What have I done wrong in the code?

Best Answer

There's an error in the code - you are only (as @mickep suggested) doing 70 time steps in total through your simulation, whatever the timestep. To get the same simulation time as before you need to change the line

for n=1:maxTime

to

for n=1:(maxTime/deltaT)

and change the plot commands to plot (0:deltaT:maxTime,...). This would not change your first plot (as deltaT=1) but the second would change.

Secondly, it makes sense generally, to preallocate your arrays t=zeros(1,maxTime/deltaT); (and similarly for s, i and r).

Thirdly (and this is very minor), as r is just the remaining population, you don't need to solve this differential equation - just plot 1-(s+i) for your recovered fraction.

Related Question