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

В результате экспериментальных исследований и испытаний предложенной виртуальной среды моделирования были выявлены пределы достижения адекватности и реалистичности изображений кадров визуализации. Использование текстур значительно (до 30%) снижает скорость рендеринга, поэтому текстуры использовались только в исключительных случаях, или в случаях низких по графической информативности виртуальных сцен.

МЕТОДЫ ПАРАЛЛЕЛЬНОГО МОДЕЛИРОВАНИЯ ЛИНЕЙНЫХ ДИНАМИЧЕСКИХ СИСТЕМ С АППРОКСИМАЦИЕЙ ПРАВОЙ ЧАСТИ

Существует целый ряд обстоятельств, обуславливающих необходимость адаптации известных численных методов решения систем обыкновенных дифференциальных уравнений в вычислительных системах с параллельной архитектурой. В первую очередь к таким обстоятельствам можно отнести тот факт, что более высокая скорость решения может достигаться за счет рациональной организации распараллеливания вычислений.

Данная статья является продолжением работ [1-3], которые посвящены разработке и исследованию параллельных алгоритмов численного решения систем ОДУ, используемых для моделирования сложных динамических систем с сосредоточенными параметрами. Предлагаемые алгоритмы ориентированы на использование в многопроцессорных вычислительных системах SIMD (single instruction stream - multiple data stream) структуры с решеткой или линейкой процессорных элементов. Набор процессоров известен до начала вычислений и не меняется в процессе счета, при этом каждый процессорный элемент может выполнить любую арифметическую операцию за один такт, временные затраты, связанные с обращением к запоминающему устройству, отсутствуют.

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

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

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

A - матрица коэффициентов системы.

В этом случае решение можно получить последовательно по шагам с помощью численных методов заданного порядка точности.

Здесь вычисление значения вектора неизвестных хп+1 на очередном шаге

требует предварительного определения значений хп . В [1] рассмотрены вопросы, связанные с возможностью параллельной реализации таких алгоритмов. В частности, если система (1) является однородной, т.е. )= 0, 1=1,2,... т, тогда, в зависимости от выбранного метода интегрирования, можно искать решение в виде

где О - оператор (матрица) переходов.

Полученный оператор перехода О, который необходимо определить один раз до начала вычислений, позволяет вычислять значения вектора неизвестных

О.А. Дмитриева

