# Plotting Incidence function of the SIR Model

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}]


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"]



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

• 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? Aug 21 at 16:56
• @MarcoB Oh! I got it.Τhanks a lot for the boost Aug 21 at 16:58
• incidence[t_] is not actually a function of t in your definition since it's just a dummy variable under Integrate. Aug 22 at 13:07

See if this works:

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

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


• 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)?
– ydd
Aug 21 at 17:06
• @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.
– Alan
Aug 21 at 23:00
• 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}}] Aug 22 at 13:37
• @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? Aug 22 at 14:03
• @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.) Aug 22 at 14:17
• 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}}]


• 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}]


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


• Plot[incident[t + 10] - incident[t] /. sol[[1]], {t, 0, 90}] Aug 22 at 0:17
• Plot[incident[t] /. sol[[1]], {t, 40, 50}]. Aug 22 at 0:26
• 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 Aug 22 at 0:45
• 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 Aug 22 at 1:33
• 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.
– Alan
Aug 22 at 18:40

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}]
`