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

В численном анализе используются методы Рунге –Кутты ( англ.: / ˈ r ʊ ŋ ə ˈ k ʊ t ɑː / ^ⓘ RUUNG -ə- KUUT -tah^[1]) представляют собой семействонеявных и явныхитерационных методов, включающееметод Эйлера, используемый привременной дискретизациидля приближенных решенийодновременных нелинейных уравнений.^[2]Эти методы были разработаны около 1900 года немецкими математикамиКарлом РунгеиВильгельмом Куттой.

Наиболее известный член семейства методов Рунге–Кутты обычно называют «RK4», «классическим методом Рунге–Кутты» или просто «методом Рунге–Кутты».

Пусть задача с начальными данными задана следующим образом:

${\displaystyle {\frac {dy}{dt}}=f(t,y),\quad y(t_{0})=y_{0}.}$

Вот неизвестная функция (скалярная или векторная) времени , которую мы хотели бы аппроксимировать; нам говорят, что , скорость, с которой изменяется, является функцией и самой себя. В начальный момент времени соответствующее значение равно . Функция и начальные условия , даны.${\displaystyle y}$${\displaystyle t}$${\displaystyle {\frac {dy}{dt}}}$${\displaystyle y}$${\displaystyle t}$${\displaystyle y}$${\displaystyle t_{0}}$${\displaystyle y}$${\displaystyle y_{0}}$${\displaystyle f}$ ${\displaystyle t_{0}}$${\displaystyle y_{0}}$

Теперь выберем размер шага h > 0 и определим:

${\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+{\frac {h}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right),\\t_{n+1}&=t_{n}+h\\\end{aligned}}}$

для n = 0, 1, 2, 3, ..., используя ^[3]

${\displaystyle {\begin{aligned}k_{1}&=\ f(t_{n},y_{n}),\\k_{2}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{1}}{2}}\right),\\k_{3}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{2}}{2}}\right),\\k_{4}&=\ f\!\left(t_{n}+h,y_{n}+hk_{3}\right).\end{aligned}}}$

( Примечание: приведенные выше уравнения имеют разные, но эквивалентные определения в разных текстах. ^[4] )

Здесь представлена ​​аппроксимация RK4 для , а следующее значение ( ) определяется текущим значением ( ) плюс средневзвешенное значение четырех приращений, где каждое приращение является произведением размера интервала h и предполагаемого наклона, заданного функцией f в правой части дифференциального уравнения.${\displaystyle y_{n+1}}$${\displaystyle y(t_{n+1})}$${\displaystyle y_{n+1}}$${\displaystyle y_{n}}$

• ${\displaystyle k_{1}}$- наклон в начале интервала, с использованием ( метода Эйлера );${\displaystyle y}$
• ${\displaystyle k_{2}}$— наклон в средней точке интервала, с использованием и ;${\displaystyle y}$${\displaystyle k_{1}}$
• ${\displaystyle k_{3}}$снова наклон в средней точке, но теперь с использованием и ;${\displaystyle y}$${\displaystyle k_{2}}$
• ${\displaystyle k_{4}}$— наклон в конце интервала с использованием и .${\displaystyle y}$${\displaystyle k_{3}}$

При усреднении четырех наклонов больший вес придается наклонам в средней точке. Если не зависит от , так что дифференциальное уравнение эквивалентно простому интегралу, то RK4 — это правило Симпсона . ^[5]${\displaystyle f}$${\displaystyle y}$

Метод RK4 является методом четвертого порядка, что означает, что локальная ошибка усечения составляет порядка , тогда как общая накопленная ошибка составляет порядка .${\displaystyle O(h^{5})}$${\displaystyle O(h^{4})}$

Во многих практических приложениях функция независима от (так называемая автономная система или система, не зависящая от времени, особенно в физике), а их приращения вообще не вычисляются и не передаются в функцию , а используется только окончательная формула для .${\displaystyle f}$${\displaystyle t}$${\displaystyle f}$${\displaystyle t_{n+1}}$

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

${\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},}$

где ^[6]

${\displaystyle {\begin{aligned}k_{1}&=f(t_{n},y_{n}),\\k_{2}&=f(t_{n}+c_{2}h,y_{n}+(a_{21}k_{1})h),\\k_{3}&=f(t_{n}+c_{3}h,y_{n}+(a_{31}k_{1}+a_{32}k_{2})h),\\&\ \ \vdots \\k_{s}&=f(t_{n}+c_{s}h,y_{n}+(a_{s1}k_{1}+a_{s2}k_{2}+\cdots +a_{s,s-1}k_{s-1})h).\end{aligned}}}$

