Схеми Рунге-Кутта четвертого порядку

Методом Рунге-Кутта можна будувати схеми різного порядку точності.

Так схема ламаних Ейлера (8.5) є схемою Рунне-Kyттa першого порядку точності. Найбільш уживані схеми четвертого порядку точності. Наведемо без виведення найчастіше використовувані з них:

,

,

                                 (8.15)

Відзначимо, що формули більш високого порядку точності практично не вживаються через громіздкість, що зростає значно швидше, ніж точність формули.

Схеми Рунге-Кутта мають ряд переваг:

1) мають досить високий ступінь точності (за винятком схеми ламаних);

2) є явними, тобто значення  обчислюється за раніше знайденими значеннями;

3) допускають використання змінного кроку, що дає можливість зменшити його там, де функція швидко змінюється, і збільшити в іншому випадку;

4) є легко застосовними, тому що для початку розрахунку досить вибрати сітку  і задати значення .

Зазначені властивості схем досить корисні при розрахунках на ЕОМ.

Оцінки похибок різних схем Рунге-Кутта пов'язані з максимумами модулів відповідних похідних функції  досить складними формулами типу (8.8). У зв'язку з цим при розв’язанні конкретної задачі виникає питання, якою з формул Рунге-Кутта доцільно користуватися і як обирати крок сітки.

Якщо f(х,у) неперервна й обмежена разом зі своїми четвертими похідними, то гарні результати дає схема четвертого порядку (8.15). Якщо права частина рівняння (8.1) не має зазначених похідних, то граничний порядок точності схеми (8.15) не може бути реалізований. Тоді доцільно користуватися схемами меншого порядку точності, що дорівнює порядкові наявних похідних, наприклад для двічі неперервно диференційованої функції f(x,y) —схемами (8.12) і (8.1З).

Крок сітки варто вибирати настільки малим, щоб забезпечити необхідну точність розрахунку. З огляду на складність виразів залишкових членів (типу 8.8) апріорною (від лат. a priori — з попередніх) оцінкою точності для вибору кроку при практичних розрахунках не користуються, а заміняють її, наприклад, розрахунками зі згущенням сітки і дають апостеріорну (від лат. a posteriori - з наступного) оцінку точності.

Зауваження. Якщо функція f(x,y) досить гладка, але швидко змінюється на , схеми Рунге-Кутта як низького, так і високого порядку точності вимагають неприйнятно малого кроку для отримання задовільного результату. Для таких задач використовуються спеціальні методи, орієнтовані на даний вузький клас задач.

Одним із найбільш простих, широко застосовуваних і досить ефективних методів оцінки похибки й уточнення отриманих результатів у наближених обчисленнях з використанням сіток є правило Рунге.

Нехай маємо наближену формулу  для обчислення величини у(х) за значеннями на рівномірній сітці  і залишковий член цієї формули має такий вигляд:

.       (8.16)

Виконаємо тепер розрахунок за тією самою наближеною формулою для тієї самої точки х, але використовуючи рівномірну сітку з іншим кроком rh, r < 1. Тоді отримане значення  пов'язане з точним значенням співвідношенням

.  (8.17)

Зауважимо, що .

Маючи два розрахунки на різних сітках, неважко оцінити величину похибки. Для цього віднімемо (8.17) з (8.16) і одержимо першу формулу Рунге

(8.18)

Перший з доданків є головним членом похибки. Таким чином, розрахунок за другою сіткою дозволяють оцінити похибку розрахунку за першою (з точністю до членів більш високого порядку).

Відзначимо, що при користуванні правилом Рунге практично досить застосувати формулу оцінки похибки у вигляді

,               (8.19)

де  - наближені значення розв’язку рівняння в одній і тій самій точці, отримані з кроком h і h/2. При цьому необхідна точність може вважатися досягнутою, якщо величина R не перевищує заданої похибки у всіх збіжних вузлах.

Застосування правила Рунге дає можливість одержувати результати досить високої точності, використовуючи обчислення за формулами низького порядку точності, тобто дозволяє уточнити результати.

