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

Метод Ейлера

Звичайним диференціальним рівнянням називають таке рівняння, що зв’язує незалежну змінну х, невідому функцію у(х) і її похідні у´, у´´, …, у⁽ⁿ⁾. Загальний вигляд звичайного диференціального рівняння F(x, y, у´, у´´, …, у⁽ⁿ⁾) = 0. Його порядок – це порядок вищої похідної, що входить до нього. Розв’язком диференціального рівняння є дозвільна функція у(х), підстановка якої у це рівняння перетворює його у тотожність.

Наприклад, у´ = cos x – звичайне диференціальне рівняння першого порядку. Воно еквівалентне твердженню, що функція у(х) є первісною cos x і отже у(х) = = sin x + C, де С – довільна стала. Цей приклад з’ясовує традиційну термінологію: звичайно розв’язування диференціального рівняння ще називають інтегруванням, його розв’язок – інтегралом, а графік розв’язку – інтегральною кривою. Ще один приклад: у´ = 2у – це теж звичайне диференціальне рівняння першого порядку, його розв’язок у(х) = Се^2х, де С – довільна стала. В цьому можна переконатись прямою підстановкою згідно з означенням розв’язку. В обох прикладах розв’язок залежить не тільки від х, але й від довільної сталої. Розв’язок виявиться єдиним, якщо він додатково задовольняє початкову умову у(х₀) = у₀. Так у першому прикладі, якщо у(0) = 0, то у(х) = sin x – єдиний розв’язок. Якщо у(0) = 1, то у(х) = е^2х – це єдиний розв’язок у другому прикладі.

Означення 1. Задача Коші для звичайного диференціального рівняння першого порядку

у´ = f (x, y); у(х₀) = у₀ (1)

полягає в тому, щоб знайти розв’язок у(х) цього рівняння, який задовольняє початкову умову у(х₀) = у₀.

Геометрично це означає, що треба знайти ту інтегральну криву у(х) рівняння у´ = f (x, y), яка проходить через точку (х₀; у₀).

Теорема (існування та єдності розв’язку задачі Коші для звичайного диференціального рівняння першого порядку). Якщо функція f (x, y) і її часткова похідна (x, y) неперервні у точці (х₀; у₀), то існує єдиний розв’язок у(х) рівняння у´ = f (x, y), який задовольняє початкову умову у(х₀) = у₀.

Отже, на відміну від квадратного рівняння диференціальне рівняння практично завжди має розв’язок і при тому єдиний. Насправді, проінтегрувати диференціальне рівняння у скінченому вигляді, “на папері” вдається дуже рідко, лише для певних спеціальних класів рівнянь. Фактично здебільшого застосовують наближене інтегрування. Мабуть найкращий шлях зрозуміти теорему існування та єдності диференціального рівняння – це розібрати найстаріший та найпростіший метод його чисельного інтегрування – метод Ейлера.

Означення 2. 1) Розв’язати задачу Коші (1) чисельно означає для заданої послідовності х₀, х₁, …, х_n значень незалежної змінної х і числа у₀ знайти числову послідовність у₀, у₁, …, у_n так, щоб у_k з заданою точністю наближав у(х_k) для всіх k, де у(х) – єдиний розв’язок задачі Коші з початковою умовою у(х₀) = у₀.

2) Різницю у_k – у(х_k) називають похибкою наближеного значення у_k в точці х_k.

3) Якщо всі точки х₀, х₁, …, х_n рівновіддалені: х_k = х₀ + kh, то величину h називають кроком інтегрування диференціального рівняння.

Спочатку зауважимо, що оскільки кутовий коефіцієнт дотичної до інтегральної кривої у(х) у точці (х_k; у_k) дорівнює у´ = f (х_k; у_k), то рівняння дотичної у (х_k; у_k) матиме вигляд у – у_k = f (х_k; у_k)(х – х_k). Геометричний зміст диференціального рівняння у´ = f (x, y) в тому, що воно визначає напрям інтегральної кривої в кожній точці площини в деякому околі початкової точки (х₀; у₀).

