c++ - C++中的复合辛普森规则

我一直在尝试使用Composite Simpson规则编写一个函数来近似积分值。

template <typename func_type>
double simp_rule(double a, double b, int n, func_type f){

    int i = 1; double area = 0;
    double n2 = n;
    double h = (b-a)/(n2-1), x=a;

    while(i <= n){

        area = area + f(x)*pow(2,i%2 + 1)*h/3;
        x+=h;
        i++;
    }
    area -= (f(a) * h/3);
    area -= (f(b) * h/3);

    return area;
    }

我要做的是用pow(2,i%2 + 1)将函数的每个值乘以2或4(和h / 3)，并减去边缘，因为这些边缘的权重应为1。

起初，我认为它工作得很好，但是，当我将其与梯形方法函数进行比较时，情况就更加不准确了，事实并非如此。

这是我以前编写的具有相同问题的代码的简化版本，我认为如果稍稍清理一下，问题就会消失，但是可惜。从另一篇文章中，我得到的想法是，我对它们进行的类型和操作正在发生某些事情，这会导致精度下降，但我只是看不到它。

编辑:

为了完整起见，我为e ^ x从1到零运行它

\\function to be approximated
double f(double x){ double a = exp(x); return a; }

int main() {

    int n = 11; //this method works best for odd values of n
    double e = exp(1);
    double exact = e-1; //value of integral of e^x from 0 to 1

    cout << simp_rule(0,1,n,f) - exact;

最佳答案

辛普森法则使用此近似值来估计定积分:

哪里

和

这样就有 n + 1 等距采样点xi。

在发布的代码中，传递给函数的参数n似乎是对函数进行采样的点数(而在上一个公式中，n是间隔数，这不是问题)。

点之间的(恒定)距离可以正确计算

double h = (b - a) / (n - 1);

由于四舍五入的原因，while循环用于求和从x = a迭代到点的所有点的加权贡献，直到其点的尾数接近b，但可能不完全是b。这意味着f的最后计算值f(x_n)可能与预期的f(b)略有不同。

但是，与以下事实相比，这没有什么好相提并论的:这些端点在循环内以 4 的起始权重求和，然后在循环后以权重 1 减去，而所有内部点都具有重量切换。实际上，这是代码计算的结果:

另外，使用

pow(2, i%2 + 1)

就效率而言，生成序列4、2、4、2，...，4是浪费的，并且可能添加(取决于实现方式)其他不必要的舍入误差。

以下算法显示了如何在不调用该库函数的情况下获得相同(固定)的结果。

template <typename func_type>
double simpson_rule(double a, double b,
                    int n, // Number of intervals
                    func_type f)
{
    double h = (b - a) / n;

    // Internal sample points, there should be n - 1 of them
    double sum_odds = 0.0;
    for (int i = 1; i < n; i += 2)
    {
        sum_odds += f(a + i * h);
    }
    double sum_evens = 0.0;
    for (int i = 2; i < n; i += 2)
    {
        sum_evens += f(a + i * h);
    }

    return (f(a) + f(b) + 2 * sum_evens + 4 * sum_odds) * h / 3;
}

请注意，此函数需要传递间隔数(例如，使用10而不是11来获得与OP函数相同的结果)，而不是要传递的点数。

可测试的here。

关于c++ - C++中的复合辛普森规则，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/60005533/