## SI model

Author: GC

Date: July 19, 2016

Duffing attractor

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')
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()


Duffing attractor

Bautin bifurcation