3D Phase portrait of SIR model in Python

Question

How to create a 3D plot of a dynamical system with Python? I found a partial solution of my problem in this question. But adding a dimension seems to be bad.. My code

import numpy as np
import matplotlib.pyplot as plt
X, Y, Z = np.meshgrid(np.linspace(-3.0, 3.0, 30), np.linspace(-3.0, 3.0, 30), np.linspace(-3.0, 3.0, 30))
u, v, w = np.zeros_like(X), np.zeros_like(X), np.zeros_like(X)
NI, NJ, NK = X.shape
for i in range(NI):
for j in range(NJ):
 for k in range(NK):        
  x, y, z = X[i,j,k], Y[i,j,k], Z[i,j,k]
  #LOTKA
       #u[i,j] =  x-x*y
       #v[i,j] = -y+x*y
       #SIR
       u[i,j,k] =  -x*y
       v[i,j,k] =  x*y-z  
       w[i,j,k] =  z         

plt.streamplot(X, Y, Z u, v, w)
plt.axis('square')#cube???
plt.axis([-3, 3, -3, 3,-3,3])
plt.show()

Answer 1

This is by no means an answer to your question, but just a suggestion: since looking at 3D phase portraits could be a little bit confusing too, could you simply apply 2D reduction of the portrait onto an invariant 2D plane and look at that? I am asking this, because in general, some SIR = Susceptible-Infected-Recovered epidemiological models, like the one you have presented, have the property that the total population number is fixed, i.e. in your system x + y + z = c for some constant c. This is a type of conservation low (conservation of total population). It allows you to express z = c - x - y and reduce your system of original equations. Here is the original system:

Sucseptible: x,  Infected: y,  Recovered: z  

(the number of susceptible individuals x decreases at a rate which is equal to a portion of all current susceptible individuals x, where this portion is determined by the number of infected individuals y: ) 
dx/dt = - x*y

(the number of infected individuals y increases at a rate which is equal to the rate at which the susceptible individuals convert into infected, i.e. x*y from the first equation, but decreases by a the number z of formerly infected individuals that have just  recovered: )      
dy/dt =   x*y - z

(the rate of change with which infected individuals recover)
dz/dt =   z

Observe that

dx/dt + dy/dt + dz/dt = - x*y + x*y - z + z = 0

which yields, after integrating with respect to t, the conservation of population law:

x + y + z = c

Express z = c - x - y and plug in the second differential equation:

reduced system, showing the conversion dynamics of susceptible to infected individuals:

dx/dt = - x*y
     
dy/dt =   x*y + x + y - c

You can even solve the system explicitly from here. But to look at the dynamics quickly:

import numpy as np
import matplotlib.pyplot as plt

X, Y  = np.meshgrid(np.linspace(-3.0, 3.0, 30), np.linspace(-3.0, 3.0, 30))

# since x + y + z = c, which is an invariant 2D plane, set up an initial c:  
c = 2

# this is simply matrix implementation of your code, since numpy is much faster at 
# calculating matrices than running for loops:
u = -X*Y
v =  X*Y + X + Y - c

plt.streamplot(X, Y, u, v)
plt.axis('square')#cube???
plt.axis([-3, 3, -3, 3])
plt.show()

and got the following phase portrait:

It might be even better, if you want to follow what the dynamics between susceptible and recovered individuals is (or alternatively between infected and recovered), to express y = c - x - z and plug it in the first equation and combine it with the third one.

dx/dt = -c*x + x*x + x*z

dz/dt = z

Since the equation dz/dt = z is decoupled from the other two, it can be solved and then the other two equations can also be solved explicitly. But to get a quick look into the dynamics from another point of view:

import numpy as np
import matplotlib.pyplot as plt

X, Z  = np.meshgrid(np.linspace(-3.0, 3.0, 30), np.linspace(-3.0, 3.0, 30))

c=2

# x + y + z = c
#y = c - x - z 

U = -c*X + X*X + X*Z
V = -Z
 
plt.streamplot(X, Y, U, V)
plt.axis('square')#cube???
plt.axis([-3, 3, -3, 3])
plt.show()

You can do this also for equations 2 and 3 and a third point of view and get a decent idea about the dynamics of all pairs of variables.