Plotting phase plane in Matlab for SIR model

I want to plot the phase plane of the SIR model,

I used the method describe in here

This is the code that I wrote,

function SIREquilibrium()
S0=0.8;
I0=0.2;
R=1-S0-I0;
beta=1;
gamma=1/10;
mu=5e-4;

f=@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2);gamma*y(2)-mu*y(3)];

y1=linspace(0,1,20);
y2=linspace(0,1,20);
[x,y]=meshgrid(y1,y2);
u=zeros(size(x));
v=zeros(size(y));
t=0;
for i=1:numel(x)
Yprime=f(t,[x(i);y(i)]);
u(i)=Yprime(1);
v(i)=Yprime(2);
end
quiver(x,y,u,v,'r')


But I get an error as

index exceeds matrix dimensions.

Error in SIREquilibrium>@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2);gamma*y(2)-mu*y(3)] (line 16)
f=@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2);gamma*y(2)-mu*y(3)];

Error in SIREquilibrium (line 25)
Yprime=f(t,[x(i);y(i)]);


I infact do not understand what Yprime=f(t,[x(i);y(i)]); does in this

Can someone please tell me how I can plot the phase plane of this SIR model

After making the changes to the code as f=@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2)]; the phase plane I get is

But, I think it should be

Can you suggest me what I am doing wrong to get a phase plane like that.

I changed the linspace to see why the spiraling in behavior couldn't be seen. And this is the phase plot that I get.

function SIREquilibrium()
if nargin==0
S0=0.8;
I0=0.2;
beta=0.3;
gamma=1/10;
mu=5e-5;
end

f=@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2)];

y1=linspace(0,1,20);
y2=linspace(0,1,20);
[x,y]=meshgrid(y1,y2);
u=zeros(size(x));
v=zeros(size(y));
t=0;
for i=1:numel(x)

Yprime=f(t,[x(i);y(i)]);
Yprime=Yprime/norm(Yprime);
u(i)=Yprime(1);
v(i)=Yprime(2);

end
quiver(x,y,u,v,'r')
axis tight equal
hold on

for y10=[0:0.2:1]
for y20=[0:0.2:1]
options=odeset('MaxStep',0.1);
[ts,ys]=ode45(f,[0,4000],[y10,y20],options);
plot(ys(:,1),ys(:,2))
end
end
hold off

• Add a test if x(i) + y(i) < 1 inside your for loop and follow Ian's suggestion to normalize with Yprime = Yprime / norm(Yprime);. It will look better. – Fabio Somenzi Mar 29 '17 at 23:47
• @FabioSomenzi I normalized and used the if condition. But still I can't see it spiraling towards the equilibrium point. Near the equilibrium point I can't see the behavior correctly. But in the rest of the field it seems to have the correct behavior. – clarkson Mar 30 '17 at 0:21
• Part of the problem is that you are using $R_0 = \beta/\gamma = 10$ instead of $3$ as in the picture. – Fabio Somenzi Mar 30 '17 at 0:32
• @Fabio Somenzi changing R0 doesn't make any difference – clarkson Mar 30 '17 at 0:53
• It does for me. The quiver plot is useful at the larger scale, but at some point you may want to draw some real trajectories. Use [0,4000] for your time span , and [S0;I0] for your initial conditions. Do you know how to set options for ode45 with odeset? You need to limit the maximum time step with something like options = odeset('MaxStep', 0.1); – Fabio Somenzi Mar 30 '17 at 1:03

You passed in a vector of length 2 into the second argument of f and then asked for its third element in the evaluation of f. You can't do that.

Given the rest of the code, my guess would be that you don't actually care about $R$ (because you're using quiver, which will give a 2D phase portrait, not quiver3). So you can just delete that last row of f (the one representing $\frac{dR}{dt}$), and then your code will work. That is, you would replace

f=@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2);gamma*y(2)-mu*y(3)];


with

f=@(t,y)[mu-beta*y(1)*y(2)-mu*y(1);beta*y(1)*y(2)-gamma*y(2)-mu*y(2)];


Note that this makes some sense because only $\frac{dR}{dt}$ involves $R$. You also are not really throwing away information because in your case the relation $S+I+R=1$ is preserved by the dynamics.

The alternative would be to make an actual 3D phase portrait, in which case you would need to pass in a value for R to f in your loop, with something like:

Yprime=f(t,[x(i);y(i);z(i)])


and then use quiver3 at the end. In this case you could still initialize everything using meshgrid; in such a case x,y,z will be tensors of order 3 (i.e. "matrices" depending on three indices).

• Thank you very much for the answer. Since S+I+R=1, can I not consider ${dR\over dt}$ and use only the other two equations. The 3d quiver works but I want to see the phase plane between S and I. Also, since all, S,I,R lie in (0,1) I choose linspace(0,1,20);. How to decide what values to be chosen for linspace. In plotting phase plane should I first calculate the equilibrium point, and then use linspace such that it includes this equilibrium point. – clarkson Mar 29 '17 at 20:33
• @clarkson So I didn't notice the "birth" term $\mu$ in $\frac{dS}{dt}$. With that there you are right that you have a conservation law for the total population so it makes sense to ignore $R$ altogether. In this case you should just delete the third row in the definition of f and then leave it be. Your portrait does not really need to explicitly include the equilibrium; the picture will fairly clearly show it anyway provided your grid is fine enough. – Ian Mar 29 '17 at 22:20
• I have edited my post to include the phase plane I get. But it seems not to be correct. It would be a great help if you an show me what I am doing wrong here – clarkson Mar 29 '17 at 23:22
• @clarkson Your figure looks qualitatively correct but it is really hard to see the small vectors towards the lower corner. You can fix that by normalizing the derivative vectors, so that you show only the direction and not the speed. – Ian Mar 29 '17 at 23:30
• @clarkson For the lines that's another matter entirely, you would need to solve the ODE for various initial conditions to get those. You might want to google pplane, it can do all this sort of stuff quite easily. – Ian Mar 29 '17 at 23:50