Інтегрування за формулою Сімпсона

Ми збираємося обчислити значення визначеного інтеграла

$$\int_a ^ b f (x) dx$$

Описаний тут розв'язок було опубліковано в одній із дисертацій Томаса Сімпсона у 1743 році.

Коли підходить цей алгоритм?

• Підінтегральну функцію $f(x)$ важко або неможливо проінтегрувати аналітично, тож потрібне чисельне наближення? (якщо первісну легко взяти аналітично — інтегруйте напряму, без чисельного методу)
• Функцію можна обчислити в довільній точці відрізка $[a, b]$, а сам $f$ достатньо гладкий (мала четверта похідна)?
• Шукаєте площу / визначений інтеграл, а не корінь чи екстремум функції? (для кореня → Метод Ньютона, для екстремуму → Тернарний пошук)

Формула Сімпсона

Нехай $n$ — деяке натуральне число. Розіб'ємо відрізок інтегрування $[a, b]$ на $2n$ рівних частин:

$$x_i = a + i h, ~~ i = 0 \ldots 2n,$$ $$h = \frac {b-a} {2n}.$$

Тепер обчислимо інтеграл окремо на кожному з відрізків $[x_ {2i-2}, x_ {2i}]$, $i = 1 \ldots n$, а потім додамо всі значення.

Отже, припустимо, що ми розглядаємо черговий відрізок $[x_ {2i-2}, x_ {2i}], i = 1 \ldots n$. Замінимо на ньому функцію $f(x)$ параболою $P(x)$, що проходить через 3 точки $(x_ {2i-2}, x_ {2i-1}, x_ {2i})$. Така парабола завжди існує і єдина; її можна знайти аналітично. Наприклад, ми могли б побудувати її за допомогою інтерполяційного многочлена Лагранжа. Єдине, що залишається зробити, — проінтегрувати цей многочлен. Якщо зробити це для загальної функції $f$, ми отримаємо напрочуд простий вираз:

$$\int_{x_ {2i-2}} ^ {x_ {2i}} f (x) ~dx \approx \int_{x_ {2i-2}} ^ {x_ {2i}} P (x) ~dx = \left(f(x_{2i-2}) + 4f(x_{2i-1})+(f(x_{2i})\right)\frac {h} {3}$$

Додавши ці значення по всіх відрізках, ми отримуємо остаточну формулу Сімпсона:

$$\int_a ^ b f (x) dx \approx \left(f (x_0) + 4 f (x_1) + 2 f (x_2) + 4f(x_3) + 2 f(x_4) + \ldots + 4 f(x_{2N-1}) + f(x_{2N}) \right)\frac {h} {3}$$

Похибка

Похибка наближення інтеграла за формулою Сімпсона дорівнює

$$-\tfrac{1}{90} \left(\tfrac{b-a}{2}\right)^5 f^{(4)}(\xi)$$

де $\xi$ — деяке число між $a$ і $b$.

Похибка асимптотично пропорційна $(b-a)^5$. Однак наведені вище викладки наводять на думку про похибку, пропорційну $(b-a)^4$. Формула Сімпсона набуває додаткового порядку, оскільки точки, у яких обчислюється підінтегральна функція, розподілені симетрично на відрізку $[a, b]$.

Реалізація

Тут $f(x)$ — деяка визначена користувачем функція.

C++

const int N = 1000 * 1000; // кількість кроків (уже помножена на 2)

double simpson_integration(double a, double b){
    double h = (b - a) / N;
    double s = f(a) + f(b); // a = x_0 і b = x_2n
    for (int i = 1; i <= N - 1; ++i) { // дивіться остаточну формулу Сімпсона
        double x = a + h * i;
        s += f(x) * ((i & 1) ? 4 : 2);
    }
    s *= h / 3;
    return s;
}

Python

N = 1000 * 1000  # кількість кроків (уже помножена на 2)

def simpson_integration(a: float, b: float) -> float:
    h = (b - a) / N
    s = f(a) + f(b)  # a = x_0 і b = x_2n
    for i in range(1, N):  # дивіться остаточну формулу Сімпсона
        x = a + h * i
        s += f(x) * (4 if i & 1 else 2)
    s *= h / 3
    return s

TypeScript

const N = 1000 * 1000; // кількість кроків (уже помножена на 2)

function simpsonIntegration(a: number, b: number): number {
    const h = (b - a) / N;
    let s = f(a) + f(b); // a = x_0 і b = x_2n
    for (let i = 1; i <= N - 1; ++i) { // дивіться остаточну формулу Сімпсона
        const x = a + h * i;
        s += f(x) * ((i & 1) ? 4 : 2);
    }
    s *= h / 3;
    return s;
}

Go

const N = 1000 * 1000 // кількість кроків (уже помножена на 2)

func simpsonIntegration(a, b float64) float64 {
    h := (b - a) / N
    s := f(a) + f(b) // a = x_0 і b = x_2n
    for i := 1; i <= N-1; i++ { // дивіться остаточну формулу Сімпсона
        x := a + h*float64(i)
        if i&1 == 1 {
            s += f(x) * 4
        } else {
            s += f(x) * 2
        }
    }
    s *= h / 3
    return s
}

Задачі для практики

• Latin American Regionals 2012 - Environment Protection

Відеоматеріали