5
$\begingroup$

I am working on the SIR model and I am trying to plot the incidence function on a specific time interval.

(*Given parameters*)\[Gamma] = 0.1;
\[Beta] = 0.41;
tMax = 90;
n = 10^6;

(*Define the force of infection (\[Lambda])*)
\[Lambda][t] := \[Beta] i[t]/n;

(*Define the system of differential equations*)
eqns = {s'[t] == -\[Lambda][t] s[t], 
   i'[t] == \[Lambda][t] s[t] - \[Gamma]*i[t], r'[t] == \[Gamma]*i[t],
    s[0] == 999999, i[0] == 1, r[0] == 0};

(*Solve the system*)
sol = NDSolve[eqns, {s, i, r}, {t, 0, tMax}];

(*Plot the solutions*)
Plot[Evaluate[{s[t], i[t], r[t]} /. sol], {t, 0, tMax}, 
 PlotLegends -> {"Susceptible", "Infected", "Recovered"}, 
 AxesLabel -> {"Time", "Population"}, 
 PlotStyle -> {Blue, Red, Green}]

enter image description here So far so simple. But when I am trying to define and plot the incidence equation

(*Define the incidence function*)
incidence[t_] = Integrate[\[Lambda][t]*s[t], {t, 40, 50}];

(*Plot the incidence*)
Plot[incidence[t], {t, 0, tMax}, AxesLabel -> {"Time", "Incidence"}, 
 PlotStyle -> Blue, PlotLabel -> "Incidence between t=40 and t=50"]

enter image description here

Could someone explain to me why I have the errors Integrate::ilim and NIntegrate::itraw and no plot?

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ Something is odd in your definition of incidence: how can incidence be a function of t, if t is also the integration variable and it spans a pre-determined range? Should perhaps the variable be a tmax or similar and the integration range be {t, 40, tmax} or something along those lines? Also, shouldn't you inject the results of NDSolve in the definition of incidence as well? $\endgroup$
    – MarcoB
    Aug 21 at 16:56
  • $\begingroup$ @MarcoB Oh! I got it.Τhanks a lot for the boost $\endgroup$
    – Athanasios Paraskevopoulos
    Aug 21 at 16:58
  • 1
    $\begingroup$ incidence[t_] is not actually a function of t in your definition since it's just a dummy variable under Integrate. $\endgroup$
    – Sjoerd Smit
    Aug 22 at 13:07

3 Answers 3

Reset to default
6
$\begingroup$

See if this works:

Clear[incidence]
incidence[tmax_] := NIntegrate[(λ[t]*s[t]) /. sol, {t, 40, tmax}]

Plot[incidence[t], {t, 40, 50}]

