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

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

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

Формула

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

$\underset{a}{\overset{b}{\int }}f(x)dx\approx \underset{a}{\overset{b}{\int }}{p}_2(x)dx=\frac{b-a}{6}(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)),$

де $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$ переважно не є відомим, для оцінки похибки використовується нерівність:

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

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

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

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

Якщо замінити функцію $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)}.$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$\underset{a}{\overset{b}{\int }}f(x)dx\approx \frac{h}{3}\cdot (\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))$

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

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

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

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

$|E(f)|\le \frac{(b-a)}{288}{h}^3\max {\mathrm{\backslash 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;
        }
    }
}

Див. також

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

Література

• Шаблон:Клепко ВМ
• Шаблон:Фіхтенгольц.укр