Відзначимо, що правило Рунге застосовується й у випадках, якщо сітки з різним числом вузлів нерівномірні, але їх можна описати функціями h(x), відношення яких . Величина відношення кроків r у правилі Рунге може бути будь-якою, але використовується вона найчастіше для цілого r, при цьому всі вузли менш докладної сітки повинні бути вузлами більш докладної. Особливо зручно згущати сітки вдвічі. У цьому випадку у вузлах, що є загальними для декількох сіток, можна уточнювати y(x) безпосередньо за правилом Рунге (8.18). Якщо розв’язок не уточнюється, а лише оцінюється його похибка, то використовується формула (8.19). Можна уточнити значення у(х,h) у всіх вузлах найдетальнішої сітки.

Використовуємо збіжніі вузли сіток для визначення виправлень до значень фікції :

.

Значення виправлень в інших вузлах знайдемо інтерполяцією. Для рівномірних сіток покладемо . Потім обчислимо уточнені значення . Цей спосіб узагальнюється на будь-яке число сіток.

Відзначимо ще раз, що виконати уточнення, як правило, простіше, ніж скласти і використовувати схему високого порядку точності.

Похибки наведених схем Рунге-Кутта визначаються максимальними значеннями відповідних похідних.

Оцінку похибки легко одержати для окремого випадку правої частини диференціального рівняння  .

У цьому випадку розв’язок рівняння може бути зведений до квадратури й усі схеми різницевого розв’язку переходять у формули чисельного інтегрування. Наприклад, схема (8.14) набирає вигляду

 ,

тобто має вигляд формули трапецій, а схема (8.15) переходить у схему

,

що являє собою формулу Симпсона з кроком .

Мажорантні оцінки похибок формул трапецій і Симпсона відомі. З них видно, що точність схем Рунге-Кутта досить висока.

Приклад. Знайти наближений розв’язок задачі Коші для звичайного диференціального рівняння (ОДУ) 1 порядку

  у'(t)=2ty, t 0 = 0, T=1, y( 0) =1 та оцінити його похибку.

Вихідні дані:

Права частина: . Початкове значення: .

Кінці відрізка: .Крок сітки: .

Число вузлів сітки: .Функція, що реалізує явний метод Ейлера, повертає вектор розв’язку:

 

 

Вхідні параметри:

f - функція правої частини;

y0 - початкове значення;

t0 - початкова точка відрізка;

h - крок сітки;

N - число вузлів сітки.

Обчислення розв’язку за методом Ейлера:

Обчислення розв’язку за методом Рунге-Кутта 4 порядку точності:

;

вхідні параметри:

y - вектор початкових значень;

t0 - початкова точка відрізка;

T - кінцева точка відрізка;

N - число вузлів сітки;

f - функція правої частини. Функція rkfixed повертає матрицю, перший стовпець якої містить вузли сітки, а другий - наближений розв’язок у цих вузлах.

Точний розв’язок: .

Точний розв’язок у вузлах сітки:

 

Розв’язок за методом Розв’язок за методом Точний розв’язок

    Ейлера                    Рунге-Кутта  

Графіки наближених і точних розв’язків

Обчислення похибки за правилом Рунге:

Обчислення наближених розв’язків із кроком h /2:

Обчислення похибок:

.

Приклад реалізації алгоритму методу Рунге-Кутта четвертого порядку з заданою точністю  на псевдокоді.

f(x,y):

//повертає значення заданої похідної при заданих x та y

end

//BeginValue – початкові умови

//YF – відповідь – масив значень функції

//h – крок

//n – кількість точок розбиття

//eps – точність розрахунку

//Метод Рунге-Кутта 4-го порядку

SolveRungeKutt(BeginValue,YF,h,n,eps):

1 repeat

2 for j:=1 to 2 do

3 Y[j,1]:=BeginValue

4 for i:=1 to n do

5 x:=a+(i-1)*h

6 k1:=f(x,Y[j,i])

7 k2:=f((x+h/2),(Y[j,i]+h*k1/2))

8 k3:=f((x+h/2),(Y[j,i]+h*k2/2))

9 k4:=f((x+h),(Y[j,i]+h*k3))

10 Y[j,i+1]:=Y[j,i]+h*(k1+2*k2+2*k3+k4)/6

11 done

12 if j=1 then

