2.4 Масштабируемость алгоритма и его реализации
Свойства реализации:
- Реализация написана под компилятор mpicxx.
Тестирование
- На вход реализации подавалась система из суперпозиций элементарных функций.
- Начальные значения выбирались случайным образом.
Список изменяющихся параметров:
В качестве изменяемых параметров выбраны количество максимальное количество процессоров, выделяемое при запуске исследуемой реализации, и параметр m, так как из следует, что объём входных данных зависит от него.
Границы изменяющихся параметров:
- Число процессоров выделяется с 256 до 2048 с шагом 128;
- Значение m находится в границах от 1000 до 10000 с шагом 1000.
График производительности:
График эффективности:
Оценки:
Основной принцип
Рассмотрим следующую проблему:
- y′знак равнож(т,y),y(т)знак равноy{\ Displaystyle у ‘= е (т, у), \ четырехъядерный у (т_ {0}) = у_ {0}}
мы будем смотреть решимость в дискретном множестве т < т 1 <… < т N . Вместо того, чтобы искать прямой метод, методы Рунге-Кутта предлагают ввести промежуточные точки , чтобы вычислить по индукции значения ( t n , y n ) с
{(тнет,я,yнет,я)}1⩽я⩽q{\ displaystyle \ {(t_ {n, i}, y_ {n, i}) \} _ {1 \ leqslant i \ leqslant q}}
- тнет,язнак равнотнет+противя⋅часнет{\ displaystyle t_ {n, i} = t_ {n} + c_ {i} \ cdot h_ {n}}
где h n = t n + 1 — t n — временной шаг, а c i находится в интервале . Для каждой промежуточной точки отметим соответствующий наклон
- пнет,язнак равнож(тнет,я,yнет,я).{\ displaystyle p_ {n, i} = f (t_ {n, i}, y_ {n, i}).}
Таким образом, для точного решения y задачи имеем
- y(тнет,я)знак равноy(тнет)+∫тнеттнет,яж(т,y(т))dтзнак равноy(тнет)+часнет∫противяж(тнет+тычаснет,y(тнет+тычаснет))dты,∀язнак равно1,…,q,{\ Displaystyle у (т_ {п, я}) = у (т_ {п}) + \ int _ {т_ {п}} ^ {т_ {п, я}} е (т, у (т)) \ mathrm {d} t = y (t_ {n}) + h_ {n} \ int _ {0} ^ {c_ {i}} f (t_ {n} + uh_ {n}, y (t_ {n} + uh_ {n})) \ mathrm {d} u, \, \ forall i = 1, \ dotsc, q,}
- y(тнет+1)знак равноy(тнет)+∫тнеттнет+1ж(т,y(т))dтзнак равноy(тнет)+часнет∫1ж(тнет+тычаснет,y(тнет+тычаснет))dты.{\ Displaystyle у (т_ {п + 1}) = у (т_ {п}) + \ int _ {т_ {п}} ^ {т_ {п + 1}} е (т, у (т)) \ mathrm {d} t = y (t_ {n}) + h_ {n} \ int _ {0} ^ {1} f (t_ {n} + uh_ {n}, y (t_ {n} + uh_ {n}) )) \ mathrm {d} u.}
Мы вычислим эти интегралы квадратурным методом, который мы можем выбрать, чтобы они были разными для двух различных значений i :
- ∫противяграмм(ты)dты≈∑kзнак равно1я-1вяkграмм(противk),∫1грамм(ты)dты≈∑kзнак равно1qбkграмм(противk),{\ displaystyle \ int _ {0} ^ {c_ {i}} g (u) \ mathrm {d} u \ приблизительно \ sum _ {k = 1} ^ {i-1} a_ {ik} g (c_ { k}), \, \ int _ {0} ^ {1} g (u) \ mathrm {d} u \ приблизительно \ sum _ {k = 1} ^ {q} b_ {k} g (c_ {k} ),}
вычислено здесь для g ( u ) = f ( t n + uh n , y ( t n + uh n )) .
Таким образом, метод Рунге-Кутты порядка q будет определяться следующим образом:
Метод часто резюмируется таблицей различных квадратурных весов, называемой таблицей Мясника :
| против1{\ displaystyle c_ {1}} | ||||||
| против2{\ displaystyle c_ {2}} | в21 год{\ displaystyle a_ {21}} | |||||
| против3{\ displaystyle c_ {3}} | в31 год{\ displaystyle a_ {31}} | в32{\ displaystyle a_ {32}} | ||||
| ⋮{\ displaystyle \ vdots} | ⋮{\ displaystyle \ vdots} | ⋱{\ displaystyle \ ddots} | ||||
| противq{\ displaystyle c_ {q}} | вq1{\ displaystyle a_ {q1}} | вq2{\ displaystyle a_ {q2}} | ⋯{\ displaystyle \ cdots} | вq,q-1{\ displaystyle a_ {q, q-1}} | ||
| б1{\ displaystyle b_ {1}} | б2{\ displaystyle b_ {2}} | ⋯{\ displaystyle \ cdots} | бq-1{\ displaystyle b_ {q-1}} | бq{\ displaystyle b_ {q}} |
Тем не менее, метод последовательный .
∀язнак равно2,…,q,∑kзнак равно1я-1вяkзнак равнопротивя{\ displaystyle \ forall i = 2, \ dotsc, q, \ sum _ {k = 1} ^ {i-1} a_ {ik} = c_ {i}}
Усовершенствованный метод Эйлера
Рассмотрим тот же самый пример: дифференциальное уравнение , частное решение, удовлетворяющее условию , промежуток и его разбиение на 10 частей ( – длина каждой части).
Цель усовершенствования состоит в том, чтобы приблизить «красные квадратики» ломаной к соответствующим «зелёным точкам» точного решения .
И идея модификации такова: отрезки должны быть параллельны касательным, которые проведены к графику функции не на левых краях, а «посерединке» интервалов разбиения. Что, естественно, улучшит качество приближения.
Алгоритм решения работает в том же русле, но формула, как нетрудно догадаться, усложняется:, где
Плясать вновь начинаем от частного решения и сразу же находим 1-й аргумент «внешней» функции:
Далее следуют уже знакомые по предыдущему параграфу вычисления , после чего можно рассчитать 2-й аргумент «внешней» функции: .
Теперь находим нашего «монстра», который на поверку оказался не таким уж и страшным – обратите внимание, что это ТА ЖЕ функция , вычисленная в другой точке:
Умножаем результат на шаг разбиения:
Таким образом:
Алгоритм заходит на второй круг, не поленюсь, распишу его подробно:
рассматриваем пару и находим 1-й аргумент «внешней» функции:
Рассчитываем 
и находим её 2-й аргумент: 
Вычислим значение:
и его произведение на шаг:
Таким образом: 
Далее рассматриваем пару и т.д.
![Список методов рунге – кутты содержание а также явные методы [ править ]](https://fuzeservers.ru/wp-content/uploads/7/f/0/7f0973064c6f83454da19032bce8c6e1.jpeg)
Вычисления разумно провести в Экселе (растиражировав формулы по той же схеме – см. видеоролик выше), а результаты свести в таблицу:
Числа целесообразно округлять до 4-5-6 знаков после запятой. Нередко в условии той или иной задачи есть прямое указание, с какой точностью следует проводить округление. Я подровнял сильно «хвостатые» значения до 6 знаков.
По результатам 2-го и 3-го столбцов (слева) построим ломаную , и для сравнения я снова приведу график точного решения :
Результат существенно улучшился! – красные квадратики практически «спрятались» за зелёными точками точного решения.
Однако нет пределов совершенству. Одна голова хорошо, а две – лучше. И снова немецкие:
1.6 Последовательная сложность алгоритма
Рассмотрим сначала одномерный случай. Наравне с обычными операциями, такими как сложение, умножение и тому подобными, в представленном выше алгоритме присутствует вызов функции f . Так как эта функция может содержать в себе любое количество простых математических операций, будем считать, что её сложность f много больше любой из них, и равна константе C. За весь алгоритм, она вызывается 4n раза. Таким образом, сложность последовательного алгоритма для одномерного случая равна 4nC.
Для случая решения системы ОДУ, когда f и y являются векторами размерности m, функция core будет вызвана m, и итоговая сложность алгоритма составит 4nCm.
1.7 Информационный граф
Информационная структура параллельного алгоритма Рунге-Кутты 4-го порядка представляет собой матрицу изображённую на рисунке ниже. Вдоль оси j располагаются j-ые параметры вектора y, соответственно, вдоль i изображены этапы вычисления y_i значения j-параметра вектора y. Таким образом, можно сказать, что вдоль оси j, изображены параллельные процессы. А ось i представляет собой шкалу времени. Чтобы соотнести её с реальным временем, надо просто посчитать сумму \sum_{i=1}^{n} t_i, где t_i — время с момента начала вычисления i-го всеми процессами до момента его завершения всеми процессами.
Можно заметить, что на каждом i-ом шаге происходит обмен информацией между процессами. Это выражено в виде линии пересекающей все j-ые процессы на любом i-ом шаге.
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=; j<sz; j++)
10 {
11 rezult.push_back(core(a,b,yj],n,fj]));
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 = ; i < n; i++)
23 {
24 k1 = f(x, y);
25 k2 = f(x + h2, y + h*k12);
26 k3 = f(x + h2, y + h*k22);
27 k4 = f(x + h2, 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 ...
Реализация алгоритма RK4 (5)
Коэффициенты, найденные Фельбергом для Формулы 1 (вывод с его параметром α2 = 1/3), приведены в таблице ниже с использованием индексации массива по основанию 1 вместо базы 0 для совместимости с большинством компьютерных языков:
| K | А (К) | B (K, L) | С (К) | CH (K) | CT (K) | ||||
|---|---|---|---|---|---|---|---|---|---|
| L = 1 | L = 2 | L = 3 | L = 4 | L = 5 | |||||
| 1 | 1/9 | 47/450 | -1/150 | ||||||
| 2 | 2/9 | 2/9 | |||||||
| 3 | 1/3 | 1/12 | 1/4 | 9/20 | 25 декабря | 3/100 | |||
| 4 | 3/4 | 69/128 | -243/128 | 135/64 | 16/45 | 32/225 | -16/75 | ||
| 5 | 1 | -17/12 | 27/4 | -27/5 | 16/15 | 1/12 | 1/30 | -1/20 | |
| 6 | 5/6 | 65/432 | -5/16 | 13/16 | 27 апреля | 5/144 | 6/25 | 6/25 |
Фельберг предлагает решение системы n дифференциальных уравнений вида:
dуяdИксзнак равножя(Икс,у1,у2,…,уп),язнак равно1,2,…,уп{\ displaystyle {\ operatorname {d} \! y_ {i} \ over \ operatorname {d} \! x} = f_ {i} (x, y_ {1}, y_ {2}, …, y_ { n}), i = 1,2, …, y_ {n}}
итеративно решить для
уя(Икс+час),язнак равно1,2,…,п{\ Displaystyle у_ {я} (х + ч), я = 1,2, …, п}
где h — адаптивный размер шага, который определяется алгоритмически:
Решение представляет собой средневзвешенное значение шести приращений, где каждое приращение является произведением размера интервала , и предполагаемого наклона, заданного функцией f в правой части дифференциального уравнения.
час{\ textstyle h}
k1знак равночас⋅ж(Икс+А(1)⋅час,у){\ Displaystyle к_ {1} = час \ cdot f (x + A (1) \ cdot h, y)}
k2знак равночас⋅ж(Икс+А(2)⋅час,у+B(2,1)⋅k1){\ Displaystyle к_ {2} = час \ cdot f (x + A (2) \ cdot h, y + B (2,1) \ cdot k_ {1})}
k3знак равночас⋅ж(Икс+А(3)⋅час,у+B(3,1)⋅k1+B(3,2)⋅k2){\ Displaystyle к_ {3} = час \ cdot f (x + A (3) \ cdot h, y + B (3,1) \ cdot k_ {1} + B (3,2) \ cdot k_ {2}) )}
k4знак равночас⋅ж(Икс+А(4)⋅час,у+B(4,1)⋅k1+B(4,2)⋅k2+B(4,3)⋅k3){\ Displaystyle к_ {4} = час \ cdot f (x + A (4) \ cdot h, y + B (4,1) \ cdot k_ {1} + B (4,2) \ cdot k_ {2}) + В (4,3) \ cdot k_ {3})}
k5знак равночас⋅ж(Икс+А(5)⋅час,у+B(5,1)⋅k1+B(5,2)⋅k2+B(5,3)⋅k3+B(5,4)⋅k4){\ Displaystyle к_ {5} = час \ cdot f (x + A (5) \ cdot h, y + B (5,1) \ cdot k_ {1} + B (5,2) \ cdot k_ {2}) + В (5,3) \ cdot k_ {3} + B (5,4) \ cdot k_ {4})}
k6знак равночас⋅ж(Икс+А(6)⋅час,у+B(6,1)⋅k1+B(6,2)⋅k2+B(6,3)⋅k3+B(6,4)⋅k4+B(6,5)⋅k5){\ Displaystyle к_ {6} = час \ cdot f (x + A (6) \ cdot h, y + B (6,1) \ cdot k_ {1} + B (6,2) \ cdot k_ {2} + B (6,3) \ cdot k_ {3} + B (6,4) \ cdot k_ {4} + B (6,5) \ cdot k_ {5})}
Тогда средневзвешенное значение:
у(Икс+час)знак равноу(Икс)+CЧАС(1)⋅k1+CЧАС(2)⋅k2+CЧАС(3)⋅k3+CЧАС(4)⋅k4+CЧАС(5)⋅k5+CЧАС(6)⋅k6{\ displaystyle y (x + h) = y (x) + CH (1) \ cdot k_ {1} + CH (2) \ cdot k_ {2} + CH (3) \ cdot k_ {3} + CH ( 4) \ cdot k_ {4} + CH (5) \ cdot k_ {5} + CH (6) \ cdot k_ {6}}
Оценка ошибки усечения:
ТEзнак равно|CТ(1)⋅k1+CТ(2)⋅k2+CТ(3)⋅k3+CТ(4)⋅k4+CТ(5)⋅k5+CТ(6)⋅k6|{\ Displaystyle TE = | CT (1) \ cdot k_ {1} + CT (2) \ cdot k_ {2} + CT (3) \ cdot k_ {3} + CT (4) \ cdot k_ {4} + CT (5) \ cdot k_ {5} + CT (6) \ cdot k_ {6} |}
По завершении шага рассчитывается новый размер шага:
часпешзнак равно0,9⋅час⋅(ϵ|ТE|)15{\ displaystyle h_ {new} = 0,9 \ cdot h \ cdot \ left ({\ frac {\ epsilon} {| TE |}} \ right) ^ {1/5}}
Если , то замените на и повторите шаг. Если , то шаг завершен. Замените на для следующего шага.
|ТE|>ϵ{\ textstyle | TE |> \ epsilon}час{\ textstyle h}часпеш{\ textstyle h_ {новое}}|ТE|⩽ϵ{\ textstyle | TE | \ leqslant \ epsilon}час{\ textstyle h}часпеш{\ textstyle h_ {новое}}
Коэффициенты, найденные Фельбергом для формулы 2 (вывод с его параметром α2 = 3/8), приведены в таблице ниже с использованием индексации массива по основанию 1 вместо базы 0 для совместимости с большинством компьютерных языков:
| K | А (К) | B (K, L) | С (К) | CH (K) | CT (K) | ||||
|---|---|---|---|---|---|---|---|---|---|
| L = 1 | L = 2 | L = 3 | L = 4 | L = 5 | |||||
| 1 | 25/216 | 16/135 | 1/360 | ||||||
| 2 | 1/4 | 1/4 | |||||||
| 3 | 3/8 | 3/32 | 9/32 | 1408/2565 | 6656/12825 | -128/4275 | |||
| 4 | 13 декабря | 1932/2197 | -7200/2197 | 7296/2197 | 2197/4104 | 28561/56430 | -2197/75240 | ||
| 5 | 1 | 439/216 | -8 | 3680/513 | -845/4104 | -1/5 | -9/50 | 1/50 | |
| 6 | 1/2 | -8/27 | 2 | -3544/2565 | 1859/4104 | -11/40 | 2/55 | 2/55 |
В другой таблице Фельберга приведены коэффициенты для RKF4 (5), полученные Д. Сарафяном:
| K | А (К) | B (K, L) | С (К) | CH (K) | CT (K) | ||||
|---|---|---|---|---|---|---|---|---|---|
| L = 1 | L = 2 | L = 3 | L = 4 | L = 5 | |||||
| 1 | 1/6 | 1/24 | -1/8 | ||||||
| 2 | 1/2 | 1/2 | |||||||
| 3 | 1/2 | 1/4 | 1/4 | 2/3 | -2/3 | ||||
| 4 | 1 | -1 | 2 | 1/6 | 5/48 | -1/16 | |||
| 5 | 2/3 | 27 июля | 27 октября | 1/27 | 27/56 | 27/56 | |||
| 6 | 1/5 | 28/625 | -1/5 | 546/625 | 54/625 | -378/625 | 125/336 | 125/336 |
2 Программная реализация алгоритма
2.2 Локальность данных и вычислений
2.2.1 Локальность реализации алгоритма
2.2.1.3 Анализ на основе теста Apex-Map
2.4 Масштабируемость алгоритма и его реализации
2.4.2 Масштабируемость реализации алгоритма
Набор и границы значений изменяемых реализации алгоритма:
- число процессоров с шагом по степеням двойки до 64, далее с шагом 10;
- размерность системы .
В результате проведённых экспериментов был получен следующий диапазон алгоритма:
- минимальная эффективность реализации 0.0000032%;
- максимальная эффективность реализации 0.0024%.
Перечислим некоторые особенности тестируемой параллельной реализации:
- для тестирования использовалась псевдослучайная правая часть системы, состоящая из суперпозиций случайно выбранных элементарных функций;
- начальные значения искомой функции задавались случайным образом;
- так как в общем случае правая часть системы неизвестна, то не представляется возможным измерять производительность в гигафлопсах; в качестве показателя производительности использовалось количество точек, в которых была вычислена функция.
На следующих рисунках приведены графики и эффективности параллельной реализации метода Рунге-Кутты в зависимости от изменяемых параметров запуска.
Рисунок 2. Изменение производительности параллельного алгоритма в зависимости от размерности системы и числа процессоров.
Рисунок 3. Эффективность работы параллельного алгоритма в зависимости от размерности системы и числа процессоров.
Построим оценки масштабируемости протестированной параллельной реализации метода Рунге — Кутты 4-го порядка:
По числу процессов -0.000000768. При увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска уменьшается, однако уменьшение достаточно медленное, что связано с крайне низкой максимальной эффективностью работы параллельной реализации алгоритма. Уменьшение эффективности при увеличении числа процессов для систем с небольшой размерностью объясняется тем, что при количестве процессов превышающем размерность системы часть из них не будет участвовать в вычислениях. Для систем с большой размерностью увеличение процессов негативно сказывается на эффективности, так как увеличивается количество пересылок между ними, что существенно затормаживает работу.
По размеру задачи 0.000000763. При увеличении размерности системы при фиксированном количестве процессов эффективность на рассмотренной области изменений параметров запуска увеличивается, причём увеличение наблюдается на почти всей рассматриваемой области. Это объясняется тем, что при фиксированном количестве процессов увеличение размерности системы приводит к увеличению вычислительной нагрузки на каждый процессор, вследствие чего увеличивается время вычислительной работы относительно времени простоя. С увеличением числа процессов скорость возрастания эффективности при увеличении размерности задачи уменьшается, так как количество пересылок становится достаточно большим.
По двум направлениям 0.00000000174. При одновременном увеличении количества процессов и размерности системы эффективность увеличивается. Скорость возрастания эффективности крайне низкая, что объясняется высокими накладными расходами.
Метод 4-го порядка является самым часто используемым из всех схем Рунге — Кутты, поэтому существует множество его программных последовательных реализаций как коммерческих, так и бесплатных. Одной из самых известных программных реализаций является ode45 в среде MATLAB, у которой имеется большое количество возможностей для настройки численного расчёта, что делает её достаточно удобной для использования. Также метод Рунге — Кутты 4-го порядка реализован в параллельной библиотеке PETSc, которая является свободно распространяемой. Несмотря на то, что в этой библиотеки большинство методов реализовано с использованием параллельных алгоритмов, метод Рунге — Кутта в ней реализован последовательным. Довольно трудно найти параллельную реализацию метода в случае рассмотрения общей задачи Коши, но существует множество научных работ, в которых описываются параллельные реализации метода Рунге — Кутты для конкретных классов систем, например, в работе А.В. Старченко описывается параллельная реализация для линейных систем.
1.2 Математическое описание алгоритма
Далее подразумевается, что \bold y, \bold y_i,\bold f, \bold k_i \in \mathbb{R}^m, а x, h, n, a, b \in \mathbb{R}^1 .
Исходные данные:
- a, b — интервал интегрирования,
- \bold y_0 — вектор значений при x = a ,
- n — количество точек в разбиении
- \bold f — функция вычисления производной
Вычисляемые данные:
- \bold y_i — вектора значений в узлах сетки, т.е. при x = a + i * \frac{b-a}{n}, i \in
Формулы метода:
- Вначале вычисляется шаг сетки:
- h = \frac{b-a}{n}.
- Приближенное значение в следующих точках вычисляется по итерационной формуле:
- \textbf{y}_{n+1} = \textbf{y}_n + {h \over 6}(\textbf{k}_1 + 2\textbf{k}_2 + 2\textbf{k}_3 + \textbf{k}_4),
- \textbf{k}_1 = \textbf{f} \left( x_n, \textbf{y}_n \right),
- \textbf{k}_2 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_1 \right),
- \textbf{k}_3 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_2 \right),
- \textbf{k}_4 = \textbf{f} \left( x_n + h, \textbf{y}_n + h\ \textbf{k}_3 \right),
- x_{n+1} = x_n + h.
Пример реализации
public double tryRungeKuteMethod(int segments)
{
double step, kOne, kTwo, kThree, kFour;
xArray = new double;
yArray = new double;
step = (endX — startX) / segments;
xArray = startX;
yArray = startY;
xArray = endX;
for (int i = 1; i < segments — 1; i++)
{
xArray = xArray + step;
}
for (int i = 1; i < segments; i++)
{
kOne = formFunction(xArray, yArray);
kTwo = formFunction((xArray + step/2),
(yArray + kOne*step/2));
kThree = formFunction((xArray + step / 2),
(yArray + kTwo * step / 2));
kFour = formFunction((xArray + step),
(yArray + kThree * step));
yArray = yArray + step/6 *(kOne+2*kTwo+2*kThree+kFour);
}
return yArray;
}
