Одношаговые методы Рунге-Кутта
Явные методы Рунге-Кутта
Рассмотрим хитрую запись с внутренними вычислениями. Сначала вычисляем \(k_i\) по порядку, потом обновляем \(y_{n+1}\).
\[\begin{split}
\begin{aligned}
&\left\{\begin{array}{l}
k_{1}=f\left(t_{n}, \, y_{n}\right) \\
k_{2}=f\left(t_{n}+\alpha_{2} \cdot \tau, \;\; y_{n}+\tau \cdot \beta_{21} \cdot k_{1}\right) \\
\vdots \\
k_{r}=f\left(t_{n}+\alpha_{r}\cdot \tau, \;\; y_{n}+\tau \cdot\left(\beta_{r 1} \cdot k_{1}+\ldots+\beta_{r, \, r-1} k_{r}\right)\right)
\end{array}\right.
\end{aligned}
\end{split}\]
\[y_{n+1}=y_{n}+\tau \cdot\left(\gamma_{1}\cdot k_{1}+\ldots+\gamma_{r}\cdot k_{n}\right)\]
Для удобной записи коэффициентов используют таблицы Бутчера:
\[\begin{split}\\ \left[\begin{array}{l|llll}
0 & & & & \\
\alpha_{2} & \beta_{21} & & & \\
\alpha_{3} & \beta_{31} & \beta_{32} & & \\
\cdots &\cdots &\cdots &\cdots &\\
\alpha_{r} & \beta_{r 1} & \beta_{r, 2} & \cdots & \beta_{r, r-1} & \\
\hline& \gamma_{1} & \gamma_{2} & \cdots & \gamma_{r-1} & \gamma_{r}
\end{array}\right]\end{split}\]
Метод, задаваемый такой таблицей, называется r-стадийным. Все методы Рунге-Кутта являются разновидностью одношаговых методов, перепеписанными в удобном виде. Общие свойства этой таблицы Бутчера опишем ниже, когда введём неявные методы.
Явными их делаем аналитическая возможность решить написанную выше систему - по сути мы по порядку вычисляем все поправки.
Приведём некоторые примеры (за точные коэффициенты здесь не ручаюсь):
Модифицированный метод Эйлера:
\[\begin{split}
\begin{aligned}
&y_{n+1}=y_{n}+h f_{2}, \\
&f_{1}=f\left(x_{n}, y_{n}\right) \\
&f_{2}=f\left(x_{n}+\frac{h}{2} ; y_{n}+\frac{h f_{n}}{2}\right)
\end{aligned}
\end{split}\]
\[\begin{split}
\begin{array}{l|ll}
0 & 0 & 0 \\
1 / 2 & 1 / 2 & 0 \\
\hline & 0 & 1
\end{array}
\end{split}\]
Метод Эйлера с пересчётом:
\[\begin{split}
\begin{aligned}
&y_{n+1}=y_{n}+h \cdot \frac{f_{1}+f_{2}}{2} \\
&f_{1}=f\left(x_{1}, y_{n}\right) \\
&f_{2}=f\left(x_{n}+h, y_{n}+h f_{1}\right)
\end{aligned}
\end{split}\]
\[\begin{split}
\begin{array}{l|ll}
0 & 0 & 0 \\
1 & 1 & 0 \\
\hline & 1 / 2 & 1 / 2
\end{array}
\end{split}\]
Метод Хойна:
\[\begin{split}
\begin{aligned}
&y_{n+1}=y_{n}+h\left(f_{1}+3 f_{3}\right) \cdot \frac{1}{4} \\
&f_{1}=f\left(x_{n}, y_{n}\right) ; \\
&f_{2}=f\left(x_{n}+\frac{h}{3} ; y_{n}+\frac{h f_{1}}{3}\right) \\
&f_{3}=f\left(x_{n}+\frac{2 h}{3} ; y_{n}+2 h \frac{f_{2}}{3}\right)
\end{aligned}
\end{split}\]
\[\begin{split}
\begin{array}{l|ccc}
0 & 0 & 0 & 0 \\
1 / 3 & 1 / 3 & 0 & 0 \\
2 / 3 & 0 & 2 / 3 & 0 \\
\hline & 1 / 4 & 0 & 1 / 4
\end{array}
\end{split}\]
Классический метод Рунге-Кутта 4-го порядка:
\[\begin{split}
\begin{array}{c|cccc}
0 & 0 & 0 & 0 & 0 \\
1 / 2 & 1 / 2 & 0 & 0 & 0 \\
1 / 2 & 0 & 1 / 2 & 0 & 0 \\
1 & 0 & 0 & 1 & 0 \\
\hline & 1 / 6 & 2 / 6 & 2 / 6 & 1 / 6
\end{array}
\end{split}\]
Метод Дормана-Принса 5(4)
При компьютерных вычисления чаще всего используют метод Дормана-Принса. В 1980-м году Дорман и Принс построили вложенный метод 5(4)-го порядка, в котором решение 5-го порядка используется в качестве начального значения для следующего шага, а решение 4-го порядка – для определения локальной погрешности выполненного шага интегрирования (с целью последующего использования в механизме управления длинной шага).
\[\begin{split}
\begin{array}{c|cccccccc}
0 & & & & & & \\
\frac{1}{5} & \frac{1}{5} & & & & & \\
\frac{3}{10} & \frac{3}{40} & \frac{9}{40} & & & & & \\
\frac{4}{5} & \frac{44}{45} & -\frac{56}{15} & \frac{32}{9} & & & & \\
\frac{8}{9} & \frac{19372}{6561} & -\frac{25360}{2187} & \frac{64448}{6561} & -\frac{212}{729} & & & \\
1 & \frac{9017}{3168} & -\frac{355}{33} & \frac{46732}{5247} & \frac{49}{176} & -\frac{5103}{18656} & & \\
1 & \frac{35}{384} & 0 & \frac{500}{1113} & \frac{125}{192} & -\frac{2187}{6784} & \frac{11}{84} & \\
\hline x_1 & \frac{35}{384} & 0 & \frac{500}{1113} & \frac{125}{192} & -\frac{2187}{6784} & \frac{11}{84} & 0 \\
\hline \hat{x}_1 & \frac{5179}{57600} & 0 & \frac{7571}{16695} & \frac{393}{640} & -\frac{92097}{339200} & \frac{187}{2100} & \frac{1}{40}
\end{array}
\end{split}\]
Первая строка коэффициентов b дает решение пятого порядка точности, а вторая строка дает альтернативное решение, которое при вычитании из первого решения дает оценку ошибки.
Поэтому его часто называют методом «4-5-го порядка».
Примечание. Несложно обобщить данный класс методов на неявный случай. Если таблица Бутчера на нижнетреугольная, то в общем случае на каждом шаге придётся решать нелинейную систему, чтобы найти все \(k_i\).
Неявные методы Рунге-Кутта
Гораздо лучшую устойчивость имеют неявные методы Рунге-Кутта. Они вполне себе пригодны для решения жестких задач.
Пусть таблица для r-стадийного неявного метода Рунге-Кутта имеет следующий __не__нижнетреугольный вид:
\[\begin{split}\\ \left[\begin{array}{l|llll}
\alpha_{1} & \beta_{11} & \beta_{12} & \cdots & \beta_{1, r-1} & \beta_{1r} \\
\alpha_{2} & \beta_{21} & \beta_{22} & \cdots & \beta_{2, r-1} & \beta_{2r} \\
\alpha_{3} & \beta_{31} & \beta_{32} & \cdots & \beta_{3, r-1} & \beta_{3r} \\
\cdots &\cdots &\cdots &\cdots & \cdots & \cdots \\
\alpha_{r} & \beta_{r, 1} & \beta_{r, 2} & \cdots & \beta_{r,\; r-1} & \beta_{r,\; r} \\
\hline & \gamma_{1} & \gamma_{2} & \cdots & \gamma_{r-1} & \gamma_{r}
\end{array}\right]\end{split}\]
Тогда соответствующая ей система (вообще говоря нелинейная), которую необходимо решать на каждом шаге:
\[\begin{split}
\left\{\begin{array}{l}
k_1=f\left(t_n+\alpha_1 \cdot \tau, \quad y_n+\tau \cdot\left(\beta_{11} \cdot k_1+\beta_{12} \cdot k_2+\cdots+\beta_{1 r} \cdot k_r\right)\right) \\
k_2=f\left(t_n+\alpha_2 \cdot \tau, \quad y_n+\tau \cdot\left(\beta_{21} \cdot k_1+\beta_{22} \cdot k_2+\cdots+\beta_{2 r} \cdot k_r\right)\right) \\
k_3=f\left(t_n+\alpha_3 \cdot \tau, \quad y_n+\tau \cdot\left(\beta_{31} \cdot k_1+\beta_{32} \cdot k_2+\cdots+\beta_{3 r} \cdot k_r\right)\right) \\
\vdots \\
k_r=f\left(t_n+\alpha_r \cdot \tau, \quad y_n+\tau \cdot\left(\beta_{r 1} \cdot k_1+\beta_{r 2} \cdot k_2+\cdots+\beta_{r, r-1} \cdot k_{r-1}+\beta_{r, r} \cdot k_r\right)\right) \\
\end{array}\right.
\end{split}\]
\[
y_{n+1}=y_n+\tau \cdot\left(\gamma_1 \cdot k_1+\gamma_2 \cdot k_2+\cdots+\gamma_r \cdot k_r\right)
\]
Как правило, эту систему решают методом простой итерации (см. прошлый семестр).
Общее свойство всех таблиц Бутчера - свойство Кутты:
\[
\sum_{j=1}^{r} \beta_{i j}=\alpha_i
\]
Для обеспечения нужного порядка аппроксимации \(p\), также необходимы определенные условия на коэффициента таблицы Бутчера. Их несложно вывести, раскладывая \([y(x)]_{n+1}\) в ряд Тейлора и зануляя коэффициенты перед ненужными степенями \(h\) методом неопределенных коэффициентов:
\[\begin{split}
\begin{aligned}
& p=1 \hookrightarrow \sum_{i=1}^s \gamma_i=1, \\
& p=2 \hookrightarrow \sum_{i=1}^s \gamma_i \alpha_i=1 / 2, \\
& p=3 \hookrightarrow \sum_{i=1}^s \gamma_i \alpha_i^2=1 / 3, \quad \sum_{i, j=1}^s b_i \beta_{i j} \alpha_j=1 / 6, \\
& p=4 \hookrightarrow \sum_{i=1}^s \gamma_i \alpha_i^3=1 / 4, \quad \sum_{i, j=1}^s b_i \alpha_i \beta_{i j} \alpha_j=1 / 8, \quad \sum_{i, j=1}^s \gamma_i \beta_{i j} \alpha_j^2=1 / 12, \quad \sum_{i, j, k=1}^s \gamma_i \beta_{i j} \beta_{j k} \alpha_k=1 / 24 .
\end{aligned}
\end{split}\]
Например, для порядка аппроксимации \(p=4\) - нужно выполнения всех строчек выше. А для \(p=3\), только первых трёх.
Методы Гаусса соответственно 2-го, 4-го, и 6-го порядков аппроксимации. Второй метод называется методом (4-го порядка) называется методом Хаммера - Холлинсворта:
\[\begin{split}\left[
\begin{array}{l|l}
1 / 2 & 1 / 2 \\
\hline & 1
\end{array}\right]
\end{split}\]
\[\begin{split}\left[
\begin{array}{c|cc}
\frac{1}{2}-\frac{\sqrt{3}}{6} & \frac{1}{4} & \quad \frac{1}{4}-\frac{\sqrt{3}}{6} \\
\frac{1}{2}+\frac{\sqrt{3}}{6} & \frac{1}{4}+\frac{\sqrt{3}}{6} & \frac{1}{4} \\
\hline &\frac{1}{2} & \frac{1}{2} \\
\end{array}\right]
\end{split}\]
\[\begin{split}\left[
\begin{array}{c|cccc}
\frac{1}{2}-\frac{\sqrt{15}}{10} & \frac{5}{36} & \frac{2}{9}-\frac{\sqrt{15}}{15} &\frac{5}{36}-\frac{\sqrt{15}}{30} \\
\frac{1}{2} & \frac{5}{36}+\frac{\sqrt{15}}{24} &\frac{2}{9} & \frac{5}{36}-\frac{\sqrt{15}}{24} \\
\frac{1}{2}+\frac{\sqrt{15}}{10} & \frac{5}{36}+\frac{\sqrt{15}}{30} &\frac{2}{9}+\frac{\sqrt{15}}{15} & \frac{5}{36} \\
\hline &\frac{5}{18} & \frac{4}{9} & \frac{5}{18}\\
\end{array}\right]
\end{split}\]
Методы Лобатто IIIА 2-го, 4-го и 6-го порядков точности:
\[\begin{split}\left[
\begin{array}{l|ll}
0 & 0 & 0 \\
1 & 1 / 2 & 1 / 2 \\
\hline & 1 / 2 & 1 / 2 \\
\end{array}\right]
\end{split}\]
\[\begin{split}\left[
\begin{array}{l|lll}
0 & 0 & 0 & 0 \\
1 / 2 & 5 / 24 & 1 / 3 & -1 / 24 \\
1 & 1 / 6 & 2 / 3 & 1 / 6 \\
\hline & 1 / 6 & 2 / 3 & 1 / 6 \\
\end{array}\right]
\end{split}\]
\[\begin{split}\left[
\begin{array}{l|llll}
0 & 0 & 0 & 0 & 0 \\
\frac{5-\sqrt{5}}{10} & \frac{11+\sqrt{5}}{120} & \frac{25-\sqrt{5}}{120} & \frac{25-13 \sqrt{5}}{120} & \frac{-1+\sqrt{5}}{120} \\
\frac{5+\sqrt{5}}{10} & \frac{11-\sqrt{5}}{120} & \frac{25+13 \sqrt{5}}{120} & \frac{25+\sqrt{5}}{120} & \frac{-1-\sqrt{5}}{120} \\
1 & 1 / 12 & 5 / 12 & 5 / 12 & 1 / 12 \\
\hline & 1 / 12 & 5 / 12 & 5 / 12 & 1 / 12 \\
\end{array}\right]
\end{split}\]
Устойчивость методов Рунге-Кутта
Теорема (о функции устойчивости методов Рунге-Кутта). Пусть имеем r-стадийый метод Рунге-Кутта, заданный таблицей Бутчера \(\begin{array}{c|c}
\vec{\alpha} & \mathbf{B} \\
\hline & \vec{\gamma}^\top
\end{array}\).
Тогда его функция устойчивости имеет вид
\[
R(z)=1+z \cdot \vec{\gamma}^\top (\mathbf{E}-z \cdot \mathbf{B})^{-1}\vec{1}
\]
или, что эквивалентно:
\[
R(z)=\cfrac{\operatorname{det}\left(\mathbf{E}-z \mathbf{B}+z \cdot \mathbf{\vec{1} \vec{\gamma}}^T\right)}{ \operatorname{det}\left(\mathbf{E}-z \mathbf{B}\right)}
\]
В случае, если система имеет порядок аппроксимации \(p\), и этот самый порядок равняется числу стадий \(p=r\), то эту формулу можно явно вычислить:
\[R_{r=p}(z)=1+z+\ldots+ \frac{z^r}{r !}\]
Доказательство.
To be written
Пример исследования на устойчивость неявного метода Рунге-Кутта
Покажите, что неявный метод Рунге-Кутты (метод Хаммера-Холлингсуорта порядка 4), заданный таблицей Бутчера:
\[\begin{split}
\begin{array}{l|ll}
1 / 2-\sqrt{3} & 1 / 4 & 1 / 4-\sqrt{3} / 6 \\
1 / 2+\sqrt{3} & 1 / 4+\sqrt{3} / 6 & 1 / 4 \\
\hline & 1 / 2 & 1 / 2
\end{array}
\end{split}\]
имеет четвертый порядок аппроксимации. Найдите функцию устойчивости \(R(z)\) и исследуйте метод на \(A\)-устойчивость и \(L\)-устойчивость.
Решение.
Вспомним условия обеспечения порядков аппроксимации методов Рунге-Кутты:
\[\begin{split}
\begin{aligned}
& p=1 \hookrightarrow \sum_{i=1}^s \gamma_i=1, \\
& p=2 \hookrightarrow \sum_{i=1}^s \gamma_i \alpha_i=1 / 2, \\
& p=3 \hookrightarrow \sum_{i=1}^s \gamma_i \alpha_i^2=1 / 3, \quad \sum_{i, j=1}^s b_i \beta_{i j} \alpha_j=1 / 6, \\
& p=4 \hookrightarrow \sum_{i=1}^s \gamma_i \alpha_i^3=1 / 4, \quad \sum_{i, j=1}^s b_i \alpha_i \beta_{i j} \alpha_j=1 / 8, \quad \sum_{i, j=1}^s \gamma_i \beta_{i j} \alpha_j^2=1 / 12, \quad \sum_{i, j, k=1}^s \gamma_i \beta_{i j} \beta_{j k} \alpha_k=1 / 24 .
\end{aligned}
\end{split}\]
Для данного метода все эти соотношения проверяются непосредственной подстановкой (соотношения пятого порядка не выполняются). Т.е. порядок аппроксимации действительно 4.
Для нахождения функции устойчивости воспользуемся общей формулой
\[\begin{split}
\begin{aligned}
& R(z)=1+z \cdot \vec{\gamma}^\top (\mathbf{E}-z \cdot \mathbf{B})^{-1}\vec{1} =1+\frac{z}{2}(1,1)\left(\begin{array}{cc}
1-\frac{z}{4} & -\frac{z}{4}+\frac{\sqrt{3} z}{6} \\
-\frac{z}{4}-\frac{\sqrt{3} z}{6} & 1-\frac{z}{4}
\end{array}\right)^{-1}\left(\begin{array}{l}
1 \\
1
\end{array}\right)= \\
& =1+\frac{z / 2}{\left(1-\frac{z}{4}\right)^2-\left(\frac{z^2}{16}-\frac{z^2}{12}\right)}(1,1)\left(\begin{array}{cc}
1-\frac{z}{4} & \frac{z}{4}-\frac{\sqrt{3} z}{6} \\
\frac{z}{4}+\frac{\sqrt{3} z}{6} & 1-\frac{z}{4}
\end{array}\right)\left(\begin{array}{l}
1 \\
1
\end{array}\right)=1+\frac{z}{\left(1-\frac{z}{2}+\frac{z^2}{12}\right)}= \\
& =\frac{\mathrm{z}^2+6 z+12}{\mathrm{z}^2-6 z+12}
\end{aligned}
\end{split}\]
Эта функция является регулярной в левой полуплоскости \(\mathbb{C}^{-}\)и по принципу максимума модуля ее модуль принимает максимальное значение на границе (вертикальная ось \(z=i y\) ).
Напомним определение регулярной функции и принцип максимума модуля.
Функция \(f(z)\) комплексного переменного \(z\), однозначно определенная в области \(G\) и имеющая в каждой её точке конечную производную, называется регулярной в этой области.
Tеорема о принципе максимума. Модуль функиии, регулярной в некоторой области, не может достигать наибольшего значения внутри этой области, если только функиия отлична от тождественной постоянной.
На вертикальной оси \(|R(z)|=\left|\frac{12-y^2+6 i y}{12-y^2-6 i y}\right| = 1 \quad \forall y\).
Таким образом, метод \(A\)-устойчив, но не \(L\)-устойчив, поскольку \(\lim _{z \rightarrow-\infty} R(z) \neq 0\).
Построим также область абсолютной устойчивости:
\[\displaystyle \left|{\frac{z^{2} + 6 z + 12}{z^{2} - 6 z + 12}}\right| < 1\]
\[\displaystyle \frac{\sqrt{x^{4} + 12 x^{3} + 2 x^{2} y^{2} + 60 x^{2} + 12 x y^{2} + 144 x + y^{4} + 12 y^{2} + 144}}{\sqrt{x^{4} - 12 x^{3} + 2 x^{2} y^{2} + 60 x^{2} - 12 x y^{2} - 144 x + y^{4} + 12 y^{2} + 144}} < 1\]
<sympy.plotting.plot.Plot at 0x7feb48d62940>
