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

2007

Доклады БГУИР

июль-сентябрь

№ 3 (19)

УДК 518.12

РЕШЕНИЕ ДВУМЕРНОГО УРАВНЕНИЯ ЭЛЛИПТИЧЕСКОГО ТИПА МЕТОДОМ БЛОЧНОЙ МАТРИЧНОЙ ПРОГОНКИ

О.И. НАРАНОВИЧ1, А.К. СИНИЦЫН2

1 Барановичский государственный университет Войкова, 21, Барановичи, 225320, Беларусь,

2Белорусский государственный университет информатики и радиоэлектроники П. Бровки, 6, Минск, 220013, Беларусь

Поступила в редакцию 29 июня 2006

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

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

Введение

Многие стационарные задачи математической физики приводят к решению краевой задачи для двумерного эллиптического дифференциального уравнения. Решение этой задачи с использованием метода сеток сводится к системе линейных алгебраических уравнений со слабо заполненной матрицей очень большого порядка. Для решения получаемой СЛАУ наиболее широко используются итерационные методы. Однако с возрастанием возможностей компьютерной техники как по хранению в оперативной памяти больших массивов действительных чисел, так и увеличению их разрядности для решения этой задачи все большее внимание привлекают прямые методы, избавленные от необходимости заботиться о выполнении требований сходимости итераций. Среди последних ранее были предложены метод на основе быстрого преобразования Фурье, метод циклической редукции, методы факторизации, метод матричной прогонки [1-3]. Однако, как отмечается в [2], эти методы имеют те или иные ограничения при применении. В статье предлагается универсальный, устойчивый алгоритм на основе использования метода блочной матричной прогонки [4, 5], в котором реализуется прямой метод Гаусса с выбором главного элемента на упакованной матрице. За счет использования выбора главного элемента снимаются ограничения на устойчивость, присущие указанным выше методам.

Факторизация эллиптического уравнения

Постановка задачи. Запишем краевую задачу для эллиптического дифференциального уравнения вида

д2и д (, ди Л ди ,ди а—- +—I Ь— 1 + с— + ё— + еи = /, (1)

дг дг V дг) дг дг

коэффициенты а,Ь,с,ё,е,/ являются функциями от г и г.

Граничные условия для прямоугольной области П = {0 < г < Ь, Ьх < г < Ь2} задаются следующим образом:

( 0,Ь ди , Р0,Ь Л 0,Ь ( 1,2 ди , „1,2 Л 1,2

\а'— + в и I =у ,!«'— + в и I = У' . (2)

V дг )0,Ь V дг )ь1,ь2

Приведение к системе обыкновенных ДУ. Выберем на интервале { Ьх < г < Ь2} равномерную сетку а>ьг ={г = ]ЬТ, кг = (Ь2 —Ь1)/ш, ]=0,..., ш} и обозначим и ={и(г, г0),..., и(г, гш)} = {%..., иш}. Для каждого г у аппроксимируем (1) центральной конечно-разностной схемой второго порядка точности:

Л.

d Uj bj+Ц2(и] +! - и} ) - bj-1/2(uj - Uj-l) du

+ J™^ ^-J' J-JZ1L + c.—L + d,Uj+1 U-1 + e,u, = f, j = 0,...,m. (3)

j dz2 hr2 j dz j 2hr J J J ^

1 2

При a^ Ф 0 для значений j=0 и j=m в (3) необходимо производить аппроксимацию производных по r с учетом второго уравнения в граничных условиях (2). При этом в правых

частях появятся пропорциональные у1'2 дополнительные члены fom = J0 m + py1'2 . Если одно

11 2 2 из граничных условий является условием первого рода (а = 0, в = 1 или а = 0, в = 1), то

1 2

значения и = у или um = у заданы и количество уравнений (3) уменьшается на единицу (или на два, если оба), и при этом дополнительный член появится при j=1 или j=m—1 соответственно

