Blog. 用 MoonBit 实现 Gauss-Kronrod 数值求积算法 [blog/gausskronrod]

Rutubet
2025-09-03

1. Gauss-Kronrod 求积算法简介 [blog/gausskronrod/intro]

数值积分是科学计算和工程领域的基础工具，主要用于解决无法解析求解的复杂积分问题，例如有限元分析、物理建模、概率统计、算法优化。

现代编程语言通过数值积分算法实现高效数值积分，显著提升了工程计算的效率（如 Scipy 使用了 Gauss-Kronrod 求积算法进行数值积分）。

与传统求积算法相比，Gauss-Kronrod 算法是比较优秀的求积算法。Newton-Cotes 公式依赖等距节点，高阶时易出现龙格现象（Runge’s phenomenon），精度受限。纯 Gauss 求积法通过正交多项式节点实现高精度（2n - 1 次代数精度），但节点数固定，无法自适应细分区间。Gauss-Kronrod 算法在 Gauss 节点基础上插入额外 Kronrod 节点，形成嵌套结构（如 7 阶 Gauss + 15 阶 Kronrod），同时计算高精度积分值并估计误差。

但是要认识到 Gauss-Kronrod 算法也有一定的局限性，它只适用于在闭区间上连续的函数进行求积，例如无法对反常积分进行求积。

2. Gauss-Kronrod 求积算法理论 [blog/gausskronrod/theory]

给出如下机械求积公式：

$$ \int^b_a \rho(x) f(x) dx \approx \sum^n_{k = 0} A_k f(x_k) \tag{1} $$

其中 $x_k$ 称为求积节点，$A_k$ 称为求积系数。

如果某个求积公式对于次数不超过 m 的多项式均能准确成立，而对于 m + 1 次的的多项式不准确成立，则称该公式具有 m 次代数精度。

具有最高代数精度 (2n + 1) 次的机械求积公式叫 Gauss 型求积公式，相应的求积节点叫 Gauss 点。

为了能够探索 Gauss 型求积公式的节点和系数的形式，引入以下理论：

用 n 次 Lagrange 插值多项式 $L_n(x)$ 近似被积函数 $f(x)$ 得到：

$$ I_n(f) = \int^b_a L_n(x) dx = \int^b_a \sum^n_{ j = 0} f(x_j) l_j(x) dx $$ $$ I_n(f) = \sum^n_{ j = 0 } f(x_j) \int^b_a l_j(x) dx \tag{2} $$

式 $(2)$ 称为插值型求积公式，其求积系数 $A_j = \int^b_a l_j(x) dx$ 。

定理 1 设有插值型求积公式 $(2)$，其求积节点 $x_j, j = 0, 1, ..., n$ 是 Gauss 点的充分必要条件是

$$ \omega_(n + 1) (x) = \prod^n_{k = 0} (x- x_k) $$

与任何不超过 $n$ 次的代数多项式在 $[a, b]$ 上带权 $\rho(x)$ 正交。即任意 $p \in \Phi_n[a, b]$，有

$$ (\omega_{n + 1} , p) = \int^b_a \rho(x) \omega_{n + 1}(x) p(x) dx = 0 $$

令 $p(x)$ 取不同值，可以用待定系数法解得求积节点和求积系数。此处计算过程略去。取积分区间为 $[-1, 1]$，权函数 $\rho(x)=1$ 根据选择的不同的节点数量，可以解得固定的节点和系数。再通过积分区间变换，即可扩展到任意闭区间上进行求解。下面需要解决积分余项的问题，即分析误差以确保达到指定精度。

通过计算 n 阶高斯求积和 2n + 1 阶高斯求积，并将差值用作误差估计是可行的。然而，这并不是最优的，因为勒让德多项式（高斯求积的节点）的零点对于不同的阶数来说永远不会相同，因此必须执行 3n + 1 函数求值。Kronrod 考虑了如何将节点交织成高斯求积的问题，以便可以重用所有先前的函数求值，同时提高可以精确积分的多项式的阶数。这允许提供后验误差估计，同时仍然保持指数收敛。Kronrod 发现，通过将 n + 1 个节点（根据 Legendre-Stieltjes 多项式的零点计算）添加到 n 阶高斯求积中，他可以对 3n+1 阶多项式进行积分。

同理，可以取积分区间为 $[-1, 1]$，权函数 $\rho(x)=1$ 根据选择的不同的节点数量，可以解得固定的 Gauss-Kronrod 求积法的节点和系数。

