kill(all)$

biotech : S,I ,R (+Diffusion) =travelling wave model,version 1.0
peter.vlasschaert@gmail.com,01/05/2020

see: biotech : S,I incomplete model,used phase analysis,version 1.0
                        peter.vlasschaert@gmail.com,16/04/2020
====================================================
→ General :compartmental modelling  (part 1)

General :compartmental modelling  (part 2)
=================================

ratprint;

change option :

ratprint:false;

load packages
=============

load(lrats)$

load(eigen)$

load(contrib_ode)$

load(vect)$

 A)   SIR ( = without diffusion)
=========================
=========================
ode =(o)rdinary (d)ifferential (e)quations

ref:youtube :
Oxford Mathematician explains SIR  Disease Model for COVID-19
→   https://www.youtube.com/watch?v=NKMHhm2Zbkw
rem :  r = β,μ = a

Mathematical model  (S,I,R), complete ?
=============================================================
?  see part 1

symbol  : S= (s)usceptible  (= not yet infected,allowed)
                 I  = (i)nfected
                 R = (r)ecoverd
S,I ,R are : compartments (= three compartments model)

S,I ,R (complete three compartmentmodel,β,μ )
====================================
β: transmission coef.
μ: recovery rate (= dead or be cured)

1e) compartment =S ,2e) compartment =I,3e) compartment =R

info : (out) compartment  = neg sign
         (in   ) compartment =  pos sign

system nonlinear differential equations. (S,I,R)
==============================

p1a:'diff(S,t) =-β*I*S;        /*f1(x,y)=f1(S,I,R),ref : S compartment*/

p1a1:rhs(p1a);

p1b:'diff(I,t) = β*I*S-μ*I;   /*f2(x,y)=f2(S,I,R),ref; I compartement*/

p1b1:rhs(p1b);

p1c:'diff(R,t) = μ*I;                /*f3(x,y)=f3(S,I,R),ref; R compartement*/

p1c1:rhs(p1c);

rem: S ,I ,R= people  > 0  →   S>0,I>0,R>0
         t=time :S=S(t)    →  ic:  S(0) =S0 
         t=time :I=I(t)        →  ic:   I(0) = I0 
         t=time :R=R(t)    →  ic:  R(0) =R0
==============================================================
ic : (i)nitial (c)ondition
==============================================================
S(t) : # number susceptable @ time t
I(t)   : # number infected         @ time t
R(t)  : # number recovery       @ time t

 B)   SIR ( + diffusion)  = trave
=================
=================
pde =  partial differential equations

ref:youtube :
Oxford Mathematician explains SIR Travelling Wave Disease Model for COVID-19
→   https://www.youtube.com/watch?v=uSLFudKBnBI
rem :  r = β,μ = a

1e) simple : diffusion equation  (pde )
==========================
laplacian = Δψ = grad(div(ψ))
ref :https://en.wikipedia.org/wiki/Laplace_operator

------------------------------------------------------------------------------
 →  ψ = ψ(space,time)  →  ' pde ', multivariable
 →  ψ = ψ (time)             → ' ode ' , onevariable
------------------------------------------------------------------------------
∂ψ/∂t = D*Δψ,  model ,  ψ(x,y,z,t),cartensian coordinates
------------------------------------------------------------------------------
ref :https://en.wikipedia.org/wiki/Diffusion

 Δ = Laplacian ( = differential operator)
==============================

(f=ψ,cartesian coordinates), ψ = ψ(x,y,z)
-------------------------------------------------------

p2a:Δψ = 'diff(ψ,x,2)+ 'diff(ψ,y,2)+ 'diff(ψ,z,2);

grad(f)
---------

depends(f,[x,y,z]);
p2a1:express(grad(f));

div(A),  A = vector=(A[x],A[y],A[Z])
-------------------------

p2a2:declare(A,nonscalar);
p2a3:express(A);

p2a4:express(div(A));

Δf =grad(div(f))
-----------------

p2a5:map("=",p2a3,p2a1);

p2a6:Δf=subst(p2a5,p2a4);

Δf ( = only space depent)
-----------------------------------

p2a7:Δf=express(laplacian(f));

check :true or false
---------------------------

p2a8:is(rhs(p2a6)=rhs(p2a7));

