Plot the graph of an integral

Question

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.

Answer 1

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

Answer 2

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

Answer 3

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[
    {
    (* your original equations *)
    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]]