Уровень алгоритма

Участник:Борис/Алгоритм Рунге-Кутты: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показано 30 промежуточных версий 3 участников)
Строка 1: Строка 1:
 +
 
{{algorithm
 
{{algorithm
 
| name              = Метод Рунге-Кутты 4 порядка
 
| name              = Метод Рунге-Кутты 4 порядка
| serial_complexity = <math>O()</math>
+
| serial_complexity = <math>4nCm</math>
| pf_height        = <math>O()</math>
+
| pf_height        = <math>O(n)</math>
| pf_width          = <math>O()</math>
+
| pf_width          = <math>O(m)</math>
 
| input_data        = <math>m + 3</math>
 
| input_data        = <math>m + 3</math>
 
| output_data      = <math>m * n</math>
 
| output_data      = <math>m * n</math>
 
}}
 
}}
  
Основные авторы описания: [[Участник:Jest|Е.П.Степанов]]([[#Общее описание алгоритма|раздел 1.1]], [[#Математическое описание алгоритма|раздел 1.2]], [[#Вычислительное ядро алгоритма|раздел 1.3]], [[#Макроструктура алгоритма|раздел 1.4]], [[#Входные и выходные данные алгоритма|раздел 1.9]]), [[Участник:Борис|Б.Савков]] ( [[#Схема реализации последовательного алгоритма|раздел 1.5]], [[#Последовательная сложность алгоритма|раздел 1.6]], [[#Информационный граф|раздел 1.7]], [[#Ресурс параллелизма алгоритма|раздел 1.8]], [[#Существующие реализации алгоритма|раздел 2.7]])
+
Основные авторы описания: [[Участник:Jest|Е.П.Степанов]](разделы [[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Вычислительное ядро алгоритма|1.3]], [[#Макроструктура алгоритма|1.4]], [[#Входные и выходные данные алгоритма|1.9]]), [[Участник:Борис|Б.Савков]] (разделы [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Информационный граф|1.7]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Свойства алгоритма|1.10]], [[#Существующие реализации алгоритма|2.7]])
  
= ЧАСТЬ. Свойства и структура алгоритмов =
+
= Свойства и структура алгоритмов =
  
 
== Общее описание алгоритма ==
 
== Общее описание алгоритма ==
Строка 31: Строка 32:
 
:<math> y(x + h) = y(x) + \int_0^hy'(x + t)dt</math>
 
:<math> y(x + h) = y(x) + \int_0^hy'(x + t)dt</math>
  
Аппроксимируя интеграл в правой части, можно получить методы Рунге-Кутты разного порядка точности. В общем виде формулы методов Рунге-Кутты можно записать следующим образом:
+
Аппроксимируя интеграл в правой части, можно получить методы Рунге-Кутты разного порядка точности <ref>Бахвалов Н. С.,  Жидков Н. П., Кобельков. Г. М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.</ref>. В общем виде формулы методов Рунге-Кутты можно записать следующим образом:
  
 
:<math>\textbf{y}_{n+1} = \textbf{y}_n + h\sum_{i=1}^s b_i \textbf{k}_i</math>
 
:<math>\textbf{y}_{n+1} = \textbf{y}_n + h\sum_{i=1}^s b_i \textbf{k}_i</math>
Строка 44: Строка 45:
 
</math>
 
</math>
  
Увеличивая число s вспомогательных точек, можно построить методы Рунге-Кутты любого порядка точности p. Однако уже при p > 5 эти методы используются довольно редко. Это объясняется как чрезмерной громоздкостью получающихся вычислительных формул, так и тем, что преимущества методов высокого порядка точности  над методами, в которых p = 4 или p = 5, проявляются либо в тех задачах, где нужна очень высокая точность и используются ЭВМ высокой разрядности, либо в тех задачах, где решение очень гладкое. Кроме того, методы Рунге-Кутты высокого порядка точности  часто оказываются менее эффективными по сравнению с методами Адамса того же порядка точности. На этой странице представлен наиболее распространенный метод Рунге-Кутты 4 порядка с ошибкой на конечном интервале интегрирования <math>O(h^4)</math>.
+
Увеличивая число s вспомогательных точек, можно построить методы Рунге-Кутты любого порядка точности p. Однако уже при p > 5 эти методы используются довольно редко. Это объясняется как чрезмерной громоздкостью получающихся вычислительных формул, так и тем, что преимущества методов высокого порядка точности  над методами, в которых p = 4 или p = 5, проявляются либо в тех задачах, где нужна очень высокая точность и используются ЭВМ высокой разрядности, либо в тех задачах, где решение очень гладкое. Кроме того, методы Рунге-Кутты высокого порядка точности  часто оказываются менее эффективными по сравнению с методами Адамса того же порядка точности <ref>Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров: Учеб. пособие. — М.: Высш. шк., 1994. — 544 с. </ref>. На этой странице представлен наиболее распространенный метод Рунге-Кутты 4 порядка с ошибкой на конечном интервале интегрирования <math>O(h^4)</math>.
  
 
Существуют как явные, так и неявные методы Рунге-Кутты. Последние имеют ряд преимуществ над явными методами, однако это достигается за счет существенного усложнения вычислительного алгоритма, так как на каждом шаге приходится считать систему нелинейных уравнений. На этой странице представлен явный метод Рунге-Кутты.
 
Существуют как явные, так и неявные методы Рунге-Кутты. Последние имеют ряд преимуществ над явными методами, однако это достигается за счет существенного усложнения вычислительного алгоритма, так как на каждом шаге приходится считать систему нелинейных уравнений. На этой странице представлен явный метод Рунге-Кутты.
Строка 55: Строка 56:
  
 
Исходные данные:  
 
Исходные данные:  
    <math> a, b </math> - интервал интегрирования,
+
:<math> a, b </math> - интервал интегрирования,
    <math> \bold y_0 </math> - вектор значений при <math> x = a </math>,  
+
:<math> \bold y_0 </math> - вектор значений при <math> x = a </math>,  
    n - количество точек в разбиении
+
:n - количество точек в разбиении
    <math>\bold f </math> - функция вычисления производной
+
:<math>\bold f </math> - функция вычисления производной
  
 
Вычисляемые данные:
 
Вычисляемые данные:
    <math> \bold y_i </math> - вектора значений в узлах сетки, т.е. при <math> x = a + i * \frac{b-a}{n}, i \in [1, n] </math>
+
:<math> \bold y_i </math> - вектора значений в узлах сетки, т.е. при <math> x = a + i * \frac{b-a}{n}, i \in [1, n] </math>
  
 
Формулы метода:   
 
Формулы метода:   
Строка 83: Строка 84:
  
 
Макроструктура алгоритма может быть представлена в следующем виде:
 
Макроструктура алгоритма может быть представлена в следующем виде:
:Вычисление вспомогательных коэффициентов <math>\textbf{k}_j, j \in [1,4]</math>;
+
 
:Вычисление приближенного значения <math> \bold y_i </math> в следующей точке;
+
1. Вычисление вспомогательных коэффициентов <math>\textbf{k}_j, j \in [1,4]</math>;
:Переход на следующую итерацию или остановка.
+
 
 +
2. Вычисление приближенного значения <math> \bold y_i </math> в следующей точке;
 +
 
 +
3. Переход на следующую итерацию или остановка.
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
Строка 99: Строка 103:
 
     std::vector<double>::size_type sz = y0.size();
 
     std::vector<double>::size_type sz = y0.size();
 
     std::vector<double> result;
 
     std::vector<double> result;
     for (unsigned i=0; i<sz; i++)  
+
     for (unsigned j=0; j<sz; j++)  
 
     {
 
     {
         rezult.push_back(core(a,b,y[i],n,f[i]));
+
         rezult.push_back(core(a,b,y[j],n,f[j]));
 
     }
 
     }
 
     return result;
 
     return result;
Строка 129: Строка 133:
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==
  
Рассмотрим сначала одномерный случай. Наравне с обычными операциями, такими как сложение, умножение и тому подобными, в представленном выше алгоритме присутствует вызов функции <math>\bold f </math>. Так как эта функция может содержать в себе любое количество простых математических операций, будем считать, что её сложность <math>\bold f </math> много больше любой из них, и равна константе <math>C</math>. За весь алгоритм, она вызывается <math>4n</math> раза. Таким образом, сложность последовательного алгоритма для одномерного случая равна <math>4nC</math>.  
+
Рассмотрим сначала одномерный случай. Наравне с обычными операциями, такими как сложение, умножение и тому подобными, в представленном выше алгоритме присутствует вызов функции <math> f </math>. Так как эта функция может содержать в себе любое количество простых математических операций, будем считать, что её сложность <math> f </math> много больше любой из них, и равна константе <math>C</math>. За весь алгоритм, она вызывается <math>4n</math> раза. Таким образом, сложность последовательного алгоритма для одномерного случая равна <math>4nC</math>.  
  
Для многомерного случая, когда <math>f</math> и <math>y</math> являются векторами размерности <math>m</math>, функция core будет вызвана <math>m</math>, и итоговая сложность алгоритма составит <math>4nCm</math>.
+
Для случая решения системы ОДУ, когда <math>f</math> и <math>y</math> являются векторами размерности <math>m</math>, функция core будет вызвана <math>m</math>, и итоговая сложность алгоритма составит <math>4nCm</math>.
  
 
== Информационный граф ==
 
== Информационный граф ==
 +
 +
Информационная структура параллельного алгоритма Рунге-Кутты 4-го порядка представляет собой матрицу изображённую на рисунке ниже. Вдоль оси <math>j</math> располагаются <math>j</math>-ые параметры вектора <math>y</math>, соответственно, вдоль <math>i</math> изображены этапы вычисления <math>y_i</math> значения j-параметра вектора <math>y</math>. Таким образом, можно сказать, что вдоль оси <math>j</math>, изображены параллельные процессы. А ось <math>i</math> представляет собой шкалу времени. Чтобы соотнести её с реальным временем, надо просто посчитать сумму <math>\sum_{i=1}^{n} t_i</math>, где <math>t_i</math> - время с момента начала вычисления <math>i</math>-го всеми процессами до момента его завершения всеми процессами.
 +
 +
[[Файл:System_graph.png]]
 +
 +
Можно заметить, что на каждом <math>i</math>-ом шаге происходит обмен информацией между процессами. Это выражено в виде линии пересекающей все <math>j</math>-ые процессы на любом <math>i</math>-ом шаге.
  
 
== Ресурс параллелизма алгоритма ==
 
== Ресурс параллелизма алгоритма ==
  
Алгоритм Рунге-Кутты 4-го порядка является итерационным алгоритмом. В нём каждый шаг цикла зависит от предыдущего. В самом шаге, каждая последующая операция также зависит от результата предыдущей. Таким образом, упразднить цикл до независимых и соответственно параллельных операций не является возможным. Однако, так как на результат функции <math> \bold f </math> влияет только <math>i</math>-ый параметр векторов <math> f </math> и <math> \bold y </math>, то для каждого <math> i </math> можно производить независимое вычисление. В итоге, если вычислять результаты выполнения функции <math> \bold f </math> одновременно для вля всех <math>\bold i</math>-ых параметров векторов <math> f </math> и <math> \bold y </math>, получим туже алгоритмическую сложность <math>4nC</math>, что и в [[#Последовательная сложность алгоритма|разделе 1.6]] для одномерного случая.
+
Алгоритм Рунге-Кутты 4-го порядка является итерационным алгоритмом. В нём каждый шаг цикла зависит от предыдущего. В самом шаге, каждая последующая операция также зависит от результата предыдущей. Таким образом, упразднить цикл до независимых и соответственно параллельных операций не является возможным. Однако, так как на результат функции <math> f </math> влияет только <math>i</math>-ый параметр вектора <math> y </math>, то для каждого <math> i </math> можно производить независимое вычисление. В итоге, если вычислять результаты выполнения функции <math> f </math> одновременно для вля всех <math> i</math>-ых параметров вектора <math> y </math>, получим ту же алгоритмическую сложность <math>4nC</math>, что и в [[#Последовательная сложность алгоритма|разделе 1.6]] для одномерного случая.
  
 
== Входные и выходные данные алгоритма ==
 
== Входные и выходные данные алгоритма ==
  
 
'''Входные данные''':  
 
'''Входные данные''':  
    <math> a, b </math> - интервал интегрирования,     
+
:<math> a, b </math> - интервал интегрирования,     
    n - количество точек в разбиении,
+
:n - количество точек в разбиении,
    <math> \bold y_0 </math> - вектор значений размерности m при <math> x = a </math>,
+
:<math> \bold y_0 </math> - вектор значений размерности m при <math> x = a </math>,
    функцию вычисления производной <math>\bold f </math> будем считать встроенной в алгоритм.
+
:функцию вычисления производной <math>\bold f </math> будем считать встроенной в алгоритм.
  
 
'''Объём входных данных''':
 
'''Объём входных данных''':
    m + 3
+
:m + 3
  
 
'''Выходные данные:'''
 
'''Выходные данные:'''
    <math> \bold y_i </math> - вектора размера m значений в n узлах сетки (т.е. при <math> x = a + i * \frac{b-a}{n}, i \in [1, n] </math>).
+
:<math> \bold y_i </math> - вектора размера m значений в n узлах сетки (т.е. при <math> x = a + i * \frac{b-a}{n}, i \in [1, n] </math>).
  
 
'''Объем выходных данных'''
 
'''Объем выходных данных'''
    n * m
+
:n * m
  
 
== Свойства алгоритма ==
 
== Свойства алгоритма ==
  
 +
Алгоритм обладает следующими свойствами:
  
= ЧАСТЬ. Программная реализация алгоритма =
+
* Точность. Как написано в [[#Общее описание алгоритма| разделе 1.1]], алгоритм имеет порядок точности <math>O(h^4)</math>.
 +
* Устойчивость. Алгоритм устойчив на любом конечном отрезке.
 +
* Гибкость. На любом этапе вычисления можно легко изменить шаг интегрирования.
 +
* Вычислительная мощность алгоритма стремится к константе. Сумма входных и выходных данных равна mn+m+3, а общее количество операций - 4nCm. Таким образом, при возрастании значений входных параметров, их отношение стремится к 4C.
 +
 
 +
= Программная реализация алгоритма =
  
 
== [[Особенности реализации последовательного алгоритма]] ==
 
== [[Особенности реализации последовательного алгоритма]] ==
Строка 170: Строка 186:
 
== Масштабируемость алгоритма и его реализации ==
 
== Масштабируемость алгоритма и его реализации ==
  
 +
Для исследования масштабируемости алгоритма Рунге-Кутты 4-го порядка была написана своя [https://github.com/D86/superpc/blob/master/rk4.c|параллельная реализация]. этого алгоритма на языке C.
 +
 +
Свойства реализации:
 +
 +
* Реализация написана под компилятор mpicxx. 
 +
* Тестирование проводилось на [http://hpc.cmc.msu.ru/bgp|IBM Blue Gene/P], расположенном на территории ВМК МГУ.
 +
 +
Исследование проводилось по [https://algowiki-project.org/ru/Scalability_methodology|рекомендованной методологии].
 +
 +
Тестирование
 +
 +
* На вход реализации подавалась система из суперпозиций элементарных функций.
 +
* Начальные значения выбирались случайным образом.
 +
 +
Список изменяющихся параметров:
 +
 +
В качестве изменяемых параметров выбраны количество максимальное количество процессоров, выделяемое при запуске исследуемой реализации, и параметр m, так как из [[#Входные и выходные данные алгоритма| раздела 1.9]] следует, что объём входных данных зависит от него.
 +
 +
Границы изменяющихся параметров:
 +
 +
* Число процессоров выделяется с 256 до 2048 с шагом 128;
 +
* Значение m находится в границах от 1000 до 10000 с шагом 1000.
 +
 +
График производительности:
 +
 +
[[Файл:Perf_b.png]]
 +
 +
График эффективности:
 +
 +
[[Файл:Effic.png]]
 +
 +
Оценки:
 +
 +
На основании полученных результатов были составлены графики [https://algowiki-project.org/ru/Глоссарий#.D0.9F.D1.80.D0.BE.D0.B8.D0.B7.D0.B2.D0.BE.D0.B4.D0.B8.D1.82.D0.B5.D0.BB.D1.8C.D0.BD.D0.BE.D1.81.D1.82.D1.8C| производительности] и [https://algowiki-project.org/ru/Глоссарий#.D0.AD.D1.84.D1.84.D0.B5.D0.BA.D1.82.D0.B8.D0.B2.D0.BD.D0.BE.D1.81.D1.82.D1.8C_.D1.80.D0.B0.D1.81.D0.BF.D0.B0.D1.80.D0.B0.D0.BB.D0.BB.D0.B5.D0.BB.D0.B8.D0.B2.D0.B0.D0.BD.D0.B8.D1.8F| эффективности] параллельной реализации, вычислены максимальная и минимальная [https://algowiki-project.org/ru/Глоссарий#.D0.AD.D1.84.D1.84.D0.B5.D0.BA.D1.82.D0.B8.D0.B2.D0.BD.D0.BE.D1.81.D1.82.D1.8C_.D1.80.D0.B5.D0.B0.D0.BB.D0.B8.D0.B7.D0.B0.D1.86.D0.B8.D0.B8| эффективность реализации], 0.000256661928486 и 1.30949963513e-05 соответственно.
  
 
== [[Динамические характеристики и эффективность реализации алгоритма]] ==
 
== [[Динамические характеристики и эффективность реализации алгоритма]] ==

Текущая версия на 23:51, 19 декабря 2016



Метод Рунге-Кутты 4 порядка
Последовательный алгоритм
Последовательная сложность [math]4nCm[/math]
Объём входных данных [math]m + 3[/math]
Объём выходных данных [math]m * n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]O(n)[/math]
Ширина ярусно-параллельной формы [math]O(m)[/math]


Основные авторы описания: Е.П.Степанов(разделы 1.1, 1.2, 1.3, 1.4, 1.9), Б.Савков (разделы 1.5, 1.6, 1.7, 1.8, 1.10, 2.7)

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

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

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

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

[math] \begin{cases} \textbf{y}' = \textbf{f}(x, \textbf{y}),\\ \textbf{y}(x_0) = \textbf{y}_0. \end{cases} [/math]

При этом [math]\textbf{y}, \textbf{f}, x[/math] отвечают следующим свойствам: [math] \textbf{y}, \textbf{f} \in \mathbb{R}^m, x \in \mathbb{R}^1 [/math].

Пусть известно значение y(x) и требуется вычислить значение y(x+h). Это можно сделать при помощи формулы:

[math] y(x + h) = y(x) + \int_0^hy'(x + t)dt[/math]

Аппроксимируя интеграл в правой части, можно получить методы Рунге-Кутты разного порядка точности [1]. В общем виде формулы методов Рунге-Кутты можно записать следующим образом:

[math]\textbf{y}_{n+1} = \textbf{y}_n + h\sum_{i=1}^s b_i \textbf{k}_i[/math]
[math] \begin{array}{ll} \textbf{k}_1 =& \textbf{f}(x_n, \textbf{y}_n),\\ \textbf{k}_2 =& \textbf{f}(x_n+c_2h, \textbf{y}_n+a_{21}h\textbf{k}_1),\\ \cdots&\\ \textbf{k}_s =& \textbf{f}(x_n+c_sh, \textbf{y}_n+a_{s1}h\textbf{k}_1+a_{s2}h\textbf{k}_2+\cdots+a_{s,s-1}h\textbf{k}_{s-1}) \end{array} [/math]

Увеличивая число s вспомогательных точек, можно построить методы Рунге-Кутты любого порядка точности p. Однако уже при p > 5 эти методы используются довольно редко. Это объясняется как чрезмерной громоздкостью получающихся вычислительных формул, так и тем, что преимущества методов высокого порядка точности над методами, в которых p = 4 или p = 5, проявляются либо в тех задачах, где нужна очень высокая точность и используются ЭВМ высокой разрядности, либо в тех задачах, где решение очень гладкое. Кроме того, методы Рунге-Кутты высокого порядка точности часто оказываются менее эффективными по сравнению с методами Адамса того же порядка точности [2]. На этой странице представлен наиболее распространенный метод Рунге-Кутты 4 порядка с ошибкой на конечном интервале интегрирования [math]O(h^4)[/math].

Существуют как явные, так и неявные методы Рунге-Кутты. Последние имеют ряд преимуществ над явными методами, однако это достигается за счет существенного усложнения вычислительного алгоритма, так как на каждом шаге приходится считать систему нелинейных уравнений. На этой странице представлен явный метод Рунге-Кутты.

Можно отметить, что в современных программах, реализующих метод Рунге-Кутты, обязательно используется некоторый алгоритм автоматического изменения шага интегрирования. Интуитивно ясно, что на участках плавного изменения решения счет можно вести с достаточно крупным шагом. В то же время на тех участках, где происходят резкие изменения поведения решения, необходимо выбирать мелкий шаг интегрирования. Обычно начальное значение шага h1 задает пользователь. Далее шаг интегрирования меняется в соответствии с величиной получаемой в ходе вычислений оценки локальной погрешности. Для упрощения описания на этой странице используется постоянный шаг интегрирования h, зависящий от размера сетки разбиения n.

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

Далее подразумевается, что [math] \bold y, \bold y_i,\bold f, \bold k_i \in \mathbb{R}^m[/math], а [math]x, h, n, a, b \in \mathbb{R}^1 [/math].

Исходные данные:

[math] a, b [/math] - интервал интегрирования,
[math] \bold y_0 [/math] - вектор значений при [math] x = a [/math],
n - количество точек в разбиении
[math]\bold f [/math] - функция вычисления производной

Вычисляемые данные:

[math] \bold y_i [/math] - вектора значений в узлах сетки, т.е. при [math] x = a + i * \frac{b-a}{n}, i \in [1, n] [/math]

Формулы метода:

Вначале вычисляется шаг сетки:
[math] h = \frac{b-a}{n}[/math].
Приближенное значение в следующих точках вычисляется по итерационной формуле:
[math]\textbf{y}_{n+1} = \textbf{y}_n + {h \over 6}(\textbf{k}_1 + 2\textbf{k}_2 + 2\textbf{k}_3 + \textbf{k}_4)[/math],
[math]\textbf{k}_1 = \textbf{f} \left( x_n, \textbf{y}_n \right),[/math]
[math]\textbf{k}_2 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_1 \right),[/math]
[math]\textbf{k}_3 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_2 \right),[/math]
[math]\textbf{k}_4 = \textbf{f} \left( x_n + h, \textbf{y}_n + h\ \textbf{k}_3 \right),[/math]
[math]x_{n+1} = x_n + h.[/math]

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

Вычислительным ядром метода Рунге-Кутты можно считать вычисление значений функции производной

[math]\textbf{f}(x, \textbf{y})[/math],

так как при вычислении каждого нового [math]\textbf{y}_n[/math] требуется 4 вычисления значения функции производной. Таким образом, сложность работы алгоритма существенно зависит от сложности вычисления функции [math]\textbf{f}[/math].

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

Макроструктура алгоритма может быть представлена в следующем виде:

1. Вычисление вспомогательных коэффициентов [math]\textbf{k}_j, j \in [1,4][/math];

2. Вычисление приближенного значения [math] \bold y_i [/math] в следующей точке;

3. Переход на следующую итерацию или остановка.

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

Пример реализации последовательного алгоритма на языке C++ (указатель на функцию f передаётся через аргументы представленной ниже функции):

 1 #include <vector>
 2 
 3 // other code ...
 4 
 5 std::vector<double> rungekutta4(double a, double b, std::vector<double> y0, int n, std::vector<double (*)(double,double)> f)
 6 {
 7     std::vector<double>::size_type sz = y0.size();
 8     std::vector<double> result;
 9     for (unsigned j=0; j<sz; j++) 
10     {
11         rezult.push_back(core(a,b,y[j],n,f[j]));
12     }
13     return result;
14 }
15 
16 double core(double a, double b, double y0, int n, double (*f)(double,double))
17 {
18     double h, k1, k2, k3, k4, x, y;
19     h = (b - a)/n;
20     y = y0;
21     x = a;
22     for (unsigned i = 0; i < n; i++)
23     {
24         k1 = f(x, y);
25         k2 = f(x + h/2, y + h*k1/2);
26         k3 = f(x + h/2, y + h*k2/2);
27         k4 = f(x + h/2, y + h*k3);
28         x += h;
29         y += h*(k1 + 2*k2 + 2*k3 + k4)/6;
30     }
31     return y;
32 }
33 
34 // other code ...

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

Рассмотрим сначала одномерный случай. Наравне с обычными операциями, такими как сложение, умножение и тому подобными, в представленном выше алгоритме присутствует вызов функции [math] f [/math]. Так как эта функция может содержать в себе любое количество простых математических операций, будем считать, что её сложность [math] f [/math] много больше любой из них, и равна константе [math]C[/math]. За весь алгоритм, она вызывается [math]4n[/math] раза. Таким образом, сложность последовательного алгоритма для одномерного случая равна [math]4nC[/math].

Для случая решения системы ОДУ, когда [math]f[/math] и [math]y[/math] являются векторами размерности [math]m[/math], функция core будет вызвана [math]m[/math], и итоговая сложность алгоритма составит [math]4nCm[/math].

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

Информационная структура параллельного алгоритма Рунге-Кутты 4-го порядка представляет собой матрицу изображённую на рисунке ниже. Вдоль оси [math]j[/math] располагаются [math]j[/math]-ые параметры вектора [math]y[/math], соответственно, вдоль [math]i[/math] изображены этапы вычисления [math]y_i[/math] значения j-параметра вектора [math]y[/math]. Таким образом, можно сказать, что вдоль оси [math]j[/math], изображены параллельные процессы. А ось [math]i[/math] представляет собой шкалу времени. Чтобы соотнести её с реальным временем, надо просто посчитать сумму [math]\sum_{i=1}^{n} t_i[/math], где [math]t_i[/math] - время с момента начала вычисления [math]i[/math]-го всеми процессами до момента его завершения всеми процессами.

System graph.png

Можно заметить, что на каждом [math]i[/math]-ом шаге происходит обмен информацией между процессами. Это выражено в виде линии пересекающей все [math]j[/math]-ые процессы на любом [math]i[/math]-ом шаге.

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

Алгоритм Рунге-Кутты 4-го порядка является итерационным алгоритмом. В нём каждый шаг цикла зависит от предыдущего. В самом шаге, каждая последующая операция также зависит от результата предыдущей. Таким образом, упразднить цикл до независимых и соответственно параллельных операций не является возможным. Однако, так как на результат функции [math] f [/math] влияет только [math]i[/math]-ый параметр вектора [math] y [/math], то для каждого [math] i [/math] можно производить независимое вычисление. В итоге, если вычислять результаты выполнения функции [math] f [/math] одновременно для вля всех [math] i[/math]-ых параметров вектора [math] y [/math], получим ту же алгоритмическую сложность [math]4nC[/math], что и в разделе 1.6 для одномерного случая.

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

Входные данные:

[math] a, b [/math] - интервал интегрирования,
n - количество точек в разбиении,
[math] \bold y_0 [/math] - вектор значений размерности m при [math] x = a [/math],
функцию вычисления производной [math]\bold f [/math] будем считать встроенной в алгоритм.

Объём входных данных:

m + 3

Выходные данные:

[math] \bold y_i [/math] - вектора размера m значений в n узлах сетки (т.е. при [math] x = a + i * \frac{b-a}{n}, i \in [1, n] [/math]).

Объем выходных данных

n * m

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

Алгоритм обладает следующими свойствами:

  • Точность. Как написано в разделе 1.1, алгоритм имеет порядок точности [math]O(h^4)[/math].
  • Устойчивость. Алгоритм устойчив на любом конечном отрезке.
  • Гибкость. На любом этапе вычисления можно легко изменить шаг интегрирования.
  • Вычислительная мощность алгоритма стремится к константе. Сумма входных и выходных данных равна mn+m+3, а общее количество операций - 4nCm. Таким образом, при возрастании значений входных параметров, их отношение стремится к 4C.

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

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

2.4 Масштабируемость алгоритма и его реализации

Для исследования масштабируемости алгоритма Рунге-Кутты 4-го порядка была написана своя реализация. этого алгоритма на языке C.

Свойства реализации:

  • Реализация написана под компилятор mpicxx.
  • Тестирование проводилось на Blue Gene/P, расположенном на территории ВМК МГУ.

Исследование проводилось по методологии.

Тестирование

  • На вход реализации подавалась система из суперпозиций элементарных функций.
  • Начальные значения выбирались случайным образом.

Список изменяющихся параметров:

В качестве изменяемых параметров выбраны количество максимальное количество процессоров, выделяемое при запуске исследуемой реализации, и параметр m, так как из раздела 1.9 следует, что объём входных данных зависит от него.

Границы изменяющихся параметров:

  • Число процессоров выделяется с 256 до 2048 с шагом 128;
  • Значение m находится в границах от 1000 до 10000 с шагом 1000.

График производительности:

Perf b.png

График эффективности:

Effic.png

Оценки:

На основании полученных результатов были составлены графики производительности и эффективности параллельной реализации, вычислены максимальная и минимальная эффективность реализации, 0.000256661928486 и 1.30949963513e-05 соответственно.

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

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

Алгоритм Рунге-Кутты 4-го порядка является одним из наиболее используемых для решения ОДУ и их систем, поэтому существует большое число его реализаций для разных языков программирования. Наиболее популярными являются следующие из них:

  • библиотека RK4, которая существует в вариантах для C++, FORTRAN90, FORTRAN77, python и MATLAB;
  • библиотека ode45 для MATLAB;
  • модуль scipy для PYTHON.

3 Литература

  1. Бахвалов Н. С., Жидков Н. П., Кобельков. Г. М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
  2. Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров: Учеб. пособие. — М.: Высш. шк., 1994. — 544 с.