maxima : no pde tools
=> need to use mathematica or maple (symbolic)

model :diffusion
========================
D = diffusion coef., ψ =ψ(x,y,z,t)

p2a:'diff(ψ,t) = D*subst([f=ψ],rhs(p2a7));

use : grad ,div,laplacian,express , you need to load(vect),packages

intermezzo :vector analysis
=====================

--------------------------------------------------------------------------------------------------

ref :https://donnishjournals.org/djerr/pdf/2014/november/Savas-Tuylu.pdf

---------------------------------------------------------------------------------------------------

2e) combination(  =  SIR +Diffusion)
============================
New model : Diffusion model 
----------------------------------------------------------
relation : differentials → partial differentials
-----------------------------------------------------------
S = S(x,y,z,t)
 I =  I(x,y,z,t)
R = R(x,y,z,t)

p3a:lhs(p1a) =rhs(p1a)+D*ΔS;
p3b:lhs(p1b) =rhs(p1b)+D*ΔI;
p3c:lhs(p1c) =rhs(p1c)+D*ΔR;

p3a1:ΔS=subst([f=S],rhs(p2a7));
p3b1:ΔI=subst([f=I],rhs(p2a7));
p3c1:ΔR=subst([f=R],rhs(p2a7));

p3a2:subst([p3a1],p3a);
p3b2:subst([p3b1],p3b);
p3c2:subst([p3c1],p3c);

dimensional model :sir +diffusion
----------------------------------------------

two variable model : x,t
---------------------------------
S = S(x,t)
 I =  I(x,t)
R = R(x,t)
--------------------------------

p3a3:lhs(p3a1)=part(rhs(p3a1),3);
p3b3:lhs(p3b1)=part(rhs(p3b1),3);
p3c3:lhs(p3c1)=part(rhs(p3c1),3);

p3d:[p3a3,p3b3,p3c3];

p3d1:makelist(rhs(p3d[i]) =lhs(p3d[i]),i,1,length(p3d));

use :  {p3a,p3b,p3c}

ref : 2d dimensional diffusion sir model = 2ddifsirmodel

p3a4:subst([p3a3],p3a);
p3b4:subst([p3b3],p3b);
p3c4:subst([p3c3],p3c);

?  How to get nondimensional model = 2dnondifsirmodel
============================================

1e)  nondimensional variables
========================
density = S0 , 'susceptible part', use ref : nondimensional variable
------------------------------------------------------------------------------------------
dot variables : nondimensional
-------------------------------------------
→ S* = Sdot 
→  I* = Idot 
→ R* = Rdot ,
→ x*  = xdot
→ t*  = tdot
-----------------------------------------

p4a:Sdot = S/S0;
p4b:Idot=I/S0;
p4c:Rdot=R/S0;
p4d:xdot =radcan(sqrt(β*S0/D)*x);
p4e:tdot =β*S0*t;

2e)  find differential of time
=====================
-------------------------------------------------------------
first order :differential:∂(S)/∂t →∂(Sdot)/∂tdot
-------------------------------------------------------------
dS/dt →S0*d(Sdot)/dt →S0*(S0*β)*d(Sdot)/dtdot

p5a:'diff(S,t)=S0*(S0*β)*'diff(Sdot,tdot);

p5b:subst([S=I,Sdot=Idot],p5a);
p5c:subst([S=R,Sdot=Rdot],p5a);

p5:[p5a,p5b,p5c];

3e)  find differential of space (part after D,)
=================================

-----------------------------------------------------------------------------------
second order : diff(ψ,x,2) =diff(diff(ψ,x,1),x,1) = diff(diff(ψ,x),x)
a)  first order :space = x
b)    a)           : space = x
------------------------------------------------------------------------------------

a) first order to x
-----------------------
dS/dx→S0*d(Sdot)/dx →S0*((sqrt(S0)*sqrt(β))/sqrt(D))*d(Sdot)/dxdot

p6a:'diff(S,x) = S0*((sqrt(S0)*sqrt(β))/sqrt(D))*'diff(Sdot,xdot);

p6b:subst([S=I,Sdot=Idot],p6a);
p6c:subst([S=R,Sdot=Rdot],p6a);

