SIR models 

If you click on !!! (the menu bar at the top) maple will execute all commands.  

You can also enter each command as you go along  reading this document. You may enter your own commands at any time and place in this document. 

> restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
restart; -1; with(LinearAlgebra); -1; with(plots); -1; with(VectorCalculus); -1; `:=`(SIR, proc (alpha, gamma) description `(`+`(P[1], `-`(`*`..." align="center" border="0">
 

We consider SIR models for a population of total size 

> `:=`(N, 100); 1
 

100 (1)
 

SIR epidemic 

We will go through three examples of epidemic SIR models  with diminishing transmission coefficients α=0.01, 0.005, 0.003  and fixed removal rate γ=0.25.  The fixed removal rate of 0.25 means that the infected individuals are infected for 4 days. 

Example 1: Consider the SIR model with transmission coefficient α=0.01 and removal rate γ=0.25 

> `:=`(F, SIR(0.1e-1, .25)); 1
 

proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(P[1], VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(0.1e-1, P[1]), P[2]))), VectorCalculus:-`+`(VectorCalculus...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(P[1], VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(0.1e-1, P[1]), P[2]))), VectorCalculus:-`+`(VectorCalculus...
(2)
 

has relative removal rate 

> rho = rho(F); 1
 

rho = 25. (3)
 

and contact number 

> sigma = evalf[3](VectorCalculus:-`*`(N, `/`(1, `*`(rho(F))))); 1
 

sigma = 4.00 (4)
 

This disease needs at least 25 susceptibles in order to spread. An invidual in this population has 4 potentially contagious contacts over 4 days. 

Here is the time plot for the population orbit through P0=<95,5,0> 

> display(timeplot(F, `<,>`(95, 5, 0), 40), implicitplot([y = 25, x = 6.81], x = 0 .. 40, y = 0 .. 100, color = magenta)); 1
display(timeplot(F, `<,>`(95, 5, 0), 40), implicitplot([y = 25, x = 6.81], x = 0 .. 40, y = 0 .. 100, color = magenta)); 1
 

Plot_2d
 

The maximal number of infected occurs when S=ρ after 6-7 time steps. 

> I[max] = max(seq(P[2], P = orbit(F, `<,>`(95, 5, 0), 20))); 1
 

I[max] = 47.31512962 (5)
 

and there remain no 

> S[infinity] = iterate(F, `<,>`(95, 5, 0), 30)[1]; 1
 

S[infinity] = .7920061990 (6)
 

susceptibles in the population after the epidemic. 

Here is the directional field in the SI-phase plane 

> display([dirfield(F, 100, 20, N), implicitplot(P[1] = rho(F), P[1] = 0 .. 100, P[2] = 0 .. 20, color = magenta)]); 1
 

Plot_2d
 

and the orbits through various points 

> display([seq(phasedia(F, P0, 20), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(k, 5))), VectorCalculus:-`*`(k, 5), 0), k = 1 .. 10)]), implicitplot(P[1] = rho(F), P...
display([seq(phasedia(F, P0, 20), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(k, 5))), VectorCalculus:-`*`(k, 5), 0), k = 1 .. 10)]), implicitplot(P[1] = rho(F), P...
 

Plot_2d
 

The maximal number of infected lies in the range 45-62 and there are no susceptibles left after the epidemic. 

> I[max] = seq(evalf[3](max(seq(P[2], P = orbit(F, P0, 20)))), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k))), VectorCalculus:-`*`(5, k), 0), k = 1 .. 10)]); 1
I[max] = seq(evalf[3](max(seq(P[2], P = orbit(F, P0, 20)))), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k))), VectorCalculus:-`*`(5, k), 0), k = 1 .. 10)]); 1
 

I[max] = (47.3, 49.2, 50.7, 51.3, 54.3, 55.3, 57.5, 59.9, 61.7, 62.5) (7)
 

> S[infinity] = seq(evalf[3](iterate(F, P0, 30)[1]), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k))), VectorCalculus:-`*`(5, k), 0), k = 1 .. 10)]); 1
 

S[infinity] = (.796, .696, .638, .572, .493, .437, .372, .328, .275, .230) (8)
 

Example 2: The SIR model with transmission coefficient α=0.005 and removal rate ρ=0.25 models an epidemic that is less contagious. 

> `:=`(F, SIR(0.5e-2, .25)); 1
 

proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(P[1], VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(0.5e-2, P[1]), P[2]))), VectorCalculus:-`+`(VectorCalculus...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(P[1], VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(0.5e-2, P[1]), P[2]))), VectorCalculus:-`+`(VectorCalculus...
(9)
 

has relative removal rate 

> rho = rho(F); 1
 

rho = 50.00000000 (10)
 

and contact number 

> sigma = VectorCalculus:-`*`(N, `/`(1, `*`(rho(F)))); 1
 

sigma = 2.000000000 (11)
 

This disease needs at least 50 susceptibles in order to spread.  An invidual in this population has 2 potentially contagious contacts over 4 days. 

Here is the time plot for the population orbit through <95,5,0> 

> display(timeplot(F, `<,>`(95, 5, 0), 40), implicitplot([y = 50, x = 11.5], x = 0 .. 40, y = 0 .. 100, color = magenta)); 1
display(timeplot(F, `<,>`(95, 5, 0), 40), implicitplot([y = 50, x = 11.5], x = 0 .. 40, y = 0 .. 100, color = magenta)); 1
 

Plot_2d
 

The maximal number of infected occurs when S=ϱ 

> I[max] = max(seq(P[2], P = orbit(F, `<,>`(95, 5, 0), 30))); 1
 

I[max] = 19.07431229 (12)
 

and there remain 

> S[infinity] = iterate(F, `<,>`(95, 5, 0), 30)[1]; 1
 

S[infinity] = 18.02571494 (13)
 

susceptibles in the population after the epidemic. 

Here are several population orbits in the SI-phase plane   

> display([dirfield(F, 100, 29, 100), seq(phasedia(F, P0, 50), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(2, k))), VectorCalculus:-`*`(2, k), 0), k = 1 .. 10)]), im...
display([dirfield(F, 100, 29, 100), seq(phasedia(F, P0, 50), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(2, k))), VectorCalculus:-`*`(2, k), 0), k = 1 .. 10)]), im...
 

Plot_2d
 

The maximal number of infected lies in the range 19-50 and there are 6-17 susceptibles left after the epidemic. 

> I[max] = seq(evalf[3](max(seq(P[2], P = orbit(F, P0, 50)))), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k))), VectorCalculus:-`*`(5, k), 0), k = 1 .. 10)]); 1
I[max] = seq(evalf[3](max(seq(P[2], P = orbit(F, P0, 50)))), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k))), VectorCalculus:-`*`(5, k), 0), k = 1 .. 10)]); 1
 

I[max] = (19.2, 21.8, 24.7, 28.1, 31.4, 34.6, 38.3, 42.0, 46.2, 50.) (14)
 

> S[infinity] = seq(evalf[3](iterate(F, P0, 50)[1]), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k))), VectorCalculus:-`*`(5, k), 0), k = 1 .. 10)]); 1
 

S[infinity] = (17.1, 15.8, 14.3, 12.7, 11.6, 10.5, 9.26, 8.25, 7.35, 6.38) (15)
 

>
 

Example 3:The SIR model with transmission coefficient α=0.003 and removal rate γ=0.25 

> `:=`(F, SIR(0.3e-2, .25)); 1
 

proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(P[1], VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(0.3e-2, P[1]), P[2]))), VectorCalculus:-`+`(VectorCalculus...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(P[1], VectorCalculus:-`-`(VectorCalculus:-`*`(VectorCalculus:-`*`(0.3e-2, P[1]), P[2]))), VectorCalculus:-`+`(VectorCalculus...
(16)
 

has relative removal rate 

> `ϱ` = rho(F); 1
 

`ϱ` = 83.33333332 (17)
 

and contact number 

> sigma = VectorCalculus:-`*`(N, `/`(1, `*`(rho(F)))); 1
 

sigma = 1.200000000 (18)
 

This disease needs at least 83 susceptibles in order to spread.  An invidual in this population has 1.2 potentially contagious contacts over 4 days. 

Here is the time plot for the population orbit through <95,5,0> 

> display(timeplot(F, `<,>`(95, 5, 0), 40), implicitplot([y = rho(F), x = 9], x = 0 .. 40, y = 0 .. 100, color = magenta)); 1
display(timeplot(F, `<,>`(95, 5, 0), 40), implicitplot([y = rho(F), x = 9], x = 0 .. 40, y = 0 .. 100, color = magenta)); 1
 

Plot_2d
 

The maximal number of infected occurs when S=ϱ 

> I[max] = max(seq(P[2], P = orbit(F, `<,>`(95, 5, 0), 30))); 1
 

I[max] = 5.838911423 (19)
 

and there remain 

> S[infinity] = iterate(F, `<,>`(95, 5, 0), 30)[1]; 1
 

S[infinity] = 61.42707396 (20)
 

susceptibles in the population after the epidemic. 

Here are several population orbits in the SI-phase plane 

> display([dirfield(F, 100, 22, 100), seq(phasedia(F, P0, 40), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(2, k))), VectorCalculus:-`*`(2, k), 0), k = 1 .. 10)]), im...
display([dirfield(F, 100, 22, 100), seq(phasedia(F, P0, 40), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(2, k))), VectorCalculus:-`*`(2, k), 0), k = 1 .. 10)]), im...
 

Plot_2d
 

The maximal number of infected lies in the range 5-50 and there remain 18-58 susceptibles in the population after the epidemic. 

> I[max] = seq(evalf[3](max(seq(P[2], P = orbit(F, P0, 40)))), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k))), VectorCalculus:-`*`(5, k), 0), k = 1 .. 10)]); 1
I[max] = seq(evalf[3](max(seq(P[2], P = orbit(F, P0, 40)))), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k))), VectorCalculus:-`*`(5, k), 0), k = 1 .. 10)]); 1
 

I[max] = (5.86, 10.4, 15., 20., 25., 30., 35., 40., 45., 50.) (21)
 

> S[infinity] = seq(evalf[3](iterate(F, P0, 40)[1]), P0 = [seq(`<,>`(VectorCalculus:-`+`(100, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k))), VectorCalculus:-`*`(5, k), 0), k = 1 .. 10)]); 1
 

S[infinity] = (58.0, 48.3, 41.9, 37.4, 33.3, 29.3, 26.0, 23.3, 20.5, 18.0) (22)
 

>
 

 

SIR endemic 

> `:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
`:=`(SIRend, proc (alpha, gamma, mu) description `(`+`(`*`(`+`(1, `-`(mu)), `*`(P[1])), `*`(mu, `*`(N)), `-`(`*`(alph..." align="center" border="0">
 

> `:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
`:=`(equi2, proc (F) local equipts, J, dm, EV; description
 

General theory 

> `:=`(F, SIRend(alpha, gamma, mu)); 1
 

proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(mu)), P[1]), VectorCalculus:-`*`(mu, N)),...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(mu)), P[1]), VectorCalculus:-`*`(mu, N)),...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(mu)), P[1]), VectorCalculus:-`*`(mu, N)),...
(23)
 

