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

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

Основные авторы описания: Е.П.Степанов(раздел 1.1, раздел 1.2, раздел 1.3, раздел 1.4, раздел 1.9), Б.Савков ( раздел 1.5, раздел 1.6, раздел 1.7, раздел 1.8, раздел 1.10, раздел 2.7)

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]

1.3 Вычислительное ядро алгоритма

Вычислительным ядром метода Рунге-Кутты можно считать вычисление значений функции производной

[math]\textbf{f}(x, \textbf{y})[/math],

так как при вычислении каждого нового [math]\textbf{y}_n[/math] требуется 4 вычисления значения функции производной. Таким образом, сложность работы алгоритма существенно зависит от сложности вычисления функции [math]\textbf{f}[/math].

1.4 Макроструктура алгоритма

Макроструктура алгоритма может быть представлена в следующем виде:

Вычисление вспомогательных коэффициентов [math]\textbf{k}_j, j \in [1,4][/math];

Вычисление приближенного значения [math] \bold y_i [/math] в следующей точке;

Переход на следующую итерацию или остановка.

1.5 Схема реализации последовательного алгоритма

Пример реализации последовательного алгоритма на языке C++ (указатель на функцию f передаётся через аргументы представленной ниже функции):

 1 #include <vector>
 2 
 3 // other code ...
 4 
 5 std::vector<double> rungekutta4(double a, double b, std::vector<double> y0, int n, std::vector<double (*)(double,double)> f)
 6 {
 7     std::vector<double>::size_type sz = y0.size();
 8     std::vector<double> result;
 9     for (unsigned j=0; j<sz; j++) 
10     {
11         rezult.push_back(core(a,b,y[j],n,f[j]));
12     }
13     return result;
14 }
15 
16 double core(double a, double b, double y0, int n, double (*f)(double,double))
17 {
18     double h, k1, k2, k3, k4, x, y;
19     h = (b - a)/n;
20     y = y0;
21     x = a;
22     for (unsigned i = 0; i < n; i++)
23     {
24         k1 = f(x, y);
25         k2 = f(x + h/2, y + h*k1/2);
26         k3 = f(x + h/2, y + h*k2/2);
27         k4 = f(x + h/2, y + h*k3);
28         x += h;
29         y += h*(k1 + 2*k2 + 2*k3 + k4)/6;
30     }
31     return y;
32 }
33 
34 // other code ...

1.6 Последовательная сложность алгоритма

Рассмотрим сначала одномерный случай. Наравне с обычными операциями, такими как сложение, умножение и тому подобными, в представленном выше алгоритме присутствует вызов функции [math] f [/math]. Так как эта функция может содержать в себе любое количество простых математических операций, будем считать, что её сложность [math] f [/math] много больше любой из них, и равна константе [math]C[/math]. За весь алгоритм, она вызывается [math]4n[/math] раза. Таким образом, сложность последовательного алгоритма для одномерного случая равна [math]4nC[/math].

Для случая решения системы ОДУ, когда [math]f[/math] и [math]y[/math] являются векторами размерности [math]m[/math], функция core будет вызвана [math]m[/math], и итоговая сложность алгоритма составит [math]4nCm[/math].

1.7 Информационный граф

Информационная структура параллельного алгоритма Рунге-Кутты 4-го порядка представляет собой матрицу изображённую на рисунке ниже. Вдоль оси [math]j[/math] располагаются [math]j[/math]-ые параметры вектора [math]y[/math], соответственно, вдоль [math]i[/math] изображены этапы вычисления [math]y_i[/math] значения j-параметра вектора [math]y[/math]. Таким образом, можно сказать, что вдоль оси [math]j[/math], изображены параллельные процессы. А ось [math]i[/math] представляет собой шкалу времени. Чтобы соотнести её с реальным временем, надо просто посчитать сумму [math]\sum_{i=1}^{n} t_i[/math], где [math]t_i[/math] - время с момента начала вычисления [math]i[/math]-го всеми процессами до момента его завершения всеми процессами.

Можно заметить, что на каждом [math]i[/math]-ом шаге происходит обмен информацией между процессами. Это выражено в виде линии пересекающей все [math]j[/math]-ые процессы на любом [math]i[/math]-ом шаге.

1.8 Ресурс параллелизма алгоритма

Алгоритм Рунге-Кутты 4-го порядка является итерационным алгоритмом. В нём каждый шаг цикла зависит от предыдущего. В самом шаге, каждая последующая операция также зависит от результата предыдущей. Таким образом, упразднить цикл до независимых и соответственно параллельных операций не является возможным. Однако, так как на результат функции [math] \bold f [/math] влияет только [math]i[/math]-ый параметр вектора [math] y [/math], то для каждого [math] i [/math] можно производить независимое вычисление. В итоге, если вычислять результаты выполнения функции [math] \bold f [/math] одновременно для вля всех [math] i[/math]-ых параметров вектора [math] y [/math], получим туже алгоритмическую сложность [math]4nC[/math], что и в разделе 1.6 для одномерного случая.

1.9 Входные и выходные данные алгоритма

Входные данные:

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

Объём входных данных:

   m + 3

Выходные данные:

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

Объем выходных данных

   n * m

1.10 Свойства алгоритма

Алгоритм обладает следующими свойствами:

• Точность. Как написано в разделе 1.1, алгоритм имеет порядок точности [math]O(h^4)[/math].
• Устойчивость. Алгоритм устойчив на любом конечном отрезке.
• Гибкость. На любом этапе вычисления можно легко изменить шаг интегрирования.

2 ЧАСТЬ. Программная реализация алгоритма

2.7 Существующие реализации алгоритма

Алгоритм Рунге-Кутты 4-го порядка является одним из наиболее используемых для решения ОДУ и их систем, поэтому существует большое число его реализаций для разных языков программирования. Наиболее популярными являются следующие из них:

• библиотека RK4, которая существует в вариантах для C++, FORTRAN90, FORTRAN77, python и MATLAB;
• библиотека ode45 для MATLAB;
• модуль scipy для PYTHON.