1

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

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。

起初，我认为它工作得很好，但是，当我将它与我的梯形方法函数进行比较时，它更加不准确，这不应该是这种情况。

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

编辑：

为了完整起见，我将 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;
4

2 回答 2

2

辛普森法则使用这个近似来估计一个定积分：

在哪里

和

这样就有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}) \正确的 ）

此外，使用

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 函数相同的结果），而不是点数。

可在这里测试。

于 2020-01-31T16:55:54.530 回答
1

上述优秀且被接受的解决方案可以受益于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;
}
于 2020-04-07T17:44:47.157 回答