用 c++ 寫個簡單的 bisecton 方法找出方程式的實數根

數學上有個著名的多項式稱為 Legendre polynomial, 該方程式的根及權重可以用來算出其它方程式的定積分值,底下簡單利用雙端夾擊法(英文是 bisection method)找出 Legendre 多項式的根,之後利用這些參數快速估算出其他方程式的定積分,例如: 底下算出 10 階 Legendre 多項式的參數(根及權重), 用這些參數推算出 ∫(2*x + 3/x)²dx  定積分值,還蠻貼近的. 可以用 python 搭配 scipy 下指令 scipy.integrate.quad(y, 1, 100) 驗證結果, 其中 y = lambda x: (2*x + 3/x)**2
// quad.cpp
#include <stdio.h>
#include <math.h>
double Lp(double x, int n) {// 快速疊代法, order n-1 Legendre polynomial
  if (n <  0) return 0; // Invalid order
  if (n == 0) return 1; // p0 = 1
  if (n == 1) return x; // p1 = x
  double pk_1 = 1;// 前一刻初始值 Lp(x, k - 1)
  double pn_1 = x;// 此刻初始值   Lp(x, k)
  int k = 1; // 此刻, begin to evaluate order n-1 Legendre polynomial
  do {  // pk_n = (pk*x*(2*n-1)-pk_1*k)/n, k=n-1 => n=k+1, 2*n-1=2*(k+1)-1=2*k+1
    double pk = (pn_1 * x * (2 * k + 1)  - pk_1 * k) / (k + 1);// 下一刻 pk 值
    pk_1 = pn_1;// 下一刻 pk_1 疊代
    pn_1 = pk;  // 下一刻 pn_1 疊代
  } while (++ k < n); // upto n => order n - 1
  return pn_1;// final interation of order n-1 Legendre polynomial
}
double d_Lp(double x, int n) { // derivative of order n-1 Legendre Polynomial
  if (n <= 0) return 0;
  if (n == 1) return 1;// LP′ₙ(x) = (− x*LPₙ(x) + LPₖ(x)) * n / (1 − x*x), k = n-1
  return (Lp(x, n - 1) - Lp(x, n) * x) * n / (1.0 - x*x);
}
double Lp_find_root(double left, double right, int nth,int iteration=1000, double eps=1e-16) {
  auto LPn = [nth](double x){ return Lp(x, nth); }; // 綁定 Lp(nth, .) 函式
  if (LPn(left)  == 0) return left ;// 先測左邊界, 若函數值等於0, 毫無懸念, 一定是根
  if (LPn(right) == 0) return right;// 再測右邊界, 若函數值等於0, 毫無懸念, 一定是根
  auto zero_cross = [LPn](double l, double r) { return LPn(l) * LPn(r); };// 再綁定上面的 LPn(.) 函式
  double product = zero_cross(left , right);
  if (product > 0) return 2; // 當左右邊界, 正負符號相同, 不可能有根
  double root = 0;   // Lp 函式的根介於 -1 到 1 之間  
  int step = 0; // 開始疊帶, to find root of Legendre Polynomials from left to right step by step
  double pre_root = 0;
  while (step ++ < iteration) {// zero cross 可以確認中間有否有交點 {0}, 當一邊是正, 另一邊是負, 乘積是負的, 必有根
    root = (left + right) / 2; // 左右兩邊夾擊取中點當作根
    double f = LPn(root);
    if (f < 0) f = -f; // 函數值取絕對值
    if (f < eps) break;// 勉強找到根了
    double temp = pre_root - root;// 前後根差異
    pre_root = root;
    if (temp < 0) temp = -temp; // 取絕對值
    if (temp < eps) {// todo: 尚可接受
      break;
    }
    product = zero_cross(left, root);// 待決定的根與左邊界, 算出零交越值
    if (product > 0) left = root; // 若與左邊界正負符號相同, 表示根在右邊,換掉左邊界, 下次再試
    else if (product < 0) right = root; //  與左邊界符號不同, 確認有零交越, 換掉右邊界, 下次再試
  }
  return root;// 介於 left 與 right 之間的根
}
double *Lp_solver(int nth = 10) {
  static double xi[1000 * 2];// todo: validate n
  for (int i = 0; i <= 2 * nth; i++) xi[i] = 0; // clear all data first
  int k = 0;
  double step = 1e-3;
  for (double x = -1.0; x <= 1.0; x += step) { // step by step
    double root = Lp_find_root(x, x + step, nth); // find root between x and x + steps
    if (root > 1) continue;
    double derivative = d_Lp(root, nth);
    xi[k] = root;// position of x
    xi[nth + k ++] = 2 / (derivative*derivative*(1 - root*root));// weight of w
  }
  return (k > 0) ? xi : nullptr;
}
double integral_quad(double f(double), double a, double b, int nth = 10) {
    static double *xi;
    static int Ln;
    if (xi == nullptr || Ln != nth) {
        Ln = nth; // reset Ln
        xi = Lp_solver(Ln);// solve once
    }
    double *wi = xi + Ln;
    double scale = (b - a) / 2; // linear transform xd: bias + scale *  1  = b , scale  = (b - a) / 2
    double bias  = (b + a) / 2; // linear transform xd: bias + scale *(-1) = a , bias = (b + a) / 2
    double sum = 0;
    for(int i = 0; i < Ln; i ++) {
        sum +=  f(xi[i] * scale + bias) * wi[i];
    }
    double area = sum * scale;
    return area;
}
int main() {
  double area = integral_quad([](double x) { // f(x) = (2*x + 3/x)²
      double y = 2*x + 3/x;// 要避開 x = 0
      return y * y;
    },
    1, 100
  );
  printf("Area = %f\n", area);
}