( Примечание: приведенные выше уравнения могут иметь разные, но эквивалентные определения в некоторых текстах. ^[4] )

Чтобы указать конкретный метод, необходимо указать целое число s (количество стадий) и коэффициенты a _ij (для 1 ≤ j < i ≤ s ), b _i (для i = 1, 2, ..., s ) и c _i (для i = 2, 3, ..., s ). Матрица [ a _ij ] называется матрицей Рунге–Кутты , в то время как b _i и c _i известны как веса и узлы . ^[7] Эти данные обычно располагаются в мнемоническом устройстве, известном как таблица Бутчера (в честь Джона К. Бутчера ):

${\displaystyle 0}$
${\displaystyle c_{2}}$	${\displaystyle a_{21}}$
${\displaystyle c_{3}}$	${\displaystyle a_{31}}$	${\displaystyle a_{32}}$
${\displaystyle \vdots }$	${\displaystyle \vdots }$		${\displaystyle \ddots }$
${\displaystyle c_{s}}$	${\displaystyle a_{s1}}$	${\displaystyle a_{s2}}$	${\displaystyle \cdots }$	${\displaystyle a_{s,s-1}}$
	${\displaystyle b_{1}}$	${\displaystyle b_{2}}$	${\displaystyle \cdots }$	${\displaystyle b_{s-1}}$	${\displaystyle b_{s}}$

Разложение в ряд Тейлора показывает , что метод Рунге–Кутты является последовательным тогда и только тогда, когда

${\displaystyle \sum _{i=1}^{s}b_{i}=1.}$

Существуют также сопутствующие требования, если требуется, чтобы метод имел определенный порядок p , что означает, что локальная ошибка усечения равна O( h ^{p +1} ). Они могут быть выведены из определения самой ошибки усечения. Например, двухэтапный метод имеет порядок 2, если b ₁ + b ₂ = 1, b ₂ c ₂ = 1/2 и b ₂ a ₂₁ = 1/2. ^[8] Обратите внимание, что популярным условием для определения коэффициентов является ^[8]

${\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ for }}i=2,\ldots ,s.}$

Однако одно это условие не является ни достаточным, ни необходимым для согласованности. ^[9]

В общем случае, если явный -этапный метод Рунге-Кутты имеет порядок , то можно доказать, что число этапов должно удовлетворять , а если , то . ^[10] Однако неизвестно, являются ли эти границы точными во всех случаях. В некоторых случаях доказано, что граница не может быть достигнута. Например, Бутчер доказал, что для , явного метода с этапами не существует . ^[11] Бутчер также доказал, что для , явного метода Рунге-Кутты с этапами не существует . ^[12] В общем случае, однако, остается открытой проблема, каково точное минимальное число этапов для явного метода Рунге-Кутты, чтобы иметь порядок . Некоторые известные значения: ^[13]${\displaystyle s}$${\displaystyle p}$${\displaystyle s\geq p}$${\displaystyle p\geq 5}$${\displaystyle s\geq p+1}$${\displaystyle p>6}$${\displaystyle s=p+1}$${\displaystyle p>7}$${\displaystyle p+2}$${\displaystyle s}$${\displaystyle p}$

${\displaystyle {\begin{array}{c|cccccccc}p&1&2&3&4&5&6&7&8\\\hline \min s&1&2&3&4&6&7&9&11\end{array}}}$

Доказуемая граница выше затем подразумевает, что мы не можем найти методы порядков , которые требуют меньше стадий, чем методы, которые мы уже знаем для этих порядков. Работа Бутчера также доказывает, что методы 7-го и 8-го порядка имеют минимум 9 и 11 стадий соответственно. ^{[11] [12]} Пример явного метода порядка 6 с 7 стадиями можно найти в ^[14] . Явные методы порядка 7 с 9 стадиями ^[11] и явные методы порядка 8 с 11 стадиями ^[15] также известны. См. ^{[16] [17]} для резюме.${\displaystyle p=1,2,\ldots ,6}$

Метод RK4 попадает в эту структуру. Его таблица ^[18]

0
1/2	1/2
1/2	0	1/2
1	0	0	1
	1/6	1/3	1/3	1/6

