4: Диференціальні рівняння

Зараз ми обговоримо чисельне рішення звичайних диференціальних рівнянь. Ми включимо початкову задачу і крайову задачу.

Задача про початкове значення

Почнемо з простого методу Ейлера, потім обговоримо більш складні методи RungEkutta, і завершуємо методом Рунге-Кутта-Фейлберга, реалізованим в оді функції MATLAB\(45 . \mathrm{m}\). Наші диференціальні рівняння призначені для\(x=x(t)\), де час\(t\) є незалежною змінною.

метод Ейлера

Метод Ейлера є найбільш простим методом інтеграції диференціального рівняння. Розглянемо диференціальне рівняння першого порядку

\[\dot{x}=f(t, x), \nonumber \]

з початковим умовою\(x(0)=x_{0}\). Визначте\(t_{n}=n \Delta t\) і\(x_{n}=x\left(t_{n}\right)\). Розширення\(x_{n+1}\) результатів серії Тейлора в

\[\begin{aligned} x_{n+1} &=x\left(t_{n}+\Delta t\right) \\ &=x\left(t_{n}\right)+\Delta t \dot{x}\left(t_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) \\ &=x\left(t_{n}\right)+\Delta t f\left(t_{n}, x_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) . \end{aligned} \nonumber \]

Таким чином, метод Ейлера пишеться як

\[x_{n+1}=x\left(t_{n}\right)+\Delta t f\left(t_{n}, x_{n}\right) . \nonumber \]

Ми говоримо, що метод Ейлера крокує вперед у часі, використовуючи тайм-крок\(\Delta t\), починаючи від початкового значення\(x_{0}=x(0)\). Локальна помилка методу Ейлера є\(\mathrm{O}\left(\Delta t^{2}\right)\). Однак глобальна помилка, що виникає при інтеграції в час\(T\), є фактором\(1 / \Delta t\) більшого і дається\(\mathrm{O}(\Delta t)\). Тому прийнято називати метод Ейлера методом першого порядку.

Модифікований метод Ейлера

Цей метод є так званим предикторно-коректорним методом. Це також перше з того, що ми побачимо, - це методи Рунге-Кутта. Як і раніше, ми хочемо вирішити (4.1). Ідея полягає в тому, щоб\(\dot{x}\) усереднити значення на початку і кінці часового кроку. Тобто ми хотіли б модифікувати метод Ейлера і написати

\[x_{n+1}=x_{n}+\frac{1}{2} \Delta t\left(f\left(t_{n}, x_{n}\right)+f\left(t_{n}+\Delta t, x_{n+1}\right)\right) . \nonumber \]

Очевидна проблема цієї формули полягає в тому, що з правого боку\(x_{n+1}\) з'являється невідоме значення. Однак ми можемо оцінити це значення в тому, що називається кроком предиктора. Для кроку предиктора ми використовуємо метод Ейлера, щоб знайти

\[x_{n+1}^{p}=x_{n}+\Delta t f\left(t_{n}, x_{n}\right) \text {. } \nonumber \]

Крок коректора тоді стає

\[x_{n+1}=x_{n}+\frac{1}{2} \Delta t\left(f\left(t_{n}, x_{n}\right)+f\left(t_{n}+\Delta t, x_{n+1}^{p}\right)\right) . \nonumber \]

Модифікований метод Ейлера може бути переписаний у такій формі, яку ми згодом ідентифікуємо як метод Рунге-Кутта

\[\begin{align}\nonumber &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{1}\right), \\ &x_{n+1}=x_{n}+\frac{1}{2}\left(k_{1}+k_{2}\right) . \nonumber \end{align} \nonumber \]

Методи Рунге-Кутта другого порядку

Тепер ми виведемо повне сімейство методів Рунге-Кутта другого порядку. Методи вищого порядку можуть бути аналогічно виведені, але вимагають значно більшої алгебри.

Ми знову розглянемо диференціальне рівняння, задане (4.1). Загальний метод Рунге-Кутта другого порядку може бути записаний у формі

\[\begin{align} \nonumber &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta k_{1}\right), \\ &x_{n+1}=x_{n}+a k_{1}+b k_{2}, \nonumber \end{align} \nonumber \]

