SIR Epidemic Model

# Population size N = 10000 # Initial infections IInit = 1 SInit = N - IInit RInit = 0 # Transmission rate beta = 0.1 # Recovery rate gamma = 0.01 R0 = beta / gamma show(r'$R_0 = %g$' % R0) # End time tMax = 1000 # Standard SIR model def ODE_RHS(t, Y): (S, I, R) = Y dS = - beta * S * I / N dI = beta * S * I / N - gamma * I dR = gamma * I return (dS, dI, dR) # Set up numerical solution of ODE solver = ode_solver(function = ODE_RHS, y_0 = (SInit, IInit, RInit), t_span = (0, tMax), algorithm = 'rk8pd') # Numerically solve solver.ode_solve(num_points = 1000) # Plot solution show( plot(solver.interpolate_solution(i = 0), 0, tMax, legend_label = 'S(t)', color = 'green') + plot(solver.interpolate_solution(i = 1), 0, tMax, legend_label = 'I(t)', color = 'red') + plot(solver.interpolate_solution(i = 2), 0, tMax, legend_label = 'R(t)', color = 'blue') ) 

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{$R_0 = 10$}

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{$R_0 = 10$}

var('beta gamma S I R N') 

(beta, gamma, S, I, R, N)

(beta, gamma, S, I, R, N)

dS = - beta * S * I / N dI = beta * S * I / N - gamma * I dR = gamma * I 

solve([dS == 0, dI == 0, dR == 0], [S, I, R]) 

[[S == r5, I == 0, R == r6]]

[[S == r5, I == 0, R == r6]]

J = jacobian([dS, dI, dR], [S, I, R]) 

J 

[       -I*beta/N        -S*beta/N                0]
[        I*beta/N S*beta/N - gamma                0]
[               0            gamma                0]

[       -I*beta/N        -S*beta/N                0]
[        I*beta/N S*beta/N - gamma                0]
[               0            gamma                0]

J0 = J.subs(S = N, I = 0, R = 0) 

J0.eigenvalues() 

[beta - gamma, 0, 0]

[beta - gamma, 0, 0]