Метод Рунге — Кутты

Ме́тоды Ру́нге — Ку́тты (в литературе встречается название методы Рунге — Кутта) — большой класс численных методов решения задачи Коши для обыкновенных дифференциальных уравнений и их систем. Первые методы данного класса были предложены около 1900 года немецкими математиками К. Рунге и М. В. Куттой.

К классу методов Рунге — Кутты относятся явный метод Эйлера и модифицированный метод Эйлера с пересчётом, которые представляют собой соответственно методы первого и второго порядка точности. Существуют стандартные явные методы третьего порядка точности, не получившие широкого распространения. Наиболее часто используется и реализован в различных математических пакетах (Maple, MathCAD, Maxima) классический метод Рунге — Кутты, имеющий четвёртый порядок точности. При выполнении расчётов с повышенной точностью всё чаще применяются методы пятого и шестого порядков точности. Построение схем более высокого порядка сопряжено с большими вычислительными трудностями.

Методы седьмого порядка должны иметь по меньшей мере девять стадий, а методы восьмого порядка — не менее 11 стадий. Для методов девятого и более высоких порядков (не имеющих, впрочем, большой практической значимости) неизвестно, сколько стадий необходимо для достижения соответствующего порядка точности.

Классический метод Рунге — Кутты четвёртого порядка

Метод Рунге — Кутты четвёртого порядка при вычислениях с постоянным шагом интегрирования столь широко распространён, что его часто называют просто методом Рунге — Кутты.

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений первого порядка (далее $𝐲,𝐟,{𝐤}_i\in {ℝ}^n$, а $x,h\in {ℝ}^1$).

${\mathbf{\text{y}}}^{\prime }=\mathbf{\text{f}}(x,\mathbf{\text{y}}),\mathbf{\text{y}}(x_0)={\mathbf{\text{y}}}_0.$

Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:

${\mathbf{\text{y}}}_{n+1}={\mathbf{\text{y}}}_n+\frac{h}{6}({\mathbf{\text{k}}}_1+2{\mathbf{\text{k}}}_2+2{\mathbf{\text{k}}}_3+{\mathbf{\text{k}}}_4)$

Вычисление нового значения проходит в четыре стадии:

${\mathbf{\text{k}}}_1=\mathbf{\text{f}}(x_n,{\mathbf{\text{y}}}_n),$

${\mathbf{\text{k}}}_2=\mathbf{\text{f}}(x_n+\frac{h}{2},{\mathbf{\text{y}}}_n+\frac{h}{2}{\mathbf{\text{k}}}_1),$

${\mathbf{\text{k}}}_3=\mathbf{\text{f}}(x_n+\frac{h}{2},{\mathbf{\text{y}}}_n+\frac{h}{2}{\mathbf{\text{k}}}_2),$

${\mathbf{\text{k}}}_4=\mathbf{\text{f}}(x_n+h,{\mathbf{\text{y}}}_n+h{\mathbf{\text{k}}}_3).$

где $h$ — величина шага сетки по $x$.

Этот метод имеет четвёртый порядок точности. Это значит, что ошибка на одном шаге имеет порядок $O(h^5)$, а суммарная ошибка на конечном интервале интегрирования имеет порядок $O(h^4)$ .

Явные методы Рунге — Кутты

Семейство явных методов Рунге — Кутты является обобщением как явного метода Эйлера, так и классического метода Рунге — Кутты четвёртого порядка. Оно задаётся формулами

${\mathbf{\text{y}}}_{n+1}={\mathbf{\text{y}}}_n+h\sum _{i=1}^s{b}_i{\mathbf{\text{k}}}_i,$

где $h$ — величина шага сетки по $x$, а вычисление нового значения проходит в $s$ этапов:

$\begin{array}{cc}{\mathbf{\text{k}}}_1=& \mathbf{\text{f}}(x_n,{\mathbf{\text{y}}}_n),\\ {\mathbf{\text{k}}}_2=& \mathbf{\text{f}}(x_n+c_{2}h,{\mathbf{\text{y}}}_n+a_{21}h{\mathbf{\text{k}}}_1),\\ \cdots & \\ {\mathbf{\text{k}}}_s=& \mathbf{\text{f}}(x_n+c_{s}h,{\mathbf{\text{y}}}_n+a_{s1}h{\mathbf{\text{k}}}_1+a_{s2}h{\mathbf{\text{k}}}_2+\cdots +a_{s,s-1}h{\mathbf{\text{k}}}_{s-1})\end{array}$

Конкретный метод определяется числом $s$ и коэффициентами $b_i,a_{ij}$ и $c_i$. Эти коэффициенты часто упорядочивают в таблицу (называемую таблицей Бутчера):