3. Gauss-Kronrod 求积算法实现 [blog/gausskronrod/implement]

如下为 MoonBit 实现的 15 (7 + 8) 点 Gauss-Kronrod 求积算法，其中 abscissae 为求积节点，weights 为求积系数，均为取积分区间为 $[-1, 1]$，权函数 $\rho(x)=1$ 计算出来的值，可以通过查表得到。为达到精度要求，使用二分法进行递归求积。

pub fn quad(func : (Double) -> Double, a : Double, b : Double) -> Double {
  fn recursiveQuad(
    func : (Double) -> Double,
    a : Double,
    b : Double,
    accuracy : Double,
  ) -> Double {
    let gaussAbscissae : FixedArray[Double] = [
      0.0, -4.058451513773971669066064120769615e-1, 4.058451513773971669066064120769615e-1,
      -7.415311855993944398638647732807884e-1, 7.415311855993944398638647732807884e-1,
      -9.491079123427585245261896840478513e-1, 9.491079123427585245261896840478513e-1,
    ]
    let gaussWeights : FixedArray[Double] = [
      4.179591836734693877551020408163265e-1, 3.818300505051189449503697754889751e-1,
      3.818300505051189449503697754889751e-1, 2.797053914892766679014677714237796e-1,
      2.797053914892766679014677714237796e-1, 1.29484966168869693270611432679082e-1,
      1.29484966168869693270611432679082e-1,
    ]
    let kronrodAbscissae : FixedArray[Double] = [
      0.0, -4.058451513773971669066064120769615e-1, 4.058451513773971669066064120769615e-1,
      -7.415311855993944398638647732807884e-1, 7.415311855993944398638647732807884e-1,
      -9.491079123427585245261896840478513e-1, 9.491079123427585245261896840478513e-1,
      -2.077849550078984676006894037732449e-1, 2.077849550078984676006894037732449e-1,
      -5.860872354676911302941448382587296e-1, 5.860872354676911302941448382587296e-1,
      -8.648644233597690727897127886409262e-1, 8.648644233597690727897127886409262e-1,
      -9.914553711208126392068546975263285e-1, 9.914553711208126392068546975263285e-1,
    ]
    let kronrodWeights : FixedArray[Double] = [
      2.094821410847278280129991748917143e-1, 1.903505780647854099132564024210137e-1,
      1.903505780647854099132564024210137e-1, 1.406532597155259187451895905102379e-1,
      1.406532597155259187451895905102379e-1, 6.309209262997855329070066318920429e-2,
      6.309209262997855329070066318920429e-2, 2.044329400752988924141619992346491e-1,
      2.044329400752988924141619992346491e-1, 1.690047266392679028265834265985503e-1,
      1.690047266392679028265834265985503e-1, 1.04790010322250183839876322541518e-1,
      1.04790010322250183839876322541518e-1, 2.293532201052922496373200805896959e-2,
      2.293532201052922496373200805896959e-2,
    ]
    let halfH = (b - a) / 2
    let mut gaussResult = 0.0
    let mut kronrodResult = 0.0
    for i = 0; i < gaussAbscissae.length(); i = i + 1 {
      let xi = halfH * gaussAbscissae[i] + a + halfH
      let yi = func(xi)
      gaussResult += gaussWeights[i] * yi
      kronrodResult += kronrodWeights[i] * yi
    }
    for i = gaussAbscissae.length(); i < kronrodAbscissae.length(); i = i + 1 {
      let xi = halfH * kronrodAbscissae[i] + a + halfH
      let yi = func(xi)
      kronrodResult += kronrodWeights[i] * yi
    }
    fn abs(x : Double) {
      if x < 0 {
        -x
      } else {
        x
      }
    }

    if abs(kronrodResult - gaussResult) * halfH < accuracy {
      kronrodResult * halfH
    } else {
      let m = a + (b - a) / 2
      let acc = accuracy / 2
      recursiveQuad(func, a, m, acc) + recursiveQuad(func, m, b, acc)
    }
  }
  
  recursiveQuad(func, a, b, 1.0e-15)
}

4. 参考文献 [blog/gausskronrod/ref]

https://cniter.github.io/posts/eee1c041.html

https://zhuanlan.zhihu.com/p/376594584

https://zhuanlan.zhihu.com/p/374569089

https://zhuanlan.zhihu.com/p/372961979

https://zhuanlan.zhihu.com/p/98431641

https://www.boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/gauss_kronrod.html