3
$\begingroup$

I am trying to solve and plot delayed SIRD model. My code is:

S4 = 
  NDSolve[
    {x'[t] == - (0.3/(80*10^6)) x[t] y[t], 
     y'[t] == (0.3/(80*10^6)) x[t] y[t] - 
                 (0.3/(80*10^6)) x[t - 14] y[t - 14] - 0.01 y[t], 
     z'[t] == (0.3/(80*10^6)) x[t - 14] y[t - 14], a'[t] == 0.01 y[t],  
     x[0] == 80*10^6, y[t /; t <= 0] == E^t, y[0] == 150, z[0] == 0, 
     a[0] == 0 }, {x, y, z, a}, {t, 0, 200}]`

In the solution I am getting the error message:

NDSolve::ndsz: At t == 164.8072392646089, step size is effectively zero; singularity or stiff system suspected.

enter image description here

And this results in values of graph that shouldnt be possible, active cases can't be negative.

What am I doing wrong?

$\endgroup$
5
$\begingroup$

I think there's a subtle inconsistency in your equations. You have to discount the recovery rate by the survival probability of infecteds, otherwise you double-count removals from y (both due to death and recovery), hence the negative y[t]. Luckily that is easy to do because survivorship is an exponential of the death rate and the infection time:

S4 = NDSolve[{
   x'[t] == -(0.3/(80*10^6)) x[t] y[t], 
   y'[t] == (0.3/(80*10^6)) x[t] y[t] - (0.3/(80*10^6)) x[t - 14] y[t - 14] E^(-0.01*14) - 0.01 y[t], 
   z'[t] == (0.3/(80*10^6)) x[t - 14] y[t - 14] E^(-0.01*14), 
   a'[t] == 0.01 y[t], x[0] == 80*10^6, y[t /; t <= 0] == E^t, 
   y[0] == 150, z[0] == 0, a[0] == 0}, {x, y, z, a}, {t, 0, 200}];

Plot[Evaluate[{x[t], y[t], z[t], a[t]} /. S4], {t, 0, 200}]

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ ah thanks.it seems was indeed double counting recoveries. $\endgroup$ Nov 30 '20 at 17:11
2
$\begingroup$

Remove y[0] == 150and set $\tau=2$ and replace $y(t \cdots)$ by $y(t-\tau \cdots)$ can remove some warm message.

τ = 2; 
S4 = 
 NDSolve[{x'[t] == -(0.3/(80*10^6)) x[t] y[t - τ], 
   y'[t] == (0.3/(80*10^6)) x[t] y[t - τ] - (0.3/(80*10^6)) x[
       t - 14] y[t - τ - 14] - 0.01 y[t - τ], 
   z'[t] == (0.3/(80*10^6)) x[t - 14] y[t - τ - 14], 
   a'[t] == 0.01 y[t - τ], x[0] == 80*10^6, 
   y[t /; t <= 0] == E^t, z[0] == 0, a[0] == 0}, {x, y, z, a}, {t, 0, 
   200}]
Plot[{x[t], y[t], z[t], a[t]} /. S4 // Evaluate, {t, 0, 200}, 
 AxesOrigin -> {0, 0}, PlotRange -> All, AspectRatio -> 1]
$\endgroup$
1
  • $\begingroup$ It still blows up, y is gets very negative... (Version 12.1.1 Mac OS X.) $\endgroup$ Nov 30 '20 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.