Quick tangent: I took a look at IteratedIntegration.jl. Am I right that it can only be used for regular multiple integrals? That is, it can’t help me if I have nested integrals where the outer integrand is a nonlinear function of the inner integral? Schematically, I’m working with something like this:
\begin{aligned}
I &= \int f(x, g(x))\, dx \, \, ,\\
g(x) &= \int h(x, y)\, dy \, ,
\end{aligned}
where f(x, g(x)) is nonlinear in the second argument, so you can’t pull it inside the inner integral. I would love to have global error estimates, but so far, I haven’t found a better approach than nesting 1D quadratures (and, as mentioned earlier, approximating g(x) with an interpolant because I need to do this for many different f with the same g).
In that case basically you are defining a new special function g(x), and like any special function that you want to evaluate many times it makes sense to define some approximations (generically, you could use Chebyshev polynomials, e.g. via ApproxFun.jl or FastChebInterp.jl … or if you are willing to do more specialized analyses you could use mixtures of approximations, e.g. Taylor + continued-fraction + …, as is done for other special functions like Bessel etc.).
Thank you all for your replies!
While the original MWE that I posted reflects the issue that I am experiencing in my case, it does not faithfully represent what I am trying to do.
In my case I have an oscillatory integral of the type:
\int f(x, y) \sin{( x y)} \ \mathrm{d}x \mathrm{d}y
Moreover, I actually want the integration points over x that the adaptive algorithm gives as an output (so that I am actually calling quadgk_segbuf) rather than the value of the integral. Later on I use them to perform some other computations.
In my approach, I perform the integral
g(x) = \int f(x, y) \sin{( x y)} \ \mathrm{d}y
over y by using a combination of a specialized method for oscillatory integrals and quadgk. The decision on when to use each method is dynamic, making the boundary a complicated curve in the 2D integration region. Then, I run quadgk_segbuf on g to get the integration points.
So, the nested quadgk calls that I use do not account for the whole integral.
Using IteratedIntegration.jl directly on the integrand would probably lead to slow executions due to the oscillatory nature of the integral. Also, since I want to do this only once, interpolating g seems to duplicate work in terms of evaluating it.
I realize that this strategy may sound a bit cumbersome, but I would really like it to work efficiently. I am unsure if there is any other package that I could use, if anyone has suggestions, I would appreciate hearing them.
Actually, I want to find the quadrature points for the original integral (also computing the corresponding weights) to later compute the integrals of a bunch of other functions that are similar to the first one.
I am doing it like that because the functions that I want to integrate are defined via a recurrence relation with parameters depending on some dynamic variables. It turns out that evaluating the functions at a set of points is easy, so I compute their values at the quadrature nodes, and together the previously computed weights, their integrals are fast to compute.
Hi, I’m the author of IteratedIntegration.jl and based on the problem you are solving involving numerous oscillatory integrals it sounds like you would be better off using a fixed 2d quadrature rule, such a product of Gauss-Legendre rules (see FastGaussQuadrature.jl) or Clenshaw-Curtis rules. You would have to check that the number of quadrature points is large enough for convergence and the standard way to do this is p-adaptive quadrature. Cubature.jl from @stevengj includes an explanation and implementation of p-adaptive integration which can be convenient to call using Integrals.jl. If you need a pure-Julia implementation, I started a p-adaptive Clenshaw-Curtis implementation as a test case in AutoSymPTR.jl, which provides a includes a generic p-adaptive algorithm implemented specifically for integrals of periodic functions.