> equi2(proc (P) options operator, arrow; `<,>`(F(P)[1], F(P)[2]) end proc); 1
 

 

 

 

 

 

`*`(The, `*`(equilibria, `*`(are)))
[[P[1] = 100, P[2] = 0], [P[1] = `/`(`*`(`+`(mu, gamma)), `*`(alpha)), P[2] = `+`(`-`(`/`(`*`(mu, `*`(`+`(mu, gamma, `-`(`*`(100, `*`(alpha)))))), `*`(alpha, `*`(`+`(mu, gamma))))))]]
`*`(The, `*`(eigenvalues, `*`(are)))
Vector[column](%id = 139336668), Vector[column](%id = 140259432)
Vector[column](%id = 139336668), Vector[column](%id = 140259432)
Vector[column](%id = 139336668), Vector[column](%id = 140259432)
Vector[column](%id = 139336668), Vector[column](%id = 140259432)
Vector[column](%id = 139336668), Vector[column](%id = 140259432)
`*`(The, `*`(absolute, `*`(values, `*`(of, `*`(the, `*`(eigenvalues, `*`(are)))))))
Vector[column](%id = 145316132), Vector[column](%id = 143569992)
Vector[column](%id = 145316132), Vector[column](%id = 143569992)
Vector[column](%id = 145316132), Vector[column](%id = 143569992)
Vector[column](%id = 145316132), Vector[column](%id = 143569992)
Vector[column](%id = 145316132), Vector[column](%id = 143569992)
(24)
 

