用 python 實現定積分 ∫ₐᵇ f(t)dt

看了一些文章後, 自己手動寫了一些簡單的程式, 實現各種定積分的方式
import numpy as np
def Lp(x, n) : # evaluate order n-1 Legendre polynomials at x
  if (n == 0) :
    return 1.0
  if (n == 1) :
    return x
  pn_1, pk_1 = x, 1.0 # 疊代初始化, 因為 n > 1, 所以至少疊代一次
  k = 1 # 此刻從 k 開始, 直到 n - 1
  while k < n :
    pn_1, pk_1 = (pn_1*x*(2 * k + 1) - pk_1*k) / (k + 1), pn_1
    k += 1 # # 因為 k=n-1 所以 n=k+1, 2*n-1 = 2*(k+1)-1 = 2*k+1
  return pn_1  # (Lp(x, k)*x*(2*n-1) - Lp(x, k-1)*k) / n, k = n - 1

def d_Lp(x, n) : # derivative of order n-1 Legendre Polynomials at x
  if (n == 0) :
    return 0.0
  if (n == 1) :
    return 1.0;# LP′ₙ(x) = (−LPₙ(x)*x + LPₖ(x)) * n / (1−x*x), k=n-1
  return (Lp(x, n - 1) - Lp(x, n) * x) * n / (1.0 - x*x)
# Multiple root finder algorithm for Legendre polynomial: https://publikacio.uni-eszterhazy.hu/3009/1/AMI_33_from3to13.pdf
def by_newton(f, a=1.0, b=2.0, n=10):
  def newton_raphson(n, eps=1e-16): # using: xₖ = xₖ - f(xₖ) / [f'(xₖ) - f(xₖ) Σₖ 1/(xₖ - xᵢ)]
    e = np.cos(np.array([np.pi*(k + 0.5)/n for k in range(n)])) #initial guess
    for k in range(n) : # find all roots by newton raphson method
      xk = e[k]
      iteration = 0
      while (iteration < 1000) :
        iteration += 1
        f = Lp(xk, n)
        temp = f if (f > 0) else -f
        if (temp < eps): # 收斂
          break
        sum_r = 0.0
        for j in range(k) :# sum_r = Σₖ 1/(xₖ - xᵢ) to remove previouse root
          delta = xk - e[j]
          temp = delta if (delta < 0) else -delta
          if (temp < eps) :
            continue# skip singular value!
          sum_r += 1.0 / delta      
        dx = f / (d_Lp(xk, n) - f * sum_r)
        temp = dx if (dx > 0) else -dx
        if (temp < eps) : # 收斂
          break
        xk -= dx # xₖ = xₖ - f(xₖ) / [f'(xₖ) - f(xₖ) Σₖ 1/(xₖ - xᵢ)]
      e[k] = xk # final root update
    xi = [e]
    derivative = d_Lp(xi[0], n)
    xi.append(2 / (derivative*derivative*(1 - xi[0]*xi[0])))
    return np.array(xi)
  xw = newton_raphson(n)
  scale = (b - a) / 2.0 # linear transform tx: bias + scale *  1  = b => scale = (b - a) / 2
  bias  = (b + a) / 2.0  # linear transform tx: bias + scale *(-1) = a => bias  = (b + a) / 2  
  tx = xw[0] * scale + bias # x -> tx
  sum = f(tx).dot(xw[1])
  return sum * scale

def by_jacobi(f, a=1.0, b=2.0, n=10):
  def Jacobi_method(n, eps=1e-16): # using: xₖ = xₖ - f(xₖ) / [f'(xₖ) - f(xₖ) Σₖ 1/(xₖ - xᵢ)]
    J = np.zeros((n, n))
    d = len(J) - 1 # to fill into the jacobian tridiagonal matrix, trace = 0, and it is a symmetry matrix
    for k in range(d): # fill
      m = k + 1
      J[k][m] = m / np.sqrt(4.0 * m * m - 1.0)
      J[m][k] = J[k][m]
    e = np.linalg.eigvals(J) # nxn jacobi matrix, solve eigenvalues     
    xi = [np.sort(e)] # eigenvalue is same as root of the order n-1 Legendre polynopmial  
    derivative = d_Lp(xi[0], n)
    xi.append(2 / (derivative*derivative*(1 - xi[0]*xi[0])))# weight relative eigenvalue
    return np.array(xi)
  xw = Jacobi_method(n)
  scale = (b - a) / 2.0 # linear transform tx: bias + scale *  1  = b => scale = (b - a) / 2
  bias  = (b + a) / 2.0  # linear transform tx: bias + scale *(-1) = a => bias  = (b + a) / 2  
  tx = xw[0] * scale + bias # x -> tx
  sum = f(tx).dot(xw[1])
  return sum * scale

def by_trapezoid(f, a=1.0, b=2.0, n=50): # integral with the trapezoid method, 使用 梯形面積=(上底 + 下底)/2 積分
  x = np.linspace(a, b, n) # total n pints include a, b
  y = f(x) # total n
  n -= 1   # split to n - 1 interval
  delta_x = (b - a) / n
  sum = y[1:n].sum() + (y[0] + y[n]) / 2.0
  return sum * delta_x

def by_simpson(f, a=1.0, b=2.0, n=50): # (b-a)/3 Σ[f(a) + 4f(a+h) + 2f(a + h) + 4f(a+2h) + 2f(a+3h)  ... + f(b)]
  if n % 2 == 1:
    n += 1
  dx = (b - a) / n
  ddx = dx + dx # double dx
  x4 = a + dx   # 2nd term * 4
  x2 = a + ddx  # 3rd term * 2
  sum = f(a) + f(b) # head + tail
  for i in range(1, n - 2, 2) : # exclude head and tail
    sum += f(x4) * 4 + f(x2) * 2
    x4 += ddx # next
    x2 += ddx # next    
  sum += f(x4) * 4 # last one
  return sum * dx / 3
 
f = lambda x: np.exp(-x**2)
a = 0.0
b = 10.0
n = 15
anser  = np.sqrt(np.pi) / 2
jacobi = by_jacobi(f, a, b, n)
newton = by_newton(f, a, b, n)
trapezoid = by_trapezoid(f, a, b, n)
simpson = by_simpson(f, a, b, n)
print(f"∫ :\tjacobi={jacobi}  \t, newton={newton} \t, trapezoid={trapezoid}\t, simpson={simpson}\t, compare to {anser}")
print(f"Δ :\t {jacobi - anser} \t, {newton - anser}\t, {trapezoid - anser}  \t, {simpson - anser}")