Numerical Integral of Complete Elliptic Integral

I'm looking to numerically evaluate an integral of the form $$I=\int_0^1rK(r)f(r)\,dr$$ where $f(r)$ is a smooth function known at a set of grid points $\{r_i\}$ with error $O(\Delta r^2)$ ($\Delta r\sim10^{-3}$) and $K$ is the complete elliptic integral of the first kind. A first attempt would be to simply use the trapezoidal rule, but $K(r)$ is singular at $r=1$ so this is poorly conditioned as the grid becomes very fine. One approach I have seen is to use an asymptotic expansion of $K(r)\approx-\log\sqrt{1-r^2}$ for $r\approx 1$ and then taking $f$ to be constant in some small interval $(1-\epsilon,1)$ (chosen to fall on a grid point) so that $$\int_{1-\epsilon}^1rK(r)f(r)\,dr\approx f(1)\int_{1-\epsilon}^1-\frac{r}{2}\log(1-r^2)\,dr=\frac{f(1)}{4}\left[1+\epsilon(2-\epsilon)\log(1-(1-\epsilon)^2)+(1-\epsilon)^2\right]$$ (if I worked out the integral correctly). This approach seems very sensitive to the value of $\epsilon$ chosen, and it further assumes $f$ is constant on this interval which does not preserve the second-order accuracy. Is there a way to approach this integral that does not require an arbitrary choice of $\epsilon$, and further can compute $I$ with $O(\Delta r^2)$ error?

Solution 1:

Basically you're doing everything right, but several improvements can be suggested.

First, the expansion for $K(r)$ near $r = 1$ is $$ K(r) = \log 4 - \frac{1}{2} \log(1-r^2) + \dots $$ The missing $\log 4$ term introduces $O(\Delta r)$ error in your computations.

I suggest splitting $K(r)$ into regular and singular parts as $$ K(r) = -\frac{1}{2} \log(1-r^2) + R(r). $$ Here $R(r)$ is smooth enough to use in numerical quadrature (though high order quadratures would loose orders of convergence due to $R'(1) = -\infty$)

Then do the splitting $$ \int_0^1 r K(r) f(r) dr = \int_0^1 r R(r) f(r) dr - \int_0^1 \frac{r}{2} \log(1-r^2) f(r) dr. $$ Note that the singular part is used for the whole domain $[0, 1]$ and not only for the $[1-\epsilon, 1]$.

The first integral can be evaluated using the trapezoidal rule. Use $R(1) = \log 4$ at the last point and $R(r) = K(r) + \frac{1}{2} \log (1-r^2)$ for the rest.

The second integral needs to be integrated analytically. To do this we need to split the $[0, 1]$ into segments $[r_m, r_{m+1}]$ in which we treat $f(r)$ as a linear function. In other words we consider $f(r)$ as a piecewise-linear function. Using $a = r_m, b = r_{m+1}$ for simplicity we need to evaluate the following integral $$ \int_a^b \frac{r}{2} \log(1 - r^2) \left[ f(a) + \frac{f(b) - f(a)}{b-a} (r - a) \right] dr = \\ = f(a) J_1(a, b) + \frac{f(b) - f(a)}{b-a} J_2(a, b),\\ J_1(a, b) = \int_a^b \frac{r}{2} \log(1 - r^2) dr\\ J_2(a, b) = \int_a^b \frac{r}{2} \log(1 - r^2) (r-a) dr. $$ The exact expressions for $J_1$ and $J_2$ are $$ J_1(a, b) = \frac{1}{4} \left( a^2 - b^2 + (1-a^2)\log(1-a^2) - (1-b^2)\log(1-b^2) \right)\\ J_2(a, b) = \frac{b-a}{36} (5 a^2+5 a b-4 b^2-12) +{} \\ {} + \frac{ (3 a - 3ab^2 + 2b^3) \log (1-b^2) - a(3-a^2) \log (1-a^2) }{12} + {} \\ {} + \frac{\operatorname{arctanh}(b) - \operatorname{arctanh}(a)}{3} $$ I've obtained them using Wolfram Mathematica with some simplifications by myself.

Below is the implementation of the method in Python:

from scipy.special import ellipk
import numpy as np


def K(r):
    return ellipk(r*r)

def f(r):
    return r**3

Catalan = 0.91596559417721901505
exact = (13 + 18 * Catalan) / 64

def regular(r):
    if not isinstance(r, np.ndarray) and r == 1:
        return np.log(4)
    return r * (K(r) + 0.5 * np.log(1 - r*r))

def J1(a, b):
    if b != 1:
        return 0.25*(a*a - b*b +
                     (1-a*a)*np.log(1-a*a) - (1-b*b)*np.log(1-b*b))
    return 0.25 * (a*a - 1 + (1-a*a)*np.log(1-a*a))

def J2(a, b):
    if b != 1:
        t1 = (b-a) * (5*a*a + 5*a*b - 4*b*b - 12) / 36
        t2 = ((3*a - 3*a*b*b + 2*b*b*b) * np.log(1-b*b) -
              a*(3-a*a) * np.log(1-a*a)) / 12
        t3 = (np.arctanh(b) - np.arctanh(a)) / 3
    else:
        t1 = (1-a) * (5*a*a + 5*a - 16) / 36
        t2 = -a*(3-a*a) * np.log(1-a*a) / 12
        t3 = np.log(2) / 3 - np.arctanh(a) / 3
    return t1 + t2 + t3

def compute(M):
    rs = np.linspace(0, 1, M+1)
    dr = rs[1]
    I1 = 0.5 * dr * (regular(0) * f(0) + regular(1) * f(1))
    I1 += dr * np.sum(regular(rs[1:-1]) * f(rs[1:-1]))

    I2 = 0
    for i in range(M):
        a = rs[i]
        b = rs[i+1]
        der = (f(b) - f(a)) / dr
        I2 += f(a) * J1(a, b) + der * J2(a, b)

    return I1 - I2

for M in (100, 200, 500, 1000):
    approx = compute(M)
    print('M =', M, 'error =', exact - approx)

# Produces output
# M = 100 error = -4.296346864496314e-05
# M = 200 error = -1.0382533520481019e-05
# M = 500 error = -1.585096305212197e-06
# M = 1000 error = -3.818538253375081e-07
```