з\(\alpha, \beta, a\) і\(b\) константами, що визначають конкретний метод Рунге-Кутта другого порядку. Ці константи повинні бути обмежені встановленням локальної помилки методу Рунге-Кутта другого порядку be\(\mathrm{O}\left(\Delta t^{3}\right)\). Інтуїтивно ми можемо здогадатися, що два обмеження будуть\(a+b=1\) і\(\alpha=\beta\).

Ми обчислюємо серії Тейлора\(x_{n+1}\) безпосередньо, і з методу Рунге-Кутта, і вимагаємо, щоб вони були однаковими для замовлення\(\Delta t^{2}\). По-перше, ми обчислюємо ряд Тейлора\(x_{n+1}\). У нас є

\[\begin{aligned} x_{n+1} &=x\left(t_{n}+\Delta t\right) \\ &=x\left(t_{n}\right)+\Delta t \dot{x}\left(t_{n}\right)+\frac{1}{2}(\Delta t)^{2} \ddot{x}\left(t_{n}\right)+\mathrm{O}\left(\Delta t^{3}\right) . \end{aligned} \nonumber \]

Тепер,

\[\dot{x}\left(t_{n}\right)=f\left(t_{n}, x_{n}\right) . \nonumber \]

Друга похідна більш складна і вимагає часткових похідних. У нас є

\[\begin{aligned} \ddot{x}\left(t_{n}\right) &\left.=\frac{d}{d t} f(t, x(t))\right]_{t=t_{n}} \\ &=f_{t}\left(t_{n}, x_{n}\right)+\dot{x}\left(t_{n}\right) f_{x}\left(t_{n}, x_{n}\right) \\ &=f_{t}\left(t_{n}, x_{n}\right)+f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right) . \end{aligned} \nonumber \]

Тому,

\[\begin{align} \nonumber x_{n+1}=x_{n}+\Delta t f\left(t_{n}, x_{n}\right) & \\ &+\frac{1}{2}(\Delta t)^{2}\left(f_{t}\left(t_{n}, x_{n}\right)+f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) . \end{align} \nonumber \]

По-друге, обчислюємо\(x_{n+1}\) з методу Рунге-Кутта, заданого (4.7). Об'єднавши (4.7) в єдиний вираз, ми маємо

\[x_{n+1}=x_{n}+a \Delta t f\left(t_{n}, x_{n}\right) \nonumber \]

\[+b \Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta \Delta t f\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) . \nonumber \]

Ми Taylor серії розширюємо за допомогою

\[\begin{aligned} f\left(t_{n}+\alpha \Delta t, x_{n}+\beta \Delta t f\left(t_{n}, x_{n}\right)\right) & \\ &=f\left(t_{n}, x_{n}\right)+\alpha \Delta t f_{t}\left(t_{n}, x_{n}\right)+\beta \Delta t f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) . \end{aligned} \nonumber \]

Тому формула Рунге-Кутта

\[\begin{aligned} x_{n+1}=x_{n}+(a+b) \Delta t f &\left(t_{n}, x_{n}\right) \\ &+(\Delta t)^{2}\left(\alpha b f_{t}\left(t_{n}, x_{n}\right)+\beta b f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) . \end{aligned} \nonumber \]

Порівнюючи (4.9) і (4.10), знаходимо

\[a+b=1, \quad \alpha b=1 / 2, \quad \beta b=1 / 2 . \nonumber \]

Оскільки існує всього три рівняння для чотирьох параметрів, існує сімейство методів Рунге-Кутта другого порядку.

Модифікований метод Ейлера, заданий (4.6), відповідає\(\alpha=\beta=1\) і\(a=\)\(b=1 / 2\). Інший метод Рунге-Кутта другого порядку, який називається методом середньої точки, відповідає\(\alpha=\beta=1 / 2, a=0\) і\(b=1\). Цей спосіб пишеться як

\[\begin{aligned} &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}\right), \\ &x_{n+1}=x_{n}+k_{2} . \end{aligned} \nonumber \]

Методи Рунге-Кутта вищого порядку

Загальний метод Рунге-Кутта другого порядку був заданий (4.7). Загальна форма методу третього порядку задається

