Plot the graph of an integral

If I have a system of differential equations, I solve it and save the result in Sol:

Clear[a, b, S0, I0];

a = .0015; b = .9; S0 = 1000; I0 = 10; tmax = 20;

t0 = 0;

Sol = NDSolve[{S'[t] == -a*II[t]*S[t],
II'[t] == a*II[t]*S[t] - b*II[t], S[0] == S0, II[0] == I0}, {S[t],
II[t]}, {t, 0, tmax}]


For II[t], the graph is correct, which demonstrates that the result is also correct.

Then, I want to plot the graph of $$\int_{t0}^{t} II[t] dt$$. The domain of the graph should be [-10, 60].

Plot[NIntegrate[Evaluate[II[t] /. Sol], {t, t0, tmax}], {t, 0, tmax},
AspectRatio -> .5, AxesLabel -> {"t", "I(t)"}, PlotRange -> All,
AxesOrigin -> {0, 0}]


The graph looks like this.

It is obviously wrong.

I was wondering what the problem is here, and how I should perform the integration.

• What does "the graph seems to be wrong" mean. It doesn't display the behavior you expect? It doesn't plot anything? It spits out some errors? Please be more specific. Also (and this might be the problem), the syntax for your NDSolve is wrong. Please look up the correct syntax in the documentation for NDSolve to make that work. Also, you need to specify all numerical parameters. What are the values of a, b, etc.? Mar 20 at 21:13
• Thanks for getting back! I edited my question. Hopefully, it clears up some confusion. The problem here doesn't seem to be with NDSolve, since I get the correct graphs for S[t] and R[t], etc. I think the problem lies in the syntax of integration. Mar 20 at 22:00
• Without your full (working!) code, we can't test out what's going wrong. I appreciate the update for how the graph is wrong, but working NDSolve code with the values for the parameters is necessary so that we can run the code on our own copies of Mathematica and tinker with it to figure out a solution to your problem. Mar 20 at 22:15
• That said: does this work? Plot[NIntegrate[Evaluate[II[tp] /. Sol], {tp, t0, t}], {t, 0, tmax}, AspectRatio -> .5, AxesLabel -> {"t", "I(t)"}, PlotRange -> All, AxesOrigin -> {0, 0}] Mar 20 at 22:16
• Try ListLinePlot[Head@Integrate[II[t] /. First@Sol, t], AspectRatio -> .5, AxesLabel -> {"t", "I(t)"}, PlotRange -> All, AxesOrigin -> {0, 0}] Mar 21 at 3:22

In this bit of code:

NIntegrate[Evaluate[II[t] /. Sol], {t, t0, tmax}]


you really want to integrate to t, not to tmax, because otherwise it's just a number. To do this, we instead use the code

Integrate[II[t] /. Sol /. t -> tp, {tp, t0, t}]


which is equivalent to $$\int_{t_0}^{t}II(t_p)dt_p$$ (note that the integration variable is different than the $$t$$ in the limits!). Because Sols is defined as II[t] -> stuff, we need to replace II[t] first, but then replace t with tp in order to do the integration (this is not strictly necessary, because Integrate scopes its integration variable, but nevermind, and anyway it's a good idea to make the Mathematica code reflect the mathematics).

(Note: I replaced NIntegrate with Integrate because InterpolatingFuntions can be integrated symbolically, since they're just made of polynomials under the hood. This isn't necessary, but it can be a nice trick elsewhere.)

With those changes, we get the following:

Clear[a, b, S0, I0];
a = .0015; b = .9; S0 = 1000; I0 = 10; tmax = 20;
t0 = 0;
Sol = NDSolve[{S'[t] == -a*II[t]*S[t], II'[t] == a*II[t]*S[t] - b*II[t], S[0] == S0, II[0] == I0},
{S[t], II[t]}, {t, 0, tmax}];

Plot[
Evaluate@Integrate[II[t] /. Sol /. t -> tp // Evaluate, {tp, 0, t}],
{t, 0, tmax},
AspectRatio -> .5,
AxesLabel -> {"t", "I(t)"},
PlotRange -> All,
AxesOrigin -> {0, 0}
]


yielding

• I truly appreciate the specific explanation, and I now understand that it should be {tp, t0, t} instead of {t, t0, tmax}. However, I think that the graph should not have the same shape as the original one, and it should keep increasing instead of decreasing, since it accumulates the former values. In my opinion, it should resemble more of a logistic curve. The value of II[t] decreases after t = 7, so the slope of the integration should decrease instead of the value. Mar 20 at 23:17
• @yixjia You're right! A simple matter of getting the order of evaluation wrong. I've fixed the issue with an Evaluate. Mar 21 at 3:12
Clear["Global*"];

a = 15/10000; b = 9/10; S0 = 1000; I0 = 10; tmax = 20;

t0 = 0;

Sol = NDSolve[{S'[t] == -a*II[t]*S[t],
II'[t] == a*II[t]*S[t] - b*II[t], S[0] == S0,
II[0] == I0}, {S[t], II[t]}, {t, t0, tmax}][[1]];

Plot[II[t] /. Sol, {t, t0, tmax}]


intII[tv_?NumericQ] := NIntegrate[II[t] /. Sol, {t, t0, tv}]

Plot[intII[t], {t, 0, tmax},
AspectRatio -> .5,
AxesLabel -> {"t", "I(t)"},
PlotRange -> All,
AxesOrigin -> {0, 0}]


I would just include the integral into your ODE, because then you avoid using numerical integration for each point of the plot. Specifically, if int[t] is your function, then int'[t] == II[t], and the initial condition is int[0] == 0. So, use:

Sol = NDSolveValue[
{
S'[t] == -a II[t] S[t],
II'[t] == a II[t] S[t] - b II[t],
S[0] == S0,
II[0] == I0,

(* equations for integral *)
int'[t] == II[t],
int[0] == 0
},
{S, II, int},
{t, 0, tmax}
];


Then, the third element is the integral:

ListLinePlot @ Sol[[3]]
`

• Definitely the more elegant way to do it! Mar 21 at 17:42