enter image description here

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ I wonder why NIntegrate[(λ[t]*s[t]) /. (First@sol), {t, 40, tmax}] looks different from NIntegrate[(λ[t]*s[t]) /. sol, {t, 40, tmax}] even though all First is doing is removing an extra curly around sol (there's only one element in sol)? $\endgroup$
    – ydd
    Aug 21 at 17:06
  • 1
    $\begingroup$ @ydd You need to use First@sol to get a correct answer, but someone else will have to explain why a list even seems to work at all. $\endgroup$
    – Alan
    Aug 21 at 23:00
  • 1
    $\begingroup$ And why if we replace t to another letter such u, then there are no such problem. Clear[incidence1, incidence2]; incidence1[tmax_] := NIntegrate[(λ[t]*s[t]) /. sol, {t, 40, tmax}]; incidence2[tmax_] := NIntegrate[(λ[t]*s[t]) /. sol[[1]], {t, 40, tmax}]; Plot[{incidence1[u], incidence2[u]}, {u, 40, 90}, PlotStyle -> {Automatic, Directive@{AbsoluteThickness[10], Opacity[.2], Automatic}}] $\endgroup$
    – cvgmt
    Aug 22 at 13:37
  • 2
    $\begingroup$ @cvgmt Ah, I thought the "looks different" referred to the extra braces on the output of incidence[] output. The reason is in part because incidence[] depends on a global parameter t, whose evaluation leaks sometimes. Better code: incidence[tmax_?NumericQ] := Block[{t}, NIntegrate[(\[Lambda][t]*s[t]) /. sol, {t, 40, tmax}]] -- But maybe it should be called a bug? $\endgroup$
    – Michael E2
    Aug 22 at 14:03
  • 1
    $\begingroup$ @Alan If you're asking why NIntegrate[{f1, f2,...},...] works, it's because it integrates each component separately. See docs, Scope > Basic Uses > "Vector- and tensor-valued functions". (The separately is not documented but is discussed somewhere on this site; but I can't find the Q&As, sorry.) $\endgroup$
    – Michael E2
    Aug 22 at 14:17
5
$\begingroup$
  • Add an equation incident'[t] == λ[t]*s[t], incident[0] == 0.
Clear["Global`*"];
γ = 0.1;
β = 0.41;
tMax = 90;
n = 10^6;
λ[t] := β i[t]/n;
eqns = {s'[t] == -λ[t] s[t], 
   i'[t] == λ[t] s[t] - γ*i[t], r'[t] == γ*i[t],
    incident'[t] == λ[t]*s[t], s[0] == 999999, i[0] == 1, 
   r[0] == 0, incident[0] == 0};
sol = NDSolve[eqns, {s, i, r, incident}, {t, 0, tMax}];
Plot[Evaluate[{s[t], i[t], r[t], incident[t]} /. sol[[1]]], {t, 0, tMax}, 
 PlotLegends -> {"Susceptible", "Infected", "Recovered", "Incidence"},
  AxesLabel -> {"Time", "Population"}, 
 PlotStyle -> {Blue, Red, Green, Directive@{Black, Dashed}}]

enter image description here

  • And calculus on specific time interval like 40<=t<=50
incident[50] - incident[40] /. sol[[1]]

569781.

  • When we fixed the length of the interval to 10 and vary the endding point b of the inverval,the plot is
Plot[incident[b] - incident[b-10] /. sol[[1]], {b, 10, 90}]

enter image description here

  • When we fixed the beginning point and endding point of the interval, the plot is
Plot[incident[t] /. sol[[1]], {t, 40, 50}]

enter image description here

Improve this answer
$\endgroup$
5
  • $\begingroup$ Plot[incident[t + 10] - incident[t] /. sol[[1]], {t, 0, 90}] $\endgroup$
    – cvgmt
    Aug 22 at 0:17
  • $\begingroup$ Plot[incident[t] /. sol[[1]], {t, 40, 50}]. $\endgroup$
    – cvgmt
    Aug 22 at 0:26
  • $\begingroup$ Could you explain to me please why the first graph" your answer" and the second graph "first comment" are different? I am a bit confused $\endgroup$
    – Athanasios Paraskevopoulos
    Aug 22 at 0:45
  • $\begingroup$ Understood that! I am confused because when I try my initial code + the following code I have different results inc[tmax_] := NIntegrate[{\[Lambda][t]*s[t]} /. sol, {t, 0, tmax}] and Plot[inc[t], {t, 40, 50}, AxesLabel -> {"Time", "Incidence"}, PlotStyle -> Orange] and inc[50] - inc[40] /. sol $\endgroup$
    – Athanasios Paraskevopoulos
    Aug 22 at 1:33
  • $\begingroup$ I think you shifted the incidence plot by 10. That is, incidence at t should be observed incidence over an interval, not projected incidence. See my earlier answer. $\endgroup$
    – Alan
    Aug 22 at 18:40
3
$\begingroup$

I suspect (?) you are after something more like

Clear[incidence]
incidence[tmax_, dt_] := 
 NIntegrate[(\[Lambda][t]*s[t]) /. First@sol, {t, tmax - dt, tmax}]
Plot[incidence[t, 10], {t, 10, 90}]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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