Подвійний перерахунок для методів Рунге – Кутта

Нехай задані послідовність х₀, х₁, …, х_n рівновіддалених точок: х_k = х₀ + kh з кроком інтегрування диференціального рівняння h і послідовність , , …, з кроком h/2 так, щоби = х_k (k = 0, 1, …, n). Розглянемо задачу Коші (1) у´ = f (x, y); у(х₀) = у() = у₀, її чисельні розв’язки у₀, у₁, …, у_n з кроком h і , , …, з кроком h/2, отримані деяким методом Рунге – Кутта порядку точності s. Нехай ε_k = у_k – у(х_k) і = – у(х_k) відповідні похибки наближених значень у_k та у деякій точці х_k = .

Метод подвійного перерахунку для методів Рунге – Кутта є узагальненням такого методу для інтегрування функцій і теж ґрунтується на двох аналогічних формулах.

1. ε = – ≈ (правило Рунге) (6)

2. у(х_k) ≈ + ε = + (аналог формули екстраполяції за Річардсоном). (7)

Практичне застосування цих формул розглянемо в наступному параграфі. Тут далі наведені обґрунтування цих формул і, як результат, їх точне формулювання.

На кожному з відрізків [х_k; х_k+1] (k = 0, 1, …, n) завдовжки h розглянемо задачу Коші (1) у´ = f (x, y); у(х_k) = у_k і заданим методом Рунге – Кутта дістанемо наближене значення розв’язку у(х_k + h) в точці х_k+1, яке за умовою дорівнює у_k+1 та похибку методу на кроці φ_k(h) = у_k+1 – у(х_k + h) = h^s+1 (0 < θ_k < 1) згідно з (4). Зокрема, на першому кроці φ₀(h) = у₁ – у(х₀ + h) = h^s+1 (0 < θ₀ < 1), на другому φ₁(h) = у₂ – у(х₁ + h) = h^s+1 (0 < θ₁ < 1). Як було зазначено та й видно безпосередньо, величина похибки методу Рунге – Кутта порядку точності s на кроці має порядок s + 1: φ_k(h) = О(h^s+1). Розглянемо перетворення T_mk(z) (m < k), яке кожній початковій умові z для задачі Коші у´ = f (x, y); у(х_m) = z ставить у відповідність значення розв’язку в точці х_k: у(х_k) = T_mk(z). Тоді під дією диференціального рівняння відрізок [у₁; у(х₀ + h)] між точним та наближеним значенням, отриманими на кроці 1, за другий крок пересунеться в [T₁₂(у₁); T₁₂(у(х₀ + h)) ]. Довжина цього відрізку ôT₁₂(у₁) – T₁₂(у(х₀ + h))ô = А₁₂ôу₁ – у(х₀ + h)ô = А₁₂ôφ₀(h)ô, де А₁₂ > 0 – стала. Оскільки ôT₁₂(у₁) – T₁₂(у(х₀ + h))ô = ô ô = ôT´₁₂(ζ)ô∙ôу₁ – у(х₀ + h)ô = ôT´₁₂(ζ)ô∙ôφ₀(h)ô, то А₁₂ = ôT´₁₂(ζ)ô, де ζ Î [у₁; у(х₀ + h)]. Сумарна похибка на другому кроці ε₂ = у₂ – у(х₂) дорівнює сумі похибки на цьому кроці φ₁(h) і похибки на першому кроці φ₀(h), перетвореній за другий крок диференціальним рівнянням, тобто перетворенням T₁₂: ε₂ = у₂ – у(х₂) = φ₁(h) + А₁₂ φ₀(h). На третьому кроці сумарна похибка дорівнює сумі похибки на цьому кроці φ₂(h), похибки на другому кроці φ₁(h), перетвореній за третій крок диференціальним рівнянням, тобто перетворенням T₂₃ та похибки на першому кроці φ₀(h), перетвореній за другий та третій кроки дією T₁₃: ε₃ = у₃ – у(х₃) = φ₂(h) + А₂₃ φ₁(h) + А₁₃ φ₀(h), де А₁₃ = ôT´₁₃(ζ₁₃)ô, ζ₁₃ Î [у₁; у(х₀ + h)], А₂₃ = ôT´₂₃(ζ)ô, ζ Î [у₂; у(х₁ + h)]. У загальному випадку сумарна похибка на кроці k

ε_k = у_k – у(х_k) = φ_k-1(h) + А_k-1k φ_k-2(h) + А_k-2k φ_k-3(h) + … + А_1k φ₀(h) = ,

де А_kk = 1, а для i < k А_ik = ôT´_ik(ζ_ik)ô, де ζ_ik Î [у_i; у(х_i-1 + h)].