Отже, рівняння дотичної у початковій точці М₀(х₀; у₀) має вигляд у – у₀ = f (х₀; у₀)(х – х₀). Зафіксуємо крок інтегрування h і на цій дотичній знайдемо точку М₁(х₁; у₁) з х₁ = х₀ + h. Звідси у₁ = у₀ + f (х₀; у₀)h. Рівняння дотичної у М₁(х₁; у₁) – це у – у₁ = f (х₁; у₁)(х – х₁). На цій дотичній точка М₂(х₂; у₂) з абсцисою х₂ = х₁ + h = х₀ + 2h має ординату у₂ = у₁ + f (х₁; у₁)h, …. На дотичній у точці М_k(х_k; у_k), де х_k = х₀ + kh обираємо точку М_k+1(х_k+1; у_k+1), де х_k+1 = х_k + h = х₀ + (k+1)h, звідки у_k+1 = у_k + f (х_k; у_k)h. В результаті дістанемо ламану, яка сполучає точки М_k(х_k; у_k) відрізками дотичних і наближено зображує інтегральну криву шуканого розв’язку y(x) задачі Коші. Її назва – ламана Ейлера (див. рисунок 1).

y

y(x)

М₂ O(h)

O(h²) М₁ М₃ М₈

М₇

М₀ М₄ – ламана Ейлера

y₀ М₆

М₅

 

x

x₀ h h h h h h h h

 

Рисунок 1

 

Отже, метод Ейлера – це ітераційний метод, рекурентна формула якого має вигляд:

x_k+1 = x_k + h, у_k+1 = у_k + f (х_k; у_k)h (k = 0, 1, 2, …, n – 1). (2)

Зазначимо, що оскільки ламана Ейлера складається з відрізків дотичних, то похибка наближеного значення у_k на кожному кроці k є величина порядку O(h²). З переходом від точки до точки ця похибка систематично зростає і на відрізку завдовжки одиниця сумарна похибка є величиною порядку O(h), оскільки ітерацій з кроком інтегрування h на такому відрізку порядку . Коли h ® 0, то ламана Ейлера з початковою умовою у(х₀) = у₀ збігається до інтегральної кривої єдиного розв’язку задачі Коші для диференціального рівняння з такою початковою умовою.

Очевидною причиною похибки на кожному кроці є те, що нахил інтегральної кривої неперервно змінюється в той час як нахил ламаної Ейлера є незмінним на кроці. Тому можна сподіватись, що похибка дещо зменшиться, якщо на кожному кроці нахил відрізку дотичної буде співпадати з нахилом інтегральної кривої не на початку, а в середині кроку. Точніше кажучи, спочатку за формулою (1) знайти точку М(х_{k + ½}; у_{k + ½}) посередині відрізку М_kМ_k+1: х_k+½ = х_k + h, у_{k + ½} = у_k + h f (х_k; у_k). А потім на [х_k; x_k+1] взяти відрізок прямої L, яка проходить через точку М_k(х_k; у_k), але ж її кутовий коефіцієнт дорівнює f (х_{k + ½}; у_{k + ½}) (див. рисунок 2).

y M_k+1

L

y_{k + ½} M_k+½

 

 

y_k

M_k

 

x

x_k x_{k + ½} x_{k + 1}

 

Рисунок 2.

 

Такий дещо модифікований алгоритм називають удосконаленим методом Ейлера. Отже, його рекурентну формулу можна записати так:

у_{k + ½} = у_k + f (х_k; у_k); у_k+1 = у_k + f (х_k + ; у_{k + ½})h (k = 0, 1, 2, …, n – 1). (3)

Чи справді таке удосконалення призведе до помітного зменшення похибки? Несподівано виявляється, що похибка наближеного значення у_k на кожному кроці k є величиною порядку O(h³) і отже на відрізку завдовжки одиниця сумарна похибка є величиною порядку O(h²) (на відміну від О(h) для методу Ейлера).