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;
}
Поделиться публикацией
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 12

    +7
    Девчонки и мальчишки, а также их родители! Пожалуйста, помните, что в реальной жизни std::function надо передавать по значению, векторы — по ссылке, а для организации двумерных массивов есть более эффективные средства, чем вектор векторов.
      0

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

        +3

        Вам в реальности не нужно будет в рантайме менять длины некоторых строк/столбцов. Самое простое, но намного более эффективное решение это std::array из std::array. Никаких лишних аллокаций, всё в памяти одним блоком.

          0

          Ну а динамических многомерных массивов в стандартной библиотеке C++, к сожалению, до сих пор нет.

            0

            Мне очень сложно представить ситуацию, когда они бы понадобились на практике. А если и понадобятся в узкоспецифическом проекте, то их не очень сложно реализовать самому.

            +1

            Изменять в рантайме размеры может и не нужно, но выделять, к примеру, матрицу размера, определённого в compile time, выглядит не всегда разумно. Вообще, как мне кажется, не нужно делать преждевременные оптимизации в ущерб читаемости кода, так как может оказаться, что проседания производительности не в лишних аллокациях. Да не спорю, круто писать сразу оптимальный во всех смыслах код, но кажется, что это может оказаться сильно дороже, чем хотелось бы.

            0
            Дать исчерпывающее перечисление всё равно не получится. Это может быть простой одномерный вектор вкупе с ручным вычислением эффективного индекса, std::array, самописный класс, обёртка над стандартным контейнером, что-то из контейнеров Boost и многое другое. И выбор зависит от немалого количества факторов. Применительно к данной статье мог бы подойти «голый» одномерный вектор или самописная обёртка, хотя первое несколько усложнит код по сравнению с «идеальным» двумерным массивом (вектор векторов в этом смысле тоже не подарочек), а второе потребует как минимум привести интерфейс и пояснения.
            +1

            Когда речь заходит об оптимизации вычислений, то передавать вообще надо не std::function, а шаблонный функтор. Ну и алгоритмические моменты: нет никакого смысла вычислять значение функции в одной и той же точке дважды.

              –1
              По алгоритмическим моментам поддержу. И повторные вычисления желательно сокращать, и возведение в квадрат и четвёртую степень лучше делать через умножение, а не дорогостоящий pow.
              А вот накладные расходы на вызов std::function будут заметны на общем фоне только для совсем уж простых вычислений. Стоит ли жертвовать ради этого выразительностью и универсальностью, которые даёт function? Зависит от обстоятельств, но в учебном проекте вряд ли стоит сражаться за считанные такты ценой снижения ясности кода.
            0

            Графики чем строили?

              0
              Что-то все больше таких статей появляется. И дело даже не в том, каково качество реализации, сколько в том, что не раскрывается тема, обозначенная в заголовке.

              Говоря о формулах Ньютона-Котеса, либо же, как их еще называют, методы Ньютона-Котеса, необходимо понимать, что это семейство методов подразумевает “аппроксимацию подинтегральной функции на выбранных промежутках многочленами”. Вот только автор как-то опускает аспект о многочленах (Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков. Численные методы. Глава 3, параграф 3.)

              Во всех методах (которые, бесспорно, являются методами Ньютона-Котеса) опущен принципиальный для численных методов вопрос — оценка погрешности вычислений.

              Здесь же можно упомянуть много, с точки зрения вычислительных методов, неточностей, которые, тем не менее, являются существенными для изложения такого рода методов.
              В качестве примера:
              “Точность приближения зависит от числа N отрезков, на которые разбивается промежуток интегрирования. Таким образом, длина промежутка: dx = (b-a)/N” (не понял, как можно вставить в коммент верстку LaTeX, простите).
              Возникают вопросы о:
              — «промежуток интегрирования»
              — «длина промежутка»
              — «отрезок интегрирования» (упомянуто немного выше)
              и как связаны между собой эти понятия.

              Далее. В разделе о методе Симпсона излагается следующее: “Метод Симпсона заключается в интегрировании интерполяционного многочлена второй степени функции f(x) с узлами интерполяции a, b и m = (a+b)/2 — параболы p(x).”
              Узлы интерполяции это a, b и m = (a + b) / 2?

              “Площадь параболы может быть найдена суммированием площадей 6 прямоугольников равной ширины. Высота первого из них должна быть равна f(a), с третьего по пятый — f(m), шестого — f(m)”
              А какие это прямоугольники? Где они обозначены на вашем рисунке?
              Формула, которую вы используете, “слегка” неправильна. Полагаю, что вы взяли формулу из этой статьи ttps://ru.wikipedia.org/wiki/%D0%A4%D0%BE%D1%80%D0%BC%D1%83%D0%BB%D0%B0_%D0%A1%D0%B8%D0%BC%D0%BF%D1%81%D0%BE%D0%BD%D0%B0 но не разобрались до конца.

              О применении метода Ромберга здесь… Хм…
              В качестве примера — “а каждое решение второго столбца R (n, 1)”. Какие столбцы? Вы это откуда-то скопировали?
              На этом хватит примеров.
              Общее впечатление от этой статьи — очень поверхностно, местами нелогично. Как пример на эту же тему могу привести следующий пост habr.com/ru/post/420867.

              Могу рекомендовать к изучению книгу Н.Н. Калиткина “Численные методы”. Думаю, что вы ее легко отыщете

                0
                У Вас в 4 формуле написано x[i]-x[i+1]. Это точно законно?

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

                Самое читаемое