Метод Галеркина. Способ параллельной реализации задач математической физики в едином вычислительном потоке

УДК 681.3

Ю.Я. ЛЕДЯНКИН

МЕТОД ГАЛЕРКИНА. СПОСОБ ПАРАЛЛЕЛЬНОЙ РЕАЛИЗАЦИИ ЗАДАЧ МАТЕМАТИЧЕСКОЙ ФИЗИКИ В ЕДИНОМ ВЫЧИСЛИТЕЛЬНОМ ПОТОКЕ

Анотація. Розвивається ідея застосування методу Гальоркіна як різновиду методу зважених не-в'язок у створенні єдиного технологічного потоку в організації обчислень. Запропоновано спосіб паралельної обробки інформації на рівні складних структур даних і описано на конкретному прикладі рішення нелінійної задачі методом Гальоркіна. Стаття представляє інтерес для спеціалістів у галузі математики і розробників методів і структур спецпроцесорів для розв'язання задач математичної фізики та інших завдань.

Ключові слова: метод Гальоркіна, паралельні обчислення, складні структури даних, скалярний помножувач, завдання математичної фізики, єдиний обчислювальний потік.

Аннотация. Развивается идея применения метода Галеркина как разновидности метода взвешенных невязок в создании единого технологического потока в организации вычислений. Предложен способ параллельной обработки информации на уровне сложных структур данных и описан на конкретном примере решения нелинейной задачи методом Галеркина. Статья представляет интерес для специалистов в области математики и разработчиков методов и структур спецпроцессоров для решения задач математической физики и других задач.

Ключевые слова: метод Галеркина, параллельные вычисления, сложные структуры данных, скалярный умножитель, задачи математической физики, единый вычислительный поток.

Abstract. The idea of using the Galerkin method develops as a form of the weighted residuals method to create a single technological flow in calculating organization. A way ofparallel information processing at the level of complex data structures is proposed and it described by the specific example of solving the nonlinear problem by the Galerkin method. This article is useful for mathematicians and developers of methods and structures of specially designed processor for solving problems of mathematical physics and other problems.

Keywords: Galerkin method, parallel computing, complex data structures, scalar multiplier, problems of mathematical physics, single computational flow.

1. Введение

Решение задач математической физики (МФ) связано с аппроксимацией с различной точностью дифференциальных и интегральных уравнений и систем, приводя их к линейным алгебраическим уравнениям (СЛАУ). Успешно решая СЛАУ с помощью средств вычислительной техники, в итоге получают некоторое псевдорешение задачи. Степень точности и суммарные затраты на получение псевдорешения определяются в большой степени применяемыми методами аппроксимации. Последние десятилетия отмечены бурным ростом (и внедрением) группы методов под общим названием “методы взвешенных невязок” (МВН). К этой группе относят и метод Г алеркина (МГ) как частный случай, в котором весовые функции совпадают с базисными.

Применение метода Г алеркина в сочетании с методом конечных элементов (МКЭ) особенно эффективно за счет кусочно-гладкой полиномиальной аппроксимации. Считают, что МГ позволяет выбрать такие пробные функции, которые позволяют с их минимально возможным набором получить решение с высокой точностью. Кроме этого, следует заметить, что МГ практически не имеет ограничений относительно типа задач, в то время как другие (метод наименьших квадратов (МНК) и др.) не подходят к использованию при решении задач с начальными и граничными условиями, задач на собственные значения. А в целом МВН имеют преимущества перед вариационным подходом решения задач МФ с помощью МКЭ.

© Ледянкин Ю.Я., 2012

ISSN 1028-9763. Математичні машини і системи, 2012, № 3

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

Следует сказать, что сделана еще одна попытка трансформировать стандартные численные методы. В целом они разработаны под конкретные задачи, но решаются на ЭВМ с последовательной числовой обработкой информации и часто не учитывают специфику технических средств, на которых впоследствии будут реализованы. А современные требования к решению задач (время решения, точность, размерность,...) предполагают переход к обработке информации в виде ССД с широким распараллеливанием вычислительного процесса. И для этого имеется современная техническая (технологическая) база.

