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

Решение задачи Коши для системы ОДУ методом Рунге-Кутты
Последовательный алгоритм
Последовательная сложность	$4 m n$
Объём входных данных	$n + 3$
Объём выходных данных	$m(n + 1)$
Параллельный алгоритм
Высота ярусно-параллельной формы	$O(m)$
Ширина ярусно-параллельной формы	$O(n)$

Основные авторы описания: Филимонова Юлия

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

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

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

Формально, методом Рунге-Кутты является модифицированный и исправленный метод Эйлера, они представляют собой схемы второго порядка точности. Существуют стандартные схемы третьего порядка, не получившие широкого распространения. Наиболее часто используется и реализована в различных математических пакетах (Maple, MathCAD, Maxima) стандартная схема четвёртого порядка. Классический метод Рунге-Кутты четвёртого порядка столь широко распространён, что его часто называют просто методом Рунге-Кутты, опуская порядок.

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

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений размерности $n$

$\dot{x} = f(t, x),\ t_0 \leqslant t \leqslant t_1,\ x(t_0) = x_0.$

Здесь $x(t), x_0 \in \mathbb{R}^n,\ t, t_0, t_1 \in \mathbb{R}, f: \mathbb{R} \rightarrow \mathbb{R}^n$.

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

$t_i = t_0 + ih,\ i = \overline{1, n},\ h = \frac{t_1 - t_0}{m},$

$x(t_i) = x_i.$

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

$x_{i+1} = x_i + \frac{h}{6} (K_1 + 2 K_2 + 2 K_3 + K_4).$

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

$K_1 = f (t_i, x_i),$

$K_2 = f (t_i + \frac{h}{2}, x_i + \frac{h}{2} K_1),$

$K_3 = f (t_i + \frac{h}{2}, x_i + \frac{h}{2} K_2),$

$K_4 = f (t_i + h, x_i + h K_3).$

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

В описанной выше вычислительной схеме наиболее трудоемкой является операция обращения к функции $f$ при вычислении коэффициентов $K_i$, эта операция является вычислительным ядром.

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

Макроструктура алгоритма представлена одним шагом итерационного процесса, описанного в пункте 1.2. Основной макрооперацией алгоритма является операция обращения к функции $f$, и основное внимание будет уделено распараллеливанию этой операции.

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

Приведем здесь псевдокод

начало;

   ввод начальных параметров (<math>x_0, t_0, t_1, m</math>);

   цикл по числу узлов сетки: <math>i = \overline{1,m-1}</math> 

      цикл по числу стадий: <math>j = \overline{1,4}</math>

         цикл по числу компонент вектора <math>x:\ k = \overline{1,n}</math>

            вычисление коэффициентов <math>K_{j, k}^{i}</math>;

         конец цикла по <math>k</math>;

      конец цикла по <math>j</math>;

      цикл по числу компонент вектора <math>x:\ k = \overline{1,n}</math>

         вычисление <math>y_{j}^{i+1}</math>;

      конец цикла по <math>k</math>;

      цикл по числу компонент вектора <math>x:\ k = \overline{1,n}</math>

         сохранение <math>y_{j}^{i} = y_{j}^{i+1}</math>;

      конец цикла по <math>k</math>;

      вывод <math>y^{i+1}</math>;

   конец цикла по <math>n</math>;

конец;

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

В данном алгоритме производится четыре операции обращения к функции $f$ (порядка $n m$ арифметических операций) и шестнадцать операций сложения векторов и умножения вектора на число (порядка $n$ арифметических операций). Поскольку операция обращения к функции является более сложной, то сложность последовательного алгоритма можно считать равной $4 n m$.

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

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

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

Последовательная схема алгоритма дает основание для организации параллельных вычислений с помощью декомпозиции. Пусть доступно $p$ процессоров (для определенности будем считать, что размерность системы $n$ кратна количество процессоров $p$: $n = qp$). На каждый из процессоров распределяется и вычисляется последовательно $q$ компонент векторов коэффициентов. Далее расчеты проводятся следующим образом:

1. На каждом процессоре вычисляется q соответствующих компонент вектора $[K_1]$ по формуле $[K_1] = [f(t_i, x_i)]$. Проводится сборка вектора $[K_1]$ целиком на каждом процессоре.

2-4. Аналогичным образом проводятся вычисления и сборка векторов $[K_2], [K_3], [K_4]$.

5. На каждом процессоре вычисляется q соответствующих компонент вектора $[x_{i+1}]$ по формуле $[x_{i+1}] = [x_i] + \frac{h}{6} ([K_1] + 2 [K_2] + 2 [K_3] + [K_4])$. Проводится сборка вектора $[x_{i+1}]$ целиком на каждом процессоре. Если вычисления не закончены, то сохраняется $[x_{i}] = [x_{i+1}]$, осуществляется переход на 1.

Алгоритм производит четыре обращения к функции $f$, шестнадцать операций сложения векторов и умножения вектора на число и четыре операции глобальной сборки векторов. Сложность алгоритма по высоте составляет $m$, а сложность по ширине равна $q = \frac{n}{p}$.

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

На вход алгоритма подаются следующие данные:

1. вектор начальных значений $x_0$ размерности $n$;
2. границы временного интервала $t_0, t_1$;
3. частота дискретизации $m$.

Общий размер входных данных $n + 3$.

На выходе получаются следующие данные:

1. вектор времени $t$ размерности $m$;
2. матрица значений $x$ ($m$ векторов длины $n$) размерности $m n$.

Общий размер выходных данных $m + m n$.

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

Метод Рунге-Кутты 4-го порядка обладает следующими свойствами:

1. Эти метод (как и все семейство методов Рунге-Кутты) является одношаговым: чтобы найти $x_{n+1}$, нужна информация об одной предыдущей точке $(t_i, x_i)$. (Это позволяет в любой момент изменить шаг интегрирования.)

2. Метод согласуются с рядом Тейлора вплоть до членов порядка $4h$. (Используя большее количество вспомогательных точек, можно увеличить точность метода.)

3. Метод не требуют вычисления производных от функции $f(t,x)$, а только требует вычисления самой функции.

Точность и устойчивость метода достаточна для широкого круга задач, метод прост в реализации - эти достоинства определили популярность метода среди большого количества исследователей.

Число операций алгоритма равно $m (n - 1)$, общее число входных и выходных данных равно $mn + m + n + 3$. Следовательно, вычислительная мощность алгоритма $\frac{mn-m}{mn + m + n + 3} \longrightarrow 1$ при увеличении размерности задачи.

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

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

Стандартная схема четвертого порядка реализована в различных математических пакетах: Maple (rkf45), MathCAD (rkfixed), Matlab (ode45).

3 Литература

[1] А. В. Старченко, Высокопроизводительные вычисления на кластерах, Издательство Томского университета, 2013, 127-130 с

[2] Л. П. Фельдман, И. А. Назарова, Параллельные алгоритмы численного решения задачи Коши для систем обыкновенных дифференциальных уравнений [1]

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

[4] Ващенко Г.В. Явные методы типа Рунге-Кутты и их параллельные аналоги // Численные методы, программные системы и комплексы программ. [2]