13 h:=h/2

14 n:=round((b-a)/h)

15 fi

16 done

17 maxr:=0

18 for i:=1 to 10 do

19 r:=abs(Y[1][((n*i) div 20)+1]-

20 -Y[2][((n*i) div 10)+1])/15;

21 if r>maxr then

22 maxr:=r

23 fi

24 done

25 if maxr<eps then

26 for i:=1 to 10 do

27 YF[i]:=Y[2][((n*i) div 10)+1]

28 done

29 fi

30 until maxr<eps

End

Багатокрокові методи

Метод прогнозу і корекції

Підправивши схему Ейлера в середній точці одержимо схему прогнозу

                            ,                       (8.20)

де наближене значення. Використовувати формулу (8.20) не можна через те, що схема прогнозу нестійка. З цієї причини використовуємо схему корекції 

              (8.21)

Для оцінки похибки корекції розкладемо в ряд Тейлора в околі точки  корекцію

та саму функцію           .

Віднявши ці два розкладання, отримаємо

              .                  (8.22)

Для оцінки похибки прогнозу розкладемо в ряд Тейлора в околі точки  прогноз

 та саму функцію   

.

Віднявши ці два розкладання, отримаємо похибку прогнозу

        .                  (8.23)

З цих оцінок зрозуміло, що прогноз наближає розв’язок з недостачею, а корекція з надлишком. Отже, точне значення розв’язку лежить між прогнозом і корекцією. На будь-якому кроці можна оцінити точність розв’язку. Якщо задається її значення, то .

Віднімаючи рівності (8.23) та (8.22), маємо

                                    .

Уточнюємо розв’язок, виходячи з формули (8.22)

                         .               (8.24)

Відтак формула (8.24) завершує схему прогнозу і корекції..

Приклад 1 Розв’язати задачу Коші для диференціального рівняння другого порядку:

Розв’язання. 1 В ведемо нові змінні:

U=y; V= . Тоді, обравши сітку , визначимо на ній сіткові функції :

2    Знаходимо  і  () за допомогою розкладання в ряд Тейлора:

  3    З диференціального рівняння знаходимо  і :

Визначаємо прогноз розв’язку в точці

.

Тимчасово покладаючи U₂=P₂;, можна отримати

Обчислюємо корекцію . Тепер можна визначити V2

І нарешті на даному кроці уточнюємо розв’язок

Похибка обчислення.

Приклад 2 Використовуючи метод прогнозу і корекції, реалізувати алгоритм розв’язку крайової задачі для звичайного диференціального рівняння з точністю  на псевдокоді.

Розв’язання

Після заміни

отримаємо (користуючись системою символічної математики Derive)

vo1:=-0.3-vo/(0.7)

vo2:=1-2*vo-(vo1*(0.7)-vo)/(0.7)^2

vo3:=-2*vo1-(vo2*(0.7)^3-2*(0.7)*(vo1*(0.7)-vo))/(0.7)^4

u1:=0.5+0.3*vo+vo1*(0.3)^2/2+vo2*(0.3)^3/6

(24018*vo+48307)/98000

v1:=vo+0.3*vo1+vo2*(0.3)^2/2+vo3*(0.3)^3/6

(2099500*vo-129339)/3430000

 

2*u1+3*v1=1.2

SOLVE(2*u1+3*v1=1.2,vo);

