C++ и Численные Методы: Приближенное интегрование по Ньютону-Котесу

Моя цель - предложение широкого ассортимента товаров и услуг на постоянно высоком качестве обслуживания по самым выгодным ценам.
Методы Ньютона-Котеса — это совокупность техник приближенного интегрирования, основанных на:

  • разбиении отрезка интегрирования на равные промежутки;
  • аппроксимации подинтегральной функции на выбранных промежутках многочленами;
  • нахождении суммарной площади полученных криволинейных трапеций.

В этой статье будут рассмотрены несколько методов Ньютона-Котеса:

  • метод трапеций;
  • метод Симпсона;
  • метод Ромберга.

Метод трапеций


Метод трапеций — простейший из рассмотренных. В качестве примера возьмем следующий интеграл:

image

Точность приближения зависит от числа N отрезков, на которые разбивается промежуток интегрирования. Таким образом, длина промежутка:

image

image

Площадь трапеции может быть вычислена по формуле:

image

Суммируя все вышесказанное, приближенное значение интеграла вычисляется по формуле:

image

Функция, вычисляющая интеграл методом трапеций должна принимать 4 параметра:

  • границы отрезка интегрирования;
  • подинтегральную функцию;
  • число N промежутков разбиения.

double trapezoidalIntegral(double a, double b, int n, const std::function<double (double)> &f) {
    const double width = (b-a)/n;

    double trapezoidal_integral = 0;
    for(int step = 0; step < n; step++) {
        const double x1 = a + step*width;
        const double x2 = a + (step+1)*width;

        trapezoidal_integral += 0.5*(x2-x1)*(f(x1) + f(x2));
    }

    return trapezoidal_integral;
}

Метод Симпсона


Метод Симпсона заключается в интегрировании интерполяционного многочлена второй степени функции f(x) с узлами интерполяции a, b и m = (a+b)/2 — параболы p(x).Для повышения точности имеет смысл разбить отрезок интегрирования на N равных промежутков(по аналогии с методом трапеций), на каждом из которых применить метод Симпсона.

image

Площадь параболы может быть найдена суммированием площадей 6 прямоугольников равной ширины. Высота первого из них должна быть равна f(a), с третьего по пятый — f(m), шестого — f(m). Таким образом, приближение методом Симпсона находим по формуле:

image

double simpsonIntegral(double a, double b, int n, const std::function<double (double)> &f) {
    const double width = (b-a)/n;

    double simpson_integral = 0;
    for(int step = 0; step < n; step++) {
        const double x1 = a + step*width;
        const double x2 = a + (step+1)*width;

        simpson_integral += (x2-x1)/6.0*(f(x1) + 4.0*f(0.5*(x1+x2)) + f(x2));
    }

    return simpson_integral;
}

Метод Ромберга


Пусть T(x) — приближение интеграла, полученное методом трапеций с шагом x. Получим 3 таких приближения, уменьшая размер шага в 2 раза при каждом вычислении.

image

Построим теперь симметричную относительно оси y параболу, проходящую через точки T(1) и T(1/2) чтоб экстраполировать полученные значения для x стремящегося к 0.

image

Следовательно, каждый член первого столбца R(n, 0) приближений Ромберга эквивалентен решениям полученным методом трапеций, а каждое решение второго столбца R (n, 1) — методом Симпсона. Таким образом, формулы для приближенного интегрирования методом Ромберга:

image

image
image

Реализация на C++:

std::vector<std::vector<double>> rombergIntegral(double a, double b, size_t n, const std::function<double (double)> &f) {
    std::vector<std::vector<double>> romberg_integral(n, std::vector<double>(n));

    romberg_integral.front().front() = trapezoidalIntegral(a, b, 1, f);

    double h = b-a;
    for(size_t step = 1; step < n; step++) {
        h *= 0.5;

        double trapezoidal_integration = 0;
        size_t stepEnd = pow(2, step - 1);
        for(size_t tzStep = 1; tzStep <= stepEnd; tzStep++) {
            const double deltaX = (2*tzStep - 1)*h;
            trapezoidal_integration += f(a + deltaX);
        }
        romberg_integral[step].front() = 0.5*romberg_integral[step - 1].front() + trapezoidal_integration*h;

        for(size_t rbStep = 1; rbStep <= step; rbStep++) {
            const double k = pow(4, rbStep);
            romberg_integral[step][rbStep] = (k*romberg_integral[step][rbStep-1] - romberg_integral[step-1][rbStep-1])/(k-1);
        }
    }

    return romberg_integral;
}
Источник: https://habr.com/ru/post/479202/


Интересные статьи

Интересные статьи

Кто бы что ни говорил, но я считаю, что изобретение велосипедов — штука полезная. Использование готовых библиотек и фреймворков, конечно, хорошо, но порой стоит их отложить и создать ...
Предыстория Когда-то у меня возникла необходимость проверять наличие неотправленных сообщений в «1С-Битрикс: Управление сайтом» (далее Битрикс) и получать уведомления об этом. Пробле...
Если Вы используете в своих проектах инфоблоки 2.0 и таблицы InnoDB, то есть шанс в один прекрасный момент столкнуться с ошибкой MySQL «SQL Error (1118): Row size too large. The maximum row si...
Реализация ORM в ядре D7 — очередная интересная, перспективная, но как обычно плохо документированная разработка от 1с-Битрикс :) Призвана она абстрагировать разработчика от механики работы с табл...
Один из самых острых вопросов при разработке на Битрикс - это миграции базы данных. Какие же способы облегчить эту задачу есть на данный момент?