I am doing a nested integral via quadrature. To give a definite example, lets say:

$$ \int_0^2dx \left[x + \int_0^x dy \, 2y\right] $$

So effectively I'm integrating $x + x^2$ from 0 to 2 (although obviously my actual problem can't be analytically solved like this). Its implemented in code something like this:

inner(x) = quad(y->2*y, 0, x, tolerance=t_inner)
quad(x + inner(x), 0, 2, tolerance=t_outer)

EDIT: Also note, in this example $x$ appears only as a limit to the inner integral, but my actual problem it appears in the integrand as well. Generally I just mean that the inner quad call depends on $x$ in some way.

The point of this example is that for small $x$, the first term is bigger than the second, whereas for larger $x$, the second term (the inner integration) dominates.

This would dictate that I could speed up my code by setting the inner quadrature error tolerance (t_inner) large when $x$ is small, since error on the inner integral at small $x$ contributes less to my overall error. Then I can increase t_inner when I get to large $x$ when this inner quadrature really matters.

My question is, is there some general prescription for figuring out the optimal way to do this kind of thing? Is there some way I can just choose my "global" tolerance, t_outer, and have t_inner picked automatically as a function of $x$?


If the two terms in your integral have substantially different behavior, the better approach is likely to split the integral into two separate integrals. This way, you can choose parameters for each integrand separately.

  • $\begingroup$ Agree and understand, though a possible point of confusion is that the original expression has up to 3 integrals. $\endgroup$ – innisfree Sep 10 '16 at 12:00
  • $\begingroup$ The double integral and the single integral $\endgroup$ – innisfree Sep 10 '16 at 12:00
  • $\begingroup$ Yes. Then you treat the single integral independently from the double integral. $\endgroup$ – Wolfgang Bangerth Sep 11 '16 at 18:45

I haven't fully thought this through, but it seems like t_inner should be changing by function of the derivative of the inner integral with respect to x. For example, lets say we have a function t_inner = A*(B - d(2y)/dx). You can pick values for A and B by defining boundary conditions at the integral limits, or by setting t_outer = t_inner at the value of x where x = inner(x). Intuitively it seems like you don't need both A and B (either A=1 or B=0) but I'm not sure which is correct.

In any case, as you're iterating through x you can then calculate a new t_inner for each x by looking at the numerical change in inner(x) over the previous iteration. As long your x spacing is fairly smooth with respect to inner(x) it seems like this should work.


Your Answer

