Чисельне розв’язання звичайних диференціальних рівнянь. Різницева апроксимація диференціальних рівнянь однокроковими методами, страница 5

Виконаємо тепер розрахунок за тією самою наближеною формулою для тієї самої точки х, але використовуючи рівномірну сітку з іншим кроком 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,  t0=0, T=1,  y(0) =1 та оцінити його похибку.

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

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

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

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

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

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

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

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

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

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

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

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

;

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

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

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

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

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

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

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

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

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