Evaluating an integral numerically at many points

Question

Given a real function $f$, how can one efficiently evaluate $\int_0^{a_i}f(x)dx$ for millions of different $a_i$?

Applying a standard quadrature method (such as Simpson's rule or Gaussian quadrature) will incur an independent cost for each $a_i$, since in general $f(x)$ will need to be evaluated at an entirely different set of points for each $a_i$.

I would prefer to pre-generate/calculate a data structure at substantial initial cost, if the integral can then be evaluated from this data structure cheaply, for different $a_i$.

One possibility I see is to evaluate the integral at a regular array of points $x_j=a_{max}j/n$ (which would allow reuse of points in Simpson's method), and then to "correct" this integral on final evaluation with a small $\int_{a_{max}j/n}^{a_i}f(x)dx$ which can be evaluated with a very low-order Gaussian quadrature because $j$ can be chosen to make this correction interval small. Although probably accurate, it seems like this method may give rise to discontinuities at the half-way points between subsequent $j$.

Another possibility may be to fit cubic splines to a regular array of points (evaluated from Simpson's method) and then rely on these to interpolate the integral.

I would rather follow an established approach if there is one.

Answer 1

The point of @WolfgangBangerth is exactly what I mentionend in my comment, so I'd always try this first. In the best case, with millions of partitions $[a_{i},a_{i+1}]_{i\in \{0,\ldots,N-1\}}$ (where $a_{i} < a_{i+1}$), the length of these is quite small such that some low-order method could be appropriate. This will require roughly $N \times (s-1)$ function evaluations, where $s$ is the number of evaluations required to integrate one subinterval (e.g., $s=2$ for the trapezoidal rule), and $N \approx 10^6$ should be no problem for a modern PC.

Alternatively, in order to map it on a standard problem, you can consider the integral $I(x) = \int_0^x f(t) dt$ in differential form, $$ \frac{\partial I(x)}{\partial x} = f(x), \quad \ I(0) = 0. $$ One can then apply Runge-Kutta or some other standard ODE solver to integrate for $I(x)$ from $x=0$ to $x=a_N$. The steps that are taken in the solution in general won't meet your specified gridpoints $a_i$, thus you need to apply a method called dense output which is basically a polynomial interpolation with an order that is appropriate for your solver. Dense output usually requires no or only a few additional function evaluations per step, and thus can save you a lot of effort.

Let's apply this approach to the idea you mentioned: one can integrate the equation using fourth order Runge-Kutta (which in this case directly corresponds to Simpson's rule) with a constant stepsize $h = a_N / M$. Here, the number of steps $M$ should be chosen such that the final solution $I(a_N)$ is accurate, what should be possible with $M \ll N$ in order to justify the overhead over the simple integration mentioned in the beginning.

Then, in each subinterval $[x_{n},\,x_{n+1}]_{n\in \{0,\ldots,N-1\}}$, where $x_{n}=nh$, $I_n=I(x_n)$, $f_n=f(x_n)$, you can apply Hermite interpolation as [see Numerical recipes, 3rd edition, formula (17.2.15)]: $$I(x_n+\theta h)=(1-\theta)I_n + \theta I_{n+1} + \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;\;\\ \theta(\theta-1)\Big((1-2\theta)(I_{n+1}-I_n)+(\theta-1)h f_n + \theta hf_{n+1}\Big)$$

This is a cubic interpolation polynomial which is appropriate for the fourth-order Runge-Kutta method -- one order less than the local order is ok, as you don't use the interpolant for taking further steps, see e.g. here. All that remains to do is to adjust the parameter $\theta$ such that $a_i = x_n+\theta h$ (you have to do that repeatedly for each $a_i \in [x_n,x_{n+1}]$).

Answer 2

Let's say you've ordered the points so that $0<a_1<a_2<\ldots<a_n$. Then you have of course that $$ \int_0^{a_2} f(x) dx = \int_0^{a_1} f(x) dx + \int_{a_1}^{a_2} f(x) dx. $$ In other words, if you have previously computed the integrals from zero to $a_1$ through some numerical method, then computing the integral from zero to $a_2$ can reuse this information, and all you have to do is now compute the integral from $a_1$ to $a_2$, which you can do independently.

The same applied for all of the other integrals, of course.

Answer 3

The standard method is double-exponential sinc quadrature. The fast convergence requires some regularity assumptions on the integrand, namely that it is the restriction to the real line of a complex differentiable function bounded on a Bernstein ellipse.