Участник:Anlesnichiy/Решение задачи Коши для системы ОДУ методом Рунге-Кутты 4 порядка: различия между версиями
(Код из пункта 1.3 перемещен в пунт 1.5 и дописан по функции) |
|||
Строка 8: | Строка 8: | ||
}} | }} | ||
− | Основные авторы описания: [[Участник:Anlesnichiy|А.А. Лесничий]] | + | Основные авторы описания: [[Участник:Anlesnichiy|А.А. Лесничий]], [[Участник:Алимов Дамир Алиевич|Д.А. Алимов]] |
== Свойства и структура алгоритма == | == Свойства и структура алгоритма == | ||
Строка 19: | Строка 19: | ||
==== Метод Рунге-Кутты для задачи Коши для ДУ первого порядка ==== | ==== Метод Рунге-Кутты для задачи Коши для ДУ первого порядка ==== | ||
− | |||
Рассмотрим задачу Коши | Рассмотрим задачу Коши | ||
:<math> | :<math> | ||
Строка 69: | Строка 68: | ||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
+ | Вычислительным яром алгоритма можно считать вычисление значение вектора <math> \bar y_i </math>, то есть вычислительным ядром является цикл по i, описанный в пункте [[#Метод Рунге-Кутты для задачи Коши для системы ДУ первого порядка|1.2.2]]. | ||
− | + | === Макроструктура алгоритма === | |
+ | Макроструктурой алгоритма является одна итерация цикла по i, описанного в пункте [[#Вычислительное ядро алгоритма|1.3]]. | ||
− | + | === Схема реализации последовательного алгоритма === | |
+ | Приведем возможную реализацию последовательного алгоритма на языке Matlab (правая часть ОДУ задана в виде функции func(x,y)) | ||
<source lang="matlab" line="1"> | <source lang="matlab" line="1"> | ||
− | for i = | + | function [x, y] = RungeKutta4(y0, a, b, n) |
− | + | h = (b-a)/n; | |
− | + | x = a : h : b | |
− | + | ||
− | + | y(:, 1) = y0; | |
− | + | for i = 2 : n | |
+ | k1 = h*func(x(i), y(:, i)) | ||
+ | k2 = h*func(x(i) + h/2, y(:, i) + k1/2) | ||
+ | k3 = h*func(x(i) + h/2, y(:, i) + k2/2) | ||
+ | k4 = h*func(x(i) + h, y(:, i) + k3) | ||
+ | y(:, i+1) = y(:, i) + (k1 + 2*k2 + 2*k3 + k4)/6 | ||
+ | end | ||
end | end | ||
</source> | </source> | ||
− | |||
− | |||
− | |||
− | |||
=== Последовательная сложность алгоритма === | === Последовательная сложность алгоритма === | ||
− | |||
Каждый шаг цикла алгоритма состоит из 4 обращений к функции <math> \bar f </math>, 11 умножений и 10 сложений. Так как обращение к функции <math> \bar f </math> является наиболее сложной операцией, то сложность линейного выполнения алгоритма можно определить как <math> 4mn </math>. | Каждый шаг цикла алгоритма состоит из 4 обращений к функции <math> \bar f </math>, 11 умножений и 10 сложений. Так как обращение к функции <math> \bar f </math> является наиболее сложной операцией, то сложность линейного выполнения алгоритма можно определить как <math> 4mn </math>. | ||
Строка 107: | Строка 110: | ||
=== Входные и выходные данные алгоритма === | === Входные и выходные данные алгоритма === | ||
− | |||
Входными данными алгоритма являются: | Входными данными алгоритма являются: | ||
# вектор <math> y^0 </math> размерности <math> m </math>; | # вектор <math> y^0 </math> размерности <math> m </math>; |
Версия 22:28, 11 октября 2016
Решение задачи Коши для системы ОДУ методом Рунге-Кутты 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 Метод Рунге-Кутты для задачи Коши для ДУ первого порядка
Рассмотрим задачу Коши
- [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 Метод Рунге-Кутты для задачи Коши для системы ДУ первого порядка
Численное решение задачи Коши для систем ОДУ 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 Масштабируемость алгоритма
2.4.2 Масштабируемость реализации алгоритма
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
3 Литература
[1] А. В. Старченко, В. Н. Берцун. Методы параллельных вычислений - Издательство Томского университета, 2013 - 173 с.
<references \>