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

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

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

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

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

Методы Рунге-Кутты (распространено неправильное название Методы Рунге-Кутта) — важное семейство численных алгоритмов решения обыкновенных дифференциальных уравнений и их систем. Данные итеративные методы явного и неявного приближённого вычисления были разработаны около 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{ K_i }[/math], следовательно вычислительным ядром является их последовательное вычисление.

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

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

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

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

начало;

ввод начальных параметров ([math]\displaystyle{ x_0, t_0, t_1, m }[/math]);

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

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

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

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

конец цикла по [math]\displaystyle{ k }[/math];

конец цикла по [math]\displaystyle{ j }[/math];

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

вычисление [math]\displaystyle{ y_{j}^{i+1} }[/math];

конец цикла по [math]\displaystyle{ k }[/math];

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

сохранение [math]\displaystyle{ y_{j}^{i} = y_{j}^{i+1} }[/math];

конец цикла по [math]\displaystyle{ k }[/math];

вывод [math]\displaystyle{ y^{i+1} }[/math];

конец цикла по [math]\displaystyle{ n }[/math];

конец;

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

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

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

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

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