# Numerical Integral of Complete Elliptic Integral

asymptoticselliptic integralsnumerical methods

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?

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
$$$$
`