b) first order(='a)')  to x
-------------------------------
d(dS/dx)/dx

p6a1:'diff(S,x,2) ='diff(rhs(p6a),xdot)*'diff(rhs(p4d),x);

p6b1:subst([S=I,Sdot=Idot],p6a1);
p6c1:subst([S=R,Sdot=Rdot],p6a1);

depends(Sdot,xdot,Idot,xdot,Rdot,xdot);

p6a2:lhs(p6a1)=ev(rhs(p6a1),diff);
p6b2:lhs(p6b1)=ev(rhs(p6b1),diff);
p6c2:lhs(p6c1)=ev(rhs(p6c1),diff);

p6c3:[p6a2,p6b2,p6c2];

4e) final differential equations : nondimensional
=====================================
use :{p3a4,p3b4,p3c4}

p7:append(p5,p6c3);

p7a:subst(p7,p3a4);
p7b:subst(p7,p3b4);
p7c:subst(p7,p3c4);

p7a1:expand(p7a)/((S0)^2*β)$
p7a2:lhs(p7a1) = expand(rhs(p7a1));
p7b1:expand(p7b)/((S0)^2*β)$
p7b2:lhs(p7b1) = expand(rhs(p7b1));
p7c1:expand(p7c)/((S0)^2*β)$
p7c2:lhs(p7c1) = expand(rhs(p7c1));

use:{p4a,p4b,p4c}

p8:[rhs(p4a)=lhs(p4a),rhs(p4b)=lhs(p4b),rhs(p4c) =lhs(p4c)];

p8a:lhs(p7a2) =expand (lratsubst(p8,rhs(p7a2)));
p8b:lhs(p7b2) =expand (lratsubst(p8,rhs(p7b2)));
p8c:lhs(p7c2) =expand (lratsubst(p8,rhs(p7c2)));

reproduction rate : RO  
------------------------------

p9a:RO = (β*S0)/μ;

p9a1:rhs(p9a) = lhs(p9a);

lambda = λ  = 1/RO
----------------------------

p9b:λ =( rhs(p9a))^(-1);

p9b1:[rhs(p9b) = lhs(p9b)];

λ  replacement ,β,μ,D
------------------------------

p10a:lratsubst(p9b1,p8a);
p10b:lratsubst(p9b1,p8b);
p10c:lratsubst(p9b1,p8c);

5e) final differential equations : nondimensional
=====================================
simplification :remove  xdot →x,.......
 → S = S(x,t)
 → I   =  I(x,t)
 → R = R(x,t)

p11:[Sdot =S,Idot=I,Rdot=R,tdot=t,xdot=x];

p11a:subst(p11,p10a);
p11b:subst(p11,p10b);
p11c:subst(p11,p10c);

Travelling wave  1D
===============
===============
I(x,t)=I(q)  , I →R→S  => q=x-c*t
--------------------------------------------
x = space
t  = time
c = wave speed
-------------------------------------------

p12:q=x-c*t;

→ differentation with respect to time
--------------------------------------------------

p12a:'diff(q,t) =diff(rhs(p12),t);

p12a1:'diff(ψ,t) = 'diff(ψ,q)*'diff(q,t);

p12a2:subst([p12a],p12a1);

→ differentation with respect to x
----------------------------------------------

p12b:'diff(q,x) =diff(rhs(p12),x);

p12b1:'diff(lhs(p12b),x) = diff(rhs(p12b),x);

p12c:'diff(ψ,x) =  'diff(ψ,q)*'diff(q,x);;

p12c1:lhs(p12c)=subst([p12b],rhs(p12c));

p12c2:'diff(ψ,x,2) ='diff(rhs(p12c1),q)*'diff(q,x);

p12c3:lhs(p12c2)=subst([p12b],rhs(p12c2));

→  ψ →{S,I,R}
--------------------

p12c4:[S,I,R];

=>  time :list:use p12a2

p13a:makelist(subst([ψ=p12c4[i]],p12a2),i,1,length(p12c4));

 => space:list use p12c3

p13b:makelist(subst([ψ=p12c4[i]],p12c3),i,1,length(p12c4));

p13c:append(p13a,p13b);

final odes : wave equations (sir +wave :model)
----------------------------------------------------------------
use :{p11a,p11b,p11c}

p14a:subst(p13c,p11a);
p14b:subst(p13c,p11b);
p14c:subst(p13c,p11c);

p14a1:rhs(p14a)-lhs(p14a)=0;
p14b1:rhs(p14b)-lhs(p14b)=0;
p14c1:rhs(p14c)-lhs(p14c) =0;

How find  relation between : c  &  λ
===========================
use: p14b1

p15:p14b1;

when : q →∞ ,then S →1 and I →0
-----------------------------------------------
use : ordinary ode  , when we set  S= 1

assume(c>0);

p15a:subst([S=1],p15);

solution :linearised equation ( = p15)
--------------------------------------------------
relation 

p15a1:assume(λ<1);

reason : how to find this , maxima ask question about values.

p15a2:assume(c^2>4*(1-λ));

p15a3:p15a2[1];

p15a31:radcan(sqrt(p15a3));

p15a32:sqrt(lhs(p15a3)) > sqrt(rhs(p15a3));

→ limiting case

p15a33:sqrt(lhs(p15a3)) = sqrt(rhs(p15a3));

covid-19: lancet : RO :2.68, def : λ= 1/RO <1,see : parameters sir model.
---------------------------------------------------------------------------------------------------

→  otherwise : will be complex (not allowed,no travelling wave (= λ>1)):

p15a34:subst([λ=1/RO],p15a33);

p15a34a:rhs(p15a34)=lhs(p15a34);

p15b:ode2(p15a,I,q);

→  use:                  p15a33 or p15a34 (and use p14a1,p14b1,p14c1):

p16a:subst([p15a34],p14a1);
p16b:subst([p15a34],p14b1);
p16c:subst([p15a34],p14c1);

rem :second order differential :reason 'diffusion part'
========================================
video : youtube ' travelling wave model',":see above"
-----------------------------------------------------
constraint :
--------------
use : only p16a(= no diffusion),p16b(= diffusion)
result :'diff(S,q,2) =0

rem :factorsum

negdistrib;

p17:['diff(S,q,2)=0];
p17a:subst(p17,p16a);
p17a1:lratsubst([p15a34a],p17a);
p17a2:factor(lhs(p17a1))=0;
p17b:p16b;
p17b1:lratsubst([p15a34a],p17b);
p17b2:factorsum(lhs(p17b1))=0;

BC:past    : t → -∞,  q=x-c*t = ∞, I = 0,S=1    (no infected,S=all),
BC:future  : t →∞  ,  q=x-c*t = -∞, I=0, S= ?   (no ifected  ,S = ?)
--------------------------------------------------------------------------------
video : q=y,above :relation 'c & RO or λ'
--------------------------------------------------------------------------------
use : BC past  q =∞, I = 0,S=1  => dI/dq = 0

→ solve :p17a2 for I

p17b3:solve(p17a2,I);

→ first part:lhs from(p7b2)

p17b4a:args(lhs(p17b2)) ;

p17b4b:subst(p17b3,p17b4a[1]);

p17b4c:[p17b4b,p17b4a[2],p17b4a[3]];

p17b4d:sum(p17b4c[i],i,1,length(p17b4c)) =rhs(p17b2);

p17b4e:expand(lhs(p17b4d)) = 0;

p17b4f:lratsubst([c*λ=φ*log(S)],integrate(lhs(p17b4e),q)) =cst;

p17b4g:subst([φ=c*λ],p17b4f);

info : I''+c*I'+c*S'-c*λ/S*S' =0
=> integration 'respect to q'
d/dq(I'+c*I+c*S)-λ*c/S*S' *dq=0
=>integration
1/S*S'*dq = log(S)

find :cst , I=0,S=1
-----------------------

p17b4h:subst([I=0,S=1],p17b4g);

p17b4h1:rhs(p17b4h)=lhs(p17b4h);

youtube :'travelling wave'
? what is severity of the epidemic 
==========================
S(end)= Send, I=0

p17b4h2:subst([I=0,p17b4h1,S=Send],p17b4g)/c;

p17b4h3:expand(lhs(p17b4h2))=rhs(p17b4h2);

p17b4h4:solve(p17b4h3,Send);

-------------------------------------------------------------------------------------------------------------------------------------------------------