65.9K
CodeProject 正在变化。 阅读更多。
Home

C++ 中的数值方法第 1 部分： 牛顿-科特斯积分

coding-thoughts
starIconstarIconstarIconstarIconstarIcon

5.00/5 (2投票s)

2019 年 4 月 17 日

CPOL

3分钟阅读

viewsIcon

7644

使用梯形和辛普森法则方法实现数值积分

欢迎来到 thoughts-on-cpp.com 的新文章。这次，我想开始一个新的系列，数值方法。但如果你期待关于 n 体问题的文章，请不要失望，我仍然计划继续“天啊，到处都是星星”系列。通过这个新系列，我想谈谈我在日常工作中非常重要的数值方法，从梯形和辛普森积分方法开始，属于 牛顿-柯特斯 类型。

牛顿-柯特斯公式是一种用于积分计算的二次数值逼近方法。其思想是用具有等距节点的某个多项式来插值将被积分的函数。该多项式可以采用对齐梯形或对齐抛物线的形式。常见的牛顿-柯特斯积分多项式规则有

通过这篇文章，我们将更深入地研究前两个积分规则，梯形和辛普森。

梯形积分近似

通过梯形的积分近似非常简单。我们所需要做的就是将示例函数

I=\int_0^{\pi/2} \frac{5}{e^\pi-2}\exp(2x)\cos(x)dx=1.0

我们想要积分成等距区域，我们可以对这些区域的精确积分进行求和。显然，近似的精度取决于细分 N 的数量。因此，每个区域的宽度定义为

\Delta x = \frac{b-a}{N}

其中 a=0 且 b=\pi/2 。

Approximation of integral of f(x) between 0 and pi/2 by trapezoids

用梯形近似 f(x) 在 0 和 pi/2 之间的积分

每个梯形的面积可以通过以下公式计算

A=\frac{a+b}{2}\cdot h

因此，我们函数的近似积分 \widetilde{ I } f(x) 可以定义为

\widetilde{ I } = \sum_{i=1}^{N+1} \frac{x_i-x_{i+1}}{2}(f(x_i)+f(x_{i+1}))

梯形积分的实现采用 4 个参数，积分函数 f(x) 的范围 [a,b]，细分的数量 n，以及函数 f(x)。

double trapezoidalIntegral(double a, double b, int n, const std::function<double (double)> &f) {
    const double width = (b-a)/n;

    double trapezoidal_integral = 0;
    for(int step = 0; step < n; step++) {
        const double x1 = a + step*width;
        const double x2 = a + (step+1)*width;

        trapezoidal_integral += 0.5*(x2-x1)*(f(x1) + f(x2));
    }

    return trapezoidal_integral;
}

辛普森积分近似

辛普森法则通过对抛物线 p(x) 在 a、b 和 m=\frac{a+b}{2} 处的节点进行精确积分来近似函数 f(x) 的积分。为了提高近似精度，该函数可以像梯形积分近似一样被细分为 N。

Approximation of integral of f(x) between 0 and pi/2 by parabola

用抛物线近似 f(x) 在 0 和 pi/2 之间的积分

精确的积分可以通过对 6 个具有等距宽度的矩形的面积求和来完成。第一个矩形的高度定义为 f(x_a) ，接下来的 4 个矩形的高度定义为 f(\frac{a+b}{2}) ，最后一个矩形的高度定义为 f(x_b) 。因此，根据辛普森法则的近似积分公式可以定义为

\widetilde{ I } = \sum_{i=1}^{N+1} \frac{x_i-x_{i+1}}{6}(f(x_i)+4f(\frac{a+b}{2})+f(x_{i+1}))

辛普森积分的实现，类似于基于梯形的解决方案，采用 4 个参数，积分函数 f(x) 的范围 [a,b]，细分的数量 n，以及函数 f(x)。

double simpsonIntegral(double a, double b, int n, const std::function<double (double)> &f) {
    const double width = (b-a)/n;

    double simpson_integral = 0;
    for(int step = 0; step < n; step++) {
        const double x1 = a + step*width;
        const double x2 = a + (step+1)*width;

        simpson_integral += (x2-x1)/6.0*(f(x1) + 4.0*f(0.5*(x1+x2)) + f(x2));
    }

    return simpson_integral;
}

龙贝格积分近似

如果我们用梯形方法积分函数 f(x)，但根据上一步的步长将步长减半，我们得到近似序列

T(1)=0.18755, \> T(\frac{1}{2})=0.724727, \> T(\frac{1}{4})=0.925565

龙贝格积分的思想是引入一个 y 轴对称的抛物线，它穿过点 T(1), \> T(\frac{1}{2}) 并外推到 x \rightarrow 0 。

Romberg parable with nodes at T(1/2) and T(1/4)

龙贝格抛物线，节点位于 T(1/2) 和 T(1/4)

因此，龙贝格积分第一列 R(n,0) 的每一项都等价于梯形积分，其中第二列 R(n,1) 的每个解都等价于辛普森法则，第三列 R(n,2) 的每个解都等价于 Boole 法则。因此，龙贝格积分的公式为

R(0,0)=h_{1}(f(a)+f(b))

R(n, 0)=\frac{1}{2} R(n-1,0)+h_{n} \sum_{k=1}^{2^{\pi-1}} f\left(a+(2 k-1) h_{n}\right)

R(n, m)=R(n, m-1)+\frac{R(n, m-1)-R(n-1, m-1)}{4^{m}-1} 

std::vector<std::vector<double>> rombergIntegral
    (double a, double b, size_t n, const std::function<double (double)> &f) {
    std::vector<std::vector<double>> romberg_integral(n, std::vector<double>(n));

    //R(0,0) Start with trapezoidal integration with N = 1
    romberg_integral.front().front() = trapezoidalIntegral(a, b, 1, f);

    double h = b-a;
    for(size_t step = 1; step < n; step++) {
        h *= 0.5;

        //R(step, 0) Improve trapezoidal integration with decreasing h
        double trapezoidal_integration = 0;
        size_t stepEnd = pow(2, step - 1);
        for(size_t tzStep = 1; tzStep <= stepEnd; tzStep++) {
            const double deltaX = (2*tzStep - 1)*h;
            trapezoidal_integration += f(a + deltaX);
        }
        romberg_integral[step].front() = 0.5*romberg_integral[step - 1].front() + 
                                         trapezoidal_integration*h;

        //R(m,n) Romberg integration with R(m,1) -> Simpson rule, R(m,2) -> Boole's rule
        for(size_t rbStep = 1; rbStep <= step; rbStep++) {
            const double k = pow(4, rbStep);
            romberg_integral[step][rbStep] = (k*romberg_integral[step][rbStep-1] - 
                                             romberg_integral[step-1][rbStep-1])/(k-1);
        }
    }

    return romberg_integral;
}

Screen capture of program execution

存储库 numericalIntegration 可以在 GitHub 上找到。顾名思义，它稍后也将包含其他数值积分方法。

您喜欢这篇文章吗？

您有什么想法吗？

请随时评论和分享这篇文章。

© . All rights reserved.