Участник:Anlesnichiy/Решение задачи Коши для системы ОДУ методом Рунге-Кутты 4 порядка
Решение задачи Коши для системы ОДУ методом Рунге-Кутты 4 порядка | |
Последовательный алгоритм | |
Последовательная сложность | [math] 4mn [/math] |
Объём входных данных | [math] m + 3 [/math] |
Объём выходных данных | [math] (m+1)n [/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | ? |
Ширина ярусно-параллельной формы | ? |
Основные авторы описания: А.А. Лесничий, Д.А. Алимов
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Метод Рунге-Кутты четвертого порядка — наиболее распространенный метод из семейства методов Метод Рунге-Кутты численных алгоритмов решения обыкновенных дифференциальных уравнений и их систем. Данные итеративные методы явного и неявного приближённого вычисления были разработаны около 1900 года немецкими математиками К. Рунге и М. В. Куттой.
Формально, методом Рунге-Кутты является модифицированный и исправленный метод Эйлера, они представляют собой схемы второго порядка точности. Существуют стандартные схемы третьего порядка, не получившие широкого распространения. Основная идея алгоритмов Рунге - Кутты состоит в замене правой части дифференциального уравнения, зависящей от искомой неизвестной функции, некоторым приближением.
1.2 Математическое описание алгоритма
1.2.1 Метод Рунге-Кутты 4-го порядка для задачи Коши для ДУ первого порядка
Рассмотрим задачу Коши
- [math] y' = f(x,y),\ a \leq x \leq b;\ y(a) = y^0 [/math]
Зададим равномерную сетку
- [math] x_i = a + ih,\ i = 1,\dots, n,\ h = \frac{b-a}{n} [/math]
Введём обозначения [math]y(x_i) = y_i[/math]. Получим вычислительную формулу:
- [math] \begin{cases} k_1 = hf(x_i,y_i)\\ k_2 = hf(x_i + h/2,y_i + k_1/2)\\ k_3 = hf(x_i + h/2,y_i + k_2/2)\\ k_4 = hf(x_i + h,y_i + k_3)\\ y_{i+1} = y_i + [ k_1 + 2k_2 + 2k_3 + k_4 ]/6 \\ \end{cases} [/math]
1.2.2 Метод Рунге-Кутты 4-го порядка для задачи Коши для системы ДУ первого порядка
Численное решение задачи Коши для систем ОДУ 1-го порядка методами Рунге-Кутты ищется по тем же формулам, что и для ОДУ первого порядка. Например, решение методом Рунге-Кутты 4 го порядка можно найти, если положить:
- [math] \begin{align} y_i \rightarrow \bar y_i\\ f(x_i,y_i) \rightarrow \bar f(x_i,\bar y_i)\\ k_l \rightarrow \bar k_l,\ l = 1, \dots, 4\\ \bar k_l = \begin{pmatrix} k^i_{l,1}\\ \vdots\\ k^i_{l,m}\\ \end{pmatrix} \end{align} [/math]
где m – размерность системы. В результате получим
- [math] \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} [/math]
1.3 Вычислительное ядро алгоритма
Вычислительным яром алгоритма можно считать вычисление значение вектора [math] \bar y_i [/math], то есть вычислительным ядром является цикл по i, описанный в пункте 1.2.2.
1.4 Макроструктура алгоритма
Макроструктурой алгоритма является одна итерация цикла по i, описанного в пункте 1.3.
1.5 Схема реализации последовательного алгоритма
Приведем возможную реализацию последовательного алгоритма на языке Matlab (правая часть ОДУ задана в виде функции func(x,y))
1 function [x, y] = RungeKutta4(y0, a, b, n)
2 h = (b-a)/n;
3 x = a : h : b
4
5 y(:, 1) = y0;
6 for i = 2 : n
7 k1 = h*func(x(i), y(:, i))
8 k2 = h*func(x(i) + h/2, y(:, i) + k1/2)
9 k3 = h*func(x(i) + h/2, y(:, i) + k2/2)
10 k4 = h*func(x(i) + h, y(:, i) + k3)
11 y(:, i+1) = y(:, i) + (k1 + 2*k2 + 2*k3 + k4)/6
12 end
13 end
1.6 Последовательная сложность алгоритма
Каждый шаг цикла алгоритма состоит из 4 обращений к функции [math] \bar f [/math], 11 умножений и 10 сложений. Так как обращение к функции [math] \bar f [/math] является наиболее сложной операцией, то сложность линейного выполнения алгоритма можно определить как [math] 4mn [/math].
1.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
Поскольку в описанной выше вычислительной схеме наиболее трудоемкой является операция расчета правых частей ОДУ при вычислении [math]k_i ( i = 1, \dots, 4) [/math], то основное внимание следует уделить распараллеливанию этой операции. Здесь будет применяться подход декомпозиции уравнений системы ОДУ на подсистемы. Поэтому для инициализации рассмотрим следующую схему декомпозиции данных по имеющимся процессорным элементам с локальной памятью: на каждый [math]\mu[/math]- ПЭ (процессорный элемент) ([math]\mu = 0, \dots, p-1 [/math]) распределяется [math] m/p [/math] дифференциальных уравнений и вектор [math] \bar y_0 [/math]. Далее расчеты производятся по следующей схеме:
- на каждом ПЭ одновременно вычисляются [math] n/p [/math] соответствующих компонент вектора [math] \bar k_1 [/math] по формуле [math] [ \bar k_1 ]_{\mu} = h[ \bar f(x_i, \bar y_i) ]_{\mu} [/math]
- для обеспечения второго расчетного этапа необходимо провести сборку вектора [math] \bar k_1 [/math] целиком на каждом ПЭ. Затем независимо выполняется вычисление компонент вектора [math] \bar k_2 [/math] по формуле [math] [ \bar k_2 ]_{\mu} = h[ \bar f(x_i + h/2,\bar y_i + \bar k_1/2)]_{\mu} [/math];
- проводится сборка вектора [math] \bar k_2 [/math] на каждом ПЭ, вычисляются компоненты вектора [math] \bar k_3:\ [ \bar k_3 ]_{\mu} = h [\bar f(x_i + h/2,\bar y_i + \bar k_2/2)]_{\mu} [/math];
- проводится сборка вектора [math] \bar k_3 [/math] на каждом ПЭ, вычисляются компоненты вектора [math] \bar k_4:\ [ \bar k_4 ]_{\mu} = h [\bar f(x_i + h,\bar y_i + \bar k_3)]_{\mu} [/math];
- рассчитываются с идеальным параллелизмом компоненты вектора [math] \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] \bar y_{i+1} [/math] на каждом ПЭ. Если необходимо продолжить вычислительный процесс, то полагается [math] i = i + 1 [/math] и осуществляется переход на п. 1
Заметим, что в данном алгоритме производится четыре операции вычисления вектора правых частей ОДУ, шестнадцать операций сложения векторов и умножения вектора на число и четыре операции глобальной сборки векторов.
1.9 Входные и выходные данные алгоритма
Входными данными алгоритма являются:
- вектор [math] y^0 [/math] размерности [math] m [/math];
- границы временного интервала [math] a [/math] и [math] b [/math];
- частота дискретизации [math] n [/math];
Всего размер входных данных: [math] m + 3 [/math]
Выходными данными являются
- [math] n [/math] векторов [math] \bar y_i [/math] размерности [math] m [/math];
- вектор [math] \bar x [/math] размерности [math] n [/math]
Всего размер выходных данных: [math] (m+1)n [/math]
1.10 Свойства алгоритма
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.2.1 Локальность реализации алгоритма
2.2.1.1 Структура обращений в память и качественная оценка локальности
2.2.1.2 Количественная оценка локальности
2.2.1.3 Анализ на основе теста Apex-Map
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
2.4.1 Масштабируемость алгоритма
Минимальная эффективность реализации 0.00032%;
Максимальная эффективность реализации 0.54%;
Построим оценки масштабируемости метода Рунге - Кутты 4-го порядка:
1. По числу процессоров -0.0000417;
2. По размеру задачи 0.00000526;
3. По двум направлениям -0.000000027;
2.4.2 Масштабируемость реализации алгоритма
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
MATLAB, PETsc.
3 Литература
[1] А. В. Старченко, В. Н. Берцун. Методы параллельных вычислений - Издательство Томского университета, 2013 - 173 с.
<references \>