МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ “ЛЬВІВСЬКА ПОЛІТЕХНІКА”
ОДНОКРОКОВІ МЕТОДИ ЧИСЕЛЬНОГО РОЗВ’ЯЗУВАННЯ ЗАДАЧІ КОШІ
ДЛЯ ЗВИЧАЙНИХ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ
МЕТОДИЧНІ ВКАЗІВКИ
з курсу "Чисельні методи"
для студентів базового напрямку
6.0802 "Прикладна математика"
Затверджено
на засіданні кафедри
“Прикладна математика”
Протокол № 7 від 31.03.2005 р.
Львів 2005
Однокрокові методи чисельного розв’язування задачі Коші для звичайних диференціальних рівнянь: Методичні вказівки з курсу «Чисельні методи» для студентів базового напрямку 6.0802 «Прикладна математика»/ Укл.: М.В.Кутнів, Я.В.Пізюр. – Львів: Видавництво Національного університету «Львівська політехніка», 2005.- 18 с.
Укладачі Кутнів М.В., канд. фіз-мат. наук, доц.,
Пізюр Я.В., канд. фіз-мат. наук, доц.
Відповідальний за випуск Мединський І.П., канд. фіз-мат. наук, доц.
Рецензенти Гнатів Б.В., канд. фіз-мат. наук, доц.,
Демків І.І., канд. фіз-мат. наук, доц.
ТЕОРЕТИЧНА ЧАСТИНА
§1. Задача Коші для звичайних диференціальних рівнянь
Нехай на відрізку необхідно знайти розв'язок задачі:
(1)
де
або в покомпонентному вигляді:
, (2)
. (3)
Зауважимо, що систему звичайних диференціальних рівнянь будь-якого порядку можна звести до системи рівнянь першого порядку (1).
Добре відомі умови, які гарантують існування та єдиність розв'язку задачі (2), (3). Нехай функції неперервні в замкнутій області З неперервності функцій випливає їх обмеженість, тобто існування константи такої, що скрізь у виконуються нерівності . Крім того, припустимо, що в функції задовольняють умову Ліпшиця, тобто
для будь-яких точок , області . Якщо умови, накладені на виконуються, то існує єдиний розв'язок системи (2), визначений при , і такий, що при приймає задані початкові умови (3).
Методи розв'язування звичайних диференціальних рівнянь можна розділити на точні, наближені і чисельні. До точних відносять методи, які дозволяють виразити розв'язок диференціального рівняння через елементарні функції, або зобразити його за допомогою квадратур від елементарних функцій. Однак класи рівнянь, для яких розроблені методи знаходження точних розв'язків, порівняно вузькі і охоплюють дуже малу частину задач, що виникають на практиці. Наближеними називають методи, в яких розв'язок одержується як границя деякої послідовності , причому елементи цієї послідовності виражаються через елементарні функції або за допомогою квадратур. Обмежуючись скінченним числом , одержимо наближений вираз для . Прикладом може бути метод розкладу розв'язку в степеневий ряд або метод послідовних наближень Пікара. Ці методи також можуть бути застосовані лише для порівняно простих задач (таких, як лінійні). Чисельні методи – це алгоритми обчислення наближених значень шуканого розв'язку на деякій вибраній сітці значень аргументу. Розв'язок при цьому одержується у вигляді таблиці. Чисельні методи не дозволяють знайти загальний розв'язок; вони можуть дати тільки який-небуть частинний розв'язок, наприклад, розв'язок задачі Коші (1), (2). Це основний недолік чисельних методів. Однак, ці методи можна застосувати до дуже широкого класу рівнянь і всіх типів задач для них. Тому з появою ЕОМ чисельні методи стали основним засобом розв'язування конкретних практичних задач.
Для розв'язування задачі Коші будемо використовувати чисельні методи. Виберемо на відрізку деяку, взагалі кажучи, нерівномірну сітку значень аргументу з кроком так, щоб для вузлів сітки виконувалися співвідношення . Будемо позначати через точний розв'зок задачі, а через – наближений розв'язок. Зазначимо, що наближений розв'язок є сітковою функцією, тобто визначений тільки в точках сітки. Чисельні методи розв'язування звичайних диференціальних рівнянь розділяють на два класи: однокрокові (методи рядів Тейлора, методи Рунге-Кутта) та багатокрокові (методи Адамса, формули диференціювання назад та ін.). Однокрокові методи для обчислення розв'язку в точці на кожному кроці використовують лише значення , тоді як багатокрокові використовують значення наближеного розв'язку в кількох попередніх точках.
§2. Метод рядів Тейлора
Для простоти викладення будемо надалі розглядати одне диференціальне рівняння
(1)
розв'язок якого задовольняє початкову умову
(2)
Якщо функція аналітична в точці , то розв'язок задачі (1), (2) можна розкласти у ряд Тейлора. Обмежуючись скінченним числом членів ряду, запишемо наближену рівність:
(3)
Для знаходження похідних, які стоять у правій частині (3), послідовно продиференціюємо рівняння (1)
(4)
Похибка наближеної рівності (3) за умови, що величина менша від радіуса збіжності ряду Тейлора, прямує до нуля при .
На практиці застосовують покроковий варіант методу рядів. Виберемо на відрізку нерівномірну сітку з кроком . Будемо вважати, що наближений розв'язок в точці вже знайдено, тобто відоме значення . Для побудови однокрокового чисельного методу знаходження розв'язку в наступній точці достатньо у формулі (3) покласти , а за взяти . Тоді одержимо
(5)
де обчислюються за формулами (4). Якщо б значення збігалося зі значенням , то похибка від заміни на мала б порядок . Використовуючи пакети програм аналітичного диференціювання функцій, за допомогою формул (5), (4) можна послідовно при заданому обчислити значення . Однак застосовувати для розрахунків формулу (5) з великим числом членів не вигідно, тому що навіть при порівняно простій правій частині вирази для похідних можуть бути громіздкими.
§3. Методи Рунге-Кутта
У найпростішому випадку, обмежуючись тільки першими двома членами розкладу (5) §2, одержимо метод Ейлера (сіткову схему Ейлера)
. (1)
Метод Ейлера дуже простий для реалізації на ЕОМ: на -ому кроці обчислюється значення , яке потім підставляється в (1).
Основне питання при використанні чисельних методів полягає в оцінці точності наближених значень . Величину називають похибкою (глобальною похибкою) чисельного методу, тоді як величину , де – чисельний розв'язок, одержаний при точному значенні називають локальною похибкою. Кажуть, що метод збіжний в точці , якщо при . Метод збіжний на відрізку , якщо він збіжний в кожній точці . Кажуть, що метод має -й порядок точності, якщо існує число таке, що при .
Підставимо в (1) і поділимо одержану рівність на , тоді будемо мати
.
Оскільки за теоремою про скінченні прирости
то
.
Функція називається похибкою апроксимації або нев'язкою сіткового рівняння на розв'язку вихідного рівняння (1) §2.
Кажуть, що чисельний метод апроксимує вихідне диференціальне рівняння, якщо або при . Метод має -й порядок апроксимації, якщо або .
Встановимо порядок апроксимації схеми Ейлера (1). Для цього розкладемо в ряд Тейлора в околі точки
Оскільки , то . Отже, метод Ейлера має перший порядок апроксимації.
Покажемо, що схема Ейлера збіжна при і має перший порядок точності. Доведення проведемо припускаючи, що для всіх
і величина кроку – стала. Тоді
тобто і наближений розв'язок збіжний до точного з першим порядком точності.
Для одержання точніших розрахункових формул обчислимо проміжне значення
, використовуючи схему Ейлера
(2)
а потім обчислимо за формулою
(3)
Формули (2), (3) називають методом прогнозу-корекції, оскільки на першому етапі (2) наближене значення розв'язку прогнозується, а на другому етапі (3) прогнозоване значення коректується. Цей метод можна реалізувати інакше. А саме, спочатку обчислимо послідовно функції
, (4)
а потім знайдемо
(5)
Така форма реалізації схем (2), (3) називається методом Рунге-Кутта. Оскільки вимагається обчислити дві проміжні функції , то даний метод називають двоступеневим. Далі ми покажемо, що цей метод має другий порядок апроксимації.
У загальному випадку явні -ступеневі методи Рунге-Кутта мають вигляд:
,
де –дійсні коефіцієнти. Компактне зображення методу Рунге-Кутта дає таблиця:
Переважно коефіцієнти задовольняють такі умови
(6)
зміст яких полягає в тому, що всі точки, в яких обчислюється , є наближеннями першого порядку до розв'язків Ці умови значно спрощують вивід методів Рунге-Кутта високого порядку.
Коефіцієнти вибираються з міркувань точності. Локальна похибка -ступеневого методу Рунге-Кутта дорівнює
(7)
де
Зазначимо, що
де –похибка апроксимації або нев'язка методу.
Зупинимося детальніше на окремих методах. При одержимо схему Ейлера, а при – множину методів:
Дослідимо локальну похибку двоступеневого методу залежно від вибору параметрів, припускаючи достатню гладкість розв'зку і функції . Для цього розкладемо всі величини, які входять до виразу локальної похибки (7), за формулою Тейлора в околі точки . Запишемо спочатку розклад за степенями значення
Використовуючи співвідношення (2.4), одержимо
(8)
Тут значення функції та її частинних похідних беруться при .
Розкладемо тепер як функцію за формулою Тейлора
,
Знайдемо похідні
.
Враховуючи, що одержимо
.
Тоді
, (9)
і
Прирівнюючи до нуля коефіцієнти при і , одержимо систему рівнянь, яку повинні задовольняти коефіцієнти двоступеневого методу для того, щоб цей метод мав другий порядок апроксимації
Зокрема, при будемо мати метод другого порядку (4), (5). Оскільки головний член (перший ненульовий член) тейлорівського розкладу локальної похибки
(10)
то звідси випливає, що двоступеневого методу третього порядку апроксимації не існує.
Щоб побудувати триступеневу схему третього порядку апроксимації, необхідно розкласти функції, які входять до (7) за формулою Тейлора до величин порядку включно. Аналогічно до , розкладемо за степенями :
.
Обчислимо
Звідси на підставі (6) будемо мати
.
Тоді
(11)
Враховуючи (8), (9), (11) локальна похибка триступеневого методу буде мати вигляд
(12)
Прирівняємо до нуля коефіцієнти при , тоді одержимо умови третього порядку апроксимації
Один з методів, який задовольняє ці умови, має вигляд
(13)
У випадку не існує формул четвертого порядку.
При не можна побудувати формул п'ятого порядку апроксимації. Наведемо таблицю коефіцієнтів найбільш вживаного чотириступеневого методу четвертого порядку
(14)
При певних припущеннях щодо гладкості правої частини рівняння (2.1), подібно до того як це було зроблено для схеми Ейлера, можна довести збіжність методів Рунге-Кутта.
§4. Практична оцінка похибки та вибір довжини кроку для методів Рунге-Кутта
Найбільш традиційним способом оцінки локальної похибки є правило Рунге (подвійний перерахунок, екстраполяція Річардсона), який полягає в повторенні обчислень із зменшеною вдвоє довжиною кроку і в порівнянні результатів.
Виходячи з точки , методом Рунге-Кутта -го порядку точності знайдемо чисельний розв'язок в точці . Тоді локальна похибка розв'язку буде мати вигляд
(1)
де виражається через коефіцієнти методу і частинні похідні правої частини (див. (10) §3), обчислені в точці . Припустимо, що в результаті двох кроків, проведених, цим самим методом, виходячи з цієї ж точки , обчислено і . Похибка розв'язку складається з двох частин: з локальної похибки першого кроку, яка має вигляд
і локальної похибки другого кроку, яка також виражається формулою (1), але з частинними похідними, обчисленими в точці . Отже,
. (2)
Якщо від (2) відняти (1), то отримаємо
.
Звідси
.
Тоді похибка (2) може бути обчислена за формулою
, (3)
а вираз
(4)
апроксимує величину з порядком .
Розглянемо алгоритм автоматичного вибору довжини кроку так, щоб локальна похибка не перевищувала допустимої точності . На основі формули (3) обчислюється похибка
,
де –масштабний множник. Для обчислення абсолютної похибки кладуть , а для відносної . Можна використовувати змішане масштабування типу
або
.
Потім величина порівнюється з заданою величиною допустимої похибки . Новий крок вибирається з умови
Звідси
.
Коефіцієнт дозволяє врахувати вплив величини . Тоді, якщо , обидва обчислених кроки вважаються успішними і розв'язування продовжується, виходячи з або , причому довжина нового кроку вибирається рівною . У протилежному випадку обидва кроки відкидаються і обчислення повторюються з новою довжиною кроку . Щоб запобігти занадто сильного збільшення кроку, максимальний коефіцієнт збільшення кроку , як правило, вибирають між 1,5 і 5. Після неуспішного кроку рекомендується не збільшувати розмір кроку.
Другий підхід до оцінки локальної похибки полягає в тому, щоб побудувати такі формули Рунге-Кутта, які містили б, крім розв'язку , і розв'язок більш високого порядку . Тобто, необхідно знайти таку таблицю коефіцієнтів
щоб величина
мала порядок , а
–порядок . Такі методи називають вкладеними. Для вкладених методів
.
Виведемо вкладені формули другого і третього порядків апроксимації, таблиця яких має вигляд
.
З (3.11) та аналогічного розкладу для випливає, що для збереження відповідних порядків апроксимації, повинні задовольнятися рівняння
Вибравши і , з перших двох рівнянь одержимо . Якщо покласти , то , . Одержаний метод наведений в таблиці
(5)
В основу ефективної програми розв'язування задачі Коші для звичайних диференціальних рівнянь лягли вкладені методи четвертого і п'ятого порядків Дорманта і Прінса (див. [23]), які визначаються коефіцієнтами:
(6)
Зазначимо, що застосування цього методу вимагає обчислення лише шести правих частин диференціального рівняння, оскільки збігається з для наступного кроку.
§5. Алгоритм чисельного розв’язування задачі Коші з заданою точністю
та автоматичним вибором величини кроку
Алгоритм обчислює чисельний розв’зок задачі Коші для звичайних диференціальних рівнянь методами Рунге-Кутта з заданою точністю на відрізку . Для практичної оцінки похибки використовується правило Рунге або вкладені методи Рунге-Кутта (див. §4). Крім того для задач з відомим точном розв’язком обчислюється норма похибки чисельного розв’язку
.
Для обчислення чисельного розв’язку в точці сітки необхідно мати на увазі, що при значеннях права частина диференціального рівняння може бути, взагалі кажучи не визначена, що може привести до зупинки (або переривання) ПЕОМ. Для запобігання такої ситуації треба поступити так. Після того як одержано чисельний розв’язок у вузлі і вибрано крок , з яким буде обчислюватися розв’язок у вузлі обчислимо величину . Якщо , де — машинна точність (машинний “шум”), то можна вважати, що і на цьому процес чисельного інтегрування завершується. При інтегрування повинно продовжуватися. Якщо при цьому , наступний крок не виводить за межі інтегрування і він виконується звичайним чином. Якщо ж , за крок слід взяти величину .
5.1. Алгоритм чисельного розв’язування задачі Коші методами Рунге-Кутта з оцінкою точності за правилом Рунге
1. Ввести значення , , , , , .
2. Ініціалізувати змінні ; ; ; .
3. If then go to 13.
4. If then .
5. ; .
6. .
7. Обчислити .
8. Обчислити
9. If then
begin ; ; ; ; go to 8 end.
10. If then
begin ; ; go to 7 end.
11. Обчислити
, .
12. If then
begin ; ; обчислити точний розв’язок ;
вивести значення ; ; ; ;
if then ;
go to 3
end;
else
begin ; ; ; go to 6 end.
13. Вивести норму похибки .
14. End.
5.2. Алгоритм чисельного розв’язування задачі Коші вкладеними
методами Рунге-Кутта
1. Ввести значення , , , , , .
2. Ініціалізувати змінні ; ; ; .
3. If then go to 10.
4. If then .
5. ; .
6. Обчислити .
7. Обчислити
8. Обчислити , .
9. If then
begin ; обчислити точний розв’язок ; вивести значення , , ; ;
if then ;
go to 3 end
else begin ; ; ; go to 7 end.
10. Вивести норму похибки .
11. End.
§7. Текст програми в середовищі Maple для знаходження точного розв‘язку задачі
> restart;
Присвоєння змінним t0 і T значень початку і кінця інтервалу.
>
> t0:=1; T:=6;
Задання виду функцій → f1, → f2, →g.
> f1:=exp(-(t+2)*t);
> f2:=exp(-3*t);
> g:=2*(t-1);
Заданняння вигляду диференціального рівняння zdr та початкових умов до нього pu.
> zdr:={diff(y(t),t)-f1*f2+g*y(t)=0};
>
> pu:={y(t0)=10};
Розв'язування диференціального рівняння.
> s:=dsolve(zdr union pu,{y(t)});
>
>
> eval(y(t),s);
>
Побудова графіка отриманого розв'язку диференціального рівняння.
>
> plot(%,t=t0..T);
>
Обчислення значень розв'язку диференціального рівняння на краях інтервалу.
> y1:=unapply(%%,t);
>
>
> evalf(y1(t0));
>
> evalf(y1(T));
>
Послідовність виконання лабораторної роботи
1. Одержати варіант завдання.
2. Вивчити відповідний лекційний матеріал.
3. Знайти точний розв’язок поставленої задачі.
4. Використовуючи будь-яку з відомих Вам мов програмування, написати та відлагодити програму чисельного розв’язування задачі Коші для ЗДР методом Рунге-Кутта з автоматичним вибором кроку (див. §4) при заданій точності та .
Зміст звіту
Постановка задачі (конкретний варіант).
Аналітичний розв’язок задачі.
Алгоритм для конкретного чисельного методу.
Текст програми.
Результати розв′язування на ЕОМ заданої задачі.
Список літератури
Бахвалов Н.С., Жидков Н.П., Кобельков Г.М.. Численные методы.-М.:Наука, 1987.
Гаврилюк І.П., Макаров В.Л. Методи обчислень. –К.:Вища школа, 1995, ч.1, ч.2.
Данилович В., Кутнів М. Чисельні методи.-Львів:Кальварія, 1998.
Калиткин Н.Н. Численные методы.-М.:Наука, 1978.
Самарский А.А., Гулин А.В. Численные методы. - М.:Наука, 1989.
Самарский А.А. Теория разностных схем.- М.:Наука, 1989.
Трифонов Н.П., Пасхин Е.Н. Практикум работы на ЭВМ.-М.: Наука, 1982.
Хайрер Э., Нерсет С., Ваннер Г. Решение обыкновенных дифференциальных уравнений.-М.:Мир, 1990.
ВАРІАНТИ ЗАВДАНЬ
Вибрати один з методів 2-го, 3-го або 4-го (див. §3, формули (4),(5); (13); (14)) порядків, для практичної оцінки похибки та вибору кроку використати правило Рунге.
Вибрати пару вкладених методів Рунге Кутта (див. §4, формули (5) або (6))та відповідний алгоритм оцінки похибки та вибору величини кроку.
Задано задачу Коші для ЗДР
.
а) , функції :
№ варіанта
1
2
3
4
5
6
7
8
9
б)
№ варіанта
10
11
12
13
14
НАВЧАЛЬНЕ ВИДАННЯ
ОДНОКРОКОВІ МЕТОДИ ЧИСЕЛЬНОГО РОЗВ’ЯЗУВАННЯ ЗАДАЧІ КОШІ ДЛЯ ЗВИЧАЙНИХ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ
МЕТОДИЧНІ ВКАЗІВКИ
з курсу “Чисельні методи”
для студентів базового напрямку
6.0802 “Прикладна математика”
Укладачі Кутнів Мирослав Володимирович
Пізюр Ярополк Володимирович
Редактор
Комп’ютерне складання
Підписано до друку
Формат 70 х 1001/16. Папір офсетний.
Друк на різографі. Умовн. друк. арк. Обл.-вид. арк.
Наклад прим. Зам.
Поліграфічний центр
Видавництва Національного університету “Львівська політехніка”
вул. Ф Кодесси, 2, 79000, Львів