How to speed up the plot of NIntegrate?

Question

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] == 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?

Answer 1

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] == 0}, y[x], {x, 0, 10}]
];
t = AbsoluteTime[];
Plot[y[t], {t, 0, 10}]
AbsoluteTime[] - t

Answer 2

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} = 
  NIntegrate`GaussKronrodRuleData[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, #][[1]]} &, 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]]

Advanced example

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} *)

Animation showing adaptivity

Following a misinterpreted suggestion by J.M. I made this animation that shows the effect the options of PPartialIntegrations.

Answer 3

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.

Answer 4

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

Answer 5

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.