Розв’язування звичайних диференціальних рівнянь
Диференціальним називається рівняння, в яке входять похідні невідомої функції.
Приклад:
(1)
(2)
Диференціальне рівняння (ДР), що містить лише одну незалежну змінну і похідні за нею, називають звичайними (ДР). Це, наприклад, рівняння (1) . ДР, що містить декілька незалежних змінних і похідні за ними, називають рівняння в частинних похідних. Ми розглянемо методи розв’язування звичайних диференціальних рівнянь.
Порядком ДР називається найвищий порядок похідної ( або диференціалу), який входить в рівняння. Звичайне ДР (ЗДР) -го порядку в загальному випадку має незалежну змінну, невідому функцію та її похідні (або диференціал) до -го порядку включно:
(3)
- незалежна змінна;
- невідома функція (залежна змінна);
- похідні цієї функції.
Диференціальне рівняння -го порядку, розв’язане відносно старшої похідної, може бути записано у вигляді:
(4)
Щоб розв’язати ЗДР, необхідно мати значення залежної змінної та (або) її похідних при деяких значення незалежної змінної.
Якщо ці значення задані при одному значенні незалежної змінної - така задача називається задачею з початковими умовами або задачею Коші.
Якщо ці значення задаються при або більше значеннях незалежної змінної - задача називається крайовою.
Значення залежної змінної та її похідних називаються ще додатковими умовами.
В задачі Коші додаткові умови називаються початковими.
В крайовій задачі - граничними.
Задача Коші
Задача Коші формулюється так :
Нехай задане ДР
(5)
з початковими умовами . Потрібно знайти функцію , що задовольняє дане рівняння, та початкову умову. Чисельний розв’язок цієї задачі одержують так. Спочатку обчислюють значення похідної, потім задаючи малий приріст , переходять до нової точки
Положення нової точки визначають за нахилом кривої, обчисленому з допомогою ДР. Таким чином, графік чисельного розв’язку являє собою послідовність коротких прямолінійних відрізків, якими апроксимується істинна крива . Сам чисельний метод визначає порядок дій при переході від даної точки кривої до наступної.
Існують дві групи методів розв’язування задачі Коші.
Однокрокові методи. В них для знаходження наступної точки на кривій потрібна інформація лише про попередній крок.
Однокроковими є метод Ейлера та методи Руте-Кутта.
Багатокрокові (або методи прогнозування та коригування). Для знаходження наступної точки кривої вимагається інформація більш ніж про одну з попередніх точок. До них належать методи Адамса, Мілна, Хеммінга.
Це чисельні методи розв’язування ДР. Вони дають розв’язок у вигляді таблиці значень.
Метод Ейлера
Однокрокові методи призначені для розв’язування диференціальних рівнянь першого порядку виду
(1)
Метод Ейлера є найпростішим методом розв’язування задачі Коші. Він дозволяє інтегрувати ДР першого порядку. Точність його не велика.
- настільки мале, що значення
функції мало відрізняється від
лінійної функції
- тангенс кута нахилу дотичної в
x
h
Тобто крива заміняється дотичними. Рух відбувається не по інтегральній кривій, а по відрізках дотичної .
Метод Ейлера базується на розкладі функції в ряд Тейлора в околі точки
(2)
Якщо мале, то, члени розкладу, що містять в собі і т.д. є малими високих порядків і ними можна знехтувати.
Тоді (3)
Похідну знаходимо з рівняння (1), підставивши в нього початкову умову. Таким чином можна знайти наближене значення залежної змінної при малому зміщенні від початкової точки. Цей процес можна продовжувати, використовуючи співвідношення.
,
роблячи як завгодно багато кроків.
Похибка методу має порядок , оскільки відкинуті члени, що містять в другій і вище степенях.
Недолік методу Ейлера - нагромадження похибок, а також збільшення об’ємів обчислень при виборі малого кроку з метою забезпечення заданої точності.
В методі Ейлера на всьому інтервалі тангенс кута нахилу дотичної приймається незмінним і рівним . Очевидно, що це призводить до похибки, оскільки кути нахилу дотичної в точках та різні. Точність методу можна суттєво підвищити, якщо покращити апроксимацію похідної.
Це можна зробити, якщо, наприклад, використати середнє значення похідної на початку та в кінці інтервалу. В т.з. модифікованому методі Ейлера (метод Ейлера з перерахунком) спочатку обчислюється значення функції в наступній точці за звичайним методом Ейлера.
(4)
Воно використовується для обчислення наближеного значення похідної в кінці інтервалу .
Обчисливши середнє між цим значенням похідної та її значенням на початку інтервалу, знайдемо більш точне значення :
(5)
Цей прийом ілюструється на рисунку.
xn h/2 xn+1
Принцип модифікованого методу можна пояснити інакше. Якщо в розкладі в ряд Тейлора зберегти член з
(6)
Замість другої похідної можна використати наближення кінцевою різницею
(7)
. Підставивши (7) в (6) одержимо
(8)
Що співпадає по формі з (5). Відмінність між (8) та (5): в (5) точне значення похідної замінимо на . Похибка при такій заміні має порядок .
Відмітимо, що за підвищення точності доводиться платити додатковими затратами машинного часу.
В обчислювальній практиці використовується також метод Ейлера-Коші з ітераціями:
знаходиться грубе початкове наближення (за звичайним методом Ейлера)
будується ітераційний процес
(9)
Ітерації продовжують до тих пір, доки два послідовні наближення не співпадуть з заданою похибкою . Якщо після декількох ітерацій співпадання нема, то потрібно зменшити крок .
Тобто в модифікованому методі Ейлера, в методі Ейлера-Коші з ітераціями спочатку (на першому етапі) знаходиться наближення для , а потім воно вже коригується за формулами (5) або (9).
Історичні відомості: Ейлер
Леонард Ейлер (1707 – 1783) народився і виховувався в Базелі (Швейцарія). Одним із його вчителів був Йоганн Бернуллі. Роки по тому Бернуллі відзивався про Ейлера, як про “найзнаменитішого і наймудрішого математика”. Ейлер писав і публікувався з великою продуктивністю: майже 600 книг і статей на протязі життя. Його вплив був настільки великим, що по крайній мірі в математиці, вісімнадцяте століття можна назвати епохою Ейлера. Ейлер зробив вклад не тільки в математику, але і в фізику, астрономію гідродинаміку, оптику, в теорію електроніки і магнетизму. В 1768 він запропонував метод рішення початкових задач. Він також ввів багато із звичних на сьогодні математичних позначень, в тому числі символ е для основи натурального логарифма, і для , ∆у для кінцевої різниці, ∑ для суми, запропонував символ f використовувати разом з дужками для позначення функції, а також придумав назву для тригонометричних функцій, які використовуються і сьогодні.
Ейлер 14 років працював в Російській академії наук, в Санкт-Петербурзі. Потім він переїхав в Берлін і працював в Берлінській академії до того часу, доки розбіжності з королем Фрідріхом Великим не змусили його повернутись в Росію, де він і провів останні роки свого життя. Він ніколи не повертався в Базель, покинутий ним у 1727 році.
Талановиті і удачливі люди дуже часто беруться за безліч проблем; Ейлер не був винятком. Він працював в області суднобудування, входив до різноманітних технічних комітетів, займався перевіркою засобів для зважування і пожежних помп, керував складанням календарних і географічних карт, слідкував за ботанічними садами. Він надто рано відмовився від першочергових намірів стати міністром, однак його інтерес до релігії і гуманітарної діяльності давав йому достатню перспективу по відношенню до інших напрямків життя. У Ейлера була чудова пам’ять. В 1771 р. в результаті хвороби він став сліпим, проте майже половина його наукової роботи написана після цієї дати, і у віці семидесяти років він міг точно відтворити перші і останні рядки кожної сторінки вергілієвської “Енеїди”, яку він вивчив на пам’ять, коли був молодим. Після смерті Ейлера на його могилі на лютеранському цвинтарі Санкт-Петербурга був встановлений масивний монумент.
Метод Рунге – Кутта четвертого порядку
Метод Рунге-Кутта об’єднує ціле сімейство методів розв’язування диференціальних рівнянь першого порядку. Найбільш часто використовується метод четвертого порядку, який просто називають “ методом Рунге-Кутта”.
В методі Рунге-Кутта значення функції , як і в методі Ейлера, визначається за формулою
(1)
Якщо розкласти функцію в ряд Тейлора і обмежитись членами до включно, то приріст можна записати у вигляді
(2)
Замість того, щоб обчислювати члени ряду за формулою (2) в методі Рунге-Кутта використовують наступні формули.
Це метод четвертого порядку точності.
Похибка на кожному кроці має порядок .Таким чином метод Рунге-Кутта забезпечує значно вищу точність ніж метод Ейлера, однак вимагає більшого об’єму обчислень ніж метод Ейлера. Це досить часто дозволяє збільшити крок .
Деколи зустрічається інша форма представлення методу Рунге-Кутта 4-го порядку точності.
Знайти розв’язок задачі Коші
Для порівняння точності різних методів розв’язування наведемо такий приклад:
Точний розв’язок
Метод Ейлера
Модифіко-ваний метод Ейлера
Метод
Рунге-Кутта IV порядку
Точний розв’язок
0,0
0,2
0,4
0,6
0,8
1,0
1,0000
1,4420
2,1041
3,1183
4,6747
7,0472
1,0000
1,4923
2,2466
3,4176
5,2288
8,0032
1,0000
1,4977
2,2783
3,5201
5,4894
8,5834
1,0000
1,4977
2,2783
3,5202
5,4895
8,5836
В більшості стандартних програм ЕОМ найчастіше використовується (схема) метод четвертого порядку (Рунге-Кутта).
Система звичайних диференціальних рівнянь
Система рівнянь першого порядку з невідомими функціями має вигляд:
(1)
де - незалежна змінна; - невідомі функції. Це система в нормальній формі, або нормальна система.
Розв’язком системи називається будь-яка система функцій , що має похідні першого порядку і перетворює рівняння системи в тотожності по .
Приклад. Розв’язком системи
з початковими умовами є функції
Справді
.
Підставляємо вирази в систему рівнянь, одержимо два тотожних рівняння
Рівняння -ого порядку
завжди можна звести до системи рівнянь першого порядку. З цією метою позначимо шукану функцію через і приймемо за нові невідомі.
Складаємо систему диференціальних рівнянь для функцій .
Маємо
Тому для функцій одержимо наступну нормальну систему диференціальних рівнянь:
причому
Приклад.
Проведемо заміни:
Таким чином, задача Коші для ЗДР -порядку з відповідними початковими умовами може бути зведена до розв’язку задачі Коші для системи диференціальних рівнянь першого порядку.
Праві частини рівнянь цієї системи в програмі, яка реалізує який -небудь чисельний метод, записуються у вигляді підпрограми. Наприклад, для системи ЗДР.
Підпрограма буде мати наступний вигляд:
Тобто:
функція
Алгоритм реалізації методу Рунге-Кутта заключається в циклічних обчисленнях значень функцій на кожному кроці по наступних формулах:
Або при іншій формі представлення методу:
Нехай дана така система звичайних диференціальних рівнянь
Методи з автоматичною зміною кроку
Застосовуються в тому випадку, якщо розв’язок потрібно одержати із заданою точністю. При високій точності (похибка ) автоматична зміна кроку забезпечує зменшення загального числа кроків в декілька разів (особливо при розв’язках у вигляді кривих, що сильно відрізняються крутизною).
Метод Рунге-Кутта з автоматичною зміною кроку
Після обчислення з кроком всі обчислення виконуються повторно з кроком . Після цього порівнюються результати, отримані в точці хn+1 з кроком і . Якщо модуль різниці менший , то обчислення продовжуються з кроком , в іншому випадку крок зменшують. Якщо нерівність дуже сильна, то крок збільшують.
Маємо
значення незалежної змінної в точці
- значення функції в точці
- значення функції в точці , обчислене з кроком
- значення функції в точці , обчислене з кроком
- значення функції , обчислене з кроком
1) Якщо
обчислення повторюються з кроком і т.д., доки не виконається умова │Y(I) –Y1(I)│≤ ε.
2) Якщо виконується ця умова, то можливі два варіанти, в залежності від значення K, де K – ознака поділу кроку.
Початкове значенняі залишається таким після першого поділу кроку на два. Надалі, якщо крок ділиться, то K приймає значення одиниці.
а) Якщо , то навіть коли виконалась умова , крок не змінюється , тобто лишається тим самим (обчислення далі проводяться з попереднім кроком).
б) Якщо і виконалась умова , тоді .
В обох випадках а) і б) результат Y(I) виводиться на друк.
Метод Рунге-Кутта-Мерсона з автоматичною зміною кроку
Метод дозволяє оцінити похибку на кожному кроці інтегрування. Похибка інтегрування має порядок . При цьому не потрібно зберігати в пам’яті обчислення значень функцій на кроці і для оцінки похибки - перевага порівняно з методом Рунге-Кутта з автоматичною зміною кроку.
Алгоритм методу
1. Задаємо число рівнянь , похибку , початковий крок інтегрування , початкові умови .
2. За допомогою п’яти циклів з керуючою змінною обчислюємо коефіцієнти
3. Знаходимо значення
та похибку
4. Перевіряємо виконання умов
Можливі випадки:
а) Якщо перша умова не виконується, тобто │Rj(n+1)│< ε, то ділимо крок на 2 та повторюємо обчислення з п.2, встановивши початкові значення .
б) Якщо виконується перша та друга умови, значення та виводяться на друк.
Якщо друга умова не виконується, крок збільшується вдвічі і тоді обчислення знову повторюється з п.2 (нема потреби обчислювати при малому кроці).
Треба відмітити, що похибка на кожному кроці методу Рунге-Кутта-Мерсона оцінюється приблизно. При розв’язуванні нелінійних ДР істинна похибка може відрізнятися в декілька разів від заданої .
, де .
- крок поділити на 2 і повернутися на початок.
для всіх рівнянь: виводимо на друк , а крок збільшуємо удвічі.
Метод Рунге-Кутта-Фельберга з автоматичною зміною кроку
Це метод четвертого порядку, дає більш точну оцінку похибки (порівняно з методом Рунге-Кутта-Мерсона) на кожному кроці і реалізується послідовним циклічним обчисленням за наступними формулами:
Похибка
Якщо
а) , крок зменшується в двічі
б) Якщо збільшується вдвічі.
Час розрахунку для однієї точки удвічі більший, ніж для методу Рунге-Кутта-Мерсона.
Загальна характеристика однокрокових методів
1) Щоб одержати інформацію в наступній точці, потрібно мати дані лише в одній попередній точці. Цю властивість деколи називають “самостартування”.
2) В основі всіх однокрокових методів лежить розклад функцій в ряд Тейлора, в якому зберігаються члени, котрі містять в степені до включно. Ціле число називають порядком методу. Похибка на кроці має порядок .
3) Всі однокрокові методи не потребують дійсного обчислення похідних - обчислюється лише сама функція. Однак можуть бути потрібні її значення в декількох проміжних точках, що викликає додаткові затрати машинного часу.
4) Властивість “самостартування” дозволяє легко змінювати значення кроку .
Жорсткі задачі
Деякі звичайні диференціальні рівняння не вирішуються ні одним із розглянутих методів.
Стала часу диференціального рівняння першого порядку - це проміжок часу, по закінченні якого величина нестаціонарної частини рішення зменшується в разів. В загальному випадку диференціального рівняння - порядку має постійних часу. Якщо будь-які дві з них сильно відрізняються по величині або якщо одна з них достатньо мала порівняно з інтервалом часу, для якого шукається розв’язок, то задача називається жорсткою і її практично неможливо розв’язати звичайними методами.
Методи прогнозу і корекції
Для обчислення положення нової точки використовується інформація про декілька раніше отриманих точок. Для цього використовуються дві так звані формули прогнозу і корекції. Схеми алгоритмів для таких методів приблизно однакові, а самі методи відрізняються лише формулами.
Оскільки в цих методах використовується інформація про декілька раніше отримані точки, то на відміну від однокрокових методів вони не мають властивості самостартування. Тому, перед тим як застосовувати метод прогнозу і корекції, необхідно обчислити вихідні дані з допомогою якого-небудь однокрокового методу.
Обчислення виконують таким чином. Спочатку за формулою прогнозу та початковими значеннями змінних визначають значення . Верхній індекс (0) означає, що прогнозоване значення є одним із послідовності значень уп +1 , розташованих в порядку зростання точності. За прогнозованим значенням з допомогою диференціального рівняння:
(1)
знаходимо похідну:
(2)
Ця похідна підставляється у формулу корекції для обчислення уточненого значення . В свою чергу використовується для одержання більш точного значення похідної:
(3)
Якщо це значення похідної недостатньо близьке до попереднього, то воно вводиться у формулу корекції і ітераційний процес продовжується. Якщо ж похідна змінюється в припустимих межах, то значення використовується для обчислення остаточного значення , яке і виводиться на друк. Після цього процес повторюється - робиться наступний крок, на якому обчислюється .
НІ
ТАК
ТАК
Метод Мілна
На етапі прогнозу використовується формула Мілна.
, (1)
а на етапі корекції - формула Сімпсона
(2)
Останні члени в обох формулах в ітераційному процесі не використовуються і служать лише для оцінки помилок. Метод Мілна відносять до методів четвертого порядку точності, оскільки в ньому відкинуті члени, які містять в 5 і більш високих степенях. Потрібно мати на увазі, що для користування формулою (1) необхідно попередньо одним із однокрокових методів визначити і значення похідних .
Похибка, внесена на будь-якому кроці, зростає експоненціально, тому методу Мілна властива нестійкість. Це недолік методу.
Метод Адамса
Має четвертий порядок точності. Формула прогнозу отримана на основі інтерполяційної формули Ньютона і має вигляд:
(1)
Формула корекції
(2)
На відміну від методу Мілна, похибка, внесена на кроці, не має тенденції до експоненціального зростання.
Метод Хемінга
Це стійкий метод четвертого порядку точності.
Формула прогнозу:
(1)
Формула корекції
(2)
Формула прогнозу зазвичай доповнюється формулою уточнення прогнозу:
(3)
(4)
Переваги методу Хемінга – простота і стійкість.
Характеристика методів прогнозу і корекції.
1) Для реалізації методів прогнозу і корекції необхідно мати інформацію про декілька попередніх точок, тобто ці методи не відносяться до числа самостартуючих. Для одержання вихідної інформації потрібно скористатись одним з однокрокових методів. Якщо в процесі розв’язування диференці-ального рівняння змінюється крок, то зазвичай потрібно тимчасово переходити на однокроковий метод.
2) Оскільки для методів прогнозу і корекції потрібні дані про попередні точки, то виникають підвищені вимоги до об’єму пам’яті ЕОМ.
3) Однокрокові методи і методи прогнозу і корекції можуть забезпечити приблизно однакову точність результатів. Однак, на відміну від перших другі дозволяють оцінити похибку на кроці. Тому при однокрокових методах значення кроку вибирають трохи меншим, ніж це потрібно, тому методи прогнозу і корекції є більш ефективними.
4) Застосовуючи метод Рунге – Кутта четвертого порядку точності, на кожному кроці потрібно обчислювати чотири значення функції. В той самий час для забезпечення збіжності методу прогнозу і корекції того ж порядку точності часто достатньо двох значень функції. Тому метод прогнозу і корекції потребує майже вдвічі менше машинного часу, ніж метод Рунге – Кутта такої ж точності.