前言
最近需要计算一下曲线长度，无法直接得到被积函数的原函数，常规的积分解法牛顿莱布尼茨公式没法使用，所以只能使用数值积分计算积分。
前言
最近需要计算一下曲线长度，无法直接得到被积函数的原函数，常规的积分解法牛顿莱布尼茨公式没法使用，所以只能使用数值积分计算积分。
下面主要介绍两种数值积分方法：龙贝格积分（Romberg Quadrature） 和 高斯-克朗罗德积分（Gauss-kronrod Quadrature）。 下面附带的代码只做简单测试，不保证正确性，语言使用 Typescript。
Romberg 篇
计算积分最简单的当然是使用复化梯形公式，即 \(I=\int_a^b{f(x)dx}=\frac{b-a}{2n}[f(a)+f(b)+2\sum\limits_{i=1}^{n-1}f(x_i)]= T_n, x_i=a+i*h, h=(b-a)/n\) ，若将 n 段每段一分为 2，可得到 \(T_{2n}=T_n/2+\frac{b-a}{2n}\sum\limits_{i=0}^{n-1}f(x_{i+1/2})\) 。考虑数列 \(T=\{T_1,T_2,T_{2^2},...,T_{2^k}\}\)，显然该数列必收敛，最后收敛为对应积分，当 \(|T_{2^k}-T_{2^{k-1}}| < ε\) （\(ε\) 为精度）时，可取 \(T_{2^k}\) 作为最后的积分结果。但是，直接利用梯形公司求解，收敛很慢，导致计算效率很差，所以需要使用理查德森（Richardson）外推法加速收敛，设 \(T_{2^k}^{(m)}\) 为 m 次加速值，当 m 为 0 时，表示没加速，为原梯形公司，则 \(T_{2^k}^{(m)} = \frac{4^m}{4^m-1}T_{2^{k+1}}^{(m-1)}-\frac{1}{4^m-1}T_{2^k}^{(m-1)}\)，当 \(|T_{2^{k+1}}^{(m)}-T_{2^k}^{(m)}| < ε\) 时，则收敛，并取其中一值作为最终的积分值。未经修饰的代码如下：
1 | function rombergIntegrator(f: (x: number) => number, a: number, b: number, tolerance = 1e-8) { |
由于采用闭型积分规则（积分上下限值参与积分计算），所以以上代码不适合计算两端点值被积函数值无限大的情况（如 1/4 圆弧积分等）。而且该方法不合适求取被积函数在积分区间内导数值变化较大（如被积函数在积分下限附近剧烈波动，在积分上限附近不变化等）的积分，因为该方法是均匀分段，这种情况将导致计算量剧增，这时就需要用到下面的自适应积分。
自适应篇
自适应积分主要包括两类：全局自适应积分和局部自适应积分，通常情况下全局自适应积分的会比局部自适应积分的表现要好，全局自适应积分一般通过二分递归实现，当满足一定条件时，递归终止，即通过二分分别计算两边的积分，若一边满足一定条件，则不继续划分，从而减少计算量。全局自适应积分中比较经典的有基于辛普森（Simpson）公式的自适应算法，普通辛普森积分公式为：\(I=\int_a^b{f(x)dx}=\frac{b-a}{6}[f(a)+4f(m)+f(b)]= S(a,b), m=(a+b)/2\)，复化辛普森公式为 \(I=\int_a^b{f(x)dx}=\frac{h}{3}[f(a)+4\sum\limits_{i=1}^{n}f(x_{2i-1})+2\sum\limits_{i=1}^{n-1}f(x_{2i})+f(b)]= S(a,b)\)，其中 \(x_i=a+i*h (i=1,2,3,...,2n),h=\frac{b-a}{2n}\)，基于辛普森公式的自适应基本原理如下：令 \(S_2(a,b) = S(a,m)+S(m,b)\)，m 为 a,b 中值，若 \(|S(a,b) - S_2(a,b)| < 15ε\)，则取 \(S_2(a,b)\) 或 \(S_2(a,b)+(S(a,b)-S_2(a,b))/15\) 作为该区间的积分值，否则，将区间二分递归，同时因为误差会累积，所以每次递归都需要将精度提高两倍，即 \(ε = ε/2\)，如此最后的精度才能达到最初的精度。具体 ts 代码如下：
1 | function adaptiveSimpsonIntegrator(f: (x: number) => number, a: number, b: number, epsilon = 1e-8) { |
在 D.V. Fedorov 写的「Introduction to Numerical Methods with examples in Javascript」一书中介绍了一种全局自适应方法，即分别使用高阶和低阶的权值分别计算积分，两者之间的差值 \(E\) 作为误差估计，设绝对精度为 \(\delta\) ，相对精度为 \(\epsilon\) ，若 \(|E|<\delta+\epsilon*Q\)，Q 为高阶权值计算的积分，则取 Q 作为积分值，否则将积分区间二分，同时使 \(\delta/\sqrt{2}\)，\(\epsilon\) 保持不变。具体 ts 代码如下：
1 | function recursiveAdaptiveIntegrator(f: (x: number) => number, a: number, b: number, accuracy = 1e-15) { |
该方法很巧妙的设计了一组插值点，使得当前计算的函数值正好可以被下次迭代所使用，从而提高性能，同时该方法可以得到 1/4 圆弧长，虽然精度只能到小数点后 8 位，至于 Shaun 写的其它测试函数，都能得到理想精度。
Gauss 篇
高斯求积法是一种多项式插值积分法，同时由于不计算被积函数在区间两个端点处的值，所以高斯积分法采用的开型积分规则，高斯积分法的衍生方法有很多种，下面主要介绍高斯-勒让德（Gauss-Legendre Quadrature）以及其迭代改良的高斯-克朗罗德法。高斯-勒让德积分法的公式为积分的原始形态，即 \(\int_a^bf(x)dx=\sum\limits_{i=1}^{∞}w_if(x_{i})\approx\sum\limits_{i=1}^{n}w_if(x_{i})\) ，只不过 \(x_i \in [-1,1]\)，并且 \(x_i\) 和 \(w_i\) 都通过勒让德多项式求出，所以其原则上只能用在积分区间为 [-1, 1] 上的积分，但是可以将积分从任意区间通过简单的变形变换到 [-1, 1] 上，即 \(\int_a^b{f(x)dx} = \frac{b-a}{2}\int_{-1}^1{f(\frac{b-a}{2}t+\frac{b+a}{2})dt}\) ，从而可以将高斯-勒让德方法扩展到任意积分上。由于每个 n 对应的 \(x_i\) 和 \(w_i\) 都可以查表可得，所以具体代码方面就很简单了，以 n = 4，即插值点个数为 4 为例，ts 代码如下：
1 | function gaussLegendreIntegrate(f: (x: number) => number, a: number, b: number, n: 4 | 8 = 4) { |
若要提高高斯-勒让德积分法的精度，可通过增加插值点或分多个区间进行积分来实现，但是由于没有误差估计，所以还是没法精确控制精度，对与某些被积函数积分精度高，但对于其它被积函数，积分精度却有限，当然可以简单的引入一些常用的误差估计法，但一般需要重新计算积分，导致效率很低，而高斯-克朗罗德法为其引入了一种基于 Kronrod 点的误差估计法，可充分利用现有计算值，从而达到有效控制精度的同时，性能没有太大的损失。设 \(G(f,n)=\sum\limits_{i=1}^{n}w_if(x_{i})\) 为具有 n 个插值点的高斯-勒让德法计算结果，\(GK(f,n) = \sum\limits_{i=1}^{n}w'_if(x_{i})+\sum\limits_{k=n+1}^{2n+1}w'_kf(x_{k})\) 为高斯-克朗罗德法的计算结果，则 \(|GK(f,n)-G(f,n)|\) 可作为误差估计，有了误差估计，最后再使用全局自适应策略，即可得到精度可控的高斯积分结果。具体 ts 代码如下，以 n = 7 为例：
1 | function gaussKronrodIntegrator(f: (x: number) => number, a: number, b: number, accuracy = 1e-15) { |
简单测试了一下，Shaun 这里写的 gaussKronrodIntegrator 方法最大精度只能到 1e-15，到 16 位就报错递归深度太大了，圆的 1/4 弧长也没法算出来，当然这些问题可通过设置最大递归深度以及处理异常值来解决，Shaun 这里就不继续写了。
后记
数值积分策略非常多，尤其是针对一些特殊的函数，可能只能使用一些特殊的策略才能计算，Shaun 这里只是介绍了一些比较基础常用的积分方法，能解决大部分积分问题，唯一需要注意一点的就是如何追求性能与精度之间的平衡，因为积分常常涉及到迭代求值，通常而言精度越高，迭代越多，求积时，同时也需要注意被积函数的异常值（如无穷大等），这时可能需要拆分或变换积分区间，并且使用开型积分规则的积分方法进行重新计算。
附录
一些常见的积分变换
\[ \int_a^bf(x) = \int_0^{b-a}f(a+t)dt \\ \int_a^b{f(x)dx} = \frac{b-a}{2}\int_{-1}^1{f(\frac{b-a}{2}t+\frac{b+a}{2})dt} \\ \int_{-∞}^{+∞}{f(x)dx} = \int_{-1}^1{f(\frac{t}{1-t^2})\frac{1+t^2}{(1-t^2)^2}dt} \\ \int_{a}^{+∞}{f(x)dx} = \int_{0}^1{f(a + \frac{t}{1-t})\frac{1}{(1-t)^2}dt} \\ \int_{-∞}^{a}{f(x)dx} = \int_{-1}^0{f(a + \frac{t}{1+t})\frac{-1}{(1+t)^2}dt} \]
参考资料
[1] 积分策略