The order of error only makes sense for the composite rule, since for the simple rule the step size $h=x_f-x_0=b-a$ is constant. Only for the composite rule do you get a variable $h=(b-a)/n$ that allows to consider the asymptotic error behavior. Thus $O(h^2)$ for the composite rule, and no asymptotic error for the simple rule.
Since the (composite) Simpson rule can be seen as Richardson extrapolation (first step of the Romberg method) of the symmetric trapezoidal rule, its error order is automatically $O(h^4)$.
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
```
Best Answer
You have already found the value of $\log(2)$. Congrats!
Hint:
$$\int_0^1\frac{\mathrm{d}x}{1+x}=\log(2)$$
So there you have your context and what you have found is the value of the integration i.e. $\log(2)$ by using numerical methods (Trapezoidal rule).