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

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

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

Формула[ред.ред. код]

Формулою Сімпсона називається інтеграл від інтерполяційного многочлена другого степеня на відрізку [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) + 4 f\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| \le \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 = \tfrac12 (b-a) (f(a)+f(b)).

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

 -\tfrac1{24} (b-a)^3 f''(a) + O((b-a)^4) \quadі\quad \tfrac1{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 h3 \cdot \left( \frac 12 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 12 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| \le \frac{(b-a)}{2880}h^4 \max\limits_{x\in[a,b]} |f^{(4)} (x)| .

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

\left| E(f) \right| \le \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;
        }
    }
}

Див. також[ред.ред. код]