# Numerically integrate a plotted function

I used Plot[NIntegrate[...]...] to plot a function of 5 different variables. It took really long. Right now I need to integrate this function one more time over the 5th variable and plot the result. Can someone tell me if there is an option to Nintegrate the plot to speed things up? Maybe there's a way to convert the plot to a list of $x$ and $y$ and then take a sum over $(y2-y1)\,(x2-x1)$ and list plot it or something similar?

I guess it's not gonna work for me. I have the following function:

f[ρ_, ϕ_, ζ_, ξ_, r_] := 2.*100.*20.*0.02*0.2*ρ*Cos[ϕ]*Sech[20*(ρ - 1.)]^2.*ζ*Exp[-ζ]*
BesselJ[0, 100.*ζ*Sqrt[r^2. + ρ^2. - 2.*r*ρ*Cos[ϕ]]]*(Sqrt[ζ^2. + 0.02^2.]*
Cosh[10.*Sqrt[ζ^2. + 0.02^2.]*(ξ + 1.)] + ζ*Sinh[10.*Sqrt[ζ^2. + 0.02^2.]
*(ξ + 1.)])/((2.*ζ^2. + 0.02^2.)*Sinh[10.*Sqrt[ζ^2. + 0.02^2.]] +
2.*ζ*Sqrt[ζ^2. + 0.02^2.]*Cosh[10.*Sqrt[ζ^2. + 0.02^2.]])


Here the function $f$ have 5 variables: $\rho,\ \phi,\ \zeta,\ \xi,\ r.$. First, I need to integrate over $\rho,\ \phi,\ \zeta,\ \xi.$ and plot $f=f(r).$ It takes about 1 hour. I use as I said before

Plot[NIntegrate[f, {ξ, -1., 0.}, {ζ, 0., ∞}, {ϕ, -3.1415, +3.1415}, {ρ, 0., ∞}], {r, 0, 10}]


The next step is I need to integrate one more time over $r$, {r,0,r1}and plot it over {r1,0,7}. I tried what you've said, but I failed. Is your method with NDsolve going to work with this kind of function?

• See the accepted answer here. – Jason B. Jan 21 '16 at 15:32
• It seems I misunderstood your question. Can you clarify why you want to tie together the plotting w.r.t r and the integration w.r.t r? Is your motivation performance? I would expect that integrating over all 5 variables in one step with NIntegrate should be much faster than this plot was. You can also speed up plotting by limiting the number of plot points. Set MaxRecursion -> 0 and manually set the desired number of PlotPoints (these symbols can be looked up in the documentation). – Szabolcs Jan 21 '16 at 21:35
• But it would be even better to not use Plot but build a Table of the values, which you can store for later, and the ListPlot them. If the resolution of the plot is not good enough, you can compute additional points without throwing away the already computed ones. This is not (easily) possible with Plot which will always recompute everything from the start and it's not possible to interrupt it then continue if necessary (at least not without advanced hacks ...) – Szabolcs Jan 21 '16 at 21:37

The trick is to use NDSolve instead of NIntegrate and thus in effect obtain a numerical antiderivative that can be evaluated fast at different points. NIntegrate will only do definite integration, so it needs to be run each time the integration bounds are changed. This is very slow, as you noticed. NDSolve will only need to be run once.

### Slow way (what you used)

We are going to integrate f:

f[x_] := Sin[x^2]

F[y_?NumericQ] := NIntegrate[f[x], {x, 0, y}]

Plot[F[x], {x, 0, 5}] // AbsoluteTiming


It took 5 seconds on my machine.

NIntegrate[F[x], {x, 0, 5}] // AbsoluteTiming
(* {0.686052, 2.63519} *)


Both integration and plotting are slow with this method. For complex functions, they can be prohibitively slow.

### Fast way (NDSolve)

iF = NDSolveValue[{F2'[x] == f[x], F2 == 0}, F2, {x, 0, 5}]


iF is an interpolating function.

Plot[iF[x], {x, 0, 5}] // AbsoluteTiming


Here's a not so widely known trick: the antiderivative of an interpolating function can be computed exactly using Integrate (not NIntegrate). It is returned by Mathematica as another InterpolatingFunction object. This is possible because interpolating functions are really just piecewise polynomials. Derivatives can be taken too.

So we can compute the second integral directly (and very quickly) as

Integrate[iF[x], {x, 0, 5}]

(* 2.63519 *)


Another trick: To get the antiderivative as a function object directly, you can also do

Derivative[-1][iF]


The disadvantage of the fast way is that there's less oversight about precision loss and numerical errors. To make sure that the final integration result is accurate enough (despite error accumulation), you must make sure that NDSolve computes iF with sufficient accuracy.

• "Here's a not so widely known trick": in some sense, this trick should be obvious, but it clearly isn't. +1 for showing the tricks! – march Jan 21 '16 at 17:33
• I added some details. Could you check it out? – LexRomah Jan 21 '16 at 19:10
• @LexRomah Right, this method won't work for you. :-( – Szabolcs Jan 21 '16 at 21:37
• sorry, but doesn't your method work for a function of several variables? – illuminates Dec 3 '17 at 16:41