Правило составного Симпсона в 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.

Сначала я подумал, что это работает просто отлично, однако, когда я сравнил это с моей функцией трапециевидного метода, это было намного более неточно, чего не должно быть.

Это более простая версия кода, который я ранее написал, которая имела ту же проблему, я думал, что, если немного ее почистить, проблема исчезнет go, но, увы. Из другого поста я понял, что с типами и операциями, которые я выполняю с ними, что-то происходит, что приводит к потере точности, но я просто не вижу этого.

Редактировать:

Для полноты я запустил его для е ^ х от 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

Приведенное выше отличное и принятое решение могло бы выиграть от свободного использования 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;
}

Answer 2

Правило Симпсона использует это приближение для оценки определенного интеграла:

Где

и

Таким образом, есть n + 1 равноотстоящие точки выборки x _i .

В размещенном коде параметр n, переданный функции, представляется как число точек, в которых выполняется выборка функции (в то время как в предыдущей формуле n - количество интервалов , это не проблема).

(Постоянное) расстояние между точками рассчитывается правильно

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

В то время как l oop используется для суммирования взвешенных вкладов всех точки повторяются от x = a до точки с асциссой, близкой к b, но, вероятно, не совсем b из-за ошибок округления. Это означает, что последнее вычисленное значение f, f(x_n), может немного отличаться от ожидаемого f(b).

Это ничто, однако, по сравнению с ошибкой, вызванной тем, что эти окончания точки суммируются внутри l oop с начальным весом 4 , а затем вычитаются после l oop с весом 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), а не количество Очки.

Тестируемый здесь .