[Scilab-users] Covid19 model

classic Classic list List threaded Threaded
1 message Options
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] Covid19 model

Hi Scilabers

A friend (Lars Risbo) published some MATLAB code in LinkedIn for simulating the infection with a company lockdown after some time. I figured I'd try to convert it to Scilab.

The code uses ODE45 and I'm not sure that I understand how to convert this MATLAB code to Scilab. I hope you can explain what I need to do to make the code work. Some of the original MATLAB code is found in the comments.

// covid19risbo.sce
//
// SEIRsim1
// Susceptible-Exposed-Infectious-Recovered (SEIR)

function dydt=odefun(t, y)
    if t<days
      R=R0;
    else
      R=0.6; // Change reproductive rate after lockdown
    end
    A=[0    0   -delta*R*y(1)/Npop  0 0 0;...
      0   -gam delta*R*y(1)/Npop  0 0 0;...
      0    gam -delta        0 0 0;...
      0    0   delta*(1-Fhosp)   0 0 0;...
      0    0   delta*Fhosp     0 -1/Thosp 0;...
      0    0   delta*R*y(1)/Npop  0 0 0];
    dydt=A*y;
end

Begin=datenum(2020,02,20,0,0,0); // begin date
Dlock=datenum(2020,03,12,0,0,0); // date of lockdown
days=Dlock-Begin;
R0=2.6;   // inital R value
gam= 1/3; // gamma, the inverse of average latent time
delta= 1/5; // inv time constant which infectious people either recover or enter hospital
Fhosp=0.16; // fraction of recovering people going to hospital
Thosp=14;  // average time of hospitalisation
Npop=6e6;  // total initial population of sensitive
y0=[Npop;50;50;0;0;0];// [S E I R Hosp ] intial cond.
t1 = 0:60;
y1 = ode(y0, 0, 60, odefun); // [t1,y1] = ode45(@odefun,[0 60],y0); // run 1st scernario
idx=[3:6];
scf();
a = gca();
plot(Begin+t1,y1(:,idx),'-'); // semilogy(Begin+t1,y1(:,idx),'-','LineWidth',3)
a.log_flags = "nln";
xlabel('Date')
ylabel('Number of cases')
xgrid();
// ax=gca;
// ax.YLim=[10 max(max(y2(:,idx)))];
xtitle({'Danish Corona lock down on 20.03.12, vs 14 days later','R goes from 2.6 to 0.6 at lockdown'});
legend('Infected','Recovered','Hospitalised','Total Cases');

P.S. Attached a simple graph of what the output should look like (file covid19risbo.jpg < 50 kb), at least partially - because I deleted an alternative case with a later lockdown date = dashed lines.

Best regards,

Claus


Virus-free. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

covid19risbo.jpg (63K) Download Attachment