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

Материал из Википедии — свободной энциклопедии
Перейти к: навигация, поиск

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

Формально, методом Рунге — Кутты является модифицированный и исправленный метод Эйлера, они представляют собой схемы второго порядка точности. Существуют стандартные схемы третьего порядка, не получившие широкого распространения. Наиболее часто используется и реализована в различных математических пакетах (Maple, MathCAD, Maxima) стандартная схема четвёртого порядка. Иногда при выполнении расчётов с повышенной точностью применяются схемы пятого и шестого порядков[1][2]. Построение схем более высокого порядка сопряжено с большими вычислительными трудностями[3]. Методы седьмого порядка должны иметь по меньшей мере девять стадий, в схему восьмого порядка входит 11 стадий. Хотя схемы девятого порядка не имеют большой практической значимости, неизвестно, сколько стадий необходимо для достижения этого порядка точности. Аналогичная задача существует для схем десятого и более высоких порядков[3].

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

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

Рассмотрим задачу Коши

\textbf{y}'=\textbf{f}(x,\textbf{y}), \quad \textbf{y}(x_0)=\textbf{y}_0.

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

 \textbf{y}_{n+1} = \textbf{y}_n + {h \over 6} (\textbf{k}_1 + 2\textbf{k}_2 + 2\textbf{k}_3 + \textbf{k}_4)

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

 \textbf{k}_1 = \textbf{f} \left( x_n, \textbf{y}_n \right),
 \textbf{k}_2 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_1 \right),
 \textbf{k}_3 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_2 \right),
 \textbf{k}_4 = \textbf{f} \left( x_n + h, \textbf{y}_n + h\ \textbf{k}_3 \right).

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

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

Прямые методы Рунге — Кутты[править | править вики-текст]

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

 \textbf{y}_{n+1} = \textbf{y}_n + \sum_{i=1}^s b_i \textbf{k}_i,

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

\begin{array}{ll}
\textbf{k}_1 =& h\textbf{f}(x_n, \textbf{y}_n),\\
\textbf{k}_2 =& h\textbf{f}(x_n+c_2h, \textbf{y}_n+a_{21}\textbf{k}_1),\\
\cdots&\\
\textbf{k}_s =& h\textbf{f}(x_n+c_sh, \textbf{y}_n+a_{s1}\textbf{k}_1+a_{s2}\textbf{k}_2+\cdots+a_{s,s-1}\textbf{k}_{s-1})
\end{array}

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

\begin{array}{c|ccccc}
  0      &&&&&\\
  c_2    & a_{21} &&&&\\
  c_3    & a_{31} & a_{32} &&&\\
  \vdots & \vdots & \vdots& \ddots&&\\
  c_s    & a_{s1} & a_{s2}& \dots & a_{ss-1}&\\
  \hline & b_1    & b_2   & \dots & b_{s-1} & b_s
\end{array}

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

\bar{\textbf{y}}(h+x_0)- {\textbf{y}}(h+x_0)=O(h^{p+1}),

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

Неявные методы Рунге-Кутты[править | править вики-текст]

Все до сих пор упомянутые методы Рунге-Кутты являются явными методами. К сожалению, явные методы Рунге-Кутты, как правило, непригодны для решения жестких уравнений, из-за малой области абсолютной устойчивости. В частности, это описано в [4]. Этот вопрос особенно важен при решении дифференциальных уравнений в частных производных.

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

 y_{n+1} = y_n + \sum_{i=1}^s b_i k_i,

где

 k_i = hf\bigl( x_n + c_i h, y_n + \sum_{j=1}^s a_{ij} k_j \bigr), \quad i = 1, \ldots, s. [5]

Явный метод характерен тем, что матрица коэффициентов a_{ij} для него имеет нижний треугольный вид, в отличие от неявного метода, где матрица имеет произвольный вид. Это также видно по таблице Бутчера.


