Euler-Maclaurin: correcting the trapezoidal rule

Adds Euler-Maclaurin endpoint corrections, built from Bernoulli numbers and odd derivatives at the endpoints, to the trapezoidal rule. Each extra term raises the convergence order by two.

import matplotlib.pyplot as plt
import numpy as np

from mathematicskit.calculus import euler_maclaurin_trapezoid

Corrections for the integral of exp over [0, 1]

exact = np.e - 1
for m in range(4):
    value = euler_maclaurin_trapezoid(np.exp, 0.0, 1.0, 8, [np.exp] * m).value
    print(f"{m} correction terms, n = 8: error = {abs(value - exact):.2e}")

0 correction terms, n = 8: error = 2.24e-03
1 correction terms, n = 8: error = 5.82e-07
2 correction terms, n = 8: error = 2.17e-10
3 correction terms, n = 8: error = 8.46e-14

Convergence with 0, 1, and 2 correction terms

ns = np.array([2, 4, 8, 16, 32, 64])
fig, ax = plt.subplots()
for m in range(3):
    errors = [abs(euler_maclaurin_trapezoid(np.exp, 0.0, 1.0, int(n), [np.exp] * m).value - exact) for n in ns]
    ax.loglog(ns, errors, "o-", label=f"{m} corrections: O(h^{2 * m + 2})")
ax.set_xlabel("n")
ax.set_ylabel("|error|")
ax.legend()
ax.set_title("Euler-Maclaurin corrections to the trapezoidal rule")

Text(0.5, 1.0, 'Euler-Maclaurin corrections to the trapezoidal rule')