f1 = f1 - Y1(b1/2 / hr2 - d1/2hr ) или fm-1 = fm-1 - г\Ьт-1/2 /hr + dm-1 /2hr ) .

После некоторых преобразований система ОДУ (3) может быть записана в векторном

виде:

E(z) ^ + Q(z) ^ + G(z)u = f(z). (4)

dz2 dz

Краевые условия при z=0, L для (4) получаются из первого уравнения граничных условий (2).

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

Конечно-разностная схема. Выберем в соответствии [4, 5] на интервале {0 < z < L} равномерную сетку a>hz = {zi = (i - 1)hz, hz = L/n, i = 1,...,n +1} . Для расчетов используем конечно-разностную центральную схему второго порядка точности:

E,Ui-1 -2U + и+ + Q,^^ + GA = f, i = 2,...,n . (5)

hi 2hz

После приведения подобных членов в (5) получим систему линейных алгебраических уравнений:

(E, -fy)U,-1 + ( + h^Gi)u, Ei + ^Qi ]U+1 = hff,,

i = 2,...,n . (6)

Используем аппроксимацию второго порядка точности для первого уравнения граничных условий (2):

а°(-3и1,- + 4ы2 ] - и^,) + 1Ь2виХ] = 2ку°, аЬ (3ил+1>, - 4^,, + и^,) + 2кгвЧ+1,, = 2Ь2уь . (7)

Представим СЛАУ для конечно-разностной схемы (6), (7) в виде Ах = й .

Матрица А размерности к*к, к=(п+1)(т+1), имеет блочно-ленточную структуру со слабо заполненными матрицами блоков [4].

г 12