$\begin{array}{cccccc}0& & & & & \\ c_2& a_{21}& & & & \\ c_3& a_{31}& a_{32}& & & \\ ⋮& ⋮& ⋮& \ddots & & \\ c_s& a_{s1}& a_{s2}& \dots & a_{ss-1}& \\ & b_1& b_2& \dots & b_{s-1}& b_s\end{array}$

Для коэффициентов метода Рунге — Кутты должны быть выполнены условия $\sum _{j=1}^{i-1}{a}_{ij}=c_i$ для $i=2,\dots ,s$. Если требуется, чтобы метод имел порядок $p$, то следует также обеспечить условие

$\overline{\mathbf{\text{y}}}(h+x_0)-\mathbf{\text{y}}(h+x_0)=O(h^{p+1}),$

где $\overline{\mathbf{\text{y}}}(h+x_0)$ — приближение, полученное по методу Рунге — Кутты. После многократного дифференцирования это условие преобразуется в систему полиномиальных уравнений относительно коэффициентов метода.

Неявные методы Рунге — Кутты

Все до сих пор упомянутые методы Рунге — Кутты являются Шаблон:Iw. К сожалению, явные методы Рунге — Кутты, как правило, непригодны для решения жёстких систем уравнений из-за малой области их абсолютной устойчивостиШаблон:Sfn. Неустойчивость явных методов Рунге — Кутты создаёт весьма серьёзные проблемы и при Шаблон:Iw.

Неустойчивость явных методов Рунге — Кутты мотивировала развитие неявных методов. Неявный метод Рунге — Кутты имеет видШаблон:SfnШаблон:Sfn:

$y_{n+1}=y_n+h\sum _{i=1}^s{b}_i{k}_i,$

где

$k_i=f(x_n+c_{i}h,y_n+h\sum _{j=1}^s{a}_{ij}{k}_j),i=1,\dots ,s.$

Явный метод характерен тем, что матрица коэффициентов $a_{ij}$ для него имеет нижний треугольный вид (включая и нулевую главную диагональ) — в отличие от неявного метода, где матрица имеет произвольный вид. Это также видно по Шаблон:Iw.

$\begin{array}{ccccc}{c}_1& a_{11}& a_{12}& \dots & a_{1s}\\ c_2& a_{21}& a_{22}& \dots & a_{2s}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ c_s& a_{s1}& a_{s2}& \dots & a_{ss}\\ & b_1& b_2& \dots & b_s\end{array}=\begin{array}{cc}𝐜& A\\ & {𝐛}^{𝐓}\end{array}$

Следствием этого различия является необходимость на каждом шагу решать систему уравнений для $k_i,i=1,2,...,s$, где $s$ — число стадий. Это увеличивает вычислительные затраты, однако при достаточно малом $h$ можно применить принцип сжимающих отображений и решать данную систему методом простой итерации. В случае одной итерации это увеличивает вычислительные затраты всего лишь в два раза.

С другой стороны, Шаблон:Iw (1961) и Шаблон:Iw (1964) показали, что при любом количестве стадий $s$ существует неявный метод Рунге — Кутты с порядком точности $p=2s$. Это значит, например, что для описанного выше явного четырёхстадийного метода четвёртого порядка существует неявный аналог с вдвое большим порядком точности.

Неявный метод Рунге — Кутты второго порядка

Простейшим неявным методом Рунге — Кутты является модифицированный метод Эйлера «с пересчётом». Он задаётся формулой:

${𝐲}_{n+1}={𝐲}_n+h\frac{𝐟(x_n,{𝐲}_n)+𝐟(x_{n+1},{𝐲}_{n+1})}{2}$.

Для его реализации на каждом шаге необходимы как минимум две итерации (и два вычисления функции).

Прогноз:

${\widetilde{𝐲}}_{n+1}={𝐲}_n+h𝐟(x_n,{𝐲}_n)$.

Коррекция:

${𝐲}_{n+1}={𝐲}_n+h\frac{𝐟(x_n,{𝐲}_n)+𝐟(x_{n+1},{\widetilde{𝐲}}_{n+1})}{2}$.

Вторая формула — это простая итерация решения системы уравнений относительно ${𝐲}_{n+1}$, записанной в форме сжимающего отображения. Для повышения точности итерацию-коррекцию можно сделать несколько раз, подставляя ${\widetilde{𝐲}}_{n+1}={𝐲}_{n+1}$. Модифицированный метод Эйлера «с пересчётом» имеет второй порядок точности.

Устойчивость

