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

Для підвищення точності замінимо функцію f(х) на відрізку [x_i-1,x_i] інтерполяційним поліномом Лагранжа другого порядку. Припустимо, що четверта похідна від функції f обмежена на відрізку [а,b]: , при всіх х, що належать відрізку [а,b].

Будемо обчислювати інтеграл: , . Позначимо

, де , R(x) – похибка інтерполяційного многочлена

, де точка [x₁,x₃].

.

Обчислимо . Тут ,

,

У свою чергу,

Інтеграли обчислюються заміною зміної. Вони обидва рівні .

обчислюється так(заміна зміної ):

обчислюється так(заміна зміної ):

Таким чином, квадратурна формула Симпсона може бути записана так: .

Похибка формули Симпсона(без виводу).

Нагадаємо позначення. Нехай ,R(x) – похибка інтерполяційного многочлена Лагранжа другого порядку, , точка [x₁,x₃] для . Позначимо .

Можна довести, що якщо четверта похідна від функції f(x) неперервна на відрізку [x₁,x₃], то існує точка t відрізка [x₁,x₃] така, то ,

Формула Симпсона точна для поліномів до третього степеня включно, формула трапецій - тільки для поліномів першого степеня.

Складена формула Симпсона

Якщо застосувати цю формулу до всього відрізка [а,b], при умові що відрізок розбитий на N однакових частин, h=(b-а) /N то одержимо складену формулу Симпсона:

Похибка складеної формули Симпсона P не перевищує суми похибок по кожному відрізку:

Останнє витікає з того, що h=(b-а) /N, а Nh=b-a.

Таким чином, складена формула Симпсона має четвертий порядок точності. Звідси видно, що формула Симпсона істотно точніше, ніж формула трапецій. На частковому відрізку вона має точність O(h⁵), а на всьому відрізку - O(h⁴).

Апостеріорна оцінка похибки методом Рунге-Ромберга.

Величина похибки чисельного інтегрування залежить як від кроку сітки h, так і від властивостей функції f(x). У оцінку погрішності, разом з h, входить величина максимуму другої або четвертої похідної M, яка може мінятися від точки до точки і, взагалі кажучи, наперед невідома. Якщо величина похибки велика, то її можна зменшити шляхом подрібнення сітки на даному відрізку [x_i-1, x_i]. Для цього перш за все треба уміти апостеріорно, тобто після проведення розрахунку, оцінювати похибку чисельного інтегрування. Апостеріорну оцінку похибки можна здійснити методом(правилом) Рунге-Ромберга, який розглянемо спочатку на прикладі формули трапецій. Нехай відрізок

[а, b] розбитий на часткові відрізки [x_i-1, x_i], i=1, 2, …, N, x₀=a, x_N=b. Відомо, що I=I_h+Сh²+O(h³), де I- точне значення інтеграла, I_h - значення інтеграла обчислене по складених наближених формулах на рівномірній сітці з кроком h, С - деяка невідома константа. Вважатимемо, що h мале і прямує до нуля. Знак рівності треба розуміти в тому значенні, що решта членів розкладу інтеграла має вищий порядок по h: h³,h⁴ . Якщо h мале, то вони не роблять практично впливу на похибку. Тому можна вважати майже точною рівність: I=I_h+Сh² . Тепер обчислимо наближено той же інтеграл з кроком h/2 по формулі трапецій: I=I_h/2+Сh²/4. Виключимо константу С з рівностей: I-I_h=Ch², I-I_h/2=Ch²/4. Звідси I-I_h=4(I-I_h/2), I=I_h/2+(I_h/2-I_h)/3.

Останню рівність треба розуміти так: наближення I_h/2+(I_h/2-I_h)/3 набагато точніше обчислює інтеграл, ніж I_h/2 і тим більше I_h. Слід звернути увагу, що точне значення С нам не знадобилося! Ці міркування носять наближений характер, бо були відкинуті доданки з h³,h⁴ ітд. Вони пропорційні третій, четвертій похідним.

Прийом Рунге- Ромберга носить евристичний характер.

Повторюючи попередні міркування для методу Симпсона також можна застосувати метод Рунге-Ромберга: I=I_h+Сh⁴+O(h⁶) , де I - точне значення інтеграла, I_h - значення інтеграла обчислене по складених наближених формулах на рівномірній сітці з кроком h, С - деяка константа. Тут в розкладанні вказано степінь h⁶, оскільки формула Симпсона симетрична щодо середини відрізка інтегрування, коєфіціент при h⁵ рівний нулю. Обчислимо той же інтеграл з кроком h/2: I=I_h/2+Сh⁴/16. Виключимо константу С з рівностей: I-I_h=Ch⁴, I-I_h/2=Ch⁴/16. Звідси I-I_h=15(I-I_h/2), I=I_h/2+(I_h/2-I_h)/15. Останню рівність треба розуміти так: наближення I_h/2+(I_h/2-I_h)/15 набагато точніше обчислює інтеграл, ніж I_h/2 і тим більше I_h . Константа у формулі Рунге-Ромберга 15=2⁴-1 виходить так: це два в степені, рівній точності складеної формули чисельного інтегрування, мінус один. Для формули трапецій вона була рівна 3=2²-1.

Блок-схема дозволяє обчислити інтеграл з заданою точністю за допомогою метода Рунге-Ромберга. Задається функція f(x), межі інтегрування а,b і точність обчислення інтеграла eps. У підпрограмі Simps реалізоване обчислення інтеграла з заданим кроком. Якщо поправка (I_h/2-I_h)/15<eps, то обчислення припиняється, інакше збільшуємо n в два рази і обчислюємо інтеграл повторно.

Для підтвердження вище викладеного розглянемо приклад : обчислити інтеграл .

Як відомо, точне значення I=2. В результаті розрахунку по приведеному алгоритму були одержані результати, приведені в таблиці. У першій колонці n- число поділів заданого відрізка, друга колонка Z2 - значення інтеграла I_h/2, в третій колонці - різниця між точним значенням інтеграла і I_h/2 . У останній колонці - різниця між точним значенням інтеграла і

I_h/2 +( I_h/2 –I_h)/15 . Видно, що точність обчислень за методом Рунге-Ромберга набагато вища, майже в 10⁵ раз в кінцевому результаті.

n Z2 Int-Z2 Int-Z2-(Z2-Z1) /15

2 2.00455975498E+0000 -4.55975498442E-0003 1.42926817616E-0003

4 2.00026916995E+0000 -2.69169948388E-0004 1.68690540144E-0005

8 2.00001659105E+0000 -1.65910479355E-0005 2.47545427969E-0007

16 2.00000103337E+0000 -1.03336941300E-0006 3.80915516729E-0009

32 2.00000006453E+0000 -6.45300019231E-0008 5.92921486648E-0011

64 2.00000000403E+0000 -4.03225740623E-0009 9.25561602187E-0013

128 2.00000000025E+0000 -2.52002534075E-0010 1.44574022892E-0014

256 2.00000000002E+0000 -1.57499503890E-0011 2.21827764491E-0016