C++ 中复合辛普森规则中的错误答案

【问题标题】：Errors in Composite Simpson's Rule in C++C++ 中复合辛普森规则中的错误
【发布时间】：2016-05-26 21:25:48
【问题描述】：

我正在编写一个小程序来使用复合辛普森规则在 C++ 中逼近积分。它似乎有效，但是当我计算不同数量的间隔的积分时（因此是for 循环），出现了一些错误，因为它似乎完全是随机的。这是我的代码。

#include <iostream>
#include <cmath>

double f(double x)
{
    double y = x;
    return y;
}

int main()
{
    for(double n=1; n<20; n+=1)
    {
        double a = 1;
        double b = 4;
        double d = (b-a)/n;
        double i = a;
        double area = 0;
        while(i<b)
        {
            double j = i + d;
            area += (d/6) * ( f(i) + 4*f((i+j)/2) + f(j) );
            i = i + d;
        }
        std::cout << "Order " << n << " : " << area << std::endl;
    }
    return 0;
}

输出如下。

Order 1 : 7.5
Order 2 : 7.5
Order 3 : 7.5
Order 4 : 7.5
Order 5 : 7.5
Order 6 : 7.5
Order 7 : 9.30612
Order 8 : 7.5
Order 9 : 7.5
Order 10 : 8.745
Order 11 : 8.6281
Order 12 : 7.5
Order 13 : 7.5
Order 14 : 7.5
Order 15 : 7.5
Order 16 : 7.5
Order 17 : 8.22145
Order 18 : 8.18056
Order 19 : 7.5

我认为这可能与变量类型或存储问题有关，尽管我无法弄清楚问题所在。任何帮助表示赞赏。

【问题讨论】：

您应该从不使用double 作为for 循环的索引，因为舍入错误。 while 也一样。一般来说，你应该测试不等式直到某个有限的 epsilon，即std::abs(x - a) < eps。表示x 与a 足够接近，可以被视为相等。
我懒得自己算算，结果应该是7.5吗？

标签： c++ numerical-integration

【解决方案1】：

这可能是浮点数比较引起的。 4.0<4.0 在这种情况下对于计算机来说可能是正确的，这会使您的循环重复多次。

要解决这个问题，您可以使用 epsilon 来比较数字：

#define EPS 1e-8
if(a<b-EPS){...} // this is analogous to < operator for ints
if(a<b+EPS){...} // this is analogous to <= operator for ints

另一种选择是动态生成i 变量，同时使用整数来控制循环：

for(int k=0;k<n;k++){
    double i=a+(b-a)*k/n;
    // ...
}

【讨论】：