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

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

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

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

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

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

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

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

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений размерности [math]\displaystyle{ n }[/math]

[math]\displaystyle{ \dot{x} = f(t, x),\ t_0 \leqslant t \leqslant t_1,\ x(t_0) = x_0. }[/math]

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

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

[math]\displaystyle{ t_i = t_0 + ih,\ i = \overline{1, n},\ h = \frac{t_1 - t_0}{m}, }[/math]

[math]\displaystyle{ x(t_i) = x_i. }[/math]

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

[math]\displaystyle{ x_{i+1} = x_i + \frac{h}{6} (K_1 + 2 K_2 + 2 K_3 + K_4). }[/math]

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

[math]\displaystyle{ K_1 = f (t_i, x_i), }[/math]

[math]\displaystyle{ K_2 = f (t_i + \frac{h}{2}, x_i + \frac{h}{2} K_1), }[/math]

[math]\displaystyle{ K_3 = f (t_i + \frac{h}{2}, x_i + \frac{h}{2} K_2), }[/math]

[math]\displaystyle{ K_4 = f (t_i + h, x_i + h K_3). }[/math]

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

В описанной выше вычислительной схеме наиболее трудоемкой является операция обращения к функции [math]\displaystyle{ f }[/math] при вычислении коэффициентов [math]\displaystyle{ K_i }[/math], эта операция является вычислительным ядром.

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

Макроструктура алгоритма представлена одним шагом итерационного процесса, описанного в пункте 1.2. Основной макрооперацией алгоритма является операция обращения к функции [math]\displaystyle{ f }[/math], и основное внимание будет уделено распараллеливанию этой операции.

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

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

начало;

   ввод начальных параметров (x[0], t0, t1, m);

   цикл по числу узлов сетки: i = 1..m-1

      цикл по числу стадий: j = 1..4

         цикл по числу компонент вектора x: k = 1..n

            вычисление коэффициентов K[j,k,i];

         конец цикла по k;

      конец цикла по j;

      цикл по числу компонент вектора x: k = 1..n

         вычисление x[k,i+1];

      конец цикла по k;

      цикл по числу компонент вектора x: k = 1..n

         сохранение x[k,i] = x[k,i+1];

      конец цикла по k;

      вывод x[i+1];

   конец цикла по n;

конец;

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

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

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

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

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

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

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

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

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

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

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

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

1. вектор начальных значений [math]\displaystyle{ x_0 }[/math] размерности [math]\displaystyle{ n }[/math];
2. границы временного интервала [math]\displaystyle{ t_0, t_1 }[/math];
3. частота дискретизации [math]\displaystyle{ m }[/math].

Общий размер входных данных [math]\displaystyle{ n + 3 }[/math].

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

1. вектор времени [math]\displaystyle{ t }[/math] размерности [math]\displaystyle{ m }[/math];
2. матрица значений [math]\displaystyle{ x }[/math] ([math]\displaystyle{ m }[/math] векторов длины [math]\displaystyle{ n }[/math]) размерности [math]\displaystyle{ m n }[/math].

Общий размер выходных данных [math]\displaystyle{ m + m n }[/math].

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

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

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

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

3. Метод не требуют вычисления производных от функции [math]\displaystyle{ f(t,x) }[/math], а только требует вычисления самой функции.

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

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

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

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

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

Литература

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

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

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

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