积分计算
Mathematics

前言

　　最近需要计算一下曲线长度，无法直接得到被积函数的原函数，常规的积分解法牛顿莱布尼茨公式没法使用，所以只能使用数值积分计算积分。

　　下面主要介绍两种数值积分方法：龙贝格积分（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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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=\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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
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\) 作为误差估计，设绝对精度为 \(\delta\) ，相对精度为 \(\epsilon\) ，若 \(|E|<\delta+\epsilon*Q\)，Q 为高阶权值计算的积分，则取 Q 作为积分值，否则将积分区间二分，同时使 \(\delta/\sqrt{2}\)，\(\epsilon\) 保持不变。具体 ts 代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
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）以及其迭代改良的高斯-克朗罗德法。高斯-勒让德积分法的公式为积分的原始形态，即 \(\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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
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 为例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
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

本文标题： 积分计算

文章作者： Shaun

发布时间： 2020年11月29日 - 18时28分

最后更新： 2021年12月18日 - 19时54分

原始链接： http://cniter.github.io/posts/eee1c041.html

许可协议： "署名-非商业性使用-相同方式共享 4.0 国际" 转载请保留原文链接及作者。