Simp(#9)

[vo=1122527/7979760]; Approx(#10)

[vo=5122/36411=0.140671]

тоді розв’яжемо систему:

//x0 – початкове значення

//xn – кінцеве значення

//h – початковий крок

//X,Y – масиви результатів

Prognose_Correction(x0,u0,v0,h,X,Y):

1 u01:=v01;

2 v01:=x0-2*u0-v0/x0

3 u02:=v01

4 v02:=1-2*u01-((v01*x0-v0)/sqr(x0))

5 u03:=v02

6 v03:=-2*u02-((v02*x0*sqr(x0)-2*v01*sqr(x0)+2*v0*x0)/sqr(sqr(x0)))

7 u1:=u0+h*u01+h*h*u02/2+h*h*h*u03/6

8 v1:=v0+h*v01+h*h*v02/2+h*h*h*v03/6

9 i:=0

10 while (x1<=xN) do

11 i++

Repeat

13 x1:=x0+i*h

14 x2:=x0+(i+1)*h

15 v11:=x1-2*u1-v1/x1

16 u11:=v1

17 p2:=u0+2*h*u11

18 p2z:=v0+2*h*v11

19 u2:=p2

20 v2:=p2z

21 u21:=p2z

22 v21:=x2-2*p2-p2z/x2

23 c2z:=v1+h*(v11+v21)/2

24 v2:=c2z+(p2z-c2z)/5

25 u21:=v2

26 c2:=u1+h*(u11+u21)/2

27 u2:=c2+(p2-c2)/5

28 if abs(p2-c2)<eps then

29 X[i]:=x2;

30 Y[i]:=u2;

31 fi

32 if abs(p2-c2)>e then

33 h:=h/2;

Fi

35 until abs(p2-c2)<eps;

36 done //while

end

Методи Адамса

 

На відміну від однокрокових методів, у яких числовий розв’язок одержують тільки з диференціального рівняння і початкової умови, алгоритми Адамса складаються з двох частин: перша з них – стартова процедура для визначення  (наближені значення точного розв’язку в точках ), а друга – багатокрокова формула для одержання наближеного значення точного розв’язку . Потім ця формула застосовується рекурсивно для того, щоб за числовим розв’язком на  послідовних кроках обчислити  і т.д.

Стартові значення можна одержати декількома способами. Дж. К. Адамс обчислював їх за допомогою розкладання точного розв’язку в ряд Тейлора. Інший спосіб полягає у використанні якого-небудь однокрокового методу, наприклад, Рунге-Кутта. Стартові значення часто також обчислюють методами Адамса низького порядку з дуже малим кроком.

Розглянемо чисельні методи розв’язання задачі Коші (8.1)-(8.2), які можуть бути задані формулою

.  (8.25)

Тут значення розв’язку  в точці  визначається через значення розв’язку в  точках, що передують . Такий метод називається  - кроковим.

З класу (8.25) виділимо багатокрокові методи вигляду

,      (8.26)

застосовувані на сітці з постійним кроком

(8.27)

Різниця між найбільшим і найменшим значеннями індексу невідомої функції у _n, що входить у рівняння (8.26), дорівнює . Тому співвідношення (8.26) є різницевим рівнянням -го порядку, загальний розв’язок якого залежить від  параметрів. Щоб виділити єдиний розв’язок цього рівняння, необхідно задати  додаткових умов на функцію у_п. Цими додатковими умовами є значення функції у _n при n = 0,1,..., -1:

            (8.28)

які передбачаються відомими.

Використовуючи значення (8.28), з рівняння (8.26) при n =0 можна знайти , потім, використовуючи значення  і покладаючи в (8.26) n =1, знайти  і т.д. Таким чином, даний метод чисельного розв’язання диференціального рівняння полягає в розв’язанні різницевої задачі Коші для різницевого рівняння (8.26) і початкових умов (8.28).

Якщо шуканий розв’язок  входить до правої частини цього рівняння, що буває, коли , то формула (8.26) визначає неявний метод. Якщо , то шуканий розв’язок до правої частини не входить і рівняння (8.26) може бути розв’язане відносно . У цьому випадку формула (8.26) визначає явний метод.

Виведемо групу явних багатокрокових формул. Для точок сітки введемо позначення  і припустимо, що нам відомі числові наближені значення  точного розв’язку  задачі (8.1)-(8.2).

З диференціального рівняння випливає

.               (8.29)

До правої частини (8.29) входить шуканий розв’язок . Але оскільки нам відомі його наближені значення , то ми маємо також і величини

,    (8.30)

а тому природно замінити функцію  в (8.29) інтерполяційним многочленом, що проходить через точки . Його можна виразити через скінченні різниці вигляду  у такий спосіб:

           (8.31)

(інтерполяційна формула Ньютона). Тоді чисельний аналог (8.29) задається формулою , або після підстановки (8.31)

,                  (8.32)

де коефіцієнти  задовольняють рівність

.                      (8.33)

Числові значення цих коефіцієнтів наведені в таблиці 8.1.