SI model


Author: GC

Date: July 19, 2016

Previous Duffing attractor

Next Bautin bifurcation


We consider the SI model (Susceptible-Infected), which is an epidemiological model, given by the following system of parabolic PDE: \[ \begin{cases} \dfrac{\partial S}{\partial t} = \Delta S - \alpha S I \\ \dfrac{\partial I}{\partial t} = \Delta I + \alpha S I - \beta I \end{cases}. \] We solve it using a splitting scheme and the Strang formula, with an initial value that leads to a traveling wave.



SI model




Source code

#!/usr/bin/env python3

"""
SI model
"""

# scientific libraries
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import odeint

# SI model
alpha = 2
beta = 0.6

def reaction(X, t):
    S, I = X
    dS = -alpha*S*I
    dI = +alpha*S*I - beta*I
    return [dS, dI]

# discretization parameters
L = 100.0
T = 100
J = 200
N = 200
dx = L/J
dt = T/N
dtau = dt/2
c = dtau/(dx**2)

# initial condition
x = np.linspace(-L,L,J)
I0 = np.zeros(J)
for i in range(80,100): I0[i] = (i-80)/50
for i in range(100,120): I0[i] = (120-i)/50
S0 = np.ones(J)-I0/2

# matrix of the scheme
diag = np.ones(J)
diagsup = np.ones(J-1)
A = np.diag(diag*(1+2*c),0)+np.diag(diagsup*(-c),1)+np.diag(diagsup*(-c),-1)

#Neumann boundary conditions
A[0][0]=1+c
A[J-1][J-1]=1+c

fig = plt.figure(figsize=plt.figaspect(0.4))
fig.canvas.set_window_title('SI model') 
ax=fig.add_subplot(111)
chrono=fig.suptitle(r'$Time$',fontsize=20)

S = S0
I = I0

invA=np.linalg.inv(A)

for i in range(100):
    ax.plot(x,S,'g')
    ax.plot(x,I,'b')
    """Strang formula"""    
    # diffusion 1/2 step
    S, I = np.dot(invA, S), np.dot(invA, I)
    # reaction 1 step
    for j in range(len(S)):
        X0 = S[j], I[j]
        orbit = odeint(reaction, X0, [0, dt])
        newpoint = orbit[-1]
        S[j], I[j] = newpoint.T
    # diffusion 1/2 step
    S, I = np.dot(invA, S), np.dot(invA, I)

    ax.legend((r"$S(t)$",r"$I(t)$"), shadow = True, loc = (0.85, 0.8))
    ax.set_ylim(0,2)
    chrono.set_text(r"$Time : \,"+str(i)+"$")
    plt.pause(0.01)
    fig.savefig("~/SI"+str(i)+".png")
    ax.clear()


Previous Duffing attractor

Next Bautin bifurcation