Чисельні методи розв'язання звичайних диференціальних рівнянь

Чисельні методи для звичайних диференціальних рівнянь — чисельні методи для знаходження наближених розв'язків звичайних диференціальних рівнянь. Їх застосування також іноді називають «чисельним інтегруванням», хоча частіше цей термін вживають для обчислення інтегралів.

Ілюстрація чисельного інтегрування для диференціального рівняння ${\displaystyle y'=y,y(0)=1.}$ Розмір кроку становить ${\displaystyle h=1.0}$.

   Синій: Метод Ейлера

   Зелений: Метод середньої точки

   Червоний: Точний розв'язок

Та сама ілюстрація для ${\displaystyle h=0.25.}$ Метод середньої точки збігається швидше, ніж метод Ейлера, оскільки ${\displaystyle h\to 0}$.

Багато диференціальних рівнянь неможливо розв'язати точно. Однак для практичних цілей – наприклад, в інженерії – часто достатньо чисельного наближення до розв'язку. Тоді стають в нагоді алгоритми для обчислення такого наближення. Альтернативним методом є використання методів математичного аналізу для отримання розв'язку у вигляді ряду.

Звичайні диференціальні рівняння зустрічаються в багатьох наукових дисциплінах, зокрема у фізиці, хімії, біології та економіці^[1]. Крім того, деякі методи перетворюють диференціальні рівняння з частинними похідними на звичайні диференціальні рівняння, які потім теж необхідно розв'язати.

Задача

Диференціальне рівняння першого порядку — це задача з початковими умовами^[en] вигляду^[2]

${\displaystyle y'(t)=f(t,y(t)),\qquad y(t_{0})=y_{0},}$

(1)

де ${\displaystyle f}$ є функцією ${\displaystyle f:[t_{0},\infty )\times \mathbb {R} ^{d}\to \mathbb {R} ^{d}}$, а початкова умова ${\displaystyle y_{0}\in \mathbb {R} ^{d}}$ — якийсь заданий вектор. Перший порядок рівняння означає, що в рівнянні з'являється лише перша похідна від y, а вищі похідні відсутні.

Без втрати загальності для систем вищого порядку, ми обмежимося диференціальними рівняннями першого порядку, оскільки звичайне диференціальне рівняння вищого порядку можна перетворити на більшу систему рівнянь першого порядку шляхом введення додаткових змінних. Наприклад, рівняння другого порядку y′′ = −y можна переписати як два рівняння першого порядку: y′ = z та z′ = −y.

У цьому розділі ми описуємо лише чисельні методи для задач з початковими умовами та зазначаємо, що крайові задачі вимагають іншого набору інструментів. У крайових задачах значення y задані в більш ніж одній точці. Через це для розв'язання крайової задачі необхідно використовувати спеціальні методи, наприклад, метод стрільби (та його варіанти) або глобальні методи, такі як метод скінченних різниць^[3], методи Гальоркіна^[4] або методи колокації.

Теорема Пікара-Лінделёфа стверджує, що існує єдиний розв'язок, за умови, що f задовільняє критерію Ліпшиця.

Методи

Чисельні методи розв'язання задачі з початковими умовами першого порядку часто поділяють на дві великі категорії^[5]: лінійні багатокрокові методи та методи Рунге-Кутти. Подальший поділ можна здійснити, розділивши методи на явні та неявні. Наприклад, неявні лінійні багатокрокові методи включають методи Адамса — Моултона та методи зворотного диференціювання^[en], тоді як неявні методи Рунге — Кутти^[6] включають діагонально неявний метод Рунге — Кутти (англ. diagonally implicit Runge–Kutta, DIRK)^[7][8], однодіагонально неявний метод Рунге-Кутти (англ. singly diagonally implicit Runge–Kutta, SDIRK)^[9] та чисельні методи Гаусса-Радау^[10] (засновані на квадратурі Гауса^[11]). Явні приклади з лінійних багатокрокових методів включають методи Адамса — Башфорта, а будь-який метод Рунге — Кутти з таблицею Бутчера нижнього трикутного виду є явним. Грубе емпіричне правило вказує, що жорсткі диференціальні рівняння вимагають використання неявних схем, тоді як нежорсткі задачі можна ефективніше розв'язувати за допомогою явних схем.

Так звані загальні лінійні методи^[en] є узагальненням двох вищезгаданих великих класів методів^[12].

Метод Ейлера

${\displaystyle y'(t)\approx {\frac {y(t+h)-y(t)}{h}},}$

(2)

Для будь-якої точки кривої можна знайти наближене значення наступної точки на кривій, перемістившися на невелику відстань вздовж лінії, дотичної до кривої.

Починаючи з диференціального рівняння (1), замінимо похідну y′ наближеною скінченною різницею

${\displaystyle y(t+h)\approx y(t)+hf(t,y(t)).}$

(3)

Після перегрупування це дає наступну формулу:

${\displaystyle y(t+h)\approx y(t)+hy'(t).}$

Використання (1) дає:

${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}$

(4)

Цю формулу зазвичай застосовують таким чином. Вибираємо крок h та будуємо послідовність ${\displaystyle t_{0},t_{1}=t_{0}+h,t_{2}=t_{0}+2h,...}$ Ми позначаємо через ${\displaystyle y_{n}}$ чисельну оцінку точного розв'язку ${\displaystyle y(t_{n})}$. Керуючись (3), ми обчислюємо ці оцінки за такою рекурсивною схемою

${\displaystyle y'(t)\approx {\frac {y(t)-y(t-h)}{h}},}$

(5)

Цей метод Ейлера (або прямий Ейлера вперед, на відміну від зворотного методу Ейлера, який буде описано нижче) названий на честь Леонарда Ейлера, який описав його в 1768 році.

Метод Ейлера є прикладом явного методу. Це означає, що нове значення y_n+1 визначається з точки зору вже відомих величин, таких як y_n.

Зворотний метод Ейлера

Якщо замість (2) використовувати наближення

${\displaystyle y'(t)\approx {\frac {y(t)-y(t-h)}{h}},}$

(5)

отримуємо зворотний метод Ейлера:

${\displaystyle y_{n+1}=e^{-Ah}y_{n}+A^{-1}(1-e^{-Ah}){\mathcal {N}}(y(t_{n}))\ .}$

(8)

Зворотний метод Ейлера є неявним методом, тобто нам потрібно розв'язати рівняння, щоб знайти y_n+1. Для досягнення цієї мети часто використовують метод простої ітерації або деяку модифікацію методу Ньютона — Рафсона.

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

Метод експоненціального інтегратора першого порядку

Експоненціальні інтегратори описують великий клас інтеграторів, які нещодавно зазнали значного розвитку^[13]. Вони існують щонайменше з 1960-х років.

Замість (1) ми припускаємо, що диференціальне рівняння або має вигляд

${\displaystyle y'(t)=-A\,y+{\mathcal {N}}(y),}$

(7)

або його локально лінеаризували відносно фонового стану для отримання лінійного члена ${\displaystyle -Ay}$ і нелінійний члену ${\displaystyle {\mathcal {N}}(y)}$.

Експоненціальні інтегратори будують шляхом множення (7) на ${\textstyle e^{At}}$ та точного інтегрування результату на певному інтервалі часу ${\displaystyle [t_{n},t_{n+1}=t_{n}+h]}$:

${\displaystyle y_{n+1}=e^{-Ah}y_{n}+\int _{0}^{h}e^{-(h-\tau )A}{\mathcal {N}}\left(y\left(t_{n}+\tau \right)\right)\,d\tau .}$

Це інтегральне рівняння є точним, але воно не визначає інтеграл.

Експоненціальний інтегратор першого порядку можна реалізувати, вважаючи ${\displaystyle {\mathcal {N}}(y(t_{n}+\tau ))}$ сталим протягом усього інтервалу часу:

${\displaystyle y_{n+1}=e^{-Ah}y_{n}+A^{-1}(1-e^{-Ah}){\mathcal {N}}(y(t_{n}))\ .}$

(8)

Узагальнення

Метод Ейлера часто недостатньо точний. Точніше кажучи, він має лише 1-й порядок (поняття порядку пояснено нижче). Це спонукало математиків шукати методи вищого порядку.

Одна з можливостей полягає в тому, щоб використовувати не лише попередньо обчислене значення y_n для визначення y_n+1, але й зробити розв'язок залежним від більшої кількості минулих значень. Це призводить до так званого багатокрокового методу. Можливо, найпростішим є метод перекроковування, який має 2-й порядок та (грубо кажучи) спирається на два минулі значення.

Майже всі практичні багатокрокові методи належать до сімейства лінійних багатокрокових методів, які мають вигляд $${\displaystyle {\begin{aligned}&{}\alpha _{k}y_{n+k}+\alpha _{k-1}y_{n+k-1}+\cdots +\alpha _{0}y_{n}\\&{}\quad =h\left[\beta _{k}f(t_{n+k},y_{n+k})+\beta _{k-1}f(t_{n+k-1},y_{n+k-1})+\cdots +\beta _{0}f(t_{n},y_{n})\right].\end{aligned}}}$$

Інша можливість полягає у використанні більшої кількості точок в інтервалі ${\displaystyle [t_{n},t_{n+1}]}$ Це призводить до появи сімейства методів Рунге-Кутти, названих на честь Карла Рунге та Мартіна Кутти. Один з їхніх методів четвертого порядку є особливо популярним.

Додаткові можливості

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

Часто неефективно постійно використовувати один і той самий розмір кроку, тому були розроблені методи зі змінним розміром кроку. Зазвичай розмір кроку вибирають таким чином, щоб (локальна) похибка на крок була нижчою за певний рівень допуску. Це означає, що методи також повинні робити оцінку локальної похибки.

Розширенням цієї ідеї є динамічний вибір між різними методами різного порядку (це називається методом змінного порядку). Методи, засновані на екстраполяції Річардсона^[14], такі як алгоритм Булірша-Штера^[15][16] часто використовують для побудови різних методів різного порядку.

Інші бажані характеристики включають:

• щільний вихід: прості для обчислень чисельні наближення для всього інтервалу інтегрування, а не лише в точках t₀, t₁, t₂, …
• локалізація подій: знаходження моментів часу, коли, наприклад, певна функція обертається на нуль. Зазвичай це вимагає чисельного розв'язання нелінійних рівняньдля знаходження коренів.
• підтримка паралельних обчислень.
• при використанні для інтегрування за часом, оборотність за часом^[en].

Альтернативні методи

Багато методів не підпадають під обговорену тут класифікацію. Деякі класи альтернативних методів:

• багатопохідні методи (англ. multiderivative methods), які використовують не лише функцію f, але й її похідні. Цей клас включає методи Ерміта — Обрешкова та Фельберга, а також такі методи, як метод Паркера — Сохацького^[en][17] або метод Бичкова — Щербакова, які рекурсивно обчислюють коефіцієнти ряду Тейлора розв'язку y.
• методи для звичайних диференціальних рівнянь другого порядку. Ми казали, що всі звичайні диференціальні рівняння вищого порядку можна перетворити на звичайні диференціальні рівняння першого порядку виду (1). Хоча це, безумовно, правда, це може бути не найкращим способом дій. Зокрема, методи Ністрема працюють безпосередньо з рівняннями другого порядку.
• Методи геометричного інтегрування^[en][18][19] спеціально розроблені для спеціальних класів звичайних диференціальних уравнень (наприклад, симплектичних інтеграторів^[en] для розв'язання рівнянь Гамільтона). Вони забезпечують дотримання чисельним розв'язком правильної геометрії цих рівнянь.
• Методи систем квантованих станів^[en] — це сімейство методів інтегрування звичайних диференціальних рівнянь, заснованих на ідеї квантування станів. Вони ефективні при моделюванні розріджених систем з частими розривами.

Паралельні в часі методи

Деякі задачі з початковими умовами вимагають інтегрування з такою високою часовою роздільною здатністю або протягом таких тривалих часових інтервалів, що застосування класичних методів з послідовними часовими кроками стає неможливим (наприклад, у задачах чисельного прогнозування погоди, моделювання плазми та молекулярної динаміки). У відповідь на ці виклики були розроблені методи паралельні у часі (англ. parallel-in-time methods, PinT), щоб скоротити час виконання моделювання за рахунок використання паралельних обчислень.

Перші паралельні у часі методи були запропоновані в 1960-х роках^[20], але спочатку вони не привертали уваги дослідників, бо необхідні для них архітектури паралельних обчислень ще не були розповсюджені. Зі збільшенням обчислювальної потужності інтерес до них відродився на початку 2000-х років з розробкою Parareal^[en], гнучкого, простого у використанні алгоритму PinT, який підходить для вирішення широкого спектру задач з початковими умовами. Поява екзафлопових обчислень призвела до того, що алгоритми PinT стали привертати дедалі більшу увагу дослідників і розроблялися таким чином, щоб вони могли використовувати найпотужніші суперкомп'ютери світу. Популярності набули такі методи, як Parareal, PFASST, ParaDiag та MGRIT^[21].

Аналіз

Чисельний аналіз — це не лише розробка чисельних методів, а й їх аналіз. Три центральні концепції в цьому аналізі:

• збіжність: чи метод наближає розв'язок,
• порядок: наскільки добре він апроксимує розв'язок, та
• стійкість: чи згладжуються похибки^[22].

Збіжність

Чисельний метод називається збіжним, якщо чисельний розв'язок наближається до точного розв'язку, коли крок h прямує до 0. Точніше, ми вимагаємо, щоб для кожного звичайного диференціального рівняння (1) з функцією f, що задовільняє критерію Ліпшица, та кожного t^* > 0,

${\displaystyle \lim _{h\to 0^{+}}\max _{n=0,1,\dots ,\lfloor t^{*}/h\rfloor }\left\|y_{n,h}-y(t_{n})\right\|=0.}$

Усі вищезгадані методи є збіжними.

Узгодженість і порядок

Припустимо, що чисельний метод має вигляд

${\displaystyle y_{n+k}=\Psi (t_{n+k};y_{n},y_{n+1},\dots ,y_{n+k-1};h).\,}$

Локальна (усічена) похибка методу — це похибка, допущена на одному кроці методу. Тобто, це різниця між результатом, отриманим методом, за припущення, що на попередніх кроках взагалі не було похибок, та точним розв'язком:

${\displaystyle \delta _{n+k}^{h}=\Psi \left(t_{n+k};y(t_{n}),y(t_{n+1}),\dots ,y(t_{n+k-1});h\right)-y(t_{n+k}).}$

Метод називається узгодженим, якщо

${\displaystyle \lim _{h\to 0}{\frac {\delta _{n+k}^{h}}{h}}=0.}$

Метод має порядок ${\displaystyle p}$, якщо

${\displaystyle \delta _{n+k}^{h}=O(h^{p+1})\quad {\mbox{as }}h\to 0.}$

Отже, метод є узгодженим, якщо його порядок перевищує 0. Прямий метод Ейлера (4) та зворотний метод Ейлера (6), представлені вище, мають порядок 1, тому вони узгоджені. Більшість методів, використовуваних на практиці, мають вищий порядок. Узгодженість є необхідною умовою збіжності, але недостатньо. Щоб метод був збіжним, він повинен бути одночасно узгодженим і стійким відносно нуля.

Пов'язаною концепцією є глобальна похибка, — сумарна похибка, що виникає на всіх кроках, необхідних для досягнення фіксованого часу ${\displaystyle t}$. Глобальна похибка в момент часу ${\displaystyle t}$ дорівнює ${\displaystyle y_{N}-y(t)}$, де ${\displaystyle N=(t-t_{0})/h}$. Глобальна похибка одногрокового методу має порядок ${\displaystyle p}$, коли вона дорівнює ${\displaystyle O(h^{p})}$; зокрема, такий метод є збіжним. Це твердження не обов'язково вірне для багатокрокових методів.

Стійкість та жорсткість

Для деяких диференціальних рівнянь застосування стандартних методів, таких як метод Ейлера, явні методи Рунге — Кутти або багатокрокові методи (наприклад, методи Адамса — Башфорта), призводить до нестійкості у розв'язках, хоча інші методи можуть давати стійкі розв'язки. Ця «складна поведінка» рівнянь (яка самі по собі не обов'язково є складними) називають жорсткістю. Часто вона спричинена наявністю різних часових масштабів в основній задачі^[23]. Наприклад, зіткнення в механічній системі, такій як ударний осцилятор, зазвичай відбувається в набагато меншому часовому масштабі, ніж час руху об'єктів; ця невідповідність призводить до дуже «різких поворотів» розв'язків рівнянь.

Проблеми жорсткості є повсюдними в хімічній кінетиці, теорії керування, механіці твердого тіла, прогнозуванні погоди, біології, фізиці плазми та електроніці. Один зі способів подолання жорсткості полягає в розширенні поняття диференціального рівняння до поняття диференціального включення^[en], яке враховує та моделює негладкість^[24][25].

Історія

Нижче наведено хронологію деяких важливих подій у цій галузі^[26][27].

• 1768 — Леонард Ейлер публікує свій метод.
• 1824 — Огюстен-Луї Коші довів збіжність методу Ейлера. У цьому доведенні Коші використовує неявний метод Ейлера.
• 1855 — Перша згадка про багатокрокові методи Джона Кауча Адамса в листі, написаному Френсісом Башфортом^[en].
• 1895 — Карл Рунге публікує перший метод Рунге — Кутти.
• 1901 — Мартін Кутта описує популярний метод Рунге — Кутти четвертого порядку.
• 1910 — Льюїс Фрай Річардсон оголошує про свій метод екстраполяції, екстраполяцію Річардсона.
• 1952 — Чарльз Кертісс та Джозеф Окленд Гіршфельдер^[en] ввели термін «жорсткі рівняння».
• 1963 — Гермунд Дальквіст^[en] вводить поняття A-стійкості методів інтегрування.

Чисельні розв'язки одновимірних крайових задач другого порядку

Крайові задачі зазвичай розв'язуються чисельно шляхом розв'язання приблизно еквівалентної матричної задачі, отриманої шляхом дискретизації вихідної крайової задачі^[28]. Найпоширеніший метод чисельного розв'язання крайової задачі в одному вимірі називається методом скінченних різниць. Цей метод використовує лінійні комбінації точкових значень для побудови коефіцієнтів скінченних різниць^[3], які описують похідні функції. Наприклад, наближення центральної різниці другого порядку до першої похідної задається так:

${\displaystyle {\frac {u_{i+1}-u_{i-1}}{2h}}=u'(x_{i})+{\mathcal {O}}(h^{2}),}$

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

${\displaystyle {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}=u''(x_{i})+{\mathcal {O}}(h^{2}).}$

В обох цих формулах, ${\displaystyle h=x_{i}-x_{i-1}}$ — відстань між сусідніми значеннями x на дискретизованій області. Потім будується лінійна система, яку можна розв'язати стандартними матричними методами. Наприклад, припустимо, що рівняння, яке потрібно розв'язати, має вигляд:

${\displaystyle {\begin{aligned}&{}{\frac {d^{2}u}{dx^{2}}}-u=0,\\&{}u(0)=0,\\&{}u(1)=1.\end{aligned}}}$

Наступним кроком буде дискретизувати задачу та використати лінійні наближення для похідних, такі як

${\displaystyle u''_{i}={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}},}$

та розв'язати отриману систему лінійних рівнянь. Це призведе до рівнянь типу

${\displaystyle {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}-u_{i}=0,\quad \forall i={1,2,3,...,n-1}.}$

На перший погляд, ця система рівнянь не містить вільних членів, і тому їй задовільняє нульовий розв'язок, але насправді це не так. При i = 1 та n − 1 існують члени, що включають граничні значення ${\displaystyle u(0)=u_{0}}$ і ${\displaystyle u(1)=u_{n}}$, і, оскільки ці два значення відомі, їх можна просто підставити в це рівняння і в результаті отримати неоднорідну систему лінійних рівнянь, яка має нетривіальні розв'язки.