Розділ 8. Числове інтегрування звичайних диференційних рівнянь та систем диференційних рівнянь

Короткі теоретичні відомості

Для ЗДР першого порядку задача Коші полягає в знаходженні такого розв’язку рівняння:

, (8.1)

що задовольняє початкову умову

, (8.2)

де – задана неперервна функція двох аргументів, – задана константа.

Для диференціального рівняння m -го порядку

(8.3)

задача Коші полягає в знаходженні функції y = y (t), що задовольняє рівняння (8.3) і початкові умови:

Розв’язуючи деякі задачі, слід знайти m функцій y ₁(t), y ₂(t),…, y_m (t), що задовольняють систему диференціальних рівнянь:

(8.4)

…

і початкові умови

(8.5)

Систему диференціальних рівнянь (8.4) називають нормальною.

Задачу Коші (8.4) з початковими умовами (8.5) можна подати у векторній формі

(8.6)

де

Систему диференціальних рівнянь

у якій – функції від t, називають лінійною.

Якщо ввести позначення

,

то систему лінійних диференціальних рівнянь записують у матричній формі:

.

Рівняння (8.3) за допомогою заміни змінних

можна звести до вигляду (8.6). Наприклад, рівняння

можна записати як систему

При чисельному розв’язку задачі Коші по змінній t вводять сітку 0< t ₁< t ₂<...< t_n та шукають значення невідомої функції в вузлах сітки. Позначимо як y (t) — точний розв’язок задачі (8.1), (8.2), а y_{i +1} — наближений розв’язок в точці t = t_{i +1}. Локальною похибкою розв’язку називають

(8.7)

При використанні наближених методів основним є питання про збіжність. Найбільше розповсюдження отримало поняття про збіжність при h_i ®0, де

h_i = t_{i +1} – t_i.

Говорять, що метод збігається в точці t_{i +1}, якщо . Метод збігається на відрізку [0, t_n ], якщо він збігається в кожній точці t Î[0, t_n ].

Говорять, що метод має p -ий порядок точності, якщо існує таке число p >0, що , при h_i ®0.

 

Метод Ейлера

Для чисельного розв’язання задачі Коші (8.1), (8.2) розіб’ємо інтервал [0, t_n ], на якому шукається розв’язок, на елементарні відрізки [ t_i, t_{i +1}], , t ₀=0. На кожному елементарному відрізку будемо шукати розв’язок у вигляді полінома першого ступеня. Нехай значення функції на границях елементарного відрізку становить y_i, y_{i+ 1}. За інтерполяційною формулою Лагранжа

.

тоді

. (8.8)

Підставляючи (8.8) в (8.1) та замінюючи функцію f (t, y) її значенням на початку елементарного відрізку, отримаємо

. (8.9)

В методі Ейлера рівняння (8.1) замінюється різницевим рівнянням (8.9). З (8.9) випливає, що наближений розв’язок задачі Коші знаходиться явно за рекурентною формулою

(8.10)

Оцінимо локальну похибку методу Ейлера. Для цього припустимо, що в точці t = t_i значення функції y (t) обчислено точно, тобто y (t_i)= y_i. Розкладаючи функцію y (t) в ряд Тейлора в околі точки t = t_i і враховуючи (8.1) отримаємо

Звідси

(8.11)

Враховуючи (8.10), вважаючи крок інтегрування h_i малим та нехтуючи малими величинами вищих порядків із (8.7) та (8.11) отримуємо

.

Таким чином, метод Ейлера має перший порядок точності. Тому в методі Ейлера зменшення кроку h_i в два рази в чотири рази зменшує локальну похибку розв’язку.

Існують інші різновиди методу. Так, підставляючи (8.8) до (8.1) і замінюючи функцію f (t, y) її значення у кінці елементарного відрізку, отримаємо

.

Звідси

. (8.12)

Для довільної функції f (t, y) значення y_{i +1} не може бути виражене з (8.12) явно. Тому метод (8.12) називають неявним методом Ейлера. При використанні неявного методу Ейлера нове значення y_{i +1} визначається на основі попереднього y_i шляхом розв’язання нелінійного рівняння (8.12). Неявний метод Ейлера також має перший порядок точності. Можна показати, що похибка неявного методу Ейлера дорівнює

