Участник:Борис/Алгоритм Рунге-Кутты

Метод Рунге-Кутты 4 порядка
Последовательный алгоритм
Последовательная сложность	[math]O()[/math]
Объём входных данных	[math][/math]
Объём выходных данных	[math][/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O()[/math]
Ширина ярусно-параллельной формы	[math]O()[/math]

Основные авторы описания: Е.П.Степанов(раздел 1.1, раздел 1.2), Б.Савков ()

1 ЧАСТЬ. Свойства и структура алгоритмов

1.1 Общее описание алгоритма

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

Метод Рунге - Кутты применяется при решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка:

[math] \begin{cases} \textbf{y}' = \textbf{f}(x, \textbf{y}),\\ \textbf{y}(x_0) = \textbf{y}_0. \end{cases} [/math]

При этом [math]\textbf{y}, \textbf{f}, x[/math] отвечают следующим свойствам: [math] \textbf{y}, \textbf{f} \in \mathbb{R}^m, x \in \mathbb{R}^1 [/math].

Пусть известно значение y(x) и требуется вычислить значение y(x+h). Это можно сделать при помощи формулы:

[math] y(x + h) = y(x) + \int_0^hy'(x + t)dt[/math]

Аппроксимируя интеграл в правой части, можно получить методы Рунге-Кутты разного порядка точности. В общем виде формулы методов Рунге-Кутты можно записать следующим образом:

[math]\textbf{y}_{n+1} = \textbf{y}_n + h\sum_{i=1}^s b_i \textbf{k}_i[/math]

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

Увеличивая число s вспомогательных точек, можно построить методы Рунге-Кутты любого порядка точности p. Однако уже при p > 5 эти методы используются довольно редко. Это объясняется как чрезмерной громоздкостью получающихся вычислительных формул, так и тем, что преимущества методов высокого порядка точности над методами, в которых p = 4 или p = 5, проявляются либо в тех задачах, где нужна очень высокая точность и используются ЭВМ высокой разрядности, либо в тех задачах, где решение очень гладкое. Кроме того, методы Рунге-Кутты высокого порядка точности часто оказываются менее эффективными по сравнению с методами Адамса того же порядка точности. На этой странице представлен наиболее распространенный метод Рунге-Кутты 4 порядка с ошибкой на конечном интервале интегрирования [math]O(h^4)[/math].

Существуют как явные, так и неявные методы Рунге-Кутты. Последние имеют ряд преимуществ над явными методами, однако это достигается за счет существенного усложнения вычислительного алгоритма, так как на каждом шаге приходится считать систему нелинейных уравнений. На этой странице представлен явный метод Рунге-Кутты.

Можно отметить, что в современных программах, реализующих метод Рунге-Кутты, обязательно используется некоторый алгоритм автоматического изменения шага интегрирования. Интуитивно ясно, что на участках плавного изменения решения счет можно вести с достаточно крупным шагом. В то же время на тех участках, где происходят резкие изменения поведения решения, необходимо выбирать мелкий шаг интегрирования. Обычно начальное значение шага h1 задает пользователь. Далее шаг интегрирования меняется в соответствии с величиной получаемой в ходе вычислений оценки локальной погрешности. Для упрощения описания на этой странице используется постоянный шаг интегрирования h, зависящий от размера сетки разбиения n.

1.2 Математическое описание алгоритма

Далее подразумевается, что [math] \bold y, \bold y_i,\bold f, \bold k_i \in \mathbb{R}^m[/math], а [math]x, h, n, a, b \in \mathbb{R}^1 [/math].

Исходные данные:

   [math] a, b [/math] - интервал интегрирования,
   [math] \bold y_0 [/math] - вектор значений при [math] x = a [/math], 
   n - количество точек в разбиении
   [math]\bold f [/math] - функция вычисления производной

Вычисляемые данные:

   [math] \bold y_i [/math] - вектора значений в узлах сетки, т.е. при [math] x = a + i * \frac{b-a}{n}, i \in [1, n] [/math]

Формулы метода:

   Вначале вычисляется шаг сетки
      [math] h = \frac{b-a}{n}[/math]

   Приближенное значение в следующих точках вычисляется по итерационной формуле:
      [math]\textbf{y}_{n+1} = \textbf{y}_n + {h \over 6}(\textbf{k}_1 + 2\textbf{k}_2 + 2\textbf{k}_3 + \textbf{k}_4)[/math],
      [math]\textbf{k}_1 = \textbf{f} \left( x_n, \textbf{y}_n \right),[/math]
      [math]\textbf{k}_2 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_1 \right),[/math]
      [math]\textbf{k}_3 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_2 \right),[/math]
      [math]\textbf{k}_4 = \textbf{f} \left( x_n + h, \textbf{y}_n + h\ \textbf{k}_3 \right),[/math]
      [math]x_{n+1} = x_n + h.[/math]