kill(all)$

biotech : S,I ,R complete model,version 1.0
peter.vlasschaert@gmail.com,16/04/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)$

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

p1cc :p1a+p1b+p1c;

def : conservative system = rhs(p1cc) = 0,
see : part 1

→  #  people
------------------

assume(S>0,I>0,R>0);

(i)nitial (c)ondition:ic
----------------------------
relative : 0 <= ..<= 1

p1dd:[I=I0,S=S0,R=R0];

total  population:N
-------------------------

p1ee:S+I+R = N;

p1ee1:rhs(p1ee)=lhs(p1ee);

parameter : μ,β>0
-------------------------

assume(β>0,μ>0);

rem :def  :RO = reproduction rate , COVID-19 = 3.0, found on youtube
→ General :compartmental modelling  (part 1)

==========================================================

end model ...............................................................................................................

simple model : time evolution
----------------------------------------
from elementary : N = constant = total population
suppose :  N =S  &  I <<<< N.( = true COVID19 virus)
reason : suppose I = small
extra : R0=0 , initial condition

simple model:time evolution
---------------------------------------

p2:N=S;

1e)  I,S
======
→ General :compartmental modelling  (part 1)

p2a:I=I0*%e^((RO-1)*t*μ);

p2b:S=S0*%e^((I0*β)/(RO*μ-μ)-(I0*β*%e^(RO*t*μ-t*μ))/(RO*μ-μ));

2e) R
=====
use : p2a,p2b

a) differential equation ,'p1c'
---------------------------------------

p2c:p1c;

p2c1:subst([p2a],p2c);

b) solved : differential equation ,'p2c1'
----------------------------------------------------

p2c2:ode2(p2c1,R,t);

c) find :%c ,R=R(t=0) =R0
------------------------------------

p2c3:subst([R=R0,t=0],p2c2);

p2c4:solve(p2c3,%c);

d) final result 
------------------
use: p2c2

p2c5:subst(p2c4,p2c2);

→  R0 = 0 ,simple model above.

p2c6:subst([R0=0],p2c5);

values : parameter,ic condition.
========================
→ General :compartmental modelling  (part 1)

→ COVID-19 =RO =3

a1) parameter
--------------------

p3a:[μ=.5,β=2,RO=3];

a2) ic condition
---------------------

p3b:[I0=0.1,S0=0.9,R0=0];

a3)  a1 & a2
------------------

p3c:flatten([p3a,p3b]);

Model time evolution : SIR model
==========================

p4a:p2a;

p4a1:rhs(p4a);

p4a2:subst(p3c,p4a1);

wxplot2d([%], [t,0,20])$

p4b:p2b;

p4b1:rhs(p4b);

p4b2:subst(p3c,p4b1);

wxplot2d([%], [t,0,20])$

p4c:p2c5;

p4c1:rhs(p4c);

p4c2:subst(p3c,p4c1);

wxplot2d([%], [t,0,20])$

combination: plot: SIR model
=======================

p5:[p4a2,p4b2,p4c2];

p5a:wxplot2d(%, [t,0,20])$

we can't see :S parameter,'see above',p5a

0 to 20 days
==========
change : y axis in the plot

p5a:wxplot2d(p5, [t,0,20],[y,0,1])$

we see all elements of the model : SIR
blue :I
red  :S
green :R

0 to 5 days
=========
change : y axis in the plot 

p5a:wxplot2d(p5, [t,0,5],[y,0,1])$

numerical differential equations
========================
with 'rk', maxima command , ( = RUNGE KUTTA  45)

syntax :
Function: rk  
    rk (ODE, var, initial, domain)  
    rk ([ODE1, …, ODEm], [v1, …, vm], [init1, …, initm], domain)

HOW TO USE:

1e)  dx/dt = 4-x^2-4*y^2     dy/dt = y^2-x^2+1
2e)  sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$
3e)  plot2d ([discrete,makelist([p[1],p[3]],p,sol)], [xlabel,"t"],[ylabel,"y"])$

1e)  dx/dt = 4-x^2-4*y^2     dy/dt = y^2-x^2+1
==================================

a) SIR (rhs: differential equations)
======
use :p1a1,p1b1,p1c1

p6a1:[p1a1,p1b1,p1c1];

→ values, p3c

p6a2:subst(p3c,p6a1);

b) ic : condition
============

p6a3:[S,I,R];

p6a4:subst(p3c,[S0,I0,R0]);

c) simulation time limits
==================
[t,t_begin,t_end]

p6a5:[t,0,20,0.1];

2e)  sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$
=================================================