Вектор правых частей й имеет вид (для граничных условий с а ' Ф 0):

й = {/,...,Л /1,0 + р/'fl,1,...,/1,т-Ъ/т + рЛ- , /п+1,0 + /п+1,1,..., Уп+1, т-1, /п+1, т +РУ2,

Методика решения Для решения системы линейных уравнений с блочноленточной матрицей (6), (7) была разработана экономичная модификация метода Гаусса с выбором главного элемента — метод блочной матричной прогонки [4]. Идея алгоритма заключается в реализации метода Гаусса на упакованном массиве, в который помещаются ненулевые блочные элементы матриц. Алгоритм метода Гаусса с выбором главного элемента обеспечивает устойчивость методов (3) и (5), даже если не выполняется условие преобладания матричного диагонального элемента, необходимое для реализации классического метода прогонки. Так как в рассматриваемом случае матрицы О и Q слабо заполнены, то наряду с упаковкой в ленточные массивы был разработан алгоритм упаковки только ненулевых элементов в массив из одно-связных динамических стеков.

Результаты тестовых расчетов

В качестве тестовых испытаний разработанных алгоритмов решены две задачи.

1. Краевая задача для волнового уравнения Гельмгольца в цилиндрической системе координат:

1 д2и 1 д ,1 ди ч ТТг2 и

■ (_ ) + ж2 — = 0.

Р дг2 Ь2 др р др

Р

(8)

Здесь искомая функция м(г,р) определена на области 0 = {0< г<Ь, 0<р = г/Ь < 1},

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

ди (г,1)

на стенке и оси волновода вида

др

= 0, и (г,0) = 0 уравнение (8) имеет точное решение

и = р-Л(хр)[с1 со$(к2г) + С2ът(к22 )Ь х = ^01 / Ь , к1 = ж 2 -Х2, % < Ьж, У01 — корень функции Бесселя Jo (х) .

Для установления произвольных коэффициентов С1, С2 зададим следующие условия на концах волновода:

и (0,р) = р^(*р), [аЬ ди + вЬи ] = / .

В соответствии с приведенной выше методикой уравнения (3) приобретают вид

(Ч = 0):

1 д 2и.

1

21.2

р, дг2 Ь2 к

р 1

7 - 2

■и,-1 -

1

1

р 1 Р 1

, -- ,+-V 2 2 У

1

и

7 ' "V+1

р 1

7+2

Ж2 и = 0, 7 = 1,...,т -1. (9)

р,

Используя граничное условие ■ди + р2ь

Vдр

= у и аппроксимацию второго порядка

_д_( 1 ды_ ^ _ J_

др^рдр) 2и;

8um-1 - Um-2 -( 7 +

в2 ( -2h2))Um - у2 ( -2h2)

+ °(h2)

получим уравнение для j=m 1

d2um

dz2 b2h2

-Um-2 + 8Um_1 - (7 + ft2 (6hr - 2h^ ) - 2b2h2W2 )Um _ (1 - 3/ hr ) Д

Запишем систему (9), (10) в векторном виде аналогично (4):

E-

dz2

- + GU _ f .

(11)

Для нее задается второе граничное условие (7). Вектор / содержит ненулевые компоненты, равные правым частям (10). Матрица Е содержит только ненулевые диагональные элементы, равные 1/р., . = 1,...,т . Матрица О размерности т*т имеет следующие ненулевые коэффициенты:

1

(

gu _-2hv

1

1

Л

r

gj, j-1 _

1

2hr2b2 Pj-1/2

р1/2 р1+1/2

1 _

'' g j, j _

W2 _ 1 1 P1 2h2b2 P1+1/2

2hr2b2

1

1

pj-1/2 pj+1/2

W2 _ 1 1

' g j, j + 1 _ ~ ,2,2 Pj 2h2b2 Pj+1/2

(12)

. = _ -1 = 8 7 + в2 (6^ - 2И2г ) - 2И?Ь2Ж2

И = 2,...,т 1 > .?т,т—2 = _ , 2,2 , ^т,т—1 = _ , 9,9 , £т,т = „,9,9 .

2н;ь2 2н;ь2 2и;ь2

Структура упакованной блочно-ленточной матрицы системы (6), (7) для рассматриваемого случая представлена на рис. 1. Квадратная матрица А размерностью кхк, где к=(п+1)т, упаковывается в ленточную матрицу размерности к х 4т . Запасной столбец необходим для реализации метода Гаусса с выбором главного элемента.

-3а 0 + 2hze0 4а0 -а0 0

E hz2G - 2 E E 0

E hz2G - 2 E E 0

E hz2G - 2 E E 0

L а -4aL 3aL + 2hzeL 0

Рис. 1. Упаковка матрицы А в ленточный массив

Так как матрица G трехдиагональная, а матрица Е диагональная, то, как видно из рис. 1, в каждой строке матрицы А имеется не более пяти ненулевых элементов. Поэтому для ее упаковки целесообразно использовать массив из односвязных динамических стеков, представленный на рис. 2. В результате выполнения прямого хода метода Гаусса с выбором главного элемента на таком массиве количество элементов в стеках увеличивается, однако не превосходит 4m. В результате эффективность метода при такой упаковке должна быть выше, чем при упаковке в ленточный массив.

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

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

1..Ш

т+1

ш+2

к-т

к-т+1..к

- 3а 0 + 2к2 в „ 0 4а 0 — а пД

"ш+1 1 "ш+1 т+1 "т+1 ш+2 "т+1 2т+1 пД

аш+2 2 "т+2 т+1 "т+2 ш+2

ак-ш к-2т-1 ак-ш к-т-2 ак-ш к-ш-1

Ь Л Ь — 4а Ъа1 + 2к2вЬ т!

"к-т к

Рис. 2. Упаковка матрицы А в массив из стеков

Пример решения краевой задачи (8) для ¿=30, W=1, Ь=3,5, 0=0, в=1, /"=0 в виде графиков и(г,гI = 0,1 /), j = 1,...,10, представлен на рис. 3.

у г^тггггтр V-

' / ' /V \\\

..Шхг.

""V

Ж

жу-У'Ж"

: ^О? ! ! !

¡м

23456789 10 11 12 13 14 15 13 17 18 1 20 21 22 23 24 25 23 27 2 29 30

7.

Рис. 3. Вариант расчета волнового уравнения

В таблице приведены времена счета задачи (8) по двум реализациям алгоритма при различных значениях параметров сетки на компьютере 8етргоп 3000. Алгоритм с упаковкой в массив из стеков для рассматриваемой задачи работает в 15-30 раз быстрее.

Временные характеристики решения краевой задачи для волнового уравнения Гельмгольца

п т Время, с Время, с

Упаковка в массив из стеков Упаковка в ленточную матрицу

1000 10 0,1 1,5

1000 20 0,5 10,7

1000 30 1,3 35,1

1000 40 2,7 82,3

1000 50 5,1 159,5

100 100 5,0 28,3

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

д2и д2и . . ди

дг2 дг2 дг

- У(г) =0, (13)

коэффициент у(г) = 4у0г(1 - г) описывает профиль скорости жидкости между пластинами.

Граничные условия для прямоугольной области О = {0 < г < Ь, 0 < г < 1}, где Ь выбирается достаточно большим, задаются следующим образом:

и(г,0) = и(г,1) = 1, и(0,г) = 0, и(Ь,г) = 1. (14)

Особенность этой задачи состоит в том, что применяемая при рассматриваемой методике центральная конечно-разностная схема (6) при У0 к2 > 1 оказывается не монотонной в том смысле, что диагональный элемент матрицы А не преобладающий, хотя имеет место аппроксимация. Поэтому итерационные алгоритмы ее решения в этом случае оказываются неустойчивыми. Как показали тестовые расчеты, предлагаемый метод является устойчивым и в этом случае. На рис. 4 для у0=50, Ь=10, кг=0,01, кг=0,05 представлены графики функции и(г, г) для различных значений г=/Ь/10, /=1,2, ..., 9 (кривые 1, 2, 3, ...). Расчет при Нг=0,01, когда условие монотонности не выполняется, дает практически тот же результат с учетом погрешности.

Проведено сравнение эффективности рассматриваемого метода с итерационным методом Зейделя [5] при оптимальном выборе параметра релаксации, которое показало, что время расчета по предлагаемому методу в 10-40 раз (в зависимости от т, п) меньше, чем по методу Зейделя.

U

3

2

............................................... 1

о

0,25

0,5 Z

0,75

1

Рис. 4. Вариант расчета уравнения теплопроводности

THE DECISION OF THE TWO-DIMENTIONAL ELLIPTIC TYPE EQUATION

A BLOK MATRIX RUN METHOD

O.I. NARANOVICH, A.K. SINITSYN

Abstract

The effective direct method of the decision finite-difference circuit approximating a boundary problem for the two-dimentional elliptic type differential equation is offered. Economic updating of Gauss method with a choice of the main element — a blok matrix run method is developed. The idea of algorithm consists in realization of Gauss method on the packed files in which zero elements of matrixes are located not. Results of comparative calculations of realization of the offered algorithms with use of packing of a matrix of system in a file from stacks and in a sheet file are resulted by the example of the decision of wave equation in cylindrical system of coordinates.

Литература

1. Калиткин Н.Н. Численные методы. М., 1978.

2. МарчукГ.И. Методы вычислительной математики. М., 1980.

3. Крылов В.И., Бобков В.В., Монастырный П.И. Начала теории вычислительных методов. Уравнения в частных производных. Минск, 1986.

4. Синицын А.К. // Докл. БГУИР. 2005. Т. 5, № 1. С. 46-49.

5. Синицын А.К. Современные информационные технологии. Проекционно-сеточные методы решения уравнений математической физики. Конспект лекций для аспирантов и магистрантов БГУИР. Минск, 2004.