\[\begin{aligned} &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta k_{1}\right), \\ &k_{3}=\Delta t f\left(t_{n}+\gamma \Delta t, x_{n}+\delta k_{1}+\epsilon k_{2}\right), \\ &x_{n+1}=x_{n}+a k_{1}+b k_{2}+c k_{3} . \end{aligned} \nonumber \]

Можна вгадати наступні обмеження на константи:\(\alpha=\beta, \gamma=\delta+\epsilon\), і\(a+b+c=1\). Решта обмеження потрібно вивести.

Метод четвертого порядку має\(k_{1}, k_{2}, k_{3}\) і\(k_{4}\). Метод п'ятого порядку вимагає, як мінімум, до\(k_{6}\). У таблиці нижче наведено мінімальне замовлення методу і кількість необхідних етапів.

замовляти	2	3	4	5	6	7	8
мінімум # етапів	2	3	4	6	7	9	11

Через стрибка кількості етапів, необхідних між методами четвертого і п'ятого порядку, метод Рунге-Кутта четвертого порядку має деяку привабливість. Загальний метод четвертого порядку з чотирма етапами має 13 констант і 11 обмежень. Особливо простий метод четвертого порядку, який широко використовувався в минулому фізиками, дається

\[\begin{aligned} &k_{1}=\Delta t f\left(t_{n}, x_{n}\right) \\ &k_{2}=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}\right) \\ &k_{3}=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}\right) \\ &k_{4}=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{3}\right) \\ &x_{n+1}=x_{n}+\frac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right) . \end{aligned} \nonumber \]

Адаптивні методи Рунге-Кутта

Як і в адаптивній інтеграції, корисно розробити інтегратор оди, який автоматично знаходить відповідний\(\Delta t\). Метод Дорманда-Принца, реалізований в оді MATLAB\(45 . \mathrm{m}\), знаходить відповідний розмір кроку шляхом порівняння результатів методу п'ятого та четвертого порядку. Він вимагає шести оцінок функцій на часовий крок і будує як метод п'ятого порядку, так і четвертого порядку з тих же оцінок функцій.

Припустимо, метод п'ятого порядку знаходить\(x_{n+1}\) з локальною помилкою\(\mathrm{O}\left(\Delta t^{6}\right)\), а метод четвертого порядку знаходить\(x_{n+1}^{\prime}\) з локальною помилкою\(\mathrm{O}\left(\Delta t^{5}\right)\). \(\varepsilon\)Дозволяти бажану похибку толерантності методу, і\(e\) нехай фактична помилка. Ми можемо оцінити\(e\) з різниці між методами п'ятого та четвертого порядку; тобто

\[e=\left|x_{n+1}-x_{n+1}^{\prime}\right| \text {. } \nonumber \]

Тепер\(e\) це з\(\mathrm{O}\left(\Delta t^{5}\right)\), де\(\Delta t\) взятий розмір кроку. \(\Delta \tau\)Дозволяти приблизний розмір кроку, необхідний для отримання бажаної помилки\(\varepsilon\). Тоді у нас є

\[e / \varepsilon=(\Delta t)^{5} /(\Delta \tau)^{5}, \nonumber \]

або рішення для\(\Delta \tau\),

\[\Delta \tau=\Delta t\left(\frac{\varepsilon}{e}\right)^{1 / 5} \nonumber \]

З одного боку, якщо фактична помилка менше бажаної помилки\(e<\varepsilon\), то приймаємо\(x_{n+1}\) і робимо наступний крок часу, використовуючи більше значення\(\Delta \tau\). З іншого боку, якщо фактична помилка більше бажаної помилки\(e>\varepsilon\), то ми відхиляємо крок інтеграції і повторюємо часовий крок, використовуючи меншу величину\(\Delta \tau\). На практиці зазвичай збільшується часовий крок трохи менше і зменшує часовий крок трохи більше, щоб запобігти марнотратству занадто багатьох невдалих кроків часу.

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

Наші числові методи можуть бути легко адаптовані для вирішення диференціальних рівнянь вищого порядку, або, еквівалентно, системи диференціальних рівнянь. Спочатку ми покажемо, як диференціальне рівняння другого порядку може бути зведено до двох рівнянь першого порядку. Розглянемо