後記: 找方程式的根有很多種方法, 例如知名的 newton raphson 法也可以迅速找到根, 只不過算出全部的根需要用點技巧. 最近玩 xAI(網址是 http://grok.com )時, 我順便問一下如何解 Legendre polynomial 的根, 它給出了我蠻意外的答案, 他用了一個 tridiagonal matrix (只有3重對角線上有值, 其餘為零)的對稱型方陣, 這個矩陣 J 是 n*n Jacobian 方陣, 在正對角線上的值全為 0, 也就是說  trace 等於 0, 換句話說全部 eigenvalues 加起來等於 0, 這是一個3重對角型對稱方陣(詳細如下列程式所述), 只要求出該矩陣的 eigenvalues 就是 n 階 Legedre polynomial 的根, 底下我用  python 語法來解 5x5 矩陣特性方程式的根(eigenvalue), 還真的與 5 階 Legedre 多項式的根一模一樣. 只不過浮點運算稍有誤差, 接近 +- 1e-17 的值應當視為 0:

    import numpy as np
    n = 5
    J = np.zeros((n,n))
    d = len(J) - 1
    for k in range(d):
        m = k + 1
        J[k][m] = m / np.sqrt(4.0 * m * m - 1.0)
        J[m][k] = J[k][m]
    e, v = np.linalg.eig(J) # 解出方陣的 eigenvalues 及 eigenvectors
    print(np.sort(e)) # 只列出 eigenvalues 根

有了根之後, 權重(weight)自然就容易算出來,  至於為何會如此呢?  The characteristic polynomial of J​, det⁡(J − λI) turns out to be the Legendre polynomial P(λ), 這句話我能理解 eigenvalue 是矩陣特性方程式的根, 但它是如何與 Legendre polynomial 產生連結的呢? 我還無法理解. 也許要上高深一點線性代數理論才有辦法明瞭, 總之要算出 500 階以上的 Legendre 多項式的根就變得容易多了. 底下是 xAI 詳細回答:

     https://grok.com/share/bGVnYWN5_c43a2f4c-6290-400d-8834-ab6b5dd5bd1e