Моделирование динамических объектов с помощью параллельных алгоритмов численного решения систем обыкновенных дифференциальных уравнений Текст научной статьи по специальности «Математика»

Л.П. Фельдман, О.А. Дмитриева

МОДЕЛИРОВАНИЕ ДИНАМИЧЕСКИХ ОБЪЕКТОВ С ПОМОЩЬЮ ПАРАЛЛЕЛЬНЫХ АЛГОРИТМОВ ЧИСЛЕННОГО РЕШЕНИЯ СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Моделирование сложных динамических систем с сосредоточенными параметрами обычно требует решения систем обыкновенных дифференциальных уравнений высокого порядка. Сокращение времени решения и, следовательно, времени моделирования, возможно при использовании мультипроцессорных систем, эффективность работы которых существенно зависит от характеристик применяемых алгоритмов. В работе предлагается подход, позволяющий распараллеливать известные последовательные методы и осуществлять их отображение на современные вычислительные структуры. Разрабатываемые алгоритмы ориентированы на использование в многопроцессорных вычислительных системах SIMD (single instruction stream/ multiple data stream) структуры типа MasPar 2000/3000 с решеткой или линейкой процессорных элементов. Набор процессоров известен до начала вычислений и не меняется в процессе счета, при этом каждый процессорный элемент может выполнить любую арифметическую операцию за один такт, временные затраты, связанные с обращением к запоминающему устройству, отсутствуют.

Если математическую модель динамической системы можно представить в виде системы ОДУ с постоянными коэффициентами и начальными условиями

dt = Ax+f(t), x(t0) = x0 = (x0,x2,...,xm)t, (1)

где x - вектор неизвестных сигналов,

f(t) - вектор воздействий, te [0,T],

A = (ay), i,j = 1,m - матрица коэффициентов,

то решение можно получить последовательно по шагам с помощью, например, следующей формулы Рунге-Кутты с контролем погрешности на шаге:

xn+' = xn + 1(kJ1+ 4k^+kn)+rn ,n = 0,1,2,...N (2)

где rn = -^(42k n+224k n+21k n-162k n-125k n)- погрешность на шаге 336 1 3 4 5 6'

k J1= t[Axn + fu], k n =t[A(xn + 0.5k J1) + fn+°'5] ; (3)

-n+0.5.

kn=t[A(xn+I(kn+k“)+f“'-] ;

kn=t

A(xn - k n + 2k n) + fn+1

k n = t[A(xn + 2y(7k ^ +1 0k n+kn) + fn+2/3] ; k n=t[A(xn + ^(28k n - 125k n+ 546 k n+ 54k n- 378k n) + fn+1 /5],

X - выбранный шаг интегрирования

Схема решения ОДУ методом Рунге-Кутты с контролем погрешности

(рис.1)

Рис. 1

Здесь вычисление значения вектора неизвестных х

.П

-П+1

требует предвари-

тельного определения значений к. , і =16. Эти векторы, как следует из формул

(3), должны вычисляться последовательно. Если переписать их для однородной системы в виде, удобном для оценки времени выполнения вычислений, то можно выразить на очередном шаге полный оператор, позволяющий определить значение вектора неизвестных на следующем шаге. Назовем его оператором (матрицей) перехода. Используя для вывода формулы систему МаґНетаґіса®, выразим ^ для і=2,3,4,5,6 через ^, подставим найденные выражения в (2) и получим для однородной системы следующие выражения:

-п+1 ,тА^ ,тА^ тА™ хАчччччч-п ...

х = (Е+тА(Е+^(Е + ^(Е+^(Е —^(Е —^))))))х , (4)

или хп+1 = ОХП, (5)

где О = Е + тА(Е + ТА(Е + Т3А(Е + ТА(Е _Х^(Е

- оператор перехода, E - единичная матрица.

Полученный оператор перехода G, который необходимо определить один раз, до начала вычислений, позволяет вычислять значения вектора неизвестных параллельно. При этом вид оператора перехода будет зависеть от точности выбранной формулы Рунге-Кутты. Так, например, для методов Рунге-Кутты 4-го порядка точности оператор перехода будет иметь вид [1]

хп+1 = (Е + хА(Е + ХА(Е + ХА(Е + >) хп .

(6)

Граф влияния полученного алгоритма, позволяющий выявить операции, которые могут выполняться параллельно, приведен на рис. 2. Первоначально вы-

полняется умножение gij * XJn во всех процессорных элементах. Затем осуществляется одиночный циклический сдвиг полученных значений и их сложение. Количество позиций, на которое производится очередной сдвиг, определяется элементами следующего ряда: 20,21,22,...,2к-1, где k = ближайшее целое сверху [ log2 т ]. Легко видеть, что на такой сетке в каждой ячейке ьой строки получается новое значение Хп+1 . Использование систолического алгоритма умножения [2] позволяет сократить количество сдвигов.

Рис. 2

При решении неоднородной системы (1) по формулам с контролем погрешности необходимо дополнительно вычислить на каждом шаге значения всех

— ---- —n —n+1/2 —n+1 —n+2/3 —n+1/5

функций fi(t) , i = 1, m в точках f ,f ,f ,f ,f . Поскольку все эти функции могут быть различными, одновременное вычисление их на SIMD-компьютере невозможно. Но для линейной системы обыкновенных дифференциальных уравнений функции f i (t) могут быть вычислены заранее.

При этом, если количество точек, в которых вычисляются промежуточные значения, не превышает число процессорных элементов, то на каждом шаге можно организовать параллельное вычисление значений f i (t). Время, которое для этого потребуется, будет определяться как m*N*8, где 8 = max(8i) - макси-

i=1,m

мальное время последовательного расчета f i (t) для одного значения аргумента.

ЛИТЕРАТУРА

1. Фельдман Л.П. Параллельные алгоритмы численного решения систем линейных обыкновенных дифференциальных уравнений// Математическое моделирование. -2000. Т. 12. № 6. - С. 15-20

2. Дмитриева О.А. Анализ параллельных алгоритмов численного решения систем обыкновенных дифференциальных уравнений методами Адамса-Башфорта и Адамса-Моултона// Математическое моделирование. 2000. Т. 12. № 5.- С. 81-86.