Интегральное приближение - правило Симпсона

Определение

Предположим, что \(f(x)\) определена на интервале \([a, b]\). Тогда правило Симпсона на всем интервале приближает определенный интеграл от \(f(x)\) на интервале по формуле

\begin{equation} \int_{a}^{b} f(x) d x \approx \frac{b-a}{6}\left(f(a)+4 f\left(\frac{a+b}{2}\right)+f(b)\right) \end{equation}

Идея состоит в том, что если \(f(x) = 1, x, x^2\), эта формула является точным равенством. Итак, правило Симпсона дает правильный интеграл от любой квадратичной функции. В общем, правило Симпсона приближает \(f(x)\) параболой через точки на графике \(f(x)\) с \(x\)-координатами \(a, \frac{a+b}{2}, b\).

Правило Симпсона обычно применяется путем разбиения интервала на \(N\) подинтервалов одинакового размера, где \(N\) - четное число, и аппроксимации интеграла по каждой паре смежных подынтервалов с использованием приведенной выше оценки.

То есть пусть \(x_0=a, x_1=a+\frac{b-a}{N}, x_2=a+2\frac{b-x}{N}, \dots, x_{N-1}=a+(N-1)\frac{b-a}{N}, x_N=b\). Тогда

\begin{equation} \int_{a}^{x_{2}} f(x) d x \approx \frac{b-a}{3 N}\left(f(a)+4 f\left(x_{1}\right)+f\left(x_{2}\right)\right) \end{equation}

\begin{equation} \int_{x_2}^{x_{4}} f(x) d x \approx \frac{b-a}{3 N}\left(f(x_2)+4 f\left(x_{3}\right)+f\left(x_{4}\right)\right) \end{equation}

и так далее. Сложение дает

\begin{equation} \int_{a}^{b} f(x) d x \approx \frac{b-a}{3 N}\left(f(a)+4 f\left(x_{1}\right)+2 f\left(x_{2}\right)+4 f\left(x_{3}\right)+2 f\left(x_{4}\right)+\cdots+4 f\left(x_{N-1}\right)+f(b)\right). \end{equation}

Использование

Представьте, что мы хотим интегрировать следующее выражение:

\begin{equation} \int_{0}^{1} (5x^3 + 2\cos{x}) dx. \end{equation}

Тогда код будет выглядеть так:

// example_integral_simpson.cpp

#include <iostream>
#include "../src/numerary.hpp"

using namespace std;
using namespace numerary;

/* Function to be integrated */
double f(double x) {
    return 5*pow(x, 3) + 2*cos(x);
}

/* The main function */
int main() {

    const double from = 0; // Lower bound of integral
    const double to = 1; // Upper bound of integral
    const string method = "simpson"; // Numerical method we will use for integration ("simpson" by default)
    const double eps = 1.e-9; // eps value for integration (1.e-9 by default)

    double *I = Numerary::integrate(f, from, to, method, eps);

    cout << "ans = " << I[0] << endl; // Value of calculated integral
    cout << "err = " << I[1] << endl; // Error of calculated integral value

    return 0;
}