Небольшая вариация метода «Рунге–Кутта» также принадлежит Кутте в 1901 году и называется правилом 3/8. ^[19] Основное преимущество этого метода в том, что почти все коэффициенты ошибок меньше, чем в популярном методе, но он требует немного больше FLOP (операций с плавающей точкой) на шаг времени. Его таблица Бутчера

0
1/3	1/3
2/3	−1/3	1
1	1	−1	1
	1/8	3/8	3/8	1/8

Однако простейшим методом Рунге–Кутты является (прямой) метод Эйлера , заданный формулой . Это единственный последовательный явный метод Рунге–Кутты с одним этапом. Соответствующая таблица имеет вид${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n})}$

0
	1

Примером метода второго порядка с двумя этапами является явный метод средней точки :

${\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {1}{2}}h,y_{n}+{\frac {1}{2}}hf(t_{n},\ y_{n})\right).}$

Соответствующая таблица:

0
1/2	1/2
	0	1

Метод средней точки — не единственный метод Рунге–Кутты второго порядка с двумя этапами; существует семейство таких методов, параметризованных α и заданных формулой ^[20]

${\displaystyle y_{n+1}=y_{n}+h{\bigl (}(1-{\tfrac {1}{2\alpha }})f(t_{n},y_{n})+{\tfrac {1}{2\alpha }}f(t_{n}+\alpha h,y_{n}+\alpha hf(t_{n},y_{n})){\bigr )}.}$

Его картина Мясника - это

0
${\displaystyle \alpha }$	${\displaystyle \alpha }$
	${\displaystyle (1-{\tfrac {1}{2\alpha }})}$	${\displaystyle {\tfrac {1}{2\alpha }}}$

В этом семействе есть метод средней точки , есть метод Хойна , ^[5] и есть метод Ралстона.${\displaystyle \alpha ={\tfrac {1}{2}}}$${\displaystyle \alpha =1}$${\displaystyle \alpha ={\tfrac {2}{3}}}$

В качестве примера рассмотрим двухэтапный метод Рунге-Кутты второго порядка с α = 2/3, также известный как метод Ральстона . Он задается таблицей

с соответствующими уравнениями

${\displaystyle {\begin{aligned}k_{1}&=f(t_{n},\ y_{n}),\\k_{2}&=f(t_{n}+{\tfrac {2}{3}}h,\ y_{n}+{\tfrac {2}{3}}hk_{1}),\\y_{n+1}&=y_{n}+h\left({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2}\right).\end{aligned}}}$

Этот метод используется для решения задачи начального значения

${\displaystyle {\frac {dy}{dt}}=\tan(y)+1,\quad y_{0}=1,\ t\in [1,1.1]}$

с размером шага h = 0,025, поэтому метод должен выполнить четыре шага.

Метод осуществляется следующим образом:

${\displaystyle t_{0}=1\colon }$
	${\displaystyle y_{0}=1}$
${\displaystyle t_{1}=1.025\colon }$
	${\displaystyle y_{0}=1}$	${\displaystyle k_{1}=2.557407725}$	${\displaystyle k_{2}=f(t_{0}+{\tfrac {2}{3}}h,\ y_{0}+{\tfrac {2}{3}}hk_{1})=2.7138981400}$
	${\displaystyle y_{1}=y_{0}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.066869388}}}$
${\displaystyle t_{2}=1.05\colon }$
	${\displaystyle y_{1}=1.066869388}$	${\displaystyle k_{1}=2.813524695}$	${\displaystyle k_{2}=f(t_{1}+{\tfrac {2}{3}}h,\ y_{1}+{\tfrac {2}{3}}hk_{1})}$
	${\displaystyle y_{2}=y_{1}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.141332181}}}$
${\displaystyle t_{3}=1.075\colon }$
	${\displaystyle y_{2}=1.141332181}$	${\displaystyle k_{1}=3.183536647}$	${\displaystyle k_{2}=f(t_{2}+{\tfrac {2}{3}}h,\ y_{2}+{\tfrac {2}{3}}hk_{1})}$
	${\displaystyle y_{3}=y_{2}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.227417567}}}$
${\displaystyle t_{4}=1.1\colon }$
	${\displaystyle y_{3}=1.227417567}$	${\displaystyle k_{1}=3.796866512}$	${\displaystyle k_{2}=f(t_{3}+{\tfrac {2}{3}}h,\ y_{3}+{\tfrac {2}{3}}hk_{1})}$
	${\displaystyle y_{4}=y_{3}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.335079087}}.}$

