Numerically integrate a plotted function

Question

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?

Answer 1

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] == 0}, F2, {x, 0, 5}]

iF is an interpolating function.

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

This took about 0.02 seconds.

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.