\begin{array}{c|cccc}
c_1    & a_{11} & a_{12}& \dots & a_{1s}\\
c_2    & a_{21} & a_{22}& \dots & a_{2s}\\
\vdots & \vdots & \vdots& \ddots& \vdots\\
c_s    & a_{s1} & a_{s2}& \dots & a_{ss} \\
\hline
       & b_1    & b_2   & \dots & b_s\\
\end{array} = 

\begin{array}{c|c}
\mathbf{c}& A\\
\hline
          & \mathbf{b^T} \\
\end{array}

Следствием этого различия является необходимость на каждом шагу решать систему алгебраических уравнений. Это существенно увеличивает вычислительные затраты, хотя существуют весьма эффективные итерационные методы решений алгебраических систем такого типа. Вместе с тем, при правильном выборе узлов и коэффициентов в неявной схеме можно обеспечить более высокий порядок точности метода (вдвое больше количества узлов: "стадий").

Устойчивость[править | править вики-текст]

Преимуществом неявных методов Рунге-Кутты в сравнении с явными является их большая устойчивость, особенно в применении к жестким уравнениям. Рассмотрим в качестве примера линейное уравнение 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+zeb^T)}{\det(I-zA)}, [6]

где e обозначает вектор единиц. Функция r называется функцией устойчивости [7] Из формулы видно, что r является отношение двух полиномов степени s, если метод имеет s стадий. Явные методы имеют строго нижнюю треугольную матрицу A , откуда следует, что det(IzA) = 1, и что функция устойчивости является многочленом.[8]

Численное решение данного примера дает чистый ноль при условии | r(z) | < 1 с z = hλ. Множество таких r называется областью абсолютной устойчивости. В частности , метод называется A-стабильным если все r с Re(z) < 0 находятся в области абсолютной стабильности. Функция устойчивости явного метода Рунге-Кутта является многочленом, поэтому явные методы Рунге-Кутты в принципе не могут быть стабильными.[8]

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

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

Произношение[править | править вики-текст]

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

Решение систем ОДУ[править | править вики-текст]

Метод Ру́нге — Ку́тты непосредственно обобщается на случай систем обыкновенных дифференциальных уравнений путём записи системы и метода в векторной форме.

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

y'' + 4y = \cos 3x,\quad y(0) = 0.8,\ y'(0) = 2,\ x \in [0,1],\ h = 0.1

производя замену y'=z и перенося 4y в правую часть, получаем систему:

\begin{cases}
y' = z = g(x,y,z)\\
z' = \cos(3x) - 4y = f(x,y,z)
\end{cases}

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

Пример решения в среде MATLAB[править | править вики-текст]

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

Имя файла и имя функции должно совпадать, но оно может быть любым неиспользуемым ранее.

Затем необходимо создать главный файл c именем, например, main.m, который будет выполнять основные вычисления. Этот главный файл будет содержать следующий текст:

Так как MATLAB ориентирован на работу с матрицами, решение по методу Рунге-Кутты очень легко выполняется для целого ряда x как, например, в приведенном примере программы. Здесь решение - график функции в пределах времён от 0 до x_fin.

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

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

Ответ: y1 = 0.98768

Ссылки[править | править вики-текст]

  1. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. — М.: Бином, 2001 — с. 363—375.
  2. Ильина В. А., Силаев П. К. Численные методы для физиков-теоретиков. т. 2. — Москва-Ижевск: Институт компьютерных исследований, 2004. — с. 16-30.
  3. 1 2 J. C. Butcher. Numerical Methods for Ordinary Differential Equations. The University of Auckland, New Zealand.
  4. Süli & Mayers 2003, pp. 349–351
  5. Iserles 1996, С. 41; Süli & Mayers 2003, pp. 351–352
  6. Hairer & Wanner 1996, pp. 40–41
  7. Hairer & Wanner 1996, С. 40
  8. 1 2 Iserles 1996, С. 60
  9. Iserles 1996, pp. 62–63
  10. Iserles 1996, С. 63
  11. Б. П. Демидович, И. А. Марон, Э. З. Шувалова. Численные методы анализа, 3-е изд. — М.: Наука, 1967.