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

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

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

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

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

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

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

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

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

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

$\frac{du}{dt} = f(t,u), u(0)=u_0$

Введем по переменному $t$ равномерную сетку с шагом $\tau \gt 0$, т.е. рассмотрим множество точек

$\begin{align} \omega_n = \left \{{t_n = n\tau, \ n = 0,\ 1,\ 2, \ \dots} \right \}. \end{align}$

Будем обозначать через $u(t)$ точное решение задачи Коши, а через $y_n = y(t_n)$ - приближенное решение.

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

$\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}$

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

1.2.2 Метод Рунге-Кутта четвертого порядка

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

$x_i = a + ih,\ i = 1,\dots, n,\ h = \frac{b-a}{n}$

Введём обозначения $y(x_i) = y_i$.

$\begin{cases} \bar k_1 = h\bar f(x_i,\bar y_i)\\ \bar k_2 = h\bar f(x_i + h/2,\bar y_i + \bar k_1/2)\\ \bar k_3 = h\bar f(x_i + h/2,\bar y_i + \bar k_2/2)\\ \bar k_4 = h\bar f(x_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}$

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

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

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

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

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

На каждой итерации алгоритма требуется 4 обращения к функции $\bar f$, 11 умножений и 10 сложений.

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

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

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

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

3 Литература

[1]

<references \>