最近需要计算一下曲线长度，无法直接得到被积函数的原函数，常规的积分解法牛顿莱布尼茨公式没法使用，所以只能使用数值积分计算积分。
下面主要介绍两种数值积分方法：龙贝格积分（Romberg Quadrature） 和 高斯-克朗罗德积分（Gauss-kronrod Quadrature）。 下面附带的代码只做简单测试，不保证正确性，语言使用 Typescript。
Romberg 篇
计算积分最简单的当然是使用复化梯形公式，即 I=∫baf(x)dx=b−a2n[f(a)+f(b)+2n−1∑i=1f(xi)]=Tn, xi=a+i∗h,h=(b−a)/n ，若将 n 段每段一分为 2，可得到 T2n=Tn/2+b−a2nn−1∑i=0f(xi+1/2) 。考虑数列 T={T1,T2,T22,...,T2k}，显然该数列必收敛，最后收敛为对应积分，当 |T2k−T2k−1|<ε （ε 为精度）时，可取 T2k 作为最后的积分结果。但是，直接利用梯形公司求解，收敛很慢，导致计算效率很差，所以需要使用理查德森（Richardson）外推法加速收敛，设 T(m)2k 为 m 次加速值，当 m 为 0 时，表示没加速，为原梯形公司，则 T(m)2k=4m4m−1T(m−1)2k+1−14m−1T(m−1)2k，当 |T(m)2k+1−T(m)2k|<ε 时，则收敛，并取其中一值作为最终的积分值。未经修饰的代码如下：
| function rombergIntegrator(f: (x: number) => number, a: number, b: number, tolerance = 1e-8) {
let h = b - a;
let n = 1;
let preTK = ((f(a) + f(b)) * h) / 2;
let tkmArray = [];
let m = 0;
tkmArray[m] = preTK;
while (true) {
m += 1;
console.log(m, tkmArray[m - 1]);
let tk = getTK(f, preTK, n, a, h);
if (Math.abs(tk - preTK) < tolerance) return tk;
let preTKM = tkmArray[0];
let preTK1M = tk;
tkmArray[0] = tk;
for (let i = 1; i <= m; i++) {
let newTKM = getTKM(preTK1M, preTKM, i);
preTKM = tkmArray[i];
preTK1M = newTKM;
tkmArray[i] = newTKM;
if (preTKM !== undefined && Math.abs(preTK1M - preTKM) < tolerance) return preTK1M;
}
preTK = tk;
n *= 2;
h /= 2;
}
function getTK(f: (x: number) => number, preTK: number, n: number, a: number, h: number) {
let sum = 0;
for (let i = 0; i < n; i++) {
let x = a + (i + 0.5) * h;
sum += f(x);
}
return (preTK + h * sum) / 2;
}
function getTKM(preTK1M: number, preTKM: number, m = 0) {
let m4 = 1 << (2 * m); // 4 ** m;
return (m4 * preTK1M - preTKM) / (m4 - 1);
}
}
|
由于采用闭型积分规则（积分上下限值参与积分计算），所以以上代码不适合计算两端点值被积函数值无限大的情况（如 1/4 圆弧积分等）。而且该方法不合适求取被积函数在积分区间内导数值变化较大（如被积函数在积分下限附近剧烈波动，在积分上限附近不变化等）的积分，因为该方法是均匀分段，这种情况将导致计算量剧增，这时就需要用到下面的自适应积分。
自适应积分主要包括两类：全局自适应积分和局部自适应积分，通常情况下全局自适应积分的会比局部自适应积分的表现要好，全局自适应积分一般通过二分递归实现，当满足一定条件时，递归终止，即通过二分分别计算两边的积分，若一边满足一定条件，则不继续划分，从而减少计算量。全局自适应积分中比较经典的有基于辛普森（Simpson）公式的自适应算法，普通辛普森积分公式为：I=∫baf(x)dx=b−a6[f(a)+4f(m)+f(b)]=S(a,b), m=(a+b)/2，复化辛普森公式为 I=∫baf(x)dx=h3[f(a)+4n∑i=1f(x2i−1)+2n−1∑i=1f(x2i)+f(b)]=S(a,b)，其中 xi=a+i∗h (i=1,2,3,...,2n),h=b−a2n，基于辛普森公式的自适应基本原理如下：令 S2(a,b)=S(a,m)+S(m,b)，m 为 a,b 中值，若 |S(a,b)−S2(a,b)|<15ε，则取 S2(a,b) 或 S2(a,b)+(S(a,b)−S2(a,b))/15 作为该区间的积分值，否则，将区间二分递归，同时因为误差会累积，所以每次递归都需要将精度提高两倍，即 ε=ε/2，如此最后的精度才能达到最初的精度。具体 ts 代码如下：
| function adaptiveSimpsonIntegrator(f: (x: number) => number, a: number, b: number, epsilon = 1e-8) {
let S = complexSimpson(f, a, b);
return adsimp(f, a, b, epsilon, S);
function adsimp(f: (x: number) => number, a: number, b: number, epsilon = 1e-8, S = 0): number {
let m = a + (b - a) / 2;
let LS = complexSimpson(f, a, m);
let RS = complexSimpson(f, m, b);
const S2 = LS + RS;
const tolerance = 15 * epsilon;
let delta = S - S2;
if (Math.abs(delta) < tolerance) return S2 + delta / 15;
else {
let doubleEPS = epsilon / 2;
return adsimp(f, a, m, doubleEPS, LS) + adsimp(f, m, b, doubleEPS, RS);
}
}
function complexSimpson(f: (x: number) => number, a: number, b: number, n = 1) {
const n2 = n * 2;
const h = (b - a) / n2;
let sum = f(a) + f(b);
for (let i = 1; i < n2; i += 2) {
sum += 4 * f(a + i * h);
}
for (let i = 2; i < n2 - 1; i += 2) {
sum += 2 * f(a + i * h);
}
return (sum * h) / 3;
}
}
|
在 D.V. Fedorov 写的「Introduction to Numerical Methods with examples in Javascript」一书中介绍了一种全局自适应方法，即分别使用高阶和低阶的权值分别计算积分，两者之间的差值 E 作为误差估计，设绝对精度为 δ ，相对精度为 ϵ ，若 |E|<δ+ϵ∗Q，Q 为高阶权值计算的积分，则取 Q 作为积分值，否则将积分区间二分，同时使 δ/√2，ϵ 保持不变。具体 ts 代码如下：
| function recursiveAdaptiveIntegrator(f: (x: number) => number, a: number, b: number, accuracy = 1e-15) {
return recursiveAdaptiveIntegrate(f, a, b, accuracy);
function recursiveAdaptiveIntegrate(
f: (x: number) => number,
a: number,
b: number,
accuracy = 1e-15,
epsilon = Number.EPSILON,
preFValue?: number[],
): number {
const abscissae = [1 / 6, 2 / 6, 4 / 6, 5 / 6];
const highOrderWeights = [2 / 6, 1 / 6, 1 / 6, 2 / 6];
const lowOrderWeights = [1 / 4, 1 / 4, 1 / 4, 1 / 4];
const isRecompute = [1, 0, 0, 1];
const h = b - a;
const fValue: number[] = [];
if (preFValue === undefined) {
abscissae.forEach((abscissa) => {
const x = a + abscissa * h;
fValue.push(f(x));
});
} else {
for (let k = 0, i = 0; i < abscissae.length; i++) {
if (isRecompute[i]) fValue.push(f(a + abscissae[i] * h));
else fValue.push(preFValue[k++]);
}
}
let highResult = 0;
let lowResult = 0;
for (let i = 0; i < highOrderWeights.length; i++) {
highResult += highOrderWeights[i] * fValue[i];
lowResult += lowOrderWeights[i] * fValue[i];
}
highResult *= h;
lowResult *= h;
const tolerance = accuracy + epsilon * Math.abs(highResult);
let errorEstimate = Math.abs(highResult - lowResult) / 3;
if (errorEstimate < tolerance) {
return highResult;
} else {
accuracy /= Math.sqrt(2);
let m = a + h / 2;
let midIndex = Math.trunc(abscissae.length / 2);
let leftFValue = fValue.slice(0, midIndex);
let rightFValue = fValue.slice(midIndex);
return (
recursiveAdaptiveIntegrate(f, a, m, accuracy, epsilon, leftFValue) +
recursiveAdaptiveIntegrate(f, m, b, accuracy, epsilon, rightFValue)
);
}
}
}
|
该方法很巧妙的设计了一组插值点，使得当前计算的函数值正好可以被下次迭代所使用，从而提高性能，同时该方法可以得到 1/4 圆弧长，虽然精度只能到小数点后 8 位，至于 Shaun 写的其它测试函数，都能得到理想精度。
Gauss 篇
高斯求积法是一种多项式插值积分法，同时由于不计算被积函数在区间两个端点处的值，所以高斯积分法采用的开型积分规则，高斯积分法的衍生方法有很多种，下面主要介绍高斯-勒让德（Gauss-Legendre Quadrature）以及其迭代改良的高斯-克朗罗德法。高斯-勒让德积分法的公式为积分的原始形态，即 ∫baf(x)dx=∞∑i=1wif(xi)≈n∑i=1wif(xi) ，只不过 xi∈[−1,1]，并且 xi 和 wi 都通过勒让德多项式求出，所以其原则上只能用在积分区间为 [-1, 1] 上的积分，但是可以将积分从任意区间通过简单的变形变换到 [-1, 1] 上，即 ∫baf(x)dx=b−a2∫1−1f(b−a2t+b+a2)dt ，从而可以将高斯-勒让德方法扩展到任意积分上。由于每个 n 对应的 xi 和 w_i 都可以查表可得，所以具体代码方面就很简单了，以 n = 4，即插值点个数为 4 为例，ts 代码如下：
| function gaussLegendreIntegrate(f: (x: number) => number, a: number, b: number, n: 4 | 8 = 4) {
const weightsAndAbscissae = getWeightsAndAbscissae(n);
const weights = weightsAndAbscissae.weights;
const abscissae = weightsAndAbscissae.abscissae;
const halfH = (b - a) / 2;
let sum = 0;
for (let i = 0; i < abscissae.length; i++) {
let xi = halfH * abscissae[i] + a + halfH;
sum += weights[i] * f(xi);
}
return sum * halfH;
function getWeightsAndAbscissae(n: 4 | 8 = 4) {
switch (n) {
case 8:
return {
weights: [
0.362683783378362,
0.362683783378362,
0.3137066458778873,
0.3137066458778873,
0.2223810344533745,
0.2223810344533745,
0.1012285362903763,
0.1012285362903763,
],
abscissae: [
-0.1834346424956498,
0.1834346424956498,
-0.525532409916329,
0.525532409916329,
-0.7966664774136267,
0.7966664774136267,
-0.9602898564975363,
0.9602898564975363,
],
};
break;
case 4:
default:
return {
weights: [0.6521451548625461, 0.6521451548625461, 0.3478548451374538, 0.3478548451374538],
abscissae: [-0.3399810435848563, 0.3399810435848563, -0.8611363115940526, 0.8611363115940526],
};
break;
}
}
}
|
若要提高高斯-勒让德积分法的精度，可通过增加插值点或分多个区间进行积分来实现，但是由于没有误差估计，所以还是没法精确控制精度，对与某些被积函数积分精度高，但对于其它被积函数，积分精度却有限，当然可以简单的引入一些常用的误差估计法，但一般需要重新计算积分，导致效率很低，而高斯-克朗罗德法为其引入了一种基于 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 为例：
| function gaussKronrodIntegrator(f: (x: number) => number, a: number, b: number, accuracy = 1e-15) {
return recursiveAdaptiveIntegrate(f, a, b, accuracy);
function recursiveAdaptiveIntegrate(f: (x: number) => number, a: number, b: number, accuracy = 1e-12): number {
const gaussAbscissae = [
0.0,
-4.058451513773971669066064120769615e-1,
4.058451513773971669066064120769615e-1,
-7.415311855993944398638647732807884e-1,
7.415311855993944398638647732807884e-1,
-9.491079123427585245261896840478513e-1,
9.491079123427585245261896840478513e-1,
];
const gaussWeights = [
4.179591836734693877551020408163265e-1,
3.818300505051189449503697754889751e-1,
3.818300505051189449503697754889751e-1,
2.797053914892766679014677714237796e-1,
2.797053914892766679014677714237796e-1,
1.29484966168869693270611432679082e-1,
1.29484966168869693270611432679082e-1,
];
const kronrodAbscissae = gaussAbscissae.concat([
-2.077849550078984676006894037732449e-1,
2.077849550078984676006894037732449e-1,
-5.860872354676911302941448382587296e-1,
5.860872354676911302941448382587296e-1,
-8.648644233597690727897127886409262e-1,
8.648644233597690727897127886409262e-1,
-9.914553711208126392068546975263285e-1,
9.914553711208126392068546975263285e-1,
]);
const kronrodWeights = [
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,
];
const halfH = (b - a) / 2;
let guassResult = 0;
let kronrodResult = 0;
for (let i = 0; i < gaussAbscissae.length; i++) {
let xi = halfH * gaussAbscissae[i] + a + halfH;
let yi = f(xi);
guassResult += gaussWeights[i] * yi;
kronrodResult += kronrodWeights[i] * yi;
}
for (let i = gaussAbscissae.length; i < kronrodAbscissae.length; i++) {
let xi = halfH * kronrodAbscissae[i] + a + halfH;
let yi = f(xi);
kronrodResult += kronrodWeights[i] * yi;
}
if (Math.abs(kronrodResult - guassResult) < accuracy / halfH) return kronrodResult * halfH;
else {
let m = a + (b - a) / 2;
accuracy /= 2;
return recursiveAdaptiveIntegrate(f, a, m, accuracy) + recursiveAdaptiveIntegrate(f, m, b, accuracy);
}
}
}
|
简单测试了一下，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] 积分策略
[2] Gaussian Quadrature Weights and Abscissae
[3] Gauss-Kronrod Quadrature
[4] Gauss-Kronrod Quadrature Nodes and Weights