На кожному відрізку [х_k; х_k+1] розташовані два відрізка: [ ; ] та [ ; ] (бо х_k = , х_k+1 = ). Аналогічно попередньому на кожному з цих відрізків завдовжки h/2 розглянемо задачу Коші (1) у´ = f (x, y); у() = і у´ = f (x, y); у() = . Даним методом Рунге – Кутта дістанемо наближене значення розв’язку у( + h/2) у точці , яке за побудовою дорівнює та наближене значення розв’язку у( + h/2) у точці , яке дорівнює . Дістанемо також похибки методу на кроці = – у( + h/2) = (h/2)^s+1 (0 < θ_2k < 1/2) і = – у( + h/2) = (h/2)^s+1 (0<θ_2k+1<1/2), які мають порядок s + 1. Розглянемо також перетворення (m < k), яке кожній початковій умові z для задачі Коші у´ = f (x, y); у() = z ставить у відповідність значення розв’язку в точці : у() = . Тоді – сумарна похибка на кроці 2k – це, як і вище, = – у(х_k) = =

= + + + … + ,

де = 1, а для i < k = ô ô, Î [ ; y( + h/2)]. Отже, одночасно

у(х_k) = у_k – ε_k = у_k – (8)

у(х_k) = – = – . (9)

На відрізку [х_i; х_i+1] розташовані два відрізка [ ; ] та [ ; ], (х_i = , х_i+1 = ), тож доданку А_ikφ_i-1(h) = А_ik h^s+1 (0 < θ_i-1 < 1) у (8) відповідають два доданки у (9): = (h/2)^s+1 (0 < θ_2i-1 < 1/2) і = (h/2)^s+1 (0 < θ_2i < 1/2). Тут А_ik = ôT´_ik(ζ_ik)ô, де ζ_ik Î [у_i; у(х_i-1+ h)], =ô ô, де Î [ ;y( + h/2)], =ô ô, де Î [ ; y( + h/2)] відрізняються між собою величиною порядку h згідно з теоремою Лагранжа, оскільки за означенням T_ik = , а відрізняється від значеннями порядку h. Так само і величини , і за означенням відрізняються між собою величиною порядку h згідно з теоремою Лагранжа. Тому спочатку знехтуємо відмінністю між величинами А_ik , , порядку h і покладемо М_i = А_ik ≈ ≈ . Тоді (8) і (9) набувають вигляду

у(х_k) = у_k – ε_k ≈ у_k – h^s+1 = у_k – Сh^s (10)

у(х_k) = – ≈ – 2 (h/2)^s+1 = – С(h/2)^s, (11)

де С = h. (Якщо відрізок [х₀, х_k] вважати порядку одиниці, то k ~ , а отже і С = h порядку одиниці). Відомими в (11) є у_k, , h і s, невідомими у(х_k) та С. Отже, це система двох лінійних рівнянь з двома невідомими. Звідси – С(h/2)^s – (у_k – Сh^s) ≈ 0; – у_k ≈ С()^s – Сh^s = С()^s(1 – 2^s) = (1 – 2^s), звідки знаходимо (6) ≈ – (правило Рунге). Тому (7): у(х_k) ≈ + (аналог формули екстраполяції за Річардсоном). Єдине неточне припущення, яке було зроблено в цьому доведенні – це рівність А_ik ≈ ≈ , насправді ці величини відрізняються між собою на величини порядку h. Тому різниці між величинами

А_ikφ_i-1(h) = А_ik h^s+1, = h^s+1 і = h^s+1 – це величини порядку h^s+2, а різниці між їх сумами по всім k = 0, 1, …, n – величини порядку h^s+1. Врахування цих величин приводить до точного варіанту формул (6) і (7):

1. = – + О(h^s+1), 2. у(х_k) = + + О(h^s+1). (12)

Отже, = + – це наближення до точного розв’язку у(х_k) порядку точності s + 1 у точці х_k. Нехай ще дана послідовність рівновіддалених точок , , …, з кроком інтегрування h/4 так, щоби = х_k (k = 0, 1, …, n). Розглянемо задачу Коші (1) у´ = f (x, y); у(х₀) = у() = у₀, її чисельні розв’язки , , …, , отримані тим же методом Рунге – Кутта порядку точності s. Тоді метод подвійного перерахунку так само можна застосувати до цієї послідовності наближень і до попередньої , , …, з кроком h/2 у точках , , …, ( = ). В результаті отримаємо оцінку похибки за правилом Рунге і уточнене значення за аналогом формули екстраполяції Річардсона:

1. = – + О(h^s+1), 2. у(х_k) = + + О(h^s+1). (13)

Тут = – – це, як і , ще одне наближення порядку точності s + 1 до точного розв’язку у(х_k) у точці х_k. Тоді можна знову застосувати правило Рунге до цих двох уточнених наближень і : = – + О(h^s+2). Тобто насправді виявиться похибкою порядку h^s+1, а уточнене наближення = + буде мати порядок уже принаймні порядок точності s + 2. Процес можна продовжити і далі, отримуючи наближення порядку s + 3, s + 4, …. Доведення такого методу кратного перерахунку принципово таке ж, хоча вимагає значних технічних зусиль. Реалізацію цього метода за допомогою Excel розглянемо у наступному параграфі.