Computing highly oscillating integrals using Python
Recently during my research I stumbled upon the following indefinite integral that needed to be computed for various values of $\omega \in (0, \infty)$. Well, maybe not precisely this integral but a really similar one. Anyway here it goes:
Naively, I just threw it into
scipy.integrate.quad routine. For small values of $\omega$ it works fine, but as $\omega$ gets larger we run into troubles - run the below code and see for yourself.
You may see the integrand for various values of $\omega$ in below plot.
Clearly this integral is convergent so why does numerical integrator fail? The reason lies in the oscillating nature of the integrand. It is often the case that highly oscillating integrals require special treatment and this is what we are going to discuss today.
Taming oscillating integrands using Python
Solution using scipy
quad can deal with our integral - we just need to feed it right parameters. The key to the solution is treating our integral as a weighted one. It turns out that if your integrand is of the form $f(x)g(\omega x)$, where $g$ is either $sin$ or $cos$,
scipy (or rather: quadpack that is used under the hood) can compute it using a quadrature specifically taiored to deal with such problems.
So, computing the above integral can be done in scipy as follows
mpmath package has a function
quadosc for the sole purpose of computing integrals of oscillating functions (docs here). The usage is simple, you need to pass
quadosc a function to integrate, integration limits and either the angular frequency (the $\omega$ in our function) or period of oscillatory term.
Which one to choose?
If we rrun both of the above examples we will see that they give similar results. So which approach to choose? There are two factors to consider: speed and accuracy.
To analyse speed of
mpmath.quadosc I computed our integral for multiple values of $\omega$ between 1 and 100, 100 times for each value using both methods. The below plot shows mean time to solution, y-axis in log scale.
As we may see using
quad is almost three orders of magnitude faster than
scipy is faster here. But what about accuracy? It is hard to measure it for our integrand, as we don’t know the exact solution. Let’s compute another integral, one which we can solve analytically, namely:
We can easily compute the above integral explicitly (nice exercise if you haven’t computed any integrals for some time!). The exact solution is $F(\omega) = \omega / (\omega^2 + 1)$. I computed the $F(\omega)$ numerically for various values of $\omega$ between 1 and 100. The below plot shows errors of numerical results with respect to the exact solution.
quadosc gives slightly smaller error, the differences are neglible, relative errors as well as absolute ones are small for boths methods.
It seems that
quad is a clear winner, offering comparable accuracy of solution in much better time.
What about mpmath.quad?
Documentation states that for exponentially or faster decaying integrals
mpmath.quad should work equally fine as
mpmath.quadosc and should be much faster. I tried it for $F(\omega)$ and it performed much worse than other solutions presented above. Below plot shows error of $F(\omega)$ values computed using
mpmath.quad with respect to exact solution.
As can be seen,
mpmath.quad doesn’t work well when oscillation frequency increases. For large values of $\omega$ the relative error reaches as high as 800%.
We presented two methods for integrating higly oscillating functions using Python scientific libraries. Both
scipy.integrate.quad gave accurate results for the tested function, however
quadosc in terms of timing.
In interpreting results of our small benchmark we should note that
mpmath.quadosc offers more flexibility when defining oscillatory terms. In particular, it offers the following possibilities that were not exploited in the benchmark:
- computing integral with oscillating term different than $\cos$ and $\sin$
- computing integrals wich uneven oscillations (for instance with oscillatoring term $\sin(x^2)$)
Keeping that in mind,
mpmath.quadosc should still be regarded as an useful tool, but for the problem discussed in this post it falls behind
The final thing to consider is that I run all of the benchmarks without tweaking any of the additional parameters of the used integrators. It is possible that better results, either in terms of accuracy or timing, can be achieved if such fine tuning is performed. The reason I didn’t dwelve into this too much is primarily because I didn’t have time to perform extensive studies. Secondly, I think that a tool is most valuable if you can use it without too much tweaking. Therefore it is good to see what numerical integrators can do with their default settings.