Численные решения соответствуют подчеркнутым значениям.

Адаптивные методы предназначены для получения оценки локальной ошибки усечения одного шага Рунге-Кутты. Это делается с помощью двух методов, одного с порядком и одного с порядком . Эти методы переплетены, т. е. имеют общие промежуточные шаги. Благодаря этому оценка ошибки имеет небольшие или пренебрежимо малые вычислительные затраты по сравнению с шагом с методом более высокого порядка.${\displaystyle p}$${\displaystyle p-1}$

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

Шаг низшего порядка определяется как

${\displaystyle y_{n+1}^{*}=y_{n}+h\sum _{i=1}^{s}b_{i}^{*}k_{i},}$

где те же, что и для метода высшего порядка. Тогда ошибка равна${\displaystyle k_{i}}$

${\displaystyle e_{n+1}=y_{n+1}-y_{n+1}^{*}=h\sum _{i=1}^{s}(b_{i}-b_{i}^{*})k_{i},}$

что равно . Таблица Бутчера для этого вида метода расширена, чтобы дать значения :${\displaystyle O(h^{p})}$${\displaystyle b_{i}^{*}}$

	0
	${\displaystyle c_{2}}$	${\displaystyle a_{21}}$
	${\displaystyle c_{3}}$	${\displaystyle a_{31}}$	${\displaystyle a_{32}}$
	${\displaystyle \vdots }$	${\displaystyle \vdots }$		${\displaystyle \ddots }$
	${\displaystyle c_{s}}$	${\displaystyle a_{s1}}$	${\displaystyle a_{s2}}$	${\displaystyle \cdots }$	${\displaystyle a_{s,s-1}}$
		${\displaystyle b_{1}}$	${\displaystyle b_{2}}$	${\displaystyle \cdots }$	${\displaystyle b_{s-1}}$	${\displaystyle b_{s}}$
		${\displaystyle b_{1}^{*}}$	${\displaystyle b_{2}^{*}}$	${\displaystyle \cdots }$	${\displaystyle b_{s-1}^{*}}$	${\displaystyle b_{s}^{*}}$

Метод Рунге –Кутты–Фельберга имеет два метода порядков 5 и 4. Его расширенная таблица Бутчера имеет вид:

Однако простейший адаптивный метод Рунге–Кутты включает в себя объединение метода Хойна , который имеет порядок 2, с методом Эйлера , который имеет порядок 1. Его расширенная таблица Бутчера имеет вид:

Другими адаптивными методами Рунге–Кутты являются метод Богацки–Шампайна (порядки 3 и 2), метод Кэша–Карпа и метод Дорманда–Принса (оба с порядками 5 и 4).

Метод Рунге–Кутты называется неконфлюэнтным ^[21], если все различны.${\displaystyle c_{i},\,i=1,2,\ldots ,s}$

Методы Рунге–Кутты–Нистрёма — это специализированные методы Рунге–Кутты, оптимизированные для дифференциальных уравнений второго порядка. ^{[22] [23]} Общий метод Рунге–Кутты–Нистрёма для системы ОДУ второго порядка

$${\displaystyle {\ddot {y}}_{i}=f_{i}(y_{1},y_{2},\ldots ,y_{n})}$$

с заказом связано с формой ${\displaystyle s}$

$${\displaystyle {\begin{cases}g_{i}=y_{m}+c_{i}h{\dot {y}}_{m}+h^{2}\sum _{j=1}^{s}a_{ij}f(g_{j}),&i=1,2,\ldots ,s\\y_{m+1}=y_{m}+h{\dot {y}}_{m}+h^{2}\sum _{j=1}^{s}{\bar {b}}_{j}f(g_{j})\\{\dot {y}}_{m+1}={\dot {y}}_{m}+h\sum _{j=1}^{s}b_{j}f(g_{j})\end{cases}}}$$

который образует стол Мясника с формой

$${\displaystyle {\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 &{\bar {b}}_{1}&{\bar {b}}_{2}&\dots &{\bar {b}}_{s}\\&b_{1}&b_{2}&\dots &b_{s}\end{array}}={\begin{array}{c|c}\mathbf {c} &\mathbf {A} \\\hline &\mathbf {\bar {b}} ^{\top }\\&\mathbf {b} ^{\top }\end{array}}}$$

Два явных метода РКН четвертого порядка задаются следующими таблицами Бутчера:

$${\displaystyle {\begin{array}{c|ccc}c_{i}&&a_{ij}&\\{\frac {3+{\sqrt {3}}}{6}}&0&0&0\\{\frac {3-{\sqrt {3}}}{6}}&{\frac {2-{\sqrt {3}}}{12}}&0&0\\{\frac {3+{\sqrt {3}}}{6}}&0&{\frac {\sqrt {3}}{6}}&0\\\hline {\overline {b_{i}}}&{\frac {5-3{\sqrt {3}}}{24}}&{\frac {3+{\sqrt {3}}}{12}}&{\frac {1+{\sqrt {3}}}{24}}\\\hline b_{i}&{\frac {3-2{\sqrt {3}}}{12}}&{\frac {1}{2}}&{\frac {3+2{\sqrt {3}}}{12}}\end{array}}}$$$${\displaystyle {\begin{array}{c|ccc}c_{i}&&a_{ij}&\\{\frac {3-{\sqrt {3}}}{6}}&0&0&0\\{\frac {3+{\sqrt {3}}}{6}}&{\frac {2+{\sqrt {3}}}{12}}&0&0\\{\frac {3-{\sqrt {3}}}{6}}&0&-{\frac {\sqrt {3}}{6}}&0\\\hline {\overline {b_{i}}}&{\frac {5+3{\sqrt {3}}}{24}}&{\frac {3-{\sqrt {3}}}{12}}&{\frac {1-{\sqrt {3}}}{24}}\\\hline b_{i}&{\frac {3+2{\sqrt {3}}}{12}}&{\frac {1}{2}}&{\frac {3-2{\sqrt {3}}}{12}}\end{array}}}$$

Эти две схемы также обладают симплектическими свойствами сохранения, когда исходное уравнение выводится из консервативной классической механической системы, т.е. когда

$${\displaystyle f_{i}(x_{1},\ldots ,x_{n})={\frac {\partial V}{\partial x_{i}}}(x_{1},\ldots ,x_{n})}$$

для некоторой скалярной функции . ^[24]${\displaystyle V}$

Все методы Рунге–Кутты, упомянутые до сих пор, являются явными методами . Явные методы Рунге–Кутты, как правило, не подходят для решения жестких уравнений , поскольку их область абсолютной устойчивости мала; в частности, она ограничена. ^[25] Этот вопрос особенно важен при решении уравнений в частных производных .

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

${\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},}$

где

${\displaystyle k_{i}=f\left(t_{n}+c_{i}h,\ y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}\right),\quad i=1,\ldots ,s.}$ ^[26]

Разница с явным методом заключается в том, что в явном методе сумма по j увеличивается только до i − 1. Это также проявляется в таблице Бутчера: матрица коэффициентов явного метода является нижней треугольной. В неявном методе сумма по j увеличивается до s , а матрица коэффициентов не является треугольной, что дает таблицу Бутчера вида ^[18]${\displaystyle a_{ij}}$

${\displaystyle {\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}\\&b_{1}^{*}&b_{2}^{*}&\dots &b_{s}^{*}\\\end{array}}={\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b^{T}} \\\end{array}}}$

Объяснение этой строки см. в разделе «Адаптивные методы Рунге-Кутты» выше.${\displaystyle b^{*}}$

Следствием этого различия является то, что на каждом шаге необходимо решить систему алгебраических уравнений. Это значительно увеличивает вычислительные затраты. Если метод с s этапами используется для решения дифференциального уравнения с m компонентами, то система алгебраических уравнений имеет ms компонентов. Это можно противопоставить неявным линейным многошаговым методам (другое большое семейство методов для ОДУ): неявный s -шаговый линейный многошаговый метод должен решить систему алгебраических уравнений только с m компонентами, поэтому размер системы не увеличивается с увеличением числа шагов. ^[27]

Простейшим примером неявного метода Рунге–Кутты является обратный метод Эйлера :

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

Таблица Мясника для этого проста:

${\displaystyle {\begin{array}{c|c}1&1\\\hline &1\\\end{array}}}$

Эта таблица Бутчера соответствует формулам

${\displaystyle k_{1}=f(t_{n}+h,\ y_{n}+hk_{1})\quad {\text{and}}\quad y_{n+1}=y_{n}+hk_{1},}$

которую можно переформулировать, чтобы получить формулу для обратного метода Эйлера, указанную выше.

Другим примером неявного метода Рунге–Кутты является правило трапеций . Его таблица Бутчера имеет вид:

${\displaystyle {\begin{array}{c|cc}0&0&0\\1&{\frac {1}{2}}&{\frac {1}{2}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\&1&0\\\end{array}}}$

Правило трапеции — это метод коллокации (как обсуждалось в этой статье). Все методы коллокации — это неявные методы Рунге–Кутты, но не все неявные методы Рунге–Кутты являются методами коллокации. ^[28]

Методы Гаусса –Лежандра образуют семейство методов коллокации, основанных на квадратуре Гаусса . Метод Гаусса–Лежандра с s этапами имеет порядок 2 s (таким образом, могут быть построены методы с произвольно высоким порядком). ^[29] Метод с двумя этапами (и, следовательно, с порядком четыре) имеет таблицу Бутчера:

${\displaystyle {\begin{array}{c|cc}{\frac {1}{2}}-{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}&{\frac {1}{4}}-{\frac {1}{6}}{\sqrt {3}}\\{\frac {1}{2}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\&{\frac {1}{2}}+{\frac {1}{2}}{\sqrt {3}}&{\frac {1}{2}}-{\frac {1}{2}}{\sqrt {3}}\end{array}}}$ ^[27]

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

${\displaystyle r(z)=1+zb^{T}(I-zA)^{-1}e={\frac {\det(I-zA+zeb^{T})}{\det(I-zA)}},}$ ^[30]

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

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

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

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

Концепция A-устойчивости для решения дифференциальных уравнений связана с линейным автономным уравнением . Далквист (1963) предложил исследование устойчивости численных схем применительно к нелинейным системам, удовлетворяющим условию монотонности. Соответствующие концепции были определены как G-устойчивость для многошаговых методов (и связанных с ними одноногих методов) и B-устойчивость (Бутчер, 1975) для методов Рунге–Кутты. Метод Рунге–Кутты, примененный к нелинейной системе , который проверяет , называется B-устойчивым , если это условие подразумевает для двух численных решений.${\displaystyle y'=\lambda y}$${\displaystyle y'=f(y)}$${\displaystyle \langle f(y)-f(z),\ y-z\rangle \leq 0}$${\displaystyle \|y_{n+1}-z_{n+1}\|\leq \|y_{n}-z_{n}\|}$

Пусть , и — три матрицы, определяемые с помощью Метод Рунге–Кутты называется алгебраически устойчивым ^[36], если матрицы и обе неотрицательно определены. Достаточным условием B-устойчивости ^[37] является: и неотрицательно определены.${\displaystyle B}$${\displaystyle M}$${\displaystyle Q}$${\displaystyle s\times s}$$${\displaystyle {\begin{aligned}B&=\operatorname {diag} (b_{1},b_{2},\ldots ,b_{s}),\\[4pt]M&=BA+A^{T}B-bb^{T},\\[4pt]Q&=BA^{-1}+A^{-T}B-A^{-T}bb^{T}A^{-1}.\end{aligned}}}$$${\displaystyle B}$${\displaystyle M}$${\displaystyle B}$${\displaystyle Q}$

В общем случае метод порядка Рунге–Кутты можно записать так:${\displaystyle s}$

${\displaystyle y_{t+h}=y_{t}+h\cdot \sum _{i=1}^{s}a_{i}k_{i}+{\mathcal {O}}(h^{s+1}),}$

где:

${\displaystyle k_{i}=\sum _{j=1}^{s}\beta _{ij}f(k_{j},\ t_{n}+\alpha _{i}h)}$

являются приращениями, полученными путем оценки производных в -ом порядке.${\displaystyle y_{t}}$${\displaystyle i}$

Мы разрабатываем вывод ^[38] для метода Рунге–Кутты четвертого порядка, используя общую формулу с оценками, как объяснено выше, в начальной точке, средней точке и конечной точке любого интервала ; таким образом, мы выбираем:${\displaystyle s=4}$${\displaystyle (t,\ t+h)}$

${\displaystyle {\begin{aligned}&\alpha _{i}&&\beta _{ij}\\\alpha _{1}&=0&\beta _{21}&={\frac {1}{2}}\\\alpha _{2}&={\frac {1}{2}}&\beta _{32}&={\frac {1}{2}}\\\alpha _{3}&={\frac {1}{2}}&\beta _{43}&=1\\\alpha _{4}&=1&&\\\end{aligned}}}$

и в противном случае. Начнем с определения следующих величин:${\displaystyle \beta _{ij}=0}$

${\displaystyle {\begin{aligned}y_{t+h}^{1}&=y_{t}+hf\left(y_{t},\ t\right)\\y_{t+h}^{2}&=y_{t}+hf\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)\\y_{t+h}^{3}&=y_{t}+hf\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)\end{aligned}}}$

где и если мы определим:${\displaystyle y_{t+h/2}^{1}={\dfrac {y_{t}+y_{t+h}^{1}}{2}}}$${\displaystyle y_{t+h/2}^{2}={\dfrac {y_{t}+y_{t+h}^{2}}{2}}.}$

${\displaystyle {\begin{aligned}k_{1}&=f(y_{t},\ t)\\k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right)\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hk_{3},\ t+h\right)\end{aligned}}}$

и для предыдущих соотношений можно показать, что справедливы следующие равенства : где: — полная производная по времени.${\displaystyle {\mathcal {O}}(h^{2})}$$${\displaystyle {\begin{aligned}k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}f\left(y_{t},\ t\right),\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t},\ t\right)+h{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\right]\end{aligned}}}$$$${\displaystyle {\frac {d}{dt}}f(y_{t},\ t)={\frac {\partial }{\partial y}}f(y_{t},\ t){\dot {y}}_{t}+{\frac {\partial }{\partial t}}f(y_{t},\ t)=f_{y}(y_{t},\ t){\dot {y}}_{t}+f_{t}(y_{t},\ t):={\ddot {y}}_{t}}$$${\displaystyle f}$

Если теперь выразить общую формулу, используя то, что мы только что вывели, то получим:$${\displaystyle {\begin{aligned}y_{t+h}={}&y_{t}+h\left\lbrace a\cdot f(y_{t},\ t)+b\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right.+\\&{}+c\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]+\\&{}+d\cdot \left[f(y_{t},\ t)+h{\frac {d}{dt}}\left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f(y_{t},\ t)+\left.{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]\right]\right\rbrace +{\mathcal {O}}(h^{5})\\={}&y_{t}+a\cdot hf_{t}+b\cdot hf_{t}+b\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+c\cdot hf_{t}+c\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+\\&{}+c\cdot {\frac {h^{3}}{4}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot hf_{t}+d\cdot h^{2}{\frac {df_{t}}{dt}}+d\cdot {\frac {h^{3}}{2}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot {\frac {h^{4}}{4}}{\frac {d^{3}f_{t}}{dt^{3}}}+{\mathcal {O}}(h^{5})\end{aligned}}}$$

и сравнивая это с рядом Тейлора примерно :${\displaystyle y_{t+h}}$${\displaystyle t}$$${\displaystyle {\begin{aligned}y_{t+h}&=y_{t}+h{\dot {y}}_{t}+{\frac {h^{2}}{2}}{\ddot {y}}_{t}+{\frac {h^{3}}{6}}y_{t}^{(3)}+{\frac {h^{4}}{24}}y_{t}^{(4)}+{\mathcal {O}}(h^{5})=\\&=y_{t}+hf(y_{t},\ t)+{\frac {h^{2}}{2}}{\frac {d}{dt}}f(y_{t},\ t)+{\frac {h^{3}}{6}}{\frac {d^{2}}{dt^{2}}}f(y_{t},\ t)+{\frac {h^{4}}{24}}{\frac {d^{3}}{dt^{3}}}f(y_{t},\ t)\end{aligned}}}$$

получаем систему ограничений на коэффициенты:

${\displaystyle {\begin{cases}&a+b+c+d=1\\[6pt]&{\frac {1}{2}}b+{\frac {1}{2}}c+d={\frac {1}{2}}\\[6pt]&{\frac {1}{4}}c+{\frac {1}{2}}d={\frac {1}{6}}\\[6pt]&{\frac {1}{4}}d={\frac {1}{24}}\end{cases}}}$

что при решении дает то, что указано выше.${\displaystyle a={\frac {1}{6}},b={\frac {1}{3}},c={\frac {1}{3}},d={\frac {1}{6}}}$

• Метод Эйлера
• Список методов Рунге–Кутты
• Численные методы для обыкновенных дифференциальных уравнений
• Метод Рунге–Кутты (SDE)
• Общие линейные методы
• Интегратор группы Ли