Преимуществом неявных методов Рунге — Кутты в сравнении с явными является их бо́льшая устойчивость, что особенно важно при решении жестких уравнений. Рассмотрим в качестве примера линейное уравнение y' = λy. Обычный метод Рунге — Кутты, применённый к этому уравнению, сведётся к итерации $y_{n+1}=r(h\lambda )y_n$, с r, равным

$r(z)=1+z{b}^T(I-zA{)}^{-1}e=\frac{\det (I-zA+ze{b}^T)}{\det (I-zA)},$

где $e$ обозначает вектор-столбец из единицШаблон:Sfn. Функция $r$ называется функцией устойчивостиШаблон:Sfn. Из формулы видно, что $r$ является отношением двух полиномов степени $s$, если метод имеет $s$ стадий. Явные методы имеют строго нижнюю треугольную матрицу $A,$ откуда следует, что $\det (I-zA)=1,$ и что функция устойчивости является многочленомШаблон:Sfn.

Численное решение данного примера сходится к нулю при условии $|r(z)|<1$ с $z=h\lambda$. Множество таких $z$ называется областью абсолютной устойчивости. В частности, метод называется A-устойчивым, если все $z$ с $\text{Re}(z)<0$ находятся в области абсолютной устойчивости. Функция устойчивости явного метода Рунге — Кутты является многочленом, поэтому явные методы Рунге — Кутты в принципе не могут быть A-устойчивымиШаблон:Sfn.

Если метод имеет порядок $p$, то функция устойчивости удовлетворяет условию $r(z)={\text{e}}^z+O(z^{p+1})$ при $z\to 0$. Таким образом, представляет интерес отношение многочленов данной степени, приближающее экспоненциальную функцию наилучшим образом. Эти отношения известны как аппроксимации Паде. Аппроксимация Паде с числителем степени $m$ и знаменателем степени $n$ А-устойчива тогда и только тогда, когда $m\le n\le m+2.$Шаблон:Sfn

$s$-стадийный метод Гаусса — Лежандра имеет порядок $2s$, поэтому его функция устойчивости является приближением Паде $m=n=s$. Отсюда следует, что метод является A-устойчивымШаблон:Sfn. Это показывает, что A-устойчивые методы Рунге — Кутты могут иметь сколь угодно высокий порядок. В отличие от этого, порядок А-устойчивости методов Адамса не может превышать два.

Произношение

Согласно грамматическим нормам русского языка, фамилия Ку́тта склоняется, поэтому говорят: «Метод Ру́нге — Ку́тты». Правила русской грамматики предписывают склонять все фамилии (в том числе и мужские), оканчивающиеся на -а, -я, которым предшествует согласный. Единственное исключение — фамилии французского происхождения с ударением на последнем слоге типа Дюма́, Золя́. Однако иногда встречается несклоняемый вариант «Метод Ру́нге — Ку́тта» (например, в книге).

Пример решения на алгоритмических языках программирования

Шаблон:EF производя замену $y^{\prime }=z$ и перенося $4y$ в правую часть, получаем систему:

Шаблон:Hider Шаблон:Hider В программе на С# используется абстрактный класс RungeKutta, в котором следует переопределить абстрактный метод F, задающий правые части уравнений.

Пример решения в среде MATLAB

Решение систем дифференциальных уравнений методом Рунге — Кутты является одним из самых распространённых численных методов решений в технике. В среде MATLAB реализована его одна из разновидностей — Шаблон:Iw. Для решения системы уравнений необходимо сначала записать функцию, вычисляющую производные, то есть функции y = g(x, y, z) и z = cos(3x) − 4y = f(x, y, z), о чём сказано выше. В одном из каталогов, к которому имеется доступ из системы MATLAB, нужно создать текстовый файл с именем (например) runge.m со следующим содержимым (для MATLAB версии 5.3): Шаблон:Hider Имя файла и имя функции должно совпадать, но оно может быть любым неиспользуемым ранее.

Затем необходимо создать главный файл c именем, например, main.m, который будет выполнять основные вычисления. Этот главный файл будет содержать следующий текст: Шаблон:Hider Так как MATLAB ориентирован на работу с матрицами, решение по методу Рунге — Кутты очень легко выполняется для целого ряда x, как, например, в приведённом примере программы. Здесь решение — график функции в пределах времён от 0 до x_fin.

Решение в среде MATLAB

Переменные x и y, полученные в результате работы функции ODE45, есть векторы значений. Очевидно, что решение конкретно заданного выше примера — второй элемент x, так как первое значение — 0, шаг интегрирования h = 0,1, а интересующее значение x = 0,1. Следующая запись в командном окне MATLAB даст искомое решение: Шаблон:Hider Ответ: y1 = 0,98768

Литература

• Шаблон:Книга
• Шаблон:Книга
• Шаблон:Книга