sir4.html

sir4.mws

This Maple worksheet derives the diffusion approximation of a bivariate Markov population process that serves as an SIR model, accounting for differential death rate. NOTE that this model uses the proper infection rate!

The file "abc.m" contains the three procedures equ, cov, denote used below. The source code for these procedures is contained in the Maple worksheet abc.mws

> restart;
with(linalg,matrix,jacobian);
read "C:/mydocu~1/ingemar/maple/desto/abc.m";
unprotect(gamma);

Reparametrize! Express old parameters in terms of new ones:

> gamma:=mu*(alpha1-delta1); beta:=mu*alpha1*R0; mu1:=mu*(delta1-1);

The transition rates (for scaled state variables x1=S/N, x2=I/N, x3=R/N) are stored in the table trans :

> trans:=table([[1,0,0]=mu,[-1,1,0]=beta*x1*x2/(x1+x2+x3),[-1,0,0]=mu*x1,[0,-1,1]=gamma*x2,[0,-1,0]=(mu+mu1)*x2,[0,0,-1]=mu*x3]);

The RHSs of the three deterministic ODEs are determined via the procedure equ :

> eq1:=equ(1,trans);
eq2a:=equ(2,trans);
eq3:=equ(3,trans);

Check the op # here!

> eq2:=op(1,eq2a)+simplify(eq2a-op(1,eq2a));

The local covariance matrix Sx is determined by the procedure cov:

> Sx:=cov(trans);

Check the op # here!

> Sx[2,2]:=op(1,Sx[2,2])+simplify(Sx[2,2]-op(1,Sx[2,2]));

Critical points:

> crit:=[solve({eq1,eq2,eq3},{x1,x2,x3})];

Denote the critical point that corresponds to endemic infection by x10, x20, x30:

> xx:=denote(crit):
x10:=xx[1]; x20:=xx[2]; x30:=xx[3];

The Jacobian:

> Bx:=jacobian([eq1,eq2,eq3],[x1,x2,x3]);

Evaluate the Jacobian at the critical point:

> B:=simplify(subs({x1=x10,x2=x20,x3=x30},evalm(Bx)));

Evaluate the local covariance matrix at the critical point:

> S:=simplify(subs(x1=x10,x2=x20,x3=x30,evalm(Sx)));

Now proceed to solve B*SIG+SIG*BT=-S:

> SIG:=matrix([[s11,s12,s13],[s12,s22,s23],[s13,s23,s33]]):

> A:=evalm(B&*SIG+SIG&*transpose(B)):

A and S are symmetric. Formulate 6 equations:

> e1:=A[1,1]=-S[1,1]: e2:=A[1,2]=-S[1,2]: e3:=A[1,3]=-S[1,3]:
e4:=A[2,2]=-S[2,2]: e5:=A[2,3]=-S[2,3]: e6:=A[3,3]=-S[3,3]:

Solve for the elements of the covariance matrix SIG:

> solve({e1,e2,e3,e4,e5,e6},{s11,s12,s13,s22,s23,s33}):

> assign(%);

Asymptotic approximations for large. Use subscript a.

> x10a:=map(simplify,convert(asympt(x10,alpha1,2),polynom));
x20a:=map(simplify,convert(asympt(x20,alpha1,3),polynom));
x30a:=map(simplify,convert(asympt(x30,alpha1,2),polynom));