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?


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?


5 Answers 5


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.

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

enter image description here

  • $\begingroup$ Thank you! That's a very smart way:) I wonder if it will also work when the integrand is complicated. $\endgroup$
    – luyuwuli
    Apr 23, 2013 at 12:37
  • $\begingroup$ @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. $\endgroup$ Apr 23, 2013 at 12:41
  • 1
    $\begingroup$ 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! $\endgroup$
    – luyuwuli
    Apr 23, 2013 at 12:55
  • $\begingroup$ Since NDSolve[] is also more cautious than NIntegrate[], this approach is also less likely to miss features of your integrand, like sharp peaks. $\endgroup$ Apr 23, 2013 at 13:54
  • 1
    $\begingroup$ @luyuwuli, your assessment does not contradict mine. It is true that more work is needed to solve a DE than to integrate a function. I talked about this at some length here. $\endgroup$ Apr 23, 2013 at 14:19

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.)


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_}] := 
   error}, {integral, 
    error} = (b - a) Total@
     MapThread[{f[#1] #2, f[#1] #3} &, {Rescale[absc, {0, 1}, {a, b}], 
       weights, errweights}];
  {integral, Abs[error]}]

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

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

ListLinePlot[PPartialIntegrations[Sin, {0, 10}, Automatic]]

enter image description here

Advanced example

Consider the function:

Plot[Sin[x] + 2 Exp[-100 (x - 5)^2], {x, 0, 10}, PlotRange -> All]

enter image description here

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.)

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

enter image description here

Animation showing adaptivity

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

enter image description here

  • $\begingroup$ 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. $\endgroup$ Jun 11, 2016 at 4:17
  • $\begingroup$ @J.M. I will upload one later today ... $\endgroup$ Jun 11, 2016 at 15:38
  • $\begingroup$ Oh my; I wasn't expecting you to go ahead and actually implement it! That's way past nice. :) $\endgroup$ Jun 11, 2016 at 15:46
  • 1
    $\begingroup$ @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.) $\endgroup$ Jun 11, 2016 at 18:38
  • $\begingroup$ @J.M. I posted the animation I made anyway... $\endgroup$ Jun 12, 2016 at 13:30

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

plot of approximate integral

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]

exact versus approximate integral

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.

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

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

enter image description here

  • 4
    $\begingroup$ Just for record,I found "SymbolicProcessing" -> 0 is faster than "SymbolicProcessing" -> False. :) $\endgroup$
    – yode
    Jun 11, 2016 at 2:48
  • 1
    $\begingroup$ Good point -- I am surprised nobody else noted that! $\endgroup$ Jun 11, 2016 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.

  • $\begingroup$ 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. $\endgroup$
    – luyuwuli
    Apr 23, 2013 at 13:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.