.

Наприклад: Нехай дано задачу Коші

.

Точним розв’язком цієї задачі є функція .

Необхідно обчислити таблицю значень y (t) на інтервалі з кроком h =0,25, тобто для t ₀=0, t ₁=0,25, t ₂=0,5, t ₃=0,75, t ₄=1.

Скористаємося явним методом Ейлера. Оскільки f (t, y)=- y ²(t +1), то згідно з формулою (8.10)

.

Тоді

Розв'яжемо цю ж задачу неявним методом Ейлера. Згідно з формулою (8.12) маємо

,

що приводить до рівняння відносно y_{i +1}:

.

Тоді

;

Наприклад: Проінтегрувати систему ЗДР з початковими умовами на відрізку (0; 0,5) з точністю до 10^-3

Розв’язання:

Використовуючи метод Ейлера з кроком h= 0,1,обчислення будемо проводити з чотирма десятинними знаками:

,

і т.д.

Отже, занесемо значення до таблиці:

x	y	z
0,1	0,01	1,1
0,2	0,0222	1,2210
0,3	0,038	1,365
0,4	0,0576	1,5353
0,5	0,0828	1,73459

Відповідь: y=[0; 0,01; 0,0222; 0,038; 0,0576; 0,0828]; z=[1; 1,1; 1,221; 1,365; 1,535; 1,735].

 

Методи Рунге – Кутта

У найзагальнішому випадку метод Рунге-Кутта порядку точності p будується за рекурентною формулою:

, (8.13)

де

, (8.14)

a ₁=0, a_j и c_j – константи, b_jl – елементи нижньотрикутної матриці, такої, що кожне k_j отримується з попередніх значень k_l.

Формули (8.13) та (8.14) містять коефіцієнтів b_jl, a_j та c, що підлягають визначенню. Для того, щоб знайти ці коефіцієнти всі функції k_j розкладають в ряд Тейлора в околі точки (t_i, y_i). Ці розклади підставляють в (8.13) і результат зрівнюють з рядом Тейлора для функції y (t_{i +1}). Оскільки локальна похибка визначена згідно з (8.7), то ставиться умова, щоб коефіцієнти при всіх в розкладах для y (t_{i +1}) и y_{i +1} були рівні. Ця вимога приводить до системи рівнянь відносно коефіцієнтів b_jl, a_j и c_j.

Як приклад розглянемо побудову методів Рунге-Кутта другого порядку точності (p =2). В цьому випадку можна отримати розв’язок для коефіцієнтів b_jl, a_j и c_j при m =2. Тоді формула (8.13) перетвориться:

, (8.15)

де , .

Розкладемо k ₁ і k ₂ в ряд Тейлора в околі точки (t_i, y_i). Для цього відмітимо, що

.

тоді

, (8.16)

(8.17)

Підставляючи (8.16) та (8.17) в (8.15), отримаємо

.

Розкладемо функцію y (t) в ряд Тейлора в околі точки (y_i, t_i) та знайдемо y (t_{i +1}). Тоді

. (8.18)

Врахуємо, що

, (8.19)

(8.20)

Підставляючи (8.9) та (8.20) в (8.18), маємо

(8.21)

Для того, щоб метод мав другий порядок точності (), вимагатимемо, щоб в різниці рівнянь (8.18) та (8.21) були відсутні доданки, які містять , и . Виходячи з цієї умови коефіцієнти a, b, c ₁ и c ₂ повинні задовольняти системі рівнянь:

(8.22)

Рівняння (8.22) складають систему з 3 рівнянь відносно 4 невідомих. Виразимо b, c ₁ и c ₂ через a ¹0. Отримаємо:

. (8.23)

Підставляючи (8.23) в (8.18) отримаємо такий узагальнений метод Рунге-Кутта другого порядку точності

, (8.24)

де , .

Широковідомим методом 2 порядку є окремий випадок методу (8.24) при :

,

де , .