sol:rk(p6a2,p6a3,p6a4,p6a5)$

3e)  plot2d ([discrete,makelist([p[1],p[3]],p,sol)], [xlabel,"t"],[ylabel,"y"])$
======================================================

wxplot2d ([discrete,makelist([p[1],p[2]],p,sol)], [xlabel,"t"],[ylabel,"S"])$

wxplot2d ([discrete,makelist([p[1],p[3]],p,sol)], [xlabel,"t"],[ylabel,"I"])$

wxplot2d ([discrete,makelist([p[1],p[4]],p,sol)], [xlabel,"t"],[ylabel,"R"])$

combination : SIR one plot   (COVID-19 ,RO =3)
=====================================

p7a:[discrete,makelist([p[1],p[2]],p,sol)]$

p7b:[discrete,makelist([p[1],p[3]],p,sol)]$

p7c:[discrete,makelist([p[1],p[4]],p,sol)]$

wxplot2d ([p7a,p7b,p7c], [xlabel,"t"])$

info:
===

=====================================
I=I(t)     ,red     →infected        ,(relative),discrete2
S=S(t)  ,blue   →susceptable ,(relative),discrete1
R=R(t)  ,green →recovered    ,(relative),discrete3
======================================

wxplot2d ([p7a,p7b,p7c], [xlabel,"t"],[ylabel,"S,I,R"])$

wxplot2d ([p7a,p7b,p7c], [xlabel,"t=time"],[ylabel,"S,I,R"], [legend, "S", "I","R"],grid2d)$

ref:see manuel 'Maxima'
==================

==========================================
(%i1) xy: [[10, .6], [20, .9], [30, 1.1], [40, 1.3], [50, 1.4]]$
(%i2) plot2d([[discrete, xy], 2*%pi*sqrt(l/980)], [l,0,50],
        [style, points, lines], [color, red, blue],
        [point_type, asterisk],
        [legend, "experiment", "theory"],
        [xlabel, "pendulum's length (cm)"],
        [ylabel, "period (s)"])$
===========================================

xlabel,ylabel,legend
===============

p8a:[xlabel,"t=time"];

p8b:[ylabel,"S,I,R"];

p8c:[legend, "S", "I","R"];

style,color
========

p8d:[style,[points,0.5,1],[lines,0.5],[points,0.5]];

p8e:[color,blue,red,green];

plot
===

wxplot2d ([p7a,p7b,p7c], p8a,p8b,p8c,p8d,p8e,grid2d)$

programming euler method in ,'Maxima'
===============================
info :see:https://rosettacode.org/wiki/Euler_method#Maxima

1e) step : euler , function of the differential
----------------------------------------------------------
use :rhs,p1a,p1b,p1c

sdot(S,I,R,β,μ):=''p1a1;

idot(S,I,R,β,μ):=''p1b1;

rdot(S,I,R,β,μ):=''p1c1;

2e) myeuler
-----------------

number of steps
=============

na:40$

algo: Euler
========

myeuler(F,G, H,a0,a1,  a2, a3,h, b1,b2) := block(
      t:a0,
      SS:a1,
      II:a2,
      RR:a3,
      f1: F(SS,II,RR,b1,b2),
      f2: G(SS,II,RR,b1,b2),       
      f3: H(SS,II,RR,b1,b2),            
     SS: SS+ f1,
      vx[i + 1]: SS,
      II: II+ f2,
      vy[i + 1]: II,
      RR:RR+ f3,
      vz[i + 1]: RR,
      t:t+h ,         
      vt[i + 1]: t)$

parameter
========

β:2;

μ:.5;

  vt: makelist(0, na + 1)$
   vx: makelist(0, na + 1)$
   vy: makelist(0, na + 1)$
   vz:makelist(0, na + 1)$
   vt[1]:0;
   vx[1]: 0.9;
   vy[1]: .1;
   vz[1]: 0;
   for i from 1 thru na  do
                                     (
   myeuler(sdot,idot, rdot,vt[i],vx[i],  vy[i], vz[i],.01, β,μ) 
                                     )$

 3e) output:simulation
-------------------------------

t :time step
=========

vt;

S:susceptible
==========

vx;

wxplot2d([discrete, vt,vx])$

I:infected
==========

vy;

wxplot2d([discrete, vt,vy]);

R:susceptible
==========

vz;

wxplot2d([discrete, vt,vz]);

Combination :plot (S,I,R)
===================

pa:[[discrete,vt,vx],[discrete,vt,vy],[discrete,vt,vz]]$

wxplot2d(pa,p8a,p8b,p8c);