\[\ddot{x}=f(t, x, \dot{x}) . \nonumber \]

Це рівняння другого порядку можна переписати як два рівняння першого порядку шляхом визначення\(u=\dot{x}\). Тоді у нас є система

\[\dot{x}=u, \quad \dot{u}=f(t, x, u) . \nonumber \]

Цей трюк також працює для рівнянь вищого порядку. Наприклад, рівняння третього порядку

\[\dddot x=f(t, x, \dot{x}, \ddot{x}), \nonumber \]

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

\[\dot{x}=u, \quad \dot{u}=v, \quad \dot{v}=f(t, x, u, v) . \nonumber \]

Можна узагальнити метод Рунге-Кутта для розв'язання системи диференціальних рівнянь. Як приклад розглянемо наступну систему двох од,

\[\dot{x}=f(t, x, y), \quad \dot{y}=g(t, x, y), \nonumber \]

з початковими умовами\(x(0)=x_{0}\) і\(y(0)=y_{0}\). Узагальнення часто використовуваного методу Рунге-Кутта четвертого порядку було б

\[\begin{aligned} k_{1} &=\Delta t f\left(t_{n}, x_{n}, y_{n}\right), \\ l_{1} &=\Delta t g\left(t_{n}, x_{n}, y_{n}\right), \\ k_{2} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}, y_{n}+\frac{1}{2} l_{1}\right), \\ l_{2} &=\Delta t g\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}, y_{n}+\frac{1}{2} l_{1}\right), \\ k_{3} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}, y_{n}+\frac{1}{2} l_{2}\right), \\ l_{3} &=\Delta t g\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}, y_{n}+\frac{1}{2} l_{2}\right), \\ k_{4} &=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{3}, y_{n}+l_{3}\right), \\ l_{4} &=\Delta t g\left(t_{n}+\Delta t, x_{n}+k_{3}, y_{n}+l_{3}\right), \\ x_{n+1} &=x_{n}+\frac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right), \\ y_{n+1} &=y_{n}+\frac{1}{6}\left(l_{1}+2 l_{2}+2 l_{3}+l_{4}\right) . \end{aligned} \nonumber \]

крайові задачі

спосіб зйомки

Розглядаємо загальну оду виду

\[\frac{d^{2} y}{d x^{2}}=f(x, y, d y / d x), \nonumber \]

з двоточковими граничними умовами\(y(0)=A\) і\(y(1)=B\). Спочатку ми сформулюємо оду як початкову ціннісну задачу. У нас є

\[\frac{d y}{d x}=z, \quad \frac{d z}{d x}=f(x, y, z) . \nonumber \]

Початкова умова\(y(0)=A\) відома, але друга початкова умова\(z(0)=b\) невідомо. Наша мета полягає в тому, щоб визначити\(b\) таке\(y(1)=B\).

По суті, це проблема пошуку коренів для відповідним чином визначеної функції. Визначимо функцію\(F=F(b)\) таку, що

\[F(b)=y(1)-B, \nonumber \]

де\(y(1)\) - числове значення, отримане в результаті інтеграції зв'язаних диференціальних рівнянь першого порядку з\(y(0)=A\) і\(z(0)=b\). Наша рутина пошуку коренів захоче вирішити\(F(b)=0\). (Метод називається зйомкою, оскільки нахил кривої розв'язку для\(y=y(x)\) at\(x=0\) задається\(b\), а рішення потрапляє\(y(1)\) на значення в\(x=1\). Це виглядає як наведення гармати і намагатися вистрілити в ціль, яка є\(B\).)

Щоб визначити значення,\(b\) що розв'язує\(F(b)=0\), ми виконуємо ітерацію за допомогою методу secant, заданого

\[b_{n+1}=b_{n}-F\left(b_{n}\right) \frac{b_{n}-b_{n-1}}{F\left(b_{n}\right)-F\left(b_{n-1}\right)} . \nonumber \]

Починати потрібно з двох початкових здогадок для\(b\), вирішуючи оду для двох відповідних значень\(y(1)\). Тоді метод secant дасть нам наступне значення,\(b\) щоб спробувати, і ми перебираємо до\(|y(1)-B|<\) tol, де tol є деяким зазначеним допуском до помилки.