Аналогічно будуються розрахункові формули методів Рунге-Кутта вищих порядків точності. При цьому відповідні системи рівнянь відносно b_jl, a_j та c_j для p =2,3,4 можуть бути розв’язані відповідно для m =2,3,4. Таким чином, для того, щоб побудувати метод порядку точності p для p =2,3 та 4 достатньо m = p викликів функції f (t, y). Проте для p =5 ця система може бути розв’язана тільки при m ³6. Таким чином, для методу Рунге-Кутта 5 порядку точності потрібно щонайменше 6 викликів функції на кожному кроці. Аналогічно для p =6 m ³7. Це, зокрема, пояснює популярність класичного методу четвертого порядку, оскільки для того, щоб отримати додатковий порядок точності додатково потрібні два обчислення функції.

Серед всього різноманіття методів Рунге-Кутта 4 порядку слід також відмітити найпоширенішу модифікацію методу з найбільш простим набором коефіцієнтів b_jl, a_j та c_j:

(8.25)

Проте, слід зазначити, що використання формул (8.25) не дає змогу побудувати досить ефективні алгоритми розв’язку задачі Коші і його слід використовувати тільки для напів-аналітичних розрахунків або на початкових етапах розробки програмного забезпечення для отримання первинних оцінок.

Наприклад:

;

.

t ₀=0;

;

Методи Рунке-Кутта можуть бути застосовані також до розвязку систем ЗДР. При цьому розвязок системи ЗДР зводиться до розвязку на кожному кроці систем нелінійних рівнянь.

Наприклад: Розв’язати попередню систему диференційних рівнянь методом Рунге-Кутта:

Розв’язання:

Так як , знаходимо :

Отже,

;

;

таким же чином знайдемо y₂^(0,1) i z₂^(0,1) в точці x₂^(0,1)= 0,2; y₂^(0,1)= 0,0252; z₂^(0,1)= 1,2451

Обчислимо значення y^(0,2) та z^(0,2) в точці x₁^(0,2)= 0,2,але з кроком h= 0,2

y₁^(0,2)= 0,0253

z₁^(0,2)= 1,2453

Похибка значень y₁^(0,2) та z₁^(0,2) дорівнює:

Продовжуючи обчислення ми маємо:

x		0,1	0,2	0,3	0,4	0,5
y		0,011	0,025	0,044	0,067	0,097
z		1,111	1,245	1,405	1,595	1,818

Відповідь: y=[0; 0,011; 0,025; 0,044; 0,067; 0,097]; z=[1; 1,111; 1,245; 1,405; 1,595; 1,818].

Аналогічні схеми можна побудувати на основі методів більш високого порядку.

Методи Ейлера та Рунге – Кутта називаються одноточковими. Це пояснюється тим, що для обчислення результату в -й точці достатньо знати значення лише в одній попередній i -й точці.

 

Багатоточкові методи

В одноточкових методах значення обчислюють на основі значення в одній попередній точці . При цьому логічно припустити, що можна підвищити точність, якщо в розрахунку значення в точці на кожному кроці бере участь не одна, а кілька попередніх точок –

В загальному випадку лінійним m -кроковим різницевим методом з постійним кроком називається система різницевих рівнянь:

; (8.26)

де a_k, b_k — числові коефіцієнти, що не залежать від i, , причому a ₀¹0, f_k = f (t_k, y_k), h = t_{i +1}- t_i, .

Метод називається лінійним, тому що кожне f_k входить до формули лінійно.

Рівняння (8.26) є рекурентним співвідношенням, за яким нове значення y_{i+ 1} виражається через раніше знайдені y_i, y_{i- 1}, y_{i- 2},..., y_{i-m+ 1}. Обчислення починається з i = m– 1, тобто з рівняння

.

Звідси видно, що для початку обчислень слід задати m початкових значень y ₀, y ₁,... y_{m -1}. Значення y ₀ задається умовою (8.26). Величини y ₁, y ₂,... y_{m -1} можна обчислити використовуючи однокрокові методи, наприклад, методи Рунге-Кутта.