2. Суть метода Г алеркина

Если заданы системы [1] уравнений

Ь

Г (и) - q = 0 и аппроксимирующая (базисная) функция

N

и = 2 ак фк ,

к=1

которая удовлетворяет граничным условиям

Г

а невязка

е=Ь

ортогонализуется по отношению к базисным функциям ф г

<е, фг >= 0, г = 1,2,..., N,

то

для линейного оператора Ь система (5) переходит в систему линейных уравнений относительно коэффициентов ак .

3. Способ параллельного решения нелинейных задач МГ

Решение нелинейных задач (НЗ) сложнее, чем линейных и, естественно, занимает больше времени и ресурсов. Но мир в целом нелинейный. Поэтому к решению НЗ требуется серьёзный подход.

Для получения псевдорешения процедуру несколько упрощают. Берут аппроксимирующие функции с линейными членами и итерационно уточняют линейное псевдореше-

ние до тех пор, пока ошибка аппроксимации станет приемлемой. Решение двух- и трехмерных задач МФ процесс трудоемкий, но в описательном плане он достаточно однотипный. Поскольку применение МКЭ при решении задач МФ (в ЕТП с вариационным подходом к решению) ранее рассмотрено в [2], ниже ограничимся рассмотрением одномерного варианта задачи.

4. Линейный случай

В уравнении, взятом из [1] в качестве примера,

Уди / Эх = БЭ 2и / Эх2 - ки .

(Замечание. Значения нижеприводимых коэффициентов также взяты из [1]. Это обеспечит корректность при сравнении численных значений, получаемых и приводимых автором в статье, с численными значениями, приведенными в [1]), где V - скорость жидкости, и -концентрация, Б - коэффициент диффузии, к, 5 - параметры затухания. Положив 5 = 2 и определив граничные условия (ГУ), и = и при х = 0, ёи / ёх =0 при х = /, аппроксимировав искомую функцию и (х) выражением

и = и

с ГУ и = и при х =0, х = /, получим

ёи / ёх | х=1 = и

откуда подстановкой в (7) значения а2 обеспечивается переход к функции с одной неизвестной:

и = и

5. Рассмотрим уравнение равновесия

Запишем его в соответствии с представлением по Галеркину:

/

| {V * Эи / Эх - Б * Э2 и / Эх 2+к * и } да * ёх = 0,

0

где ди - произвольное приращение функции и(х).

После подстановки в (10) уравнения (9), а также вычисления нами значений первой и второй производной, приведения подобных членов в подынтегральном выражении, уравнение (10) принимает следующий вид:

1

| {(V а * I + к * I + а * Б )* х/1-(V * а * I + а * Б + к * I)* х2/212 +

0

+

6. Обработка исходных данных в виде ССД

В соответствии с правилами записи в виде РМП, приведенными в [2, 3], и с учетом, что матрица для интегрирования [А] равна

[ ¿І

0 0 0 0 0

1 0 0 0 0

0 1/2 0 0 0

0 0 1/3 0 0

0 0 0 1/4 0

0 0 0 0 1/5

а также с учетом записи матрицы, определенной (там же) в виде /*(х), интегральное уравнение (11) принимает вид

"0

(V * а * / + к * / + а * Б)х / /

[А] *

[-

0

0

[(у * а * / + к * / + а * Б)х2 / /]/2/

[-(у * а * / + а * Б + к * /) х3 / 2/2 + (к * а * / - V * а) х3 / / ]/3 [у * а - к * а * /)х4/2/2 - к * ах4/2/]/4 [(к * а * х5 / 4/2 / 5

=0* / •

Замечание. Следует помнить, что количество столбцов (т) матрицы [А] должно быть равно количеству строк (т) вектора столбца /*(х), а результирующая матрица (вектор столбец) - размерность (тК ) на единицу больше (тК = т +1).

В уравнении (13) вектор-столбец в соответствии с правилом записи в виде РМП [2,

3, 4] равен первой строке верхней треугольной матрицы. Подставив в вычисленное значение интеграла (13) пределы интегрирования, выполним приведение подобных членов и получим:

( V * а * / + к * / + а * Б * / )/2-( V * а * / + а * Б * / + к • / )/6 / +( к * а * / -

- V * а ) * / 3/3 / +( V * а - к * а * / ) * / /8 / 2- к * а * / /2 / + к * / /20 * / = = V * а * / (1/2-1/6-1/3+1/8)+ а * Б * / (1/2-1/6)+1/2 к * / 2 -1/6 к * / 2 +

+ к * а * / 3 (1/3-1/8+1/20)=0.

В результате находим

V * а * /2/8+ а * Б * //3+к * /2/3 +2к * а * /3/15 = 0,

а =40 к * //(15 V * / +40 Б +16 к * /).

Замечание. В (14) каждая составляющая левой части уравнения вычисляется в спецпроцессоре параллельного типа с помощью СУ, заложенного в структуру каждого ПЭ. Следует обратить внимание, что:

1

0

- это полином, то есть каждая группа коэффициентов является множителем (что видно из уравнения (13)) при неизвестном значении х в соответствующей степени;

- пределы интегрирования (“0”-“1”) в значительной степени “комфортны“ для применения, так как сокращают и упрощают вычислительный процесс в целом.

Для решения уравнения (10) в записи РМП операцию интегрирования выполняем путем умножения матрицы [А] на вектор-столбец (23). Но эту же операцию можно расписать и выполнить в виде набора сумм произведений матрицы [ А] на составляющие компоненты глобального вектора столбца (13):

[А] *

0

0

V * а * / к * /

- V * а * / /2/2 - а * Б/2/2

= [А] *

V * а / 2/2 к * а/4/2

0 '

V * а * / / /

- V* а/2/2

V * а/2/2 к * а /4/2

к * а * / / 2/

0

а* Б //

- к * //2/2 - к * а /2/ 0

0

0

00 к * а * / / / - V * а / /

х

х1

х2 ёх + [А] *

3

х

+[ А ] *

4

0

к * а * / / / 0 0

0

к * / / /

- а * Б/2/2

- к * а * //2/2 0

х0

х1

х2 ёх + [А] *

3

х

0 0

.0

х1

х2 ёх + [А] *

х3

х4

0

0

2 ёх =

0

а* Б //

- к * / / 2/ - к * а /2/ 0

2 ёх +

х

4

х

V * а / /

х

х

ёх = 0

Замечания.

1. Интегрирование по частям (16) даст в итоге результат, аналогичный с (13)-(14). Смысл в записи вектора из (13) в виде (16) состоит в том, что реализация матричного произведения в спецпроцессоре, построенном на ЕТП, в матричном варианте позволит одновременно задействовать 5 линеек из СУ. Такое распараллеливание повысит производительность по сравнению с вариантом, когда вычисление будет выполнять один СУ.

2. Если значения коэффициентов V, Б, к и проч. известны заранее или же не требуется повторного вычисления уравнения (когда они могут быть рассчитаны и записаны в виде одного числа каждой из строк (13)), смысл в записи (16) теряется, так как она преобразуется в вектор-столбец с одним числом в строке.

3. Операции по вычислению сумм парных произведений в работе не описываются, поскольку являются известными [2], для чего в общем случае они и предназначены. Поэтому описание этих операций опускаем.

4. Вычисление дифференциалов и интегралов от функций вносит в вычисления численные значения показателей степени неизвестных функций, а это целые числа и они больше единицы. Для устройств с плавающей запятой (ПЗ) это не является проблемой. Для ПЭ со скалярными умножителями, в которых числа представлены в форме фиксированной запятой (ФЗ) , числа в знаменателе могут быть подняты в числитель в виде обратного дробного значения. Эта процедура требует согласованных операций в ПЭ. Важно иметь в виду, что

0

3

4

0

0

3

х

4

0

0

0

х

0

2

х

4

х

основной поток по обработке ССД будет выполняться в СУ с широким распараллеливанием вычислительного процесса.

7. Нелинейный случай

Получив линейное приближение, приступают [1] к решению нелинейного уравнения.

При рассмотрении (6) с целью упрощения принимается

к = V =Б= 1 = и =1.

Тогда, положив в уравнении (15) а(1) = а, вычисляется

а

Замечание. Ниже с целью упрощения записи и чтения индекс (1) над а(1) будет опущен в

местах, где его наличие или отсутствие не вызывает разночтения, а практически до его за-

мены на (2) в конце статьи, когда будет введено новое значение (а(2)).

После вычисления а(1) подставляем его значение в уравнение равновесия (6):

1

| { йи / йх - й2 и / йх2 + и 2}ди * йх =0.

0

При аппроксимации и из выражения (19) определяется функция у, характеризующая состояние равновесия,

1

у (1) = | {а(1)(1 - х)+ а(1) +[1+а(1)(х-х2/2)]}(х-х2/2)йх =

0

=-0,333+0,725 а

Далее, поскольку для состояния равновесия должно выполняться у = 0 , следуя методу Ньютона-Рафсона, применяют рекурсию:

й у = -у

или

где й у /й а =0,725-0,114а.

После подстановки в (22) значения а из а(1) имеют величину

Полученное значение ( й у(1) / й а ) позволяет вычислить

А а

и затем уточнить величину а(1):

а

в итоге вычислить функцию у :

у

Величина у(2) уже показывает, что вычисления близки к решению.

Покажем, как этот процесс можно продолжить с помощью предложенного способа.

8. Вычисления с записью полиномов в виде РМП

Представим интеграл (20) следующим образом:

у(1) = ^ {[ / 1( х )+ / 2 ( х )]( / 3( х )} йх ,

где

/ ( х )=[ а (1-х )+ а ]=2 а - ах =[ / 1 ( х )]*

2а - а 0

2а - а 2а

Вычислим / ^ ( х )=[ / 2 ( х )]2, помня, что

{ т = Шу + Ш2 +1 —2+2+1—5 },

/ 2 ( х )={1+а ( х - х 2 /2)} 2 =(1+а х - а х 2/2) * (1+а х - а х 2 /2)*[ / 2 ( х )]2 *

1 а - а/2 1 а - а/2

* 1 а * 1 а

1 1

1

а 1

- а/2 а 1

0 - а/2 а 1

0 0 - а/2 а 1 0

1

а

- а/2 0

1

- а/2 + а2 - а / 2

2

а /4

х

х

х2 / 2( х)).

х

х

Запишем /3 (х):

"0 1 -1/2" "0"

У3(х) =( х - х 2 /2)* 0 1 1 *[/ 3(х)],

0 -1/2

а также произведение /1 (х) * ^ (х) в виде РМП:

[/1(х)] * [/э(х)]*

2а 0 0

"2а - а 0 " "0 1 -1/2" - а 2а 1 2а

2а - а * 0 1 = 0 -а 2а * -1/2 = - 2а

2а 0

0 0 - а 2а 0 а / 2

х

[/ 2(х)]г * / 3(х)*

0

х

х

2

х

2а 1

(а2 - а) 2а 1

2 (а2 - а) 2а

а

2

1

а2/4 0 0

-а

2

2

(а2 - а) 2а

1

а2 /4 0

а

2

(а2 - а) 2а

1

а2 /4

0 1

-а

(а2 - а) 2а 1

0

1

-1/2 0 0 0 0

Помня, что матрица для интегрирования определена в (12), вычисляем интеграл (27) по частям следующим образом:

Ю (/ г( х))2 » / з(х)Л =[ А ]* [/ г(х )2]Г * [/ 3(/)]ах *

1

2а 1

* [ А]

(а2 - а) 2а 1

- а2 (а2 - а) 2а 1

а2/4 - а2 (а2 - а) 2а 1

0 а2/4 - а2 (а2 - а) 2а 1

0 0 а2/4 - а2 (а2 - а) 2а

0 0 0 а2/4 - а2 (а2

1

- а) 2а 1

0 ' "0 ' "0 " х0

1 1 0 х1

-1/2 (2а -1/2) 1/2 х2

0 йх =[ А ] * ((а2 - а) - а) йх = (а2 - 2а)/3 х3

0 (-а2 - (а2 - а /2)) (-3а2/2+а/2)/4 х 4

0 (а2/4 + а2/2) (3а2 / 4) / 5 х5

0 (а2/8) (-а /8) / 6 х6

0

=(1/2)х2 +1/3(а2-2а)х3+3/8(а2х4)+1/8(а х4)+3/20(а х5)-1/48(ах6) = =1/2+ а (-13 /24)+ а2 (402/480)=

=0,5+0,5417 а +0,8375 а2

2

0

х

х

*

*

1

*

0

и, учитывая ранее полученное произведение (31) функций, вычислим интеграл

Л, / 1(х) * / 3( * ) дх = [ А ] * {[ / 1( х )] * [ / 3( х )] Г } дх *

1

2а Xі 0 х1

*[ А ]‘

0 ' х0 "0

2а х1 0

- 2а * х II а

а/2 х3 - а2/3

0 х4 а/8

=( ах2- а х3 2/3+а х4 1/8)= а 11/24=0,4583 а.

Л ( f 2( х ))2 * / 3( х ) ^ *[ А ] * [ / 2 ( х )]2[ / 3( х )] ^ *

1

* [ А]

"0 ' "0 ' х0

1 0 х1

(2а -1/2) 1/2 х2

(а2 - 2а) дх = (2а -1/2)3 х3

(-3а2/2 + а/2) (а2 - 2а)/4 х 4

(3а2/4) (-3а2/2 + а/2)/5 х5

(-а2/8) 3а2/ 4) /6 х6

0 (-а2 / 8) / 7 х7

0

=1 х 2/2+(1/3)(2а-1/2)х 3+1/4(а2-2а)х 4+1/5(-а23/2+ а/2)х 5 +

+1/6( а2 3/4) х б+1/7(- а2/8) х 7 =

=(1/2-1/6)+ а (2/3-1/2+1/10)+ а 2 (1/4-3/10+3/24)=

=0,333+0,266 а +0,057 а 2.

Суммирование результатов вычисления интегралов

у(1) =0,4583 а =0,333+0,266 а +0,057 а2 =

=0,333+0,725 а(2) +0,057( а(2} )2 = 0,001

дает результат, вычисленный в [1]. Но результат, вычисленный по предлагаемому способу, будет получен в более короткое время. А оно сокращается за счет распараллеливания всего процесса обработки сложных структур данных в СУ взамен последовательного вычисления каждой составляющей подынтегрального выражения (с обращением к памяти за исходными данными и отправкой полученных промежуточных результатов в память) за счет распараллеливания вычислительного процесса (заложенного в алгоритм).

9. Анализ задач

Анализ, проведенный в этом разделе, показывает, что вычисление укрупненных значений функций под интегралом не всегда будет оптимальным. Количество столбцов в матрице типа (17), а, следовательно, и (34), (35), колеблется от одной задачи к другой. Поэтому вполне рациональным может быть некоторое «дробление» подынтегральных функций ти-

3

х

4

0

х

па Л(х), Л (х) , представив /3 (х), например, в виде / (х) = х и /5 (х) = (- х2 / 2), а значение функции /2 (х) следует вычислять с помощью скорректированного набора функций, записанных в виде РМП:

/з(х)=[/ 1(х)]+[/ з(х)]*

1 о 0 1 "0 0 -1/2"

0 1 + 0 0

0 0

Тогда

г'11 =[ А ]» {[/1(х)]»[/1(х)]}ах +[ А ]» {[/1(х)]» [/,(х)]}ах+

+ [ А ] * { /2

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

можно разбить на 2 части, например, так:

(1+2 ах + х2( а2-1))+(- а х3 + х4/4)=

=

а затем каждое из них умножать на /3

[А]* [/ б(х)]г * [/ 3(х)]*[ А ] =

1

2а 1

- а 2а 1

+[ А ]’

0

00 (а2 -1) 0 0

1/ 4 (а2 -1) 0 0

0

1

-1/2

0

0

1

1/2

ёх=.

ёх +

и т. д.

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

10. Дальнейшее вычисление уравнения

Их можно производить так. Дифференцируя выражение

[ С ]=ё ¥ / ёа ) А а ^

* [A] *

“0.333 0.725 0.057" a0 “0.333" a0 “0.725" a0

0.333 0.725 »1 a [ A] * 0.725 »1 a и 0.114 a1

0.333 a2 0.057 a2 0 a2

где [А] - матрица дифференцирования, и для а =-0,564, полученного в (18), имеем

Это позволяет вычислить из (34) Да:

А а =- ¥

После получения матрицы [ С ]-1, обратной от [С] или путем вычисления обратного к 0,661 значения (1 /0,661), вычисление уточненных значений а (уравнения (25) и (26)) выполняются просто, поскольку СУ каждого ПЭ включает сумматор накопительного типа.

Из сравнения результатов, полученных по МГ, с точными значениями функции, приведенными в [1], видно, что он дает близкие приближения.

Замечание. СУ параллельной структуры имеет п -разрядные сумматоры накопительного типа. В них значения результата и переноса в каждом разряде и на каждом такте запоминаются элементами типа триггер. На входе сумматора (до схемы конъюнкции) хранятся новые значения переменных (одно из которых при умножении хранит текущий разряд множимого, а другой - множителя), а при сложении хранятся значения обоих слагаемых.

11. Выводы

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

- записывать их компактно и в виде, удобном для интегрирования или дифференцирования как в аналитическом, так и численном варианте;

- при умножении матриц использовать транспонированную от верхней треугольной матрицы в РМП записи в качестве множимого (запись слева) на вектор-столбец, представляющий первую строку матрицы циркулянта в качестве множителя (запись справа), что дает сокращение в n -раз ( n -разрядность сумматора СУ) количества операций умножения;

- в работе [2] показано, что произведение матрицы на вектор-столбец можно свести к произведению двух векторов, что в m раз (m - порядок результирующего сплайна) сократит количество оборудования и выполняемых арифметических операций типа умножения;

- позволяет от последовательной обработки чисел перейти к параллельной обработке ССД (сплайнов), что резко повышает производительность вычислительного устройства;

- решение дифференциальных и интегральных уравнений и систем выполнять в рамках идеологии единого вычислительного (технологического) потока решения задач МФ, что в целом дополнительно повышает производительность и снижает себестоимость устройств, ориентированных на их решение.

Считаю приятным долгом выразить благодарность д.ф.-м.н. Химичу А.Н. и д.т.н. Резнику А.М. за полезные замечания по тексту статьи, сделанные в процессе её обсуждения.

СПИСОК ЛИТЕРАТУРЫ

1. Коннор Дж. Метод конечных элементов в механике жидкости / Дж. Коннор, К. Бреббиа; пер. с англ. - Л.: Судостроение, 1979. - 264 с.

2. Единый технологический поток в организации вычислений - способ повышения производительности параллельных структур на процессорных элементах транспьютерного типа / Ледянкин Ю.Я.

- Киев, 1989. - 20 с. - (Препр. / АН УССР. Ин-т кибернетики им. В.М. Глушкова; 89-57).

3. Ледянкин Ю.Я. К вопросу преобразования и параллельного ввода граничных условий при решении краевых задач в едином вычислительном потоке / Ю.Я. Ледянкин // Математичні машини і системи. - 2012. - № 1. - С. 28 - 35.

4. Ледянкин Ю.Я. Методы взвешенных невязок, коллокаций, моментов. Способ параллельной реализации в едином вычислительном потоке решения задач математической физики / Ю.Я. Ледянкин // Математичні машини і системи. - 2012. - № 2. - С. 17 - 28.

Стаття надійшла до редакції 05.10.2011