> rhoend(F); 1
 

`/`(`*`(`+`(mu, gamma)), `*`(alpha)) (25)
 

Example 1 

> `:=`(F, SIRend(0.1e-1, .25, 0.5e-1)); 1
 

proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(0.5e-1)), P[1]), VectorCalculus:-`*`(0.5e...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(0.5e-1)), P[1]), VectorCalculus:-`*`(0.5e...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(0.5e-1)), P[1]), VectorCalculus:-`*`(0.5e...
(26)
 

> equi2(proc (P) options operator, arrow; `<,>`(F(P)[1], F(P)[2]) end proc); 1
 

 

 

 

 

 

`*`(The, `*`(equilibria, `*`(are)))
[[P[1] = 100., P[2] = 0.], [P[1] = 30., P[2] = 11.66666667]]
`*`(The, `*`(eigenvalues, `*`(are)))
Vector[column](%id = 144968452), Vector[column](%id = 140304240)
`*`(The, `*`(absolute, `*`(values, `*`(of, `*`(the, `*`(eigenvalues, `*`(are)))))))
Vector[column](%id = 138424844), Vector[column](%id = 138651172) (27)
 

> display([timeplot(F, `<,>`(95, 5, 0), 40), plot(rhoend(F), 0 .. 40, color = magenta)]); 1
 

Plot_2d
 

> display([seq(phasedia(F, P0, 50), P0 = {seq(`<,>`(VectorCalculus:-`*`(5, k), 20, VectorCalculus:-`+`(N, VectorCalculus:-`-`(VectorCalculus:-`+`(20, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k)))))), ...
 

Plot_2d
 

> equi2(proc (P) options operator, arrow; `<,>`(F(P)[1], F(P)[2]) end proc); 1
 

 

 

 

 

 

`*`(The, `*`(equilibria, `*`(are)))
[[P[1] = 100., P[2] = 0.], [P[1] = 30., P[2] = 11.66666667]]
`*`(The, `*`(eigenvalues, `*`(are)))
Vector[column](%id = 139887032), Vector[column](%id = 143279472)
`*`(The, `*`(absolute, `*`(values, `*`(of, `*`(the, `*`(eigenvalues, `*`(are)))))))
Vector[column](%id = 137338928), Vector[column](%id = 137453256) (28)
 

> timeplot(F, `<,>`(10, 90, 0), 40); 1
 

Plot_2d
 

Example 2 

> `:=`(F, SIRend(0.5e-2, .25, 0.5e-1)); 1
 

proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(0.5e-1)), P[1]), VectorCalculus:-`*`(0.5e...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(0.5e-1)), P[1]), VectorCalculus:-`*`(0.5e...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(0.5e-1)), P[1]), VectorCalculus:-`*`(0.5e...
(29)
 

> display([timeplot(F, `<,>`(95, 5, 0), 70), plot(rhoend(F), 0 .. 70, color = magenta)]); 1
 

Plot_2d
 

> display([seq(phasedia(F, P0, 100), P0 = {seq(`<,>`(VectorCalculus:-`*`(5, k), 20, VectorCalculus:-`+`(N, VectorCalculus:-`-`(VectorCalculus:-`+`(20, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k)))))),...
 

Plot_2d
 

> equi2(proc (P) options operator, arrow; `<,>`(F(P)[1], F(P)[2]) end proc); 1
 

 

 

 

 

 

`*`(The, `*`(equilibria, `*`(are)))
[[P[1] = 100., P[2] = 0.], [P[1] = 60., P[2] = 6.666666667]]
`*`(The, `*`(eigenvalues, `*`(are)))
Vector[column](%id = 145272056), Vector[column](%id = 137599716)
`*`(The, `*`(absolute, `*`(values, `*`(of, `*`(the, `*`(eigenvalues, `*`(are)))))))
Vector[column](%id = 143534640), Vector[column](%id = 143517256) (30)
 

Example 3 

> `:=`(F, SIRend(0.35e-2, .25, 0.5e-1)); 1
 

proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(0.5e-1)), P[1]), VectorCalculus:-`*`(0.5e...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(0.5e-1)), P[1]), VectorCalculus:-`*`(0.5e...
proc (P) options operator, arrow; VectorCalculus:-`<,>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(0.5e-1)), P[1]), VectorCalculus:-`*`(0.5e...
(31)
 

> display([timeplot(F, `<,>`(95, 5, 0), 70), plot(rhoend(F), 0 .. 70, color = magenta)]); 1
 

Plot_2d
 

> display([seq(phasedia(F, P0, 100), P0 = {seq(`<,>`(VectorCalculus:-`*`(5, k), 30, VectorCalculus:-`+`(N, VectorCalculus:-`-`(VectorCalculus:-`+`(30, VectorCalculus:-`-`(VectorCalculus:-`*`(5, k)))))),...
 

Plot_2d
 

> equi2(proc (P) options operator, arrow; `<,>`(F(P)[1], F(P)[2]) end proc); 1
 

 

 

 

 

 

`*`(The, `*`(equilibria, `*`(are)))
[[P[1] = 100., P[2] = 0.], [P[1] = 85.71428571, P[2] = 2.380952381]]
`*`(The, `*`(eigenvalues, `*`(are)))
Vector[column](%id = 145385568), Vector[column](%id = 141060764)
`*`(The, `*`(absolute, `*`(values, `*`(of, `*`(the, `*`(eigenvalues, `*`(are)))))))
Vector[column](%id = 145081968), Vector[column](%id = 145065736) (32)
 

 

>
 

>