dx = Ax + f(t), x(t( )= x0 = (x0, x02,...,x0mf, dt

(1)

(2)

параллельно [2,3]. Для методов Рунге-Кутты, например, этот оператор может быть представлен, в зависимости от точности метода, как

в = E+тА(Е+ ^А(Е+ ^А(Е+ ^А)),, (3)

который обеспечивает точность 4-го порядка, или, например, как

в = Е+тА(Е+^А(е+ Ц-(Е+^А(Е ~(Е ~)))), (4)

точность которого оценивается 6-м порядком.

При решении неоднородной системы необходимо дополнительно вычислить на каждом шаге значения всех функций /.(г), . =1, т в нескольких промежуточных точках. Поскольку все эти функции могут быть различными, одновременное вычисление их на 81МБ-компьютере невозможно.

В связи с этим можно предложить два основных подхода, первый из которых состоит в том, что все промежуточные значения функций fi(t) могут быть

вычислены заранее. При этом, если количество промежуточных точек метода определяется как г, то можно оценить количество тактов расчета на топологических структурах линейке из т процессоров и решетке из т*т процессоров (размерность процессорного поля выбрана совпадающей с размерностью системы уравнений исключительно для удобства изложения). Определим для этого трудоемкость реализации правых частей системы как о-- и при расчете будем оперировать с

максимальным значением, которое обозначим как Оf = тах{О^ }. Если расчет

осуществляется для общего количества узлов, равного Ы, то общее число тактов работы на линейке процессоров составит для одного уравнения ближайшее целое

т*вг* N"

--------- +1. Т огда вся система может

r*Оf * N

сверху соотношения --------- или

т

быть рассчитана за г*0-*М + т тактов. На решетке процессоров одно уравнение

г*0-* N

будет решаться за ближайшее целое сверху к ---—--- тактов, а время, которое

т

■ *0-* N

+ 1 . По каждому

потребуется для расчета всей системы, составит

т

уравнению системы придется хранить двумерные массивы размерностью г*Ы и использовать коэффициенты с нужными индексами при расчете.

Второй подход, который позволит избежать последовательных участков при параллельной реализации системы, основывается на предварительном интерполировании правых частей (1). В вычислительной практике с таким подходом сталкиваются, если приходиться заменять одну функцию /(I) (известную, неизвестную или частично известную) некоторой функцией <р(1), близкой к /($ и обладающей определенными свойствами, позволяющими производить над нею те или иные аналитические или вычислительные операции. Такую замену называют приближением функции(). Тогда при решении задачи вместо функции/($ оперируют с функцией <р(1), а задача построения функции (р(.1) называется задачей приближения. Исходя из проблематики задачи, т.е. принимая во внимание большое количество узлов, которые будут участвовать в расчете, и задаваясь требуемой точностью приближения функций, можно утверждать, что наиболее перспективным является случай, когда используется кусочно-полиномиальная аппроксимация, или сплайны, так как при этом интерполяционный многочлен строится не на весь интервал

решения задачи, а на подынтервалах, что позволит избежать накопления ошибок приближения.

Основная идея такого подхода заключается в следующем: исходный отрезок решения для (1) [0,Т] разбивается на несколько подынтервалов V с шагом, определяющимся из соотношения точности методов численного интерполирования и интегрирования, а затем на каждом таком интервале строится интерполяционный многочлен. Поскольку в качестве интерполяционной функции обычно выбирают многочлены степени не выше 3-4-й, что соответственным образом влияет на точность интерполяции, то необходимо предварительно согласовать порядки точности методов численного интегрирования и предварительного интерполирования.

Рис. 1 Пример схемы расчета правых частей (1) с промежуточными точками

Если порядок метода численного интегрирования (1) определяется как 0(ту ), а порядок сплайна как 0(Н4 ), то между шагами двух решающихся задач

должно выполняться соотношение ту = Н4 . Если порядок точности метода интегрирования у=4 или более, т.е. между количеством узлов задач интерполирования и интегрирования выполняется соотношение У^~, то использование интерполирования для восстановления значений правых частей является нерациональным. Проще заранее вычислить значения правых частей на промежутке [0,Т]. Если же речь идет о методах интегрирования, которые имеют порядок погрешности ниже 4, то тогда V<N и при этом значения V и N можно связать с помощью некоторого коэффициента в, т.е.

N =b*V , где Ь >>1. (5)

При этом желательно выбирать множитель Ь целым, т.к. предпочтительнее, чтобы на одном такте расчета использовались коэффициенты одного интервала сплайна, что значительно упростит алгоритм вычисления и выбор нужного интервала по заданному аргументу.

Если исходить из (5), то узлов интерполирования будет в Ь раз меньше, чем узлов интегрирования. К тому же оценку погрешности для сплайна порядка

0(Н4 ) можно считать завышенной. Тогда исходная задача может быть сведена к двум подзадачам, каждая из которых легко распараллеливается.

Первая подзадача будет заключаться в нахождении коэффициентов сплайна. Замену исходных функций f■^(t) в (1) на сплайн-функции будем осуществлять в виде

2 з 4 • --- ---

аг/ + V + сиг + йг/г + еиг , г = 1>т>1 = 1,у. (6)

Тогда, во время реализации второй подзадачи, вместо разнородных операций (которые на 81ЫБ-структурах выполняются последовательно), все правые части системы (1) будут считаться параллельно по одним и тем же аргументам г, но с разными коэффициентами сплайн-функций а*/, Ъц, с ц, й */, ец, г = 1,т,1 = IV. Для нахождения неизвестных коэффициентов по каждому уравнению системы (1) придется решать систему линейных алгебраических уравнений размерностью 4Vx4V. Формирование системы линейных уравнений будет исходить из принципов совпадения значений функции, ее первых, вторых и третьих производных на соседних подынтервалах слева и справа от узла интерполяции.

Пусть {р1 } - множество узлов интерполяции, / = 0,У. Причем, для любой правой части системы (1) это множество узлов будет постоянным.

Между узлами р/ _1 и р/ будем представлять функцию /'- го уравнения системы в виде

$(Р)= аг/ + Ъг/ (Р _ Р1 _1)+ сг/ (Р _ Р1 _1 / + йг 1(Р _ Р1 _1 / + еП(Р _ Р1 _1 )4 , (7) Р1 _1 < р < Р1, г = 1т, I =IV.

Исходя из постановки задачи интерполирования, в узлах интерполяции значения исходной функции и интерполяционного многочлена должны совпадать, т.е.

$(р/) = 1п г = 1т, / = 0,У . (8)

С другой стороны, если аргумент сплайна совпадает с узлом интерполяции, тогда из (7) и (8)

(р/_1) = аа = /а_1, г = 1,т, / = 1,У. (9)

Поскольку интерполяционный многочлен строится для равностоящих узлов,

т.е.

р/ _ р/_1 = Н , (10)

то, используя полученные соотношения, можно переписать формулу для сплайна в

1-м узле в следующем виде:

2 3 4

!й = !й_1 + Ъг/Н + сг/Н + йг/Н + ег/Н . (11)

Для построения оставшихся соотношений воспользуемся соглашением о совпадении 1-й, 2-й и 3-й производных справа и слева от узлов интерполяции, т.е. для первых производных

(р)= Ъг/ + 2сг/(р _ р/_1 )+ 3йг/(р _ р/_1 ) + 4ег/(р _ р/_1 ?,

р _1 < р < р/, (12) $ (р)= Ъг/+1 + 2сг/+1 (р _ р/)+ 3йг/+1 (р _ р/) + 4ег/ (р _ р/?, р/ < р < р/+1■

Для каждого уравнения системы (1) будем приравнивать (12) в точке р/,

тогда

2 3

Ъг/+1 = Ъг/ + 2сг/Н + 3йг/Н + 4ег/Н . (13)

Аналогично приравнивая значения полученных выражений для 2-й и 3-й производных в узлах интерполяции и добавив граничные условия, получаем т систем линейных уравнений.

Сформированные системы относительно вектора неизвестных коэффициентов можно представить в общем виде как

О У г = gi, г = 1,т . (14)

Вектора неизвестных будут иметь вид

Уг = (аг1 ,Ъг 1 ,сг 1 ,йг 1,ег 1,...’аг V,Ъгу,сг ¥,йг¥,ег¥ ), г =1,т . (15)

Особенностью полученных систем является ленточный вид матрицы Q, т.к. каждое уравнение системы (за исключением первого и последнего) будет содержать только 4 неизвестных. В этом случае систему можно преобразовать так, чтобы ее можно было решать методом встречной прогонки [4]. Трудоемкость решения таких систем на параллельных 81ЫБ структурах линейно зависит от размерности решаемой системы и для системы размерностью к оценивается как 0(k). Тогда для нашего случая трудоемкость нахождения коэффициентов сплайн-функции для одного уравнения системы будет оцениваться на уровне 0(4У). Тогда для исходной задачи (1) число операций приближенно будет оцениваться на уровне 4*У*т.

Еще один возможный подход к решению полученных систем (14), которые имеют большую размерность и разреженную матрицу коэффициентов, заключается в приведении ее к блочно-диагональной форме с обрамлением [5,6] и формировании вспомогательной системы значительно меньшей размерности, которая определит вектор определяющих величин, или переменных связи. Трудоемкость реализации такого подхода на параллельных вычислительных структурах будет, как и в предыдущем случае, линейно зависеть от размерности системы.

Кроме того, для интерполирования необходимо предварительно вычислить значения правых частей в V точках, которые будут использоваться в качестве исходных данных для построения интерполяционного многочлена, тогда для системы из т уравнений потребуется т* (4У + У*0^ ) тактов. Также возникает

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

2 3 4

ставит для многочлена общего вида а + Ъг + сг + йг + ег на 81ЫБ-структуре умножения

1-й такт -

2-й такт - г2 * г, г2 * г2 - параллельно

3-й такт - Ъ*г,с*г2, й*г3 ,е*г4 - параллельно; сложения

4-й такт - а + Ъ*г,с*г2 + й*г3

2 3 4

5-й такт - а + Ъ*г + с*г + й*г + е*г .

Таким образом, на определение одного значения правой части необходимо 5 временных тактов. Если учесть, что число точек, в которых необходимо восстанавливать значение правой части каждого уравнения определяется как г*N, то всего на восстановление значений функции по интерполяционному многочлену для системы потребуется 5*m*г*N временных тактов.

Сведем полученные приближенные результаты для 2-х описанных способов реализации правых частей в следующую таблицу

Таблица 1

Предварительный расчет Интерполирование

Общее число операций г*0у *М*т т*(4У + У*0у + 5*г*М )

Число операций на линейке процессоров [г*0у*^]+ 1 4У + У*0у + 5*г*М

Число операций на решетке процессоров Гг*0у *м 1 у +1 т (4У + У*0у + 5*г*М )/т

Поскольку изначально предполагалось, что трудоемкости вычисления правых частей 0у являются высокими [7], то оценка трудоемкости всего алгоритма

может осуществляться относительно этих значений. Т огда очевидно, что в случае выполнения соотношения (5), предпочтительнее интерполировать правые части, хотя этот подход и сопряжен с алгоритмическими сложностями, но имеет безусловные преимущества.

Таким образом, в представляемой работе исследованы возможности распараллеливания известных последовательных алгоритмов численного решения систем обыкновенных дифференциальных уравнений и переноса их на параллельные вычислительные структуры с целью получения максимальной реальной производительности. Предложены подходы, позволяющие избегать последовательных участков работы многопроцессорных вычислительных SIMD-систем, один из этих подходов основан на предварительном вычислении правых частей неоднородной линейной системы ОДУ, который предпочтительнее использовать, если точность метода интегрирования, с помощью которого решается задача, высока. Разработан также подход, исключающий последовательные вычисления, который основывается на предварительном интерполировании правых частей системы (1) с помощью сплайнов. Выявлены соотношения между порядками погрешностей методов численного интегрирования системы (1) и методами интерполирования, которые позволяют определить оптимальный выбор метода параллельной реализации правых частей.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

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

2. Фельдман Л.П., Дмитриева О.А. Эффективные методы распараллеливания численного решения задачи Коши для обыкновенных дифференциальных уравнений. //Математическое моделирование. Том 13. 2001. № 7, - С.66-72.

3. Дмитриева О.А. Параллельное моделирование динамических объектов с сосредоточенными параметрами. Тезисы докладов XII Юбилейной международной конференции по вычислительной механике и современным прикладным программным средствам.- М.:МГИУ, 2003.

4. Гаранжа В.А., Коньшин В.Н. Прикладные аспекты параллельных высокоточных алгоритмов решения задач вычислительной гидродинамики.// Тезисы докладов Всероссийской научной конференции «Фундаментальные и прикладные аспекты разработки больших распределенных программных комплексов». - М.: Изд-во МГУ. - 1998. - С. 34-38.

5. Джорт А., Лю Д. Численное решение больших разреженных систем уравнений. - М.: Мир, 1984. - 333 с.

6. Дмитриева О.А. Анализ параллельных алгоритмов численного решения систем линейных уравнений итерационными методами. //Научн. труды ДонГТУ. Серия: Про-

блемы моделирования и автоматизации проектирования динамических систем. Вып. 10: - Донецк:, 2000. - С. 15-22.

7. Feldman L., Dmitrieva O., Gerber S. Abbildung der blockartigen Algoritmen auf die Paral-lelrechnerarchitekture. 16 Symposium Simulationstechnik ASIM 2002, Rostock, 10.09 bis 13.09 2002. - Erlangen: Gruner Druck, 2002. - P. 359-364.

Ю.И.Карпенко, Е.Б.Орлова

ПОСТРОЕНИЕ МОДЕЛЕЙ НАПРЯЖЕННО-ДЕФОРМИРОВАННЫХ СОСТОЯНИЙ В ТЕОРИИ ОРТОТРОПНЫХ ЦИЛИНДРИЧЕСКИХ

ОБОЛОЧЕК

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

Разрешающее уравнение в теории ортотропных цилиндрических оболочек, основанной на гипотезе Тимошенко, может быть записано в следующем виде:

. д1()Ф д10Ф д10Ф д10Ф . д10Ф

A]0 0 -77Г + Ав 2 - -— + Аб 4- - + Ая 6 --;-Т- + А

l10,0—JC + А8,2 8 2 + А6-4^ 0W + А,6 -л 4^ об + А28 л 2^а8 '

da da dp da dp4 da dp da dp

. d10Ф . d8(& . d8Ф . d8Ф . d8Ф

+ А0101^Т0 + А8.°1~8 + ,2 ?,„,б^п2 + А4 ^4^4 + А2б^2^б +

dp10 ’ da6 ’ da0dpJ ’ da4dp4 ’ daJdpc

dбФ

+ А4 2----:---7Г + А2

„ д8Ф „ д6Ф , д6Ф „ д6Ф „ д6Ф

+А8 р+ да6 + 42 ар+ ар+ др6+

, д4Ф д4Ф л д4Ф п

+ А4,0Т~Т + А2’2 -ч 2лп2 + А°4^4 = ° да4 да др др4

Решение разрешающего уравнения ищется в виде разложения функции Ф в тригонометрические ряды

Ф(а;Р )= '^фп1(а)cosnр + Фn2(a)sin ПР\ (2)

П =°

Таким образом, каждая функция Фп удовлетворяет обыкновенному дифференциальному уравнению десятого порядка с постоянными коэффициентами. Соответствующее ему характеристическое уравнение будет иметь вид

А1°,°к1° + {А8,° - А8,2П )к8 + (а6° - А6,2п2 + А6,4п4)К +

+ {А4,° + А4,4п4 - А4,2п‘2 - А4,6П )К + (А2,8п8 - А2,6п6 + А2,4п4 - А2,2п2 )К + (3)

+ (а°,4п4 - А),6п6 + А),8п8 - А°,1°п1° )= °=

где величина п принимает только неотрицательные целые значения.

Коэффициенты характеристического уравнения (3) зависят от геометрических и физических параметров оболочки, следовательно,