Метод (8.26) називається явним, якщо b ₀=0, і, відповідно, шукане значення y_{i +1} виражається через попередні y_i, y_i-1,..., y_{i-m +1}

(8.27)

Якщо b ₀¹0, то метод називається неявним. Тоді для знаходження y_{i +1} доводиться розв’язувати нелінійне рівняння

. (8.28)

Зазвичай це рівняння розв’язують методом Ньютона і задають початкове наближення y_{i +1}⁽⁰⁾, що дорівнює y_i.

Слід відзначити, що коефіцієнти рівняння (8.26) задані з точністю до множника. Щоб усунути цю неоднозначність, вимагають виконання умови нормування

. (8.29)

Можна показати, що для того, аби різницевий метод (8.26) мав порядок точності p, необхідно щоб були виконані умови:

; (8.30)

. (8.31)

Разом з умовою нормування (8.29) рівняння (8.30) і (8.31) становлять систему з p +2 лінійних алгебраїчних рівнянь відносно 2 m +2 невідомих .

Для того, щоб існував розв’язок цієї системи, необхідно, щоб кількість рівнянь була не більшою від кількості невідомих, тому p £2 m. Ця вимога означає, що лінійні m -крокові різницеві методі мають порядок точності не більше, ніж 2 m. Тому порядок точності неявних m -крокових методів дорівнює 2 m, коли a_i ¹0, b_i ¹0, , а явних — 2 m -1, оскільки b ₀=0 і невідомих в (8.29), (8.30), (8.31) на одну менше.

У практиці обчислень найбільше поширення отримали методи Адамса, які становлять частковий випадок багатокрокових методів (8.26), коли похідна y ¢ апроксимується тільки по двох точках t_{i +1}, t_i, тобто

.

Таким чином, методи Адамса мають вигляд

. (8.32)

У випадку b ₀=0 методи Адамса називають явними, а при b ₀¹0 — неявними.

Для методів Адамса умови (8.29), (8.31) приймають вигляд

; (8.33)

. (10.34)

Звідси видно, що найвищий порядок точності m -крокового неявного методу Адамса дорівнює m +1 (кількість рівнянь р, а невідомих — m +1), а найвищий порядок точності явного методу дорівнює m (b ₀=0).

Знайдемо коефіцієнти b_k для кількох m. Розглянемо явні методи Адамса. Коефіцієнти знаходять із (10.33), (10.34) при p = m. При m =1 отримуємо метод Ейлера

.

При m =2 (10.33) приймає вигляд

;

.

Звідси .

Таким чином, при m =2, отримаємо явний метод Адамса другого порядку

.

Звідси

.

Розв’язки для інших m наведено у таблиці 8.1.

Коефіцієнти неявного методу Адамса знаходяться із системи (8.33), (8.34) з р = m +1. Наведемо формули неявних -точкових методів Адамса для деяких значень m:

 

Таблиця 8.1. Коефіцієнти явного методу Адамса порядку m

m	b 1	b 2	b 3	b 4	b 5	b 6	b 7	b 8	b 9	b 10
	3/2	-1/2
	23/12	-4/3	5/12
	55/24	-59/24	37/24	-3/8
	1901/720	-1387/360	109/30	-637/360	251/720
	4277/1440	-2641/480	4991/720	-3649/720	959/480	-95/288
	28392/	-18637/	235183/	-10754/	135713/	-5603/	6054/
	16083/	-2281287/	242653/	-24605/	28060547/	-2384104/	32863/	-5257/
	2787217/	-294002/	1929453/	-7436324/	135961/	-11884/	261790/	50246/	-994/
	173234/	-17722/	1971119/	-4119436/	10176373/	-6963789/	843419/	-6835/	122585/	-130/

 

Для обчислювальна схема матиме вигляд:

Наприклад: Нехай дано задачу Коші

Необхідно обчислити таблицю значень y (t) на інтервалі , якщо h =0,25.

Застосуємо трьохточковий явний метод Адамса. Використовуючи таблицю 8.1 отримаємо загальну формулу для цього методу:

.

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

.

Після цього продовжимо розв’язок заданим методом: