C++ 中的复合辛普森法则

Question

我一直在尝试编写一个函数来使用复合辛普森法则来近似计算积分值。

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;
    }

我所做的是将函数的每个值乘以 2 或 4（和 h/3）和 pow(2,i%2 + 1)，然后减去边缘，因为它们的权重应该仅为 1。

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

这是我以前写的代码的简化版本，它有同样的问题，我认为如果我稍微清理一下问题就会消失，但是唉。从另一个 post，我了解到类型和我正在对它们执行的操作会导致精度损失，但我只是没有看到它。

编辑：

为了完整起见，我是运行 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;

Answer 1

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

$\int_{a}^{b}{f(x)} \approx \frac{\Delta x}{3}\left ( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + ... + 4f(x_{n-1}) + f(x_n) \right )$

在哪里

$\Delta x = \frac{b - a}{n}$

和

$x_i = a + i \cdot \Delta x$

所以有n + 1个等距样本点x_i.

在发布的代码中，传递给函数的参数n似乎是函数采样的点数（而在前面的公式中n是间隔数，这不是问题。

正确计算点之间的（恒定）距离

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

用于对所有点的加权贡献求和的while循环从x = a迭代到横坐标接近b的点，但由于舍入误差，可能不完全 b。这意味着 f、f(x_n) 的最后计算值可能与预期的 f(b).

略有不同

不过，与那些端点在循环内以 4 的起始权重相加然后在循环后减去这一事实所导致的错误相比，这不算什么权重为 1，而所有内部点的权重都已切换。事实上，这就是代码的计算结果：

$\frac{\Delta x}{3}\left ( 3f(x_0)+ 2f(x_1) + 4f(x_2) + ... + 2f(x_{n-1}) + 3f(x_{n}) \right )$

另外，使用

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。

Answer 2

上述出色且公认的解决方案可以受益于 std::fma() 的自由使用和浮点类型的模板化。 https://en.cppreference.com/w/cpp/numeric/math/fma

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

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

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

C++ 中的复合辛普森法则

Composite Simpson's Rule in C++

c++

numerical-integration

simpsons-rule