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

Метод Сімпсона є одним із методів чисельного інтегрування. Названий на честь британського математика Томаса Сімпсона (1710—1761).

Метод Сімпсона одержується за допомогою інтерполяції функції f(x) (блакитний колір) квадратичним многочленом P(x) (червоний колір).

Формула

Формулою Сімпсона називається інтеграл від інтерполяційного многочлена другого степеня на відрізку $[a,b]$:

${\int \limits _{a}^{b}f(x)dx}\approx {\int \limits _{a}^{b}{p_{2}(x)}dx}={\frac {b-a}{6}}{\left(f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right)},$

де $f(a)$, $f((a+b)/2)$ і $f(b)$ — значення функції у відповідних точках .

Похибка

При умові, що функція $f(x)$ на відрізку $[a,b]$ має похідну четвертого порядку, похибка $E(f)$, дорівнює:

$E(f)=-{\frac {(b-a)^{5}}{2880}}{f^{(4)}(\zeta )},\ \ \ \zeta \in [a,b].$

Зважаючи, що значення $\zeta$ переважно не є відомим, для оцінки похибки використовується нерівність:

$\left|E(f)\right|\leq {\frac {(b-a)^{5}}{2880}}\max \limits _{x\in [a,b]}{\left|f^{(4)}(x)\right|}.$

Виведення формули

Формула Сімпсона може бути виведена за допомогою багатьох різних способів.

Квадратична інтерполяція

Якщо замінити функцію $f(x)$ квадратичним поліномом $P(x)$ що приймає ті ж значення що й $f(x)$ у точках a,b і m = (a+b) / 2. використавши інтерполяційну формулу Лагранжа, то одержимо формулу:

$P(x)=f(a){\frac {(x-m)(x-b)}{(a-m)(a-b)}}+f(m){\frac {(x-a)(x-b)}{(m-a)(m-b)}}+f(b){\frac {(x-a)(x-m)}{(b-a)(b-m)}}.$

Після необхідних обчислень одержуємо:

$\int _{a}^{b}P(x)\,dx={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right].$

Використання методів прямокутників і трапецій

У цьому способі виведення використовуються метод прямокутників:

$M=(b-a)f\left({\frac {a+b}{2}}\right)$

і метод трапецій:

$T={\tfrac {1}{2}}(b-a)(f(a)+f(b)).$

Похибки цих наближень дорівнюють

$-{\tfrac {1}{24}}(b-a)^{3}f''(a)+O((b-a)^{4})\quad$і$\quad {\tfrac {1}{12}}(b-a)^{3}f''(a)+O((b-a)^{4}),$

відповідно. Звідси випливає, що аби позбутися третього степеня слід взяти для наближення величину

${\frac {2M+T}{3}}.$

Однак таким чином одержується формула Сімпсона.

Метод невизначених коефіцієнтів

Запишемо в загальному виді:

${\frac {1}{b-a}}\int _{a}^{b}f(x)\,dx\approx \alpha f(a)+\beta f\left({\frac {a+b}{2}}\right)+\gamma f(b).$

Коефіцієнти α, β і γ можуть бути знайдені з вимоги, що дане наближення є точним для всіх многочленів другого степеня. Таким чином знову ж одержується метод Сімпсона.

Ітераційна формула

Для точнішого обчислення інтеграла проміжок $[a,b]$ розбивають на $N$ відрізків однакової довжини і застосовують формулу Сімпсона на кожному з них. Значення інтеграла є сумою для всіх відрізків.

${\int \limits _{a}^{b}f(x)dx}\approx {\frac {h}{3}}\cdot \left({\frac {1}{2}}f(x_{0})+\sum _{k=1}^{N-1}f(x_{k})+2\sum _{k=1}^{N}f\left({\frac {x_{k-1}+x_{k}}{2}}\right)+{\frac {1}{2}}f(x_{N})\right)$

де $h={\frac {b-a}{N}}$ величина кроку, а $x_{k}=a+k\cdot h$ межі відрізків.

Загальну похибку $E(f)$ при інтегруванні на відрізку $[a,b]$ з кроком $x_{i}-x_{i-1}=h$ визначають за формулою:

$\left|E(f)\right|\leq {\frac {(b-a)}{2880}}h^{4}\max \limits _{x\in [a,b]}|f^{(4)}(x)|$.

При неможливості оцінити похибку за допомогою четвертої похідної можна використати слабшу оцінку:

$\left|E(f)\right|\leq {\frac {(b-a)}{288}}h^{3}\max \limits _{x\in [a,b]}|f^{(3)}(x)|$.

Приклади реалізації

Реалізація на C#:

using System;

namespace NumericIntgeration
{
    internal class Program
    {
        private delegate double Func(double x);

        private static void Main()
        {
            const int n = 10000;
            double result = SimpsonMethod(0.0, 2.0, n, x => x * Math.Exp(Math.Sqrt(x)));
            Console.WriteLine("x = {0}", result);

            Console.ReadKey();
        }
        private static double SimpsonMethod(double a, double b, int n, Func func)
        {
            double h = (b - a) / n;
            double s = (func(a) + func(b)) * 0.5;
            for (int i = 1; i <= n - 1; i++)
            {
                double xk = a + h * i; //xk
                double xk1 = a + h * (i - 1); //Xk-1
                s += func(xk) + 2 * func((xk1 + xk) / 2);
            }
            var x = a + h * n; //xk
            var x1 = a + h * (n - 1); //Xk-1
            s += 2 * func((x1 + x) / 2);

            return s * h / 3.0;
        }
    }
}

Див. також

• Чисельне інтегрування
• Метод прямокутників
• Метод трапецій

Література

• Формула парабол (формула Сімпсона) // Вища математика в прикладах і задачах / Клепко В.Ю., Голець В.Л.. — 2-ге видання. — К. : Центр учбової літератури, 2009. — С. 450. — 594 с.
• Григорій Михайлович Фіхтенгольц. Курс диференціального та інтегрального числення. — 2024. — 2200+ с.(укр.)