Черкаси 2009 6 страница

, (8.10)

Також крок можемо вибирати за методом трьох зон:

(8.11)

Остаточно вибирається найменший крок серед всіх значень hⁱ. У формулі (8.11) часто використовують значення m = 2. При жорстких вимогах до точності значення m вибирають більшим.

 

Приклад 8.2: Розв’яжемо систему звичайних диференційних рівнянь, задану формулою 8.7.

1. Задаємо початкові дані х_поч , х_кін , e_доп: х_{1, 0} = х_{2, 0} = 0, e_доп = 0, 1.

2. Розрахуємо крок інтегрування за формулою 8.10 чи 8.11. Для початкових ітерацій виберемо його таким як і в попередньому прикладі: h = 0, 3.

3. Розв’яжемо систему відносно х_к+1.

4. Виведемо на екран результати.

5. Розрахуємо значення невідомих на перших трьох ітераціях.

6. Після 3-ї ітерації обчислимо (формула 8.9). Якщо умови точності виконуються та , перейдемо до наступної ітерації з тим же кроком, тобто до п.3. Інакше, якщо та , розрахуємо новий крок – перейдемо до п.2. Якщо , припиняємо ітерації.

 

Неявний метод Ейлера може застосовуватися для розв’язання більшості систем звичайних диференційних рівнянь як з малим, так і з великим розсіюванням постійних часу, але трудомісткість кожного кроку значно вище ніж у попередньому методі. Він абсолютно стійкий в області уявної вісі, але дає більш сильне затухання коливальних процесів у порівнянні з реальними коливаннями, тобто має велику похибку розв‘язків [3].

 

8.4. Модифікації методу Ейлера

 

Існують різноманітні модифікації методу Ейлера, які спрямовані на уточнення напрямку переходу з точки в точку . Розглянемо декілька простих модифікацій явного методу Ейлера [17]. Всі наведені нижче методи є методами другого порядку точності, тобто локальна похибка на кожному кроці має порядок . Умовою припинення ітерацій для всіх наведених у даному пункті методів є виконання нерівності:

(8.12)

Виправлений метод Ейлера. Припустимо, що наближене значення є розв’язком задачі Коші й необхідно обчислити , де . Розкладемо розв’язок за формулою Тейлора p -го порядку, приймаючи за базову точку і покладемо в цьому розкладі . Отримаємо

(8.12)

Якщо обмежитись двома першими елементами в правій частині розкладу (8.12), то отримаємо звичайний метод Ейлера (8.6).

Нижче розглянемо, що дає врахування третього доданку. При із (8.12) отримаємо рівняння:

(8.13)

Значення першої похідної в точці , приблизно відоме:

(8.14)

Диференціюючи рівняння за формулою повної похідної

знаходимо наближене значення другої похідної:

(8.15)

Підставляючи наближені вирази , та в рів­ність (8.13), отримуємо формулу для обчислення значення функції у наступній точці при :

(8.16)

Це формула виправленого методу Ейлера.

Оскільки при формули (8.14) і (8.15) є точними, а , згідно початкової умови задачі Коші, то на першому кроці обчислення за формулою (8.16) буде виникати похибка, пов’язана тільки з усіканням ряду Тейлора. Тобто, локальна похибка методу (8.16) складає величину , таким чином виправлений метод Ейлера належить до методів другого порядку.

 

Удосконалений метод Ейлера називають методом середньої точки, оскільки в ітераційних формулах використовується проміжна точка , яка ділить відрізок між і -ю та (і+1)-ю точкою навпіл. Цим методом розв’язок задачі Коші обчислюється в два етапи за формулами:

(8.17)

 

Метод Ейлера-Коші (метод Хойна) є одним з варіантів спільного використання методу Ейлера та неявного методу трапецій. В даному методі розрахунки проводяться спочатку за формулою явного методу Ейлера (8.18), а на другому етапі значення функцій усереднюється (формула 8.19):

(8.18)

(8.19)

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

 

Удосконалений метод Ейлера-Коші з ітераційною обробкою. Можна досягти більшої точності, якщо, виходячи з того ж самого початкового наближення

(8.20)

зробити не одну, а кілька ітерацій (як правило k = 3 ¸ 6) за методом трапецій:

(8.21)

 

