Numerical evaluation of an elliptic integral in python

Question

Goal:

I need to evaluate numerically an integral of the following form:

$$ \int_0^\infty \frac{dx}{(a^2+x)\sqrt{(a^2+x)(b^2+x)(c^2+x)}} $$ where $a,b,c \in \mathbb{R}$ are in the interval $(1,1000)$.

This needs to be done in Python.

Problem:

The integration algorithm does not converge, but I know the integral is convergent. To my knowledge there is only one package available that handles improper integrals: scipy.integrate.quad. This method does not converge for the following values: $a=240,b=110,c=100.$ I receive the warning ``IntegrationWarning: The integral is probably divergent, or slowly convergent.'' The integral evaluates to $1.36 \times 10^{-9}$. Computing the same integral in Matlab and Mathematica results in a value of $1.13 \times 10^{-7}$ in each case, with no warnings or errors. This strongly suggests that the integral converges just fine, and that the problem lies in scipy's quadrature routine.

Attempts at a Solution:

The paper from which I took this integral indicates that it is elliptic. There exist several methods to integrate such functions numerically; however, I cannot find any standard elliptic integrals of this form (checking, for example, the discussion on mathworld and also posts such as this one on these forums).

Others have had problems with this package also; for example, this post. The accepted solution in this case is to convert the the integral into an ODE. I have seen, however, that this can be quite dangerous.

I have also considered a change of variables in order to transform this into an elliptic form; I cannot find any. For example, an obvious one, $x \to y^2$, results in an undesirable factor of $y$ in the numerator.

What I think I need:

I believe I will need to implement some sort of quadrature routine manually. Are there any that might be well-suited to this (seemingly benign) integral?

Answer 1

The problem is almost definitely with how QUADPACK (which is the backend used by scipy.integrate.quad) handles numerically small integrands. Essentially the integrand is so small (at $x=0$ it is $6.58\times 10^{-12}$), that it is comparable to the absolute error tolerance. I'm not sure quite why it says it's "probably divergent", but that doesn't matter. I suspect that when the error estimate is comparable with the integral's estimated value, it concludes that the most common explanation is that the integral diverges.

There are several things to try. First is to try and specify the parameter epsabs=0 (by default it is something like $10^{-8}$) to tell it to ignore absolute error completely, and use only relative error. This seems to fix the problem, and is by far the best approach. By comparing the result with the output from mpmath, the result of quad with all other default settings is correct to 13 digits.

The second is to rewrite the integral as $$ \frac{1}{a^3bc} \int \frac{1}{(1+x/a^2)\sqrt{(1+x/a^2)(1+x/b^2)(1+x/c^2)}}. $$ Then the integrand is $O(1)$, which doesn't confuse the integration routine. This also fixes the problem, giving the same result as with epsabs=0.

The third option is to unleash the full machinery of elliptic integrals on this. It is, as far as I know, somewhat tricky to reduce a general-form elliptic integral to the standard elliptic integrals. It can certainly be done, though.

The integral is, in the notation of DLMF 19 $$ \frac{2}{3} R_D(c^2, b^2, a^2). $$ The function $R_D(x, y, z)$ satisfies the recurrence relation DLMF 19.26.E20 $$ R_D(x, y, z) = 2R_D(x+\lambda, y+\lambda, z+\lambda) + \frac{3}{\sqrt{z}(z+\lambda)} $$ $$ \lambda = \sqrt{x}\sqrt{y} + \sqrt{y}\sqrt{z} + \sqrt{z}\sqrt{x}, $$ and in the limit $x,y,z\to+\infty$, $R_D(x,y,z)$ approaches 0. It is sufficient to apply this duplication formula about $30$ times, and take the single remaining elliptic integral $R_D$ with very large arguments as approximately zero, for your values of $x=c^2$, $y=b^2$, $z=a^2$, before it converges to machine precision to $$ 1.1322299840381455\times10^{-7}. $$ This will work more slowly for smaller values of $a,b,c$.