Numerical integral with symbolic integral in exponent

Question

Many times in fourier approximation we come across integrals such as

$$\int_0^1 e^{-\gamma\int_0^xu_0(\eta)d\eta}dx$$ where $\gamma$ is a constant and the data for $u_0$ is provided as a discretely sampled vector.

How does one perform integration for such numerical integration?

Answer 1

I would either interpolate the samples of $u_0$ or do a curve-fit, so as to obtain an approximate $u_0(\eta)$ that is easy to integrate analytically. The choice of interpolation or fitting methods depends on what you know about $u_0(\eta)$. For example, if it is a band-limited signal, and the sample frequency is greater than the Nyquist rate, then interpolation with sinc functions would be ideal, though the outer integral would have to be done numerically. If nothing much is known about $u_0$, I would use splines, of order 1 or more, depending on the smoothness of the samples.

For any linear interpolation method, for example, you could write, $$u_0(\eta) = \sum_i a_i \phi_i(\eta),$$ and then compute $$I = \int_0^1 \exp\left(-\gamma \sum_i a_i \psi_i(x)\right)dx,$$ where $$\psi_i(x) = \int_0^x\phi_i(\eta)d\eta.$$ For an order-1 spline, i.e., a piecewise linear approximation, it should be possible to evaluate $I$ using the error function. Otherwise, you would need some quadrature rule, and since you are allowed to chose the quadrature nodes, a rule such as Gauss-Kronrod or Clenshaw-Curtis might be the best choice.

Answer 2

I think the easiest way may simply be an outer loop over the outer integral, while interatively increasing the inner one. So you basically do a midpoint rule for the outer integral, and remember the value of your inner integral from the last step.

pseudocode:

double dx = 1e-12;
double inner_int =0.0;
double outer_int =0.0;

for(double x = 0.0; x<=1.0;x+=dx){
  inner_int += u_0(x)*dx;
  outer_int += exp(gamma*inner_int)*dx;
}

Think of it as doing a sweep from left to right on your axis, remembering the value of the last inner integral. There might be improvements to do this more elegantly (trapezoidal rule or similar). This is just a hunch, so please verify it with some testcase:-)

You have not posted the function $u_0$. Make absolutely sure that there is no analytical expression available for the inner integral. There are huge tables for integral expressions available and you might be lucky to find your particular one. (e.g. in Bronstein/Semandjajew or similar collections)