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

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

一开始，我认为它工作得很好，但是，当我把它与梯形法函数比较时，它是非常不准确的--不应该是这样的。

这是我之前编写的代码的一个简单版本，它有同样的问题，我想如果我稍微清理一下它，问题就会消失，但是唉。从另一篇文章中，我得到了这样的想法:我正在对它们所做的类型和操作发生了一些变化，导致了精度的下降，但我没有看到它。

编辑：

为了完整起见，我将e^x从1运行到0。

\\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

上述优秀且被接受的解决方案可以从浮点类型的std::fma()和templatize的广泛使用中获益。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;
}