How to speed up the plot of NIntegrate?

Here is a toy example:

f[t_] := NIntegrate[Sin[x], {x, 0, t}];
Plot[f[t], {t, 0, 10}] // Timing


Even such a simple example will take 2.8 seconds on my computer.

Since many of the plot family functions have the attribute HoldAll, Mathematica evaluates f[t] only after assigning specific numerical values to t, thus causing a lot of repeated evaluations of NIntegrate. The integration to a smaller upper-limit is forgotten when integrating to a bigger upper-limit. On the other hand, I can’t benefit from wrapping the function in Evaluate, because of the numerical nature of NIntegrate. Something stuck here. So the questions are:

• Can NIntegrate remember or make full use of the result of a smaller upper-limit integral? Or generally, how to speed up the plot involving NIntegrate or is there any principle to do it?

• Is it possible to realize ParametricNIntegrate just as ParametricNDSolve in version 9?

Edit

Following the way of changing the integral to NDSolve (Thanks to Mark McClure), I've come up with one possible way to realize ParametricNIntegrate. Here is the code for a more general way compared to the toy model:

sol = ParametricNDSolveValue[{f'[x] == Sin[a x], f == 0}, f, {x, 0, 10}, {a}];
Manipulate[Plot[{sol[a][t], 1/a (1 - Cos[a t])}, {t, 0, 10}], {a, 1, 5}]


Now I can change the parameter and plot the integral in real-time. Any better ideas?

You can express your integral in terms of a differential equation and use NDSolve. Since NDSolve builds up the solution as it goes, this is typically much faster.

Clear[y];
y[x_] = y[x] /. First[
NDSolve[{y'[x] == Sin[x], y == 0}, y[x], {x, 0, 10}]
];
t = AbsoluteTime[];
Plot[y[t], {t, 0, 10}]
AbsoluteTime[] - t • Thank you! That's a very smart way:) I wonder if it will also work when the integrand is complicated. – luyuwuli Apr 23 '13 at 12:37
• @luyuwuli It should work whenever NDSolve works. My guess is that it is not quite a broadly applicable as the integrate technique, but it should be nearly so. – Mark McClure Apr 23 '13 at 12:41
• I just came up with the idea that if the function has parameters in the integrand(no just the upper-limit),following the same way, ParametricNDSolve will do the trick:) I think this may be called ParametricNIntegrate. Thank you again! – luyuwuli Apr 23 '13 at 12:55
• Since NDSolve[] is also more cautious than NIntegrate[], this approach is also less likely to miss features of your integrand, like sharp peaks. – J. M. is away Apr 23 '13 at 13:54
• @J.M. Why is NDSolve more cautious than NIntegrate? IMHO, solving the integral is easier than the PDEs. – luyuwuli Apr 23 '13 at 14:11

An alternative approach is to form an approximate interpolating function from your actual function, which is (often) cheaper to integrate. To wit,

nsin = FunctionInterpolation[Sin[x], {x, 0, 10}];
ni = Derivative[-1][nsin];

Plot[ni[t], {t, 0, 10}] In this case, we know that the integral of sine can be expressed simply, so we can compare the exact function with the approximate one:

Plot[ni[t] - 2 Haversine[t], {t, 0, 10}, Frame -> True, PlotRange -> All] Pretty darn good, if I do say so myself! If you need a bit more accuracy, FunctionInterpolation[] has a number of options you can tweak as needed; see the docs for details.

• Thanks. I've added a solution to do something like ParametricNIntegrate. I'm just not sure whether it's good. – luyuwuli Apr 24 '13 at 3:44

Can NIntegrate remember or make full use of the result of a smaller upper-limit integral? Or generally, how to speed up the plot involving NIntegrate or is there any principle to do it?

Below is a solution in the spirit of this request.

In short, we find adaptively sampled points, compute integral estimates over the intervals having those points as boundaries, compute cumulative sums, and list plot.

General procedure

First, note that NIntegrate is similar to Plot -- both use sampling points adaptively to produce their output. (E.g. see Max(Plot)Points and MaxRecursion.)

1. Use NIntegrate over the whole range the integral values to be plotted within. Gather the sampling points.

2. Partitition the sampling points into disjoint or overlapping intervals. Sort them accordingly to their boundaries.

3. For each interval compute an integral estimate using an integration rule. (Preserving the order of the intervals.)

4. Compute the cumulative sums of the integral estimates and pair them with the interval boundaries.

5. (List)Plot the obtained pairs.

Note that because of step 1 we will have a reasonably adaptive behavior.

How to obtain integration rules and integral estimates with them is shown in NIntegrate's advanced documentation. (It is also given below.)

Implementation

This implementation also allows sampling points specification instead of using NIntegrate's sampling points in step 1.

{absc, weights, errweights} =
NIntegrateGaussKronrodRuleData[5, MachinePrecision];

IRuleEstimate[f_, {a_, b_}] :=
Module[{integral,
error}, {integral,
error} = (b - a) Total@
MapThread[{f[#1] #2, f[#1] #3} &, {Rescale[absc, {0, 1}, {a, b}],
weights, errweights}];
{integral, Abs[error]}]

Clear[PPartialIntegrations]
Options[PPartialIntegrations] = {MinRecursion -> 1, PrecisionGoal -> 3, "PointsPerInterval" -> 2};
PPartialIntegrations[f_, {a_?NumericQ, b_?NumericQ},
opts : OptionsPattern[]] :=
Block[{minrec, pgoal, npint, res, x},
minrec = OptionValue[MinRecursion];
pgoal = OptionValue[PrecisionGoal];
If[TrueQ[pgoal === Automatic], pgoal = 3];
npint = OptionValue["PointsPerInterval"];
If[npint === Automatic, npint = 2];
res = Reap@
NIntegrate[f[x], {x, a, b}, MinRecursion -> minrec,
PrecisionGoal -> pgoal, EvaluationMonitor :> Sow[x]];
PPartialIntegrations[f, {a, b}, Sort[res[[2, 1]]][[1 ;; -1 ;; npint - 1]]]
];

PPartialIntegrations[f_, {a_?NumericQ, b_?NumericQ}, stepArg_] :=
Block[{step = stepArg},
If[step === Automatic, step = (b - a)/200];
PPartialIntegrations[f, {a, b}, Range[a, b, step]]
] /; AtomQ[stepArg];

PPartialIntegrations[f_, {a_?NumericQ, b_?NumericQ},
range : {_?NumericQ ..}] :=
Block[{intervals, integrals, cumSums},
res = Map[{#, IRuleEstimate[f, #][]} &, Partition[range, 2, 1]];
intervals = res[[All, 1]];
integrals = res[[All, 2]];
cumSums = Prepend[Accumulate[integrals], 0];
Transpose[{Prepend[intervals[[All, 2]], intervals[[1, 1]]], cumSums}]
];


Basic examples

We can see that the implementation is reasonably fast.

AbsoluteTiming[PPartialIntegrations[Sin, {0, 10}, Automatic];]
(* {0.009395, Null} *)

AbsoluteTiming[
PPartialIntegrations[Sin, {0, 10}, MinRecursion -> 2, PrecisionGoal -> 5];]
(* {0.006568, Null} *)

ListLinePlot[PPartialIntegrations[Sin, {0, 10}, Automatic]] Consider the function:

Plot[Sin[x] + 2 Exp[-100 (x - 5)^2], {x, 0, 10}, PlotRange -> All] We can see that we get reasonably good list plot using the options MinRecursion and PrecisionGoal given to NIntegrate in step 1. (The grid lines show the sampling points.)

AbsoluteTiming[
pres = PPartialIntegrations[Sin[#] + 2 Exp[-100 (# - 5)^2] &, {0, 10},
MinRecursion -> 1, PrecisionGoal -> 10];
]
ListLinePlot[pres, ImageSize -> Large, PlotRange -> All,
GridLines -> {pres[[All, 1]], None},
GridLinesStyle -> {Directive[GrayLevel[0.9]], None}]

(* {0.062035, Null} *) Following a misinterpreted suggestion by J.M. I made this animation that shows the effect the options of PPartialIntegrations. • This is nice! I guess a logical extension, at least for univariate integrals, would be an animation showing how NIntegrate[] samples and splits its integrand, at each recursion level. – J. M. is away Jun 11 '16 at 4:17
• @J.M. I will upload one later today ... – Anton Antonov Jun 11 '16 at 15:38
• Oh my; I wasn't expecting you to go ahead and actually implement it! That's way past nice. :) – J. M. is away Jun 11 '16 at 15:46
• @J.M. It seems I had in mind different thing. I think the NIntegrate Explorer provides something very close to what you suggested. (It was definitely designed and programmed so people can get better idea of how numerical integration algorithms operate.) – Anton Antonov Jun 11 '16 at 18:38
• @J.M. I posted the animation I made anyway... – Anton Antonov Jun 12 '16 at 13:30

This answer is just a response to the following part of the question:

……generally, how to speed up the plot involving NIntegrate or is there any principle to do it?

When the symbolic processing of NIntegrate isn't necessary for getting the correct result of integration, "SymbolicProcessing" -> 0 is also a general way for speeding up:

f[t_] := NIntegrate[Sin[x], {x, 0, t}, Method -> {Automatic, "SymbolicProcessing" -> 0}];
Plot[f[t], {t, 0, 10}] // AbsoluteTiming • Just for record,I found "SymbolicProcessing" -> 0 is faster than "SymbolicProcessing" -> False. :) – yode Jun 11 '16 at 2:48
• Good point -- I am surprised nobody else noted that! – Anton Antonov Jun 11 '16 at 20:48

Another possibility is to sieze control over the evaluation of f[t]. For instance:

ListLinePlot[Table[f[t], {t, 0, 10, 0.1}]]


evaluates f[t] at exactly the points you choose. In this case, it gives substantially the same plot as the original, and is about 10 times faster.

• This is a discrete version. The drawback is that I have to manually decided how many and which points to calculate. For a complicated function, it would very difficult to do so. And in the aspect of saving calculation, I think that this method will repeat the integral 101 times. If the number of selected points is as much as the one in Plot`, this will also drag down the speed. – luyuwuli Apr 23 '13 at 13:52