Уточнений метод Ейлера. Для отримання наступної модифікації методу Ейлера, проінтегруємо рівняння на відрізку :

,

звідки отримаємо:

(8.22)

Застосовуючи до останнього інтегралу одноточкову квадратурну формулу середніх прямокутників і замінюючи значення і відомими наближеними значеннями та відповідно, із (8.22) отримаємо формулу для обчислення наближеного значення :

(8.23)

яку будемо називати уточненим методом Ейлера.

Звернемо увагу на одну принципову відмінність методу (8.23) від всіх інших розглянутих до цього моменту методів: уточнений метод є двокроковим. Тут для обчислення значення використовуються два попередніх значення та . Двокроковість накладає певні обмеження, принаймні, на початок чисельного процесу: значення не може бути знайдене безпосередньо за цим методом. Тому недостатню для процесу (8.23) другу початкову точку треба отримувати іншим методом, наприклад, явним методом Ейлера. Так, необхідне значення можна знайти за формулами (8.6) чи (8.17), а потім продовжувати обрахунки за формулою (8.23).

 

8.5. Метод Гюна

 

Цей метод використовує фундаментальну теорему аналізу, тобто інтегрування функції y¢ (x) на інтервалі [ x_k; x_k+1 ]:

(8.24)

Якщо з формули (8.24) виразити y(x_k+1), то отримаємо формулу (8.22):

Для наближеного обчислення визначеного інтегралу використаємо метод трапецій з кроком h = x₁ – x₀:

(8.25)

Оскільки рівняння неявне – у правій і лівій частинах є значення – для перетворення його у явний вигляд скористаємося формулою Ейлера для , що знаходиться у правій частині рівняння:

(8.26)

Для обчислення послідовності точок, що є розв‘язком задачі, використовуємо перше наближення за методом Ейлера, позначивши вираз для заміненого через :

(8.27)

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

(8.28)

Якщо ми обрахуємо на кожному кроці лише цю похибку, то накопичена похибка на m кроках:

(8.29)

Тобто при зменшенні кроку вдвічі, похибка повинна зменшитися в 4 рази.

8.6. Метод рядів Тейлора

 

Цей метод є еталонним методом, за яким порівнюють точність різних чисельних методів при розв‘язанні задачі Коші.

Теорема Тейлора: Припустимо, що і що функцію y(x) можна розкласти в ряд Тейлора порядку n у околі фіксованої точки х = х_k Î [x₀, x_n]:

(8.30)

де інтеграл наближено представлений рядом Тейлора:

(8.31)

У формулі 8.31 є (j-1)-ю похідною по (х_k) функції f. Отримаємо формулу для похідної за допомогою рекурентної процедури:

(8.32)

де є частинними похідними по відповідним змінним.

У загальному виді:

(8.33)

де 3 – диференційний оператор. При розв‘язанні задачі Коші ітерація Тейлора має вигляд:

(8.34)

де , . Остаточна похибка методу Тейлора має порядок . Точність методу Тейлора залежить від кількості елементів ряду. Для n = 4 і для кроків h та h/2 похибка обчислень:

(8.35)

8.7. Методи Рунге-Кутта

 

Ці методи виводяться з відповідних методів Тейлора таким чином, щоб остаточна загальна помилка методів не перевищувала порядку О(h ^N). Методи були виведені для того, щоб на кожному кроці виключити необхідність підрахунку похідних вищих порядків.

Методи Рунге-Кутта будуються для будь-яких порядків N, що обумовлені точністю. Найчастіше застосовується метод Рунге-Кутта 4-го порядку (N=4). Він дає досить малу помилку при невеликій кількості розрахунків і не складний для програмування. Метод оснований на обчисленні наступного значення за попереднім по формулі:

(8.36)

де k₁, k₂, k₃, k₄ – коефіцієнти, що розраховуються за формулами:

(8.37)

Рунге і Кутт отримали систему рівнянь, що має локальну похибку відсікання порядку О(h⁵):

(8.38)

Отримані 11 рівнянь вміщують 13 змінних, тому для їх роз­в’язання потрібні 2 додаткові умови. Такими умовами найчастіше є:

(8.39)

За таких умов розв’язками системи рівнянь будуть:

Підставляючи значення (8.38) та (8.39) в початкове рівняння (8.36) отримаємо формулу для стандартного методу Рунге-Кутта 4-го порядку:

, (8.40)

де

Рис. 8.3. Наближення інтегралу за методом Рунге-Кутта 4-го порядку.

 

Накопичена похибка в методі Рунге-Кутта складається з елементарних похибок, що утворюються на кожному кроці. Дана похибка методу розраховується за формулою:

. (8.41)

Приклад 8.3. Задана функція y¢ = (x-y)/2 на інтервалі [ 0; 3 ] з початковою умовою у(0)=1. Виберемо крок h = 0, 25. На першій ітерації підрахуємо значення коефіцієнтів f₁, f₂, f₃, f₄ і чергове значення y₁:

Всі інші значення на інтервалі [ 0; 3 ] розраховуються аналогічно, не застосовуючи операцію диференціювання.

 

Контрольні питання

1. Що є розв‘язком звичайного диференційного рівняння?

2. Що називається задачею Коші?

3. За яких умов задача Коші має єдиний розв‘язок?

4. Які основні вимоги пред‘являються до чисельних методів розв‘язання диференційних рівнянь?

5. Від чого залежить точність розв‘язання диференційного рівняння? Порівняйте ці складові з переліком похибок чисельного методу з п. 1.3.

6. Від яких критеріїв залежить розрахунок кроку інтегрування?

7. В яких випадках застосовують адаптивний вибір кроку інтегрування?

8. В чому полягає суть явного методу Ейлера?

9. В яких випадках застосовується неявний метод Ейлера і чому?

10. З якою метою розроблені модифікації явного методу Ейлера?

11. За рахунок чого модифікації методу Ейлера дозволяють отримати більш точний результат, ніж при застосуванні явного методу Ейлера?

12. Який метод розв‘язання задачі Коші є базовим для оцінки похибки і чому?

13. В чому полягає суть методу Гюна?

14. В чому полягає перевага методів Рунге-Кутта при розв‘язанні звичайних диференційних рівнянь?

15. В чому полягає суть алгоритму розв‘язання диференційних рівнянь методом Рунге-Кутта 4-го порядку?

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

 

Методи Ейлера, Гюна, Тейлора, Рунге-Кута називаються однокроковими методами, тому що в них використовується лише інформація про значення диференційного рівняння в одній попередній точці для того, щоб обрахувати значення в наступній, тобто лише значення в початковій точці (t ₀; y ₀) використовується для обчислення значення диференційного рівняння в точці (t ₁; y ₁), і загалом, щоб обрахувати у_{к +1} необхідне лише у_к. Але після знаходження значень диференційного рівняння в декількох точках можна їх використовувати для отримання більш точних результатів у наступній вузловій точці. Наприклад, для чотирикрокового методу Адамса-Бешфорса-Маултона необхідні значення в точках у_{к -3}, у_{к -2}, у_{к -1}, та у_к для обчислення у_{к +1}. У цьому методі значення диференційного рівняння в чотирьох початкових точках, (t ₀, y ₀), (t ₁, y ₁), (t ₂, y ₂) та (t ₃, y ₃), розраховуються спочатку для подальшої генерації точок {(t _k, y _k): k 4} з використанням багатокрокового методу.

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

 

9.1. Метод Адамса-Бешфорса-Маултона

Метод прогнозу-корекції Адамса-Бешфорса-Маултона – це багатокроковий метод, виведений із фундаментальної теореми аналізу:

. (9.1)

Прогноз використовує наближення поліномом Лагранжа для функції f (t, y (t)), що побудована по точках (t_{k -3}, y_{k -3}), (t_{k -2}, y_{k -2}), (t_{k -1}, y_{k -1}) i (t_k, y_k). Для функції f (t, y (t)) з формули (9.1), яка інтегрується на інтервалі [ t_k, t_{k +1}] прогноз Адамса-Бешфорса, має вигляд

. (9.2)

Як будь-який прогноз, він потребує корекції. Коректор отримується аналогічно. А саме, як тільки значення p_{k +1} обраховано, його використовують для побудови наступного поліному Лагранжа для функції f (t, y (t)), який будується по точкам (t_{k -2}; y_{k -2}), (t_{k -1}; y_{k -1}) (t_k; y_k) і новій точці (t_{k +1}; y_{k +1}), в якій замість y_{k +1} використовується прогноз p_{k +1}, тобто функція в точці t_{k +1}; обраховується як f (t_{k +1}, p_{k +1}). Після розрахунку поліному шляхом інтегрування на інтервалі [ t_k; t_{k +1}] отримаємо коректор Адамса-Маултона:

. (9.3)

 

Оцінка помилки і корекція. Залишковий елемент формули чисельного інтегрування використовується, щоб отримати і прогноз, і коректор порядку О(h ⁵). Локальна помилка відсікання для формул (9.2) та (9.3) має вигляд:

, (локальна похибка відсікання

для прогнозу) (9.4)

. (локальна похибка відсікання

для коректора) (9.5)

Припустимо, що h мале і y ⁽⁵⁾ (t) – майже усталена на інтервалі. Тоді можна виключити члени, що містять похідну 5-го порядку у формулах (9.4) та (9.5). В результаті отримаємо:

. (9.6)

Переваги методу прогнозу-корекції Адамса-Бешфорса-Маултона у тому, що формула (9.6) дає наближену оцінку помилки, яку отримуємо при обчисленні значень p_{k +1} i y_{k +1}, не використовуючи у ⁽⁵⁾(t).

Коректор (9.3) використовує наближення при обчисленні у_{к +1}. Оскільки у_{к +1} також є оцінкою для у (t_{k +1}), її можна використовувати у коректорі (9.3) для генерування нового наближення для f_{k +1}, яке, в свою чергу, буде генерувати нове значення для у_{к +1}. Окрім того, якщо продовжити цю ітерацію у коректорі, вона буде збігатись до фіксованої точки у формулі (9.3) швидше, ніж при класичному розв‘язку диференційного рівняння. Якщо необхідно підвищити точність, то процедура корекції більш ефективна, ніж зменшення довжини кроку.

 

Формулу (9.6.) можна використовувати для покращення значень прогнозу. Якщо припустити, що різниця між значеннями прогнозу і корекції на кожному кроці змінюється повільно, то можна підставити і замість і в (9.6) і отримати наступний управляючий параметр:

(9.7)

Його значення використовується замість на кроці корекції, і формула (9.3) приймає вигляд:

 

(9.8)

 

Формулу (9.6) також можна використовувати для визначення адаптивної довжини кроку. Можливо зменшити довжину кроку, наприклад, до h /2 і збільшити до 2 h. Припустимо ε = 5× 10^-6 – критерій відносної похибки, і нехай kε = 10^-5

Якщо , тоді h = h /2 (9.9)

Якщо , тоді h =2 h (9.10)

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

Зменшення довжини кроку потребує чотирьох нових початкових значень. Для заміщення значень, яких не вистачає і які ділять інтервали [ t_{k -2}, t_{k -1}] і [ t_{k -1}, t_k ] навпіл, використаємо інтерполяцію функції f (t, y (t)) поліномом четвертого ступеню. Нові чотири вузлові точки, t_{k -3/2}, t_{k -1}, t_{k -1/2} і t_k, використовуються в наступних розрахунках.

Інтерполяційна формула, необхідна для отримання нових початкових значень (рис. 9.1) для кроку h/2 має вигляд:

(9.11)

Рис. 9.1. Зменшення довжини кроку до h /2.

 

Набагато простіше збільшити довжину кроку, тоді для розрахунків за одно-кроковим методом потрібно розрахувати сім попередніх точок (рис. 9.2), щоб подвоїти крок.

Рис. 9.2. Збільшення довжини кроку до 2 ∙ h.

 

 

9.2. Метод Мілна-Сімпсона

Інший розповсюджений метод прогнозу-корекції це метод Мілна-Сімпсона. Його прогноз базується на інтегруванні функції f (t, y (t)) на інтервалі [ t_{k -3}, t_{k +1}]:

(9.12)

Прогноз використовує наближення поліномом Лагранжа для f(t, y(t)), побудоване по точкам (t_{k -3}, f_{k -3}), (t_{k -2}, f_{k -2}), (t_{k -1}, f_{k -1}) та (t_k, f_k). Інтегруємо його по інтервалу [ t_{k -3}, f_{k +1}] і отримуємо прогноз Мілна: