Участник:Alexboboshko/Решение начальной задачи Коши для системы ОДУ методом Рунге-Кутта 4-го порядка

РК
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ 4mn }[/math]
Объём входных данных	[math]\displaystyle{ m + 3 }[/math]
Объём выходных данных	[math]\displaystyle{ (m+1)n }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	?
Ширина ярусно-параллельной формы	?

Основные авторы описания:

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

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

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

Формально, методом Рунге — Кутты является модифицированный и исправленный метод Эйлера, они представляют собой схемы второго порядка точности. Существуют стандартные схемы третьего порядка, не получившие широкого распространения.

Запишем точную формулу для получения [math]\displaystyle{ y_{i+1} }[/math] :

[math]\displaystyle{ y_{i+1} = y_i + \int\limits_{t_i}^{t_{i+1}}f(t,y)dt }[/math]

Основная идея алгоритмов Рунге-Кутты состоит в замене функции [math]\displaystyle{ f(t,y) }[/math], которая зависит от неизвестной функции [math]\displaystyle{ y(t) }[/math], некоторым приближением. Чем точнее будет приближенное значение подынтегральной функции, тем точнее будет посчитан интеграл и точнее будет определено значение [math]\displaystyle{ y_{i+1} }[/math].

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

Общая формулировка методов

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

[math]\displaystyle{ \frac{du}{dt} = f(t,u), u(0)=u_0 }[/math]

Введем по переменному [math]\displaystyle{ t }[/math] равномерную сетку с шагом [math]\displaystyle{ \tau \gt 0 }[/math], т.е. рассмотрим множество точек

[math]\displaystyle{ \begin{align} \omega_n = \left \{{t_n = n\tau, \ n = 0,\ 1,\ 2, \ \dots} \right \}. \end{align} }[/math]

Будем обозначать через [math]\displaystyle{ u(t) }[/math] точное решение задачи Коши, а через [math]\displaystyle{ y_n = y(t_n) }[/math] - приближенное решение.

Явный m-этапный метод Рунге-Кутта состоит в следующем. Пусть решение [math]\displaystyle{ y_n = y(t_n) }[/math] уже известно. Задаются числовые коэффициенты [math]\displaystyle{ a_i, \ b_{ij},\ i = 2, 3, \dots, m,\ j = 1, 2, \dots, m-1,\ \sigma_i,\ i = 1, 2, \dots, m, }[/math] и последовательно вычисляются функции

[math]\displaystyle{ \begin{array}{l} k_1 = f(t_n, y_n), \ k_2 = f(t_n+a_2\tau,\ y_n+b_{21}\tau k_1) \\ k_3 = f(t_n+a_3 \tau, \ y_n+b_{31}\tau k_1+b_{32}\tau k_2), \ \dots, \\ k_m = f(t_n+a_m \tau, \ y_n+b_{m1}\tau k_1+b_{m2}\tau k_2+\dots+b_{m,m-1}\tau k_{m-1}) \end{array}{l} }[/math]

Затем из формулы [math]\displaystyle{ \frac{y_{n+1}-y_n}{\tau} = \sum_{i = 1}^{m} \sigma_i k_i }[/math] находится новое значение [math]\displaystyle{ y_{n+1} = y(t_{n+1}) }[/math]

Метод Рунге-Кутта четвертого порядка для системы ДУ

Запишем задачу Коши в векторном виде

[math]\displaystyle{ \frac{d \bar y}{dt} = \bar f(t, \bar y); \ \bar y(t_0) = \bar y_0 = \begin{pmatrix} y_{10}\\ \vdots\\ y_{m0}\\ \end{pmatrix} }[/math]

[math]\displaystyle{ a \leq t \leq b }[/math]

Зададим равномерную сетку

[math]\displaystyle{ t_i = a + ih,\ i = 1,\dots, n,\ h = \frac{b-a}{n} }[/math]

Введём обозначения [math]\displaystyle{ y(x_i) = y_i }[/math].

[math]\displaystyle{ \begin{cases} \bar k_1 = h\bar f(t_i,\bar y_i)\\ \bar k_2 = h\bar f(t_i + h/2,\bar y_i + \bar k_1/2)\\ \bar k_3 = h\bar f(t_i + h/2,\bar y_i + \bar k_2/2)\\ \bar k_4 = h\bar f(t_i + h,\bar y_i + \bar k_3)\\ \bar y_{i+1} = \bar y_i + [ \bar k_1 + 2\bar k_2 + 2\bar k_3 + \bar k_4 ]/6 \\ \end{cases} }[/math]

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

В описанной выше вычислительной схеме наиболее трудоемкой является операция расчета правых частей ОДУ при вычислении [math]\displaystyle{ k_i ( i = 1, \dots, 4) , }[/math] то есть основное внимание следует уделить распараллеливанию этой операции.

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

Пример программы на языке C

   #include <stdio.h>
   void main(){
      float t=0, dt = 0.05;
      float vx = 1, vy = 2;
      float x=0, y=0;
      vector U(x, y, vx, vy);
      vector k1, k2, k3, k4;
      while(U.y >= 0){
          k1 = F(U, t)*dt;
          k2 = F(U + 0.5*k1, t+0.5*dt)*dt;
          k3 = F(U + 0.5*k2, t+0.5*dt)*dt;
          k4 = F(U + k3, t+dt)*dt;
          U = U + 1.0/6.0 * (k1 + 2*k2 + 2*k3 + k4);
          t += dt;       
      }
   }

http://www.alexeypetrov.narod.ru/C/runge_about.html

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

Для каждого [math]\displaystyle{ i }[/math] последовательно вычисляются функции [math]\displaystyle{ k_i \ ( i = 1, \dots, 4) }[/math]. После этого вычисляется значение [math]\displaystyle{ y_{i+1} }[/math].

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

На каждой итерации алгоритма требуется 4 обращения к функции [math]\displaystyle{ \bar f }[/math], 11 умножений и 10 сложений.

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

Задача распараллеливания в такой постановке сводится к отысканию максимального независимого множества вершин орграфа, причем вершинам графа сопоставляются выполняемые операции и вершины соединяются дугами тогда и только тогда, когда результат выполнения одной операции влияет на результат смежной. На рисунке выше приведен граф влияния для вычисления вектора значений [math]\displaystyle{ i }[/math]-го шагового коэффициента в случае, если [math]\displaystyle{ m=p }[/math].

В общем случае для каждого шагового коэффициента [math]\displaystyle{ k_i , i = \overline{1, s} }[/math] , соответственно, [math]\displaystyle{ l }[/math] компонент вектора [math]\displaystyle{ g_i , i = \overline{1, s} }[/math] могут быть вычислены параллельно

1) [math]\displaystyle{ l=m }[/math] компонент при [math]\displaystyle{ m=p }[/math];

2) [math]\displaystyle{ l=[m/p] }[/math] компонент при [math]\displaystyle{ m\gt p }[/math].

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

Поскольку в описанной выше вычислительной схеме наиболее трудоемкой является операция расчета правых частей ОДУ при вычислении [math]\displaystyle{ k_i ( i = 1, \dots, 4) }[/math], то основное внимание следует уделить распараллеливанию этой операции. Здесь будет применяться подход декомпозиции уравнений системы ОДУ на подсистемы. Поэтому для инициализации рассмотрим следующую схему декомпозиции данных по имеющимся процессорным элементам с локальной памятью: на каждый [math]\displaystyle{ \mu }[/math]- ПЭ (процессорный элемент) ([math]\displaystyle{ \mu = 0, \dots, p-1 }[/math]) распределяется [math]\displaystyle{ m/p }[/math] дифференциальных уравнений и вектор [math]\displaystyle{ \bar y_0 }[/math]. Далее расчеты производятся по следующей схеме:

1. на каждом ПЭ одновременно вычисляются [math]\displaystyle{ m/p }[/math] соответствующих компонент вектора [math]\displaystyle{ \bar k_1 }[/math] по формуле [math]\displaystyle{ [ \bar k_1 ]_{\mu} = h[ \bar f(x_i, \bar y_i) ]_{\mu} }[/math]
2. для обеспечения второго расчетного этапа необходимо провести сборку вектора [math]\displaystyle{ \bar k_1 }[/math] целиком на каждом ПЭ. Затем независимо выполняется вычисление компонент вектора [math]\displaystyle{ \bar k_2 }[/math] по формуле [math]\displaystyle{ [ \bar k_2 ]_{\mu} = h[ \bar f(x_i + h/2,\bar y_i + 1/2 \bar k_1)]_{\mu} }[/math];
3. проводится сборка вектора [math]\displaystyle{ \bar k_2 }[/math] на каждом ПЭ, вычисляются компоненты вектора [math]\displaystyle{ \bar k_3:\ [ \bar k_3 ]_{\mu} = h [\bar f(x_i + h/2,\bar y_i + 1/2 \bar k_2)]_{\mu} }[/math];
4. проводится сборка вектора [math]\displaystyle{ \bar k_3 }[/math] на каждом ПЭ, вычисляются компоненты вектора [math]\displaystyle{ \bar k_4:\ [ \bar k_4 ]_{\mu} = h [\bar f(x_i + h,\bar y_i + \bar k_3)]_{\mu} }[/math];
5. рассчитываются с идеальным параллелизмом компоненты вектора [math]\displaystyle{ \bar y_{i+1}:\ [\bar y_{i+1}]_{\mu} = [\bar y_{i}]_{\mu} + ([ \bar k_1 ]_{\mu} + 2[ \bar k_2 ]_{\mu} + 2[ \bar k_3 ]_{\mu} + [ \bar k_4 ]_{\mu})/6\ }[/math] и производится сборка вектора [math]\displaystyle{ \bar y_{i+1} }[/math] на каждом ПЭ. Если необходимо продолжить вычислительный процесс, то полагается [math]\displaystyle{ i = i + 1 }[/math] и осуществляется переход на п. 1

Заметим, что в данном алгоритме производится четыре операции вычисления вектора правых частей ОДУ, шестнадцать операций сложения векторов и умножения вектора на число и четыре операции глобальной сборки векторов.

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

Вход:

1. [math]\displaystyle{ \bar y_{0} \in \mathbb{R}^m }[/math] - начальное значение
2. [math]\displaystyle{ a,\ b }[/math] - границы отрезка
3. [math]\displaystyle{ n }[/math] - число итераций

Выход:

1. [math]\displaystyle{ \bar y_{1},\dots, \bar y_n \in \mathbb{R}^m }[/math]
2. [math]\displaystyle{ \bar x \in \mathbb{R}^n }[/math]

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

Методы Рунге-Кутты имеют несколько достоинств, определивших их популярность среди значительного числа исследователей. Эти методы легко программируются. Они обладают достаточными для широкого круга задач свойствами точности и устойчивости. Эти методы (как и все одношаговые методы) являются самостартующими и позволяют на любом этапе вычислений легко изменять шаг интегрирования.

Если метод Рунге-Кутта аппроксимирует исходное уравнение, то он сходится при диаметре разбиения [math]\displaystyle{ h }[/math] → 0, причем порядок точности совпадает с порядком аппроксимации. В случае методов 4-го порядка это значит, что ошибка на одном шаге имеет порядок [math]\displaystyle{ O(h^5) }[/math], а суммарная ошибка на конечном интервале интегрирования имеет порядок [math]\displaystyle{ O(h^4) }[/math] .

Увеличивая число [math]\displaystyle{ m }[/math] вспомогательных точек, можно построить методы Рунге-Кутты любого порядка точности [math]\displaystyle{ p }[/math]. Однако уже при [math]\displaystyle{ p\gt 5 }[/math] эти методы используются довольно редко. Это объясняется как чрезмерной громоздкостью получающихся вычислительных формул, так и том, что преимущества методов высокого порядка точности [math]\displaystyle{ p }[/math] над методами, в которых [math]\displaystyle{ p=4 }[/math] и [math]\displaystyle{ p=5 }[/math], проявляются либо в тех задачах, где нужна очень высокая точность и используются ЭВМ высокой разрядности, либо в тех задачах, где решение очень гладкое. Кроме того, методы Рунге-Кутты высокого порядка точности часто оказываются менее эффективными по сравнению с методами Адамса того же порядка точности.

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

Программная реализация алгоритма

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

Стандартная схема 4-го порядка реализована в различных математических пакетах (Maple, MathCAD, Maxima), в библиотеке для языка программирования Python - SciPy.

Литература

[1] Тыртышников Е. Е. Методы численного анализа — М., Академия, 2007.

[2] Бахвалов Н. С., Жидков Н. П., Кобельков. Г. М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008.

[3] Самарский А. А., Гулин А. В. Численные методы. М.: Наука, 1989. <references \>