О КОНЕЧНЫХ МЕТОДАХ РЕШЕНИЯ УРАВНЕНИЯ ПУАССОНА НА ПРЯМОУГОЛЬНИКЕ С КРАЕВЫМ УСЛОВИЕМ ДИРИХЛЕ Текст научной статьи по специальности «Математика»

МАТЕМАТИКА

УДК 519.6:517.958

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

С КРАЕВЫМ УСЛОВИЕМ ДИРИХЛЕ

(Московский государственный технический университет им. Н.Э. Баумана); д-р физ.-мат. наук, проф. К.А. ВОЛОСОВ (Российский университет транспорта, Москва); канд. физ.-мат. наук А.К. ВОЛОСОВА (ООО «Трамплин», Москва);

канд. физ.-мат. наук, доц. Д.Ф. ПАСТУХОВ, канд. физ.-мат. наук, доц. Ю.Ф. ПАСТУХОВ (Полоцкий государственный университет)

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

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

Введение. Матрицы и матричные уравнения специального типа применяются во многих разделах прикладной математики. В квантовой механике динамика частиц со спином определяется матрицами кватернионов (полукватернионов) [1; 2]. Уравнение Пуассона на прямоугольнике (параллелепипеде) можно решить методом прогонки [3; 4; 5; 6; 10; 12; 13; 19]. Алгебраический метод прогонки совместно с формулой простой итерации [5] является приближенным методом, так как число итераций не ограничено, но имея формулу аппроксимации уравнения Пуассона с шестым порядком погрешности, можно значительно снизить погрешность и время вычислений [5]. В данной работе рассмотрен метод прогонки в матричной форме для численного решения уравнения Пуассона за конечное число арифметических операций. Идея работы частично основана на идее статьи [10], а также на идее модификации краевых столбцов и строк в матрице правой части уравнения Пуассона с шестым порядком аппроксимации [5]. Однако в работе [10] и в данной работе возможно обобщение задачи, то есть решение уравнения Пуассона на прямоугольной сетке щ х Щ. с квадратными ячейками 1\ = И2 = И, но матрицы коэффициентов А, В при этом

по-прежнему квадратные щ х щ . Этот эффект мы назвали эффектом прямоугольной шахты, в которой перемещается квадратная кабина лифта (квадратные матрицы А, В щ х щ ) в направлении щ , минимальное перемещение И = И = И (перемещение поперек шахты не разрешается). Возможны ситуации щ < п2, щ > «2 - длина шахты как больше размера кабины, так и меньше. Получены достаточные условия корректности предложенного алгоритма, теоремы 1, 2, 3. Метод можно использовать в численных задачах математической физики [15; 16; 17], а также в двумерных задачах гидродинамики, система уравнений которых содержит уравнение Пуассона от функции тока, где правая часть - функция вихря Ду = —ю [20].

Постановка задачи. Рассмотрим краевую задачу Дирихле для уравнения Пуассона на прямоугольнике с неизвестной функцией и(х, у) с неоднородными граничными условиями у е[с, ё], х е [а, Ь] (^ (у),

Ф2М, Фэ(х), фДх), х е [а,Ь], у е[с,й~]):

Н.К. ВОЛОСОВА

Дифференциальной задаче (1) сопоставим разностную задачу (2), в которой два параметра Н\, -шаги равномерной сетки по координатным осям х, у соответственно [4, с. 102]:

Ц»,»2 иЬфг = (V ^2 )6 : \ = Хп = а + ^П = 1 П1 - 1 У»2 = Ут = С + тК т = 1 П2 -1

V иИ,н\Т =Ф»,»2,(Х*,' У»2 )бГ*,»2 : \ = X = а + ^ П = 0' ^ >»2 = Ут = с + т»2. т = 0, Щ,» = -, »2 =

Ь-а , ( -с (2)

где и^ - неизвестная сеточная функция, решение задачи (2).

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

и симметрии шаблона. Можно получить шестой порядок аппроксимации задачи (1) задачей (2) на девятиточечном шаблоне с равными шагами сетки [5; 6] \ = \ = И. Краевой разностный оператор /АА = 1

в задаче Дирихле (2) - единичный. В работе [5, с. 73, формула (32] нами получено разностное уравнение (2) с шестым порядком погрешности для девятиточечного шаблона с центральным узлом иш п:

1 Г-10 2/ ч \( Л

I ^ ит," + 3 (ит-1,и + ит+1,п + ит,п-1 + ит,п+1 ) + ^ (ит-1,и-1 + ит+1,п-1 + ит-1,п+1 + ит+1,и+1 ) I "

= /т,п + (/х + /> ) + И ^ (/<4)+ />(4)) + ^ /4>У )+ О (» ) .

(3)

п = 1, п -1, т = 1, п -1.

Покажем, что разностное уравнение (3) можно получить из рекуррентного разностного матричного уравнения второго порядка. Определим матрицы А, В:

А = К А =

10 2 0

3 3

2 10 2 0....

3 3 3

0 2 3 10 3 20... 3

0 0... 2 3 10 3

0 0... 0 2 3

т=п т= 1, п -1, п = 1, п

2 3 10

з.

10

3

—, т = п +1V т = п -1 3

0, т > п + 2 V т < п - 2

5 = Ь =

т,п

2 1 0

з 6

1 2 1 0.... 0

6 3 6

0 1 2 10... 0

. 6 3 6

0 0... 1

6

0 0... 0

ь„„ =

—, т = п; т = 1, п, -1, п = 1, п, -1

3= = 1 =1 ,

1, т = п +1V т = п -1 6

0, т > п + 2 V т < п - 2

(4)

где А, В - квадратные симметрические трехдиагональные матрицы Теплица ранга щ -1. Перепишем уравнение (3) в эквивалентном виде:

-10 2/

3 ит,п + ^ (ит-1,п + ит+1,п + ит,п-1 + ит,п+1 (ит-1,п-1 + ит+1,п-1 + ит-1,п+1 + ит+1,п+1

1 ) + 1 (ит-1,п-1 '

) =

р = /И2 + — (/ + / ) + И6 (—(/(4)+ /(4)) + — /(4)1 + О (А8), т = 1,п - 1,п = 1,п -1. (5)

т,п ¿т,п УУ ' 1 36() V У ' дд-/хх>> I V / 2

а = <

т.п

Если ыт_ [, ит, ит+1 вектор-строки матрицы решения ранга (щ — 1) х (щ — 1), то уравнение (5) эквивалентно уравнению (6):

щ —1 щ — 1 щ — 1

BuTm-1 + AuTm + BuTm+l = FTm ö £bsy km-1 + £ k+ £k,m+i =FТ, s = 2,щ - 2,m = 2, щ - 2 ö

k=1 k=1 k =1

щ-1 Щ-1 Щ-1

ö Z bs,kUm-1,k + £ as,kUm,k + £ bs,kUm+1,k =Fm,s ö bs,s-1Um-1,s-1 + bs,sUm-1,s + bs,s+1Um-1,s+1 +

+0,*-1Um,s-1 + as,Um,s + +1 + bs,s-1Um+1,s-1 + + +1Um+1,s+1 = Fm,s , S = 2 «1 " 2 m = 2 «2 " 2 • (6)

Учитывая (4), при s = w из (6) получим

12 12 2 1 2 1 -тг

TUm-1,n-1 + T Um-1,n + T Um-1,n+1 + T Um,n-1 T" Um,n + T Um,n+1 + T Um+1,n-1 + T Um + 1,n + T Um + 1,n+1 = Fm,n •

6 36 3336 36

Группируя слагаемые с одинаковыми коэффициентами в последнем уравнении, имеем

10 2/ ч Ь ч

~ Um,n + (Um-1,n + Um+1,n + Um,n-1 + Um,n+1 (Um-1,n-1 + Um-1,n+1 + Um+1,n-1 + Um+1,n+1 ) = Fm,n ■.

2 tí

П = 2, «1 - 2, m = 2, n2 - 2 , ГДе Fm,n = fm,nh

Fm,n = /„.h + (/- + fy )+ h Í ^ fy<4)) + ^ /- )+ 0 (h8 )

í=X , У=Ум

Последнее уравнение эквивалентно уравнению (5) на внутренних узлах. Уравнение But , + Aut + BuT^, = FT при m = 1. n -1 имеет смысл, если положить AU + B«T = FT, Bmt 9+ Amt , = FT ,.

m-1 m m+1 m ^ 2 * 12 1? n2-2 n2-1 n2 -1

В результате запишем систему матричных уравнений, в которой матрицы А, В определены формулой (4):

AUT + BUT = f Т

BUt ,+ AUt + BUt , = FT,m = 2,n - 2 . (7)

m-1 m m+1 m ' ' 2 v '

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

и „ =ф , и =ф , ип =ф , и =ф , т = 0, щ, щ = 0, щ . (8)

С учетом системы уравнений (7), краевых условий (8) формула (5) имеет теперь смысл при всех

т = 1, щ — 1, и = 1, щ — 1. Рассмотрим первое уравнение системы (7) в угловом узле равномерной сетки с индексами (т = 5 = 1, п = 1):

-10 2 2 1 _-10 2 / ч 1 _ —

a1,1U1,1 + a1,2U1,2 + b1,1U2,1 + b1.2U2,2 = F1,1 U1>1 +3 U1>2 ^3 U2,1 ^ U2,2 = ^ U1,1 ^ (U2,1 + U1,2 U2,2 = F1,1 . (9)

В правой части уравнения (9) стоит черта, так как правая часть (9) отличается от уравнения (5):

1 и1,1 + 2 (и2,1 + и1,2 + и1,0 + и0,1) + 1 (и2,2 + и0,2 + и2,0 + и0,0 ) = ^1,1 + ^ (и1,0 + и0,1 ) + ^ (и0,2 + и2,0 + и0,0 ) = ^1,1 ^

_ 2 1

° ^¡,1 = Р — 2 (и1,0 + и0,1 ) — 1 (и0,2 + и2,0 + и0,0 ) .

k=1

k =1

k=1

Аналогично, для трех остальных узлов имеем -10 2/ ч 1/ ч „

— и 1 ,«,- 1 + 3 (и2,п- 1 + и 1 ,п-2 + и 1 + и0,п- 1 ) + ~ (и2,п-2 + "0,^-2 + "г.п, + "0,« ) = Р1 ,«- 1

р1 ,«,- 1 = рп- 1 - | (и 1 + и0,„1- 1 )- ^ (и0,п-2 + и2,п + и0,п) -10 2 1

— ип2-1,1 + ип2-2,1 + "п, -1,2 + и„!-1,0 + ^,1 )^(и»2-2,2 + ипг,2 + и„! - 2,0 + и ,0 ) = Р„2-1,1

р , ,= р ,,-2(и ,„ + и ,)- — (и .+ и .„ + и „)

п2 -1,1 п2 -1,1 п2 -1,0 п2,1 п2,2 п2 -2,0 п2,0

ч 1/ ч

3 ит^-1 + 3 \ит^-1 + "«2-2 + "»2-1,»2 + "«2,«,-1 )+ 6 "-2,^-2 + "«2«-2 + "»2^ + ) ^^-1 ,

р , ,= р , ,-2 (и , + и ,)-1 (и 0 + и _ + и )

«2 -1 -1 -1,«1 -1 ^ V «2 — 1,п1 «2 -1 / ^ V «2 -2 п2- 2,«! «2 /

("2,1 + "1,2 + "1,0 + "0,1 ) + 1 ("2,2 + "0,2 + "2,0 + "0,0 ) = р1,1

-10 2

— "1,1 + "("2,1 + "1,2 + "1,0 + "0,1 —("2,2 + "0,2 + "2,0 + И0,0 ) р1,1 = р1,1 - | ("1,0 + "0,1 ) - ^ ("0,2 + "2,0 + "0,0 )

где р„,„ = /т,2 + И4 / + /У ) + И6 Й14)+/У(4)) + ¿^ХУ ^ О (И8 )

Рассмотрим второе уравнение системы (7) для узла с координатами (т = 2, п = 5 = 1), из (6) имеем: Ь и , ,+ Ь и , + Ь л . . + а и ,+ а и + а м ^ + Ь м^ ,+

5,5-1 т-1,5-1 5,5 т-1,5 5,5+1 ТО-1,5+1 5,5-1 т,5-1 5,5 т,5 5,5+1 ТО,5+1 5,5-1 ТО+1,5-1

+Ь5,5ит+1,5 + Ь5,5+1"т+1,5+1 = рт,5 » Ь1,1и1,1 + Ь1,2и1,2 + а1,1и2,1 + а1,2и2,2 + Ь1,1и3,1 + Ь1,2и3,2 = р2,1 «

2 1 10 2 2 1 — 10 2/ ч 1/ ч

-и1,1 + -и1,2 ~ и2,1 +-и2,2 + Т и3,1 + Т и3,2 = р2,1 « ~ и2,1 + т(и1,1 + и2,2 + и3,1 )^(и1,2 + и3,2 )

363336 3 3 6

Добавляя тождественно несколько слагаемых к левой и правой части последнего уравнения до полноты формулы (5), получим

10 2, Ч Ь Ч — 2 1/ ч „

-У и2,1 + 3 (и1,1 + и2,2 + и3,1 + и2,0 )^(и1,0 + и3,2 + и1,2 + и3,0 ) = р2,1 + "и2,0 + ^ (и1,0 + и3,0 ) = р2,1 . (11)

Аналогично формуле (11) для узла (т = 2, п = 5 = 1) для всех граничных (неугловых) узлов (т = 1, п = 1, ..., п1 - 1), (т = п2 - 1, п = 1, ..., п1 - 1), (т = 1, ..., п2 - 1, п = 1), (т = 1, ..., п2 - 1, п = п1 - 1) имеем

10 2/ ч 1/ ч ^ --г

--Т и1п +-:(и1,«-1 + и2,п + и1,«+1 + и0п )^(и2,„-1 + и2,«+1 + и0,«-1 + и0,«+1 )= р1,„, « = ^ «1 - 2

3 3 6

--2 1/ ч -

р1,„ = р1,п--и0п -7(и0,„-1 + и0,«+1),п = 2,«1 - 2

36

(и , ,+ и 0 + и , ,+ и )+— (и 0 ,+ и 0 ,+ и ,+ и ,)= р , ,п = 2,п, -2

п -1,п-1 п -2,п п -1,п+1 п ,п п -2,п-1 п -2,п+1 п ,п-1 п ,п+1 п -1,п 1

2 _

ип -1,п +^"(и«2-1,п-1 + ипг-2,п + ип -1,«+1 + ипг ,п-2,п-1 1 "щ-2,п+1 1 "щ ,п-1 1 "щ ,п+1 )= р п2-1,

--2 1/ ч -

Рщ-1,« = -1,« - ^ "«2,« ^ (и«2 ,«-1 + и«2 ,«+1 > « = 2, «1 - 2

10 ^ ч 1/ ч ^ ----(12)

— ит,1 ^(ит-1,1 + ит,2 + "т+1,1 + "т,0 )^(ит-1,2 + "т+1,2 + "т-1,0 + "т+1,0 )= рт,1, т = 2, «2 - 2

3 3 6

--2 1/ ч -

рт,1 = рт,1 ~ ит,0 "К-1,0 + "т+1,0 ), т = 2, «2 - 2

' т,1 = ' т,1 , "т,0 , (ит-1,0 + "т+1,0 ), т = 2, «2

36

10 2/ \ 1/ \ ^ ---

--и ,Н—(и , ,+ и 0 + и , ,+ и )+—(и , 0 + и , 0 + и , + и , )= г ,, т = 2, п - 2 ;

т,п -1 т-1,п -1 т,п -2 т+1,п -1 т,п т-1,п -2 т+1,п -2 т-1,п т+1,п т,п -1 2

-- 2 1 / \ _-

рт,щ-1 = рт,щ-1 "т,« ^~(ит-1,« + "т+1,« т = 2, «2 - 2

рт,„ = рт,„, Ут £ 2, «2 - 2,« £ 2, «1 - 2

р = / Й2 + (/ + / ) + Й6 Г—1— (/(4) + /(4)) + —/(4)1-^ О (» ^

где

Таким образом, запись уравнений системы (7) с матрицами А, В, определенных формулой (4) во внутренних узлах (т = 2, ..., п2 - 2, п = 2, ..., п\ - 2), не отличается от формулы (5) - аппроксимации уравнения Пуассона на прямоугольнике с шестым порядком погрешности. В граничных узлах следует использовать формулы (10) - в четырех угловых граничных узлах или формулы (12) - в граничных узлах на четырех отрезках.

Будем искать решение рекуррентно заданной системы матричных уравнений (7) в виде

К = ХУ™^ + Ут , т = 1 «2 - 2 • (13)

К-1 =У%-1' т = «2 -1

Как следует из теории размерностей [7] в формуле (7) ит, у т - векторы ранга « -1, а коэффициенты прогонки Хт' т = 1,п2 -2 представляют собой квадратные матрицы того же ранга п1 -1 • В (\3) иТт обозначает вектор-столбец - транспонированную строку ии с номером т матрицы решения (п2 - \)(п\ - \). Из первого уравнения системы (7) имеем

иТ = -А~1БиТ2 + А'1 В? о \ = -А_1Б,у1 = А'1 ВТ , (\4)

где ВТ обозначает вектор-столбец - транспонированную строку с номером т = 1 матрицы правой части Е уравнения Пуассона.

Как следует из (13) иК = Хт_К + Ут_1, преобразуем второе уравнение системы (\3):

БиТ ,+ АиТ + Би\, = ¥Т о Б(Х К + у ,) + АиТ + Би\, = ¥Т,(БХ ,+ А)иТ = -Би\, + ¥Т -Бу ,о

т-1 т т+1 т у т-1 т т-1 / т т+1 т > \ т-1 / т т+1 т т-1

ит = -(БХт-1 + А)-1 + + (БХт-1 + А)-1 (1 -БУт-1) о

Хт = - (БХт-1 + А)-1 Б'у„ = (БХ^ + А)-1 (1 -Бут-1),т = 2'«2 - 2 • (15)

Используя последнее уравнение системы (5), форму решения (13), найдем последнюю строку мат-

Т

рицы решения и -1:

Би«2-2 + А<-1 = В«2-1' и«2-2 = Х«-2и«2-1 + У«-2 0 Б (Х«2-2и«2-1 + У«-2 ) + Аи«2-1 = <-1 о

(БХ«-2 + А)и\-1 = - Бу«-2 о иТ-1 =(БХ«-2 + А)-1 (^ - Бу«-2) • (16)

В силу второго уравнения (13) получим

У «2-1 = <-1 • (\7)

Формулы (14), (15) назовем матричными формулами прогонки вперед, а формулы (17), (13) матричными формулами прогонки назад:

иТт =Хтит+1 +Ут ' т = «2 - 2,1 • (18)

\ По формуле ЕТт« = + / + 4) + А6 ^/>+ /?>) + ±у;«;. ] + 0(А8)

Опишем построенный нами алгоритм решения уравнения Пуассона за конечное число арифметических операций матричным методом прогонки (в добрых традициях мехмата МГУ [3, с. 587-588])^

~А I „ „ \ 1 / „1л\ 1 „1л\ I „/.о\

вычис-

, . У = Ут

лить правую часть уравнения Пуассона во всех внутренних узлах равномерной сетки прямоугольника (т = \, ..., П2 - \; п = \, ..., п\ - 1). Краевые условия задачи (1) определяются формулой (8) и „ = ф, , и = ф- , к = <ф , и = ф, , т = 0,«, п = 0,п •

2^ Поправить (модифицировать) правые части системы уравнений (7) по формулам (10), (12) в узлах прямоугольного контура, соседнего с граничным контуром, то есть найти В по величинам В пункта 1.

3^ Найти матричные коэффициенты прогонки вперед по формулам (14), (15) т = 1, п2 - 2 • 4^ Найти вектор-строку иТп -1 по формуле (16).

5^ Найти остальные строки матрицы - решения итт, т = п2 - 2,1 по формулам (18).

Решим пример на классе элементарных функций, методы решения аналогичных примеров описаны в учебнике по уравнениям математической физики [9]:

мхс + = з1п(2х) $т(3у), (х, у) е (0,2л) х (0,3л). и(0, у) = и(2л, у) = зт(у), у е [0,3л] и(х, 0) = и(х, 3л) = $т(х), х е [0,2л]

(19)

Проведем линейную редукцию примера (19), положив и( х, у) = щ (х, у) + и2 (х, у) + и3 (х, у), где ихх (х, у) + и!уу (х, у) = зш(2х) 8т(3у), (х, у) е (0,2л) х (0,3л).

1) <и,(0, у) = и, (2л, у) = 0, у е [0,3л]

щ (х, 0) = щ (х, 3л) = 0, х е [0,2л] Решение последней системы ищем в виде щ (х, у) = А зт(2х) зт(3у). При этом краевые условия -второе и третье уравнения системы 1), выполняются автоматически. Подставим щ (х, у) в первое уравнение системы:

Цхх (х, у) + (х, у) = 8ш(2л) 8ш(3у) о - А (4 + 9) эт(2х) эт(3 у) = эт(2х) 5т(3 у) о А = - -1,

тогда ц,( х,у) = - ¿а*2 х)8Ш(3 у).

ихх (х, у) + и2уу (х, у) = 0, (х, у) е (0,2л) х (0,3л).

2) [ и2 (0, у) = и2 (2л, у) = 8т(у), у е [0,3л]

и2 (х, 0) = и2 (х, 3л) = 0, х е [0,2л] Ищем решение 2) в виде щ (х, у) = /(х) зт(у), автоматически удовлетворяющее краевому условию - третьему уравнению в системе 2). Подставим и2 (х, у) в первое и второе уравнения системы 2):

о

/ ■( х) - / (х) = 0

| ( /( х) - / ( х) ) 81П( у) = 0

|и2 (0, у) = щ (2л, у) = Б1П(у) о /(0) Б1П(у) = /(2л) Б1п(у) = Б1П(у) [/(0) = /(2л) = 1.

Общее решение однородного обыкновенного дифференциального уравнения в краевой задаче ищем

1 - сЬ(2л)

в виде /(х) = А бЬ(х) + В сЬ(х), /(0) = 1 о А • 0 + В • 1 = 1, В = 1; /(2 л) = 1 о А • 8Ь(2 л) + сЬ(2 л) = 1, А = ■

8Ь(2л)

где вЬ( х) = ■

е + е

, сЬ(х) =-Ух е (-да, +да).

Тогда щ (х, у) = / (х) бш( у) = (А зЬ(х) + В сЬ(х)) зт(у) =

1- сЬ(2л) зЬ(2л)

зЬ(х) + сЬ( х)

зт(у) .

и3хх (х, у) + и3уу (х, у) = 0, (х, у) е (0,2л) х (0,3л). 3) < и3(0, у) = и3(2л, у) = 0, у е [0,3л]

и2 (х, 0) = и2 (х, 3л) = 5ш(х), х е [0,2л] Ищем решение 3) в виде щ (х, у) = / (у)з1п( х), автоматически удовлетворяющее краевому условию - второму уравнению в системе 3). Подставим щ (х, у) в первое и третье уравнения системы 3):

о

/(у) - / (у) = 0 .

| ( /(у) - / (у)) ®1п( х) = 0

[щ (х, 0) = щ (х, 3л) = БШ(х) о /(0) БШ(х) = /(3л) БШ(х) = БШ(х) [/(0) = /(3л) = 1

Но последние две краевые задачи отличаются друг от друга заменой переменных х о у, поэтому и решения задач 2) и 3) отличаются только заменой переменных х о у и краевыми условиями / (2л) = 1 о /(3л) = 1:

3 (х, у) = / (у) зт(х) = (А зЬ( у) + В сЬ( у)) зт(х) =

1- сЬ(3л) зЬ(3л)

8Ь( у) + сЬ( у)

зт(х) .

ех -е

и

Тогда решение исходного примера (19) имеет вид

u( x, y) = —^1sm(2x)sm(3 y) +

1 - ch(2rc) sh(2^)

sh( x) + ch(x)

sin( y) +

1 - ch(3rc) sh(3:rc)

sh( y) + ch( y)

sin(x) . (20)

Программа, написанная нами на языке FORTRAN [20] с использованием алгоритма (8), (10), (12), (14), (15), (16), (18), возвращает численное решение umjI (numerical), которое сравнивается с точным решением u(x, y)(exact) в таблице.

Таблица. - Сравнение численного решения umjI (numerical), полученного алгоритмом (8), (10), (12), (14), (15), (16), (18) с точным решением u( x, ym )(exact).

V

ym = c + ¿2 •m x = a + h • n u(x, ym )(exact) u (numerical)

1 2 3 4

1.25663706143 0.000000000000E+000 0.951056516295 0.951056516295

2.51327412287 0.000000000000E+000 0.587785252292 0.587785252292

3.76991118430 0.000000000000E+000 -0.587785252292 -0.587785252292

5.02654824574 0.000000000000E+000 -0.951056516295 -0.951056516295

6.28318530717 0.000000000000E+000 -2.4492127076E-016 -0.44921270764E-016

7.53982236861 0.000000000000E+000 0.951056516295 0.951056516295

8.79645943005 0.000000000000E+000 0.587785252292 0.587785252292

0.000000000000E+000 1.25663706143 0.951056516295 0.951056516295

1.25663706143 1.25663706143 0.573907707085 0.573907707086

2.51327412287 1.25663706143 0.205804929549 0.205804929546

3.76991118430 1.25663706143 -0.102572841154 -0.102572841152

5.02654824574 1.25663706143 -0.285044176197 -0.285044176199

6.28318530717 1.25663706143 4.28714624319E-002 4.287146243183E-002

7.53982236861 1.25663706143 0.447878357634 0.447878357636

8.79645943005 1.25663706143 0.635305917351 0.635305917348

0.000000000000E+000 2.51327412287 0.587785252292 0.587785252292

1.25663706143 2.51327412287 0.223220476438 0.223220476435

2.51327412287 2.51327412287 0.178820392472 0.178820392476

3.76991118430 2.51327412287 -0.115019453114 -0.115019453119

5.02654824574 2.51327412287 -4.4692996239E-002 -4.46929962365E-002

6.28318530717 2.51327412287 2.64960209303E-002 2.64960209302E-002

7.539822368615 2.513274122871 0.1453300548980 0.1453300548951

8.79645943005 2.51327412287 0.444266601135 0.444266601139

0.000000000000E+000 3.76991118430 -0.587785252292 -0.587785252292

1.25663706143 3.76991118430 -2.5661795207E-002 -2.5661795204E-002

2.51327412287 3.76991118430 -5.6722412698E-002 -5.6722412703E-002

3.76991118430 3.76991118430 -7.0785266586E-003 -7.0785266539E-003

5.02654824574 3.76991118430 -0.152865684991 -0.152865684994

6.28318530717 3.76991118430 -2.6496020930E-002 -2.6496020930E-002

7.53982236861 3.76991118430 5.22286263330E-002 5.22286263358E-002

8.79645943005 3.76991118430 -0.322168621361 -0.322168621366

0.000000000000E+000 5.02654824574 -0.951056516295 -0.951056516295

1.25663706143 5.02654824574 -2.1099959942E-002 -2.1099959944E-002

2.51327412287 5.02654824574 0.135849047429 0.135849047432

3.76991118430 5.02654824574 -0.239081135823 -0.239081135826

5.02654824574 5.02654824574 -0.267763570945 -0.267763570943

6.28318530717 5.02654824574 -4.2871462431E-002 -4.2871462431E-002

7.53982236861 5.02654824574 0.104929389507 0.104929389505

8.79645943005 5.02654824574 -0.293651940372 -0.293651940369

0.000000000000E+000 6.28318530717 -2.4492127076E-016 0.000000000000E+000

1.25663706143 6.28318530717 0.951056516295 0.951056516295

2.51327412287 6.28318530717 0.587785252292 0.587785252292

3.76991118430 6.28318530717 -0.587785252292 -0.587785252292

5.02654824574 6.28318530717 -0.951056516295 -0.951056516295

6.28318530717 6.28318530717 -2.5596176402E-016 -2.44921270764E-016

7.53982236861 6.28318530717 0.951056516295 0.951056516295

8.79645943005 6.28318530717 0.587785252292 0.587785252292

Замечание 1. Сравнение третьего и четвертого столбцов таблицы 1 показывает совпадение точного аналитического u(xn, ym )(exact) и численного решения ии л(numerical) в одиннадцати значащих цифрах. Программа вычисляет норму Чебышева разности между точным и численным решениями \um n(numerical)-(u(exact))mnIP"2 = max \umn(numerical)— (u(exact^))mn 1 = 5.3980153680299-10—12, с параметрами щ = 250, n2 = 375 (число интервалов деления стороны квадрата b — a = 2я, d — c = 3я) равномерной сетки. Известно также, что все другие функциональные нормы в численных методах не превышают нормы Чебышева в пространстве равномерно непрерывных функций [3; 4; 8].

Отметим, что в нашей задаче появляется проблема трансляции аналитического решения в численный массив. Решение (20) с учетом относительной погрешности 10~16 —10~15 в компиляторе с двойной точностью [18] имеет абсолютную ошибку (ch(2rc) + ch(3^})-10—15 « 6.436-10—12 > 5.398-10—12.

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

Замечание 2. Математики Новосибирской математической школы придерживаются мнения, что тестовые примеры в численных алгоритмах должны быть в некотором смысле экстремальными [14], поэтому в примере (19) выбрана быстро осциллирующая правая часть уравнения f (x, y) = sin(2x) sin(3y).

Замечание 3. Определим порядок погрешности алгоритма (8), (10), (12), (14), (15), (16), (18) совместно с формулой (5). Программа возвращает значение нормы Чебышева

ЦП, = 50,^ = 75 _g

Не

nn, =100,П2 =150 о

Р 2 _ 1 О/СС О 1 П—9

~ " 100'n2 =150 = 8.145-10-8/1.2658-10"9

ИС II m,n 4 \ 4 /-* m,n

~ 64,35 ~ 26 ^ p = 6. То есть порядок погрешности алгоритма не может превышать порядка погрешности аппроксимации уравнения Пуассона формулой (3), равного шести.

Замечание 4. В работе [5] методами простой итерации и скалярной прогонки с числом итераций m = 1800 (n = n2 = 111) с быстро осциллирующей правой частью f (x, y) = sin(2x) sin(3y) программа возвращает норму Чебышева (с краевыми условиями u(0,y) = и(л,y) = sin(yXy e [0,л],u(x,0) = u(x, л) = sin(xX x e [0, л])

II / \ 11«! =111,n2 =111 11 1 ro Л

для уравнения Пуассона равной ||um п (numerical) — (u(exact))тп || = 2.428 -10 за время At = 58.2 с, в то

/" II / \ ||n =111,n2 =111 _i 1

время как в данной работе норма невязки и время равны ||um п (numerical) — (u(exact ))m и || = 1.06 -10 ,

At = 0,87 с соответственно.

Отметим, что, как частный случай, из предложенного нами алгоритма можно получить матричный алгоритм, описанный в общих чертах в [3, с. 584-585]. Авторы Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. ввели термин - «метод прогонки в блочной форме». Там уравнение Пуассона решается на прямоугольнике с пятиточечным шаблоном «крест» со вторым порядком погрешности. Заметим, что в [3, с. 584-585] уравнение Пуассона решается на квадрате с одинаковым числом интервалов сетки. Мы обобщили их алгоритм для решения уравнения Пуассона на прямоугольник (n2 — 1)х(щ —1). Матрицы A, B имеют вид

II /Л II«I =50,^ =75 о

||umn(numerical) — (u(exact))mn^c = 8.145-10 , а при числе интервалов вдвое больше имеем

II / \ in =100,n2 =150 _„

||umn(numerical) — (u(exact))ти|| = 1.2658-10 , тогда порядок погрешности алгоритма p = 6 :

\um n(numerical) —(u(exact ))mn ||n 2 /| |um)! (numerical) — (u(exact))mn Ц"1 = Л = = 8.145-10~8/1.2658-10"

-4, m = n; m = 1, n¡ -1,n = 1, П -1 r _ _

|l,m = n;m = 1, n, -1,n = 1,n, -1 1, m = n +1 v m = n -1 , 6ти =] 1 1 ■ (21)

i, I 0, m ^ n

0, |m - n| > 2 1

То есть матрица B представляет собой единичную матрицу ранга n1 - 1. Вместо разностного уравнения (5) в [3] приводится уравнение

Um-1,n + Um+1,n + Um,n-1 + Um,n+1 - 4Um,n = Fm,n = fm,nh . (22)

Формулы (7) по-прежнему остаются в силе:

Au\ + Bul = FT

Bui_i + AuTm + BuTm+1 = F„T,m = 2,n2 — 2 . (23)

Bu1 , + AuT . = F1 .

n-> —2 n7 —i n-> —i

В (23) матрицы A, B определяются формулами (21). Граничные поправки (10), (12) переходят в

u 1+ u, , + u, + u„ ,- 4u, ,= h2 f ,= F ,

2,n -1 -2 1,n 0,n -1 -1 ''I, n -1 -1

^-1 - F1,n-1 - (u1,n + u0,n-1)

\-2,1 + -1,2 + uni-1,0 + uni,1 4un-1,1 = h fn-1,1 - Fn2-1,1

Fn2-1,1 - Fn1 -1,1 (un-1,0 + un2,1 )

u r, , + u , . + u , + u

n2-2,«!-1 n2-1,n-2 n2-1,n n2,n1 -

Fn2-1,n1-1 — Fn2-1,n1-1 -(un2-1,n1 + un2,n1-1 ) u2,1 + u1,2 + u1,0 + u0,1 - 4u1,1 = h f1,1 — F1,1 F1,1 - F1,1 -(u1,0 + u0,1)

- 4u = h2 f - F

1 ^ии2 -1,n,-1 " J n2 -1,n,-1 _ 1 n

n2 -1, n1 -1

(24)

u1,n-1 + u2,n + u1,n+1 + u0,n - 4u1,n = h fin - F1,n > n = 2 n1 - 2

F1,n = F1,n - u0,n > n = 2 n1 - 2

•'n -1,n-1 + un-2,n

+ u , ,,+ u -4u , = h f , - F ,,n = 2,n, -2

Fn,-1,n = Fn, -1,n - un, ,n > n = 2 n1 - 2

um-1,1 + um,2 + um+1,1 + u„,0 - 4um,1 = h fm,1 - Fm,1> m = 2 n2 - 2

Fm,1 = Fm,1 - u«,0> m = 2 n2 - 2

u , ,+ u . + u ,, ,+ u - 4u ,= h f , - F ,, m = 2, n - 2

m-1,n1 -1 m,n1 -2 m+1,n1 -1 m,n1 m,n1 -1 J m,n1 -1 m,n1 -1 > > 2

Fmn -1 = Fmn -1 - um,n > m = 2 n2 - 2

Fm,n = Fm,n > Vm 6 2 n2 - 2 n 6 2 n1 - 2

(25)

Формулы (14), (15), (16), (17), (18) по-прежнему остаются в силе для матриц А, В, определенных формулой (22):

= - A-1, V1 = A-1 F? , Xm = - (X„-1 + A)"1 , Vm = (X^ + A) ' (F,m - Vm-1 )

m = 2, n - 2 ;

(26)

Т

- =

1 =(X^-2 + A) 1 (Fl-1 -V^-2 ) > V^-1

(27)

um =Xmum+1 + Vm , m = "i - 2,1 .

(28)

Программа с использованием алгоритма (22)-(29) и условий примера (19) на квадрате со сторо-

/—II / \ in = 111,n2 = 111 _4

ной л возвращает норму невязки Чебышева pm „(numerical)-(u(exact ))m„|| = 1.027 •Ю за время

At = 0,312 с.

Сравнивая полученные программой результаты (Замечание 3) в алгоритмах (22)-(29) и в (8), (10), (12), (14), (15), (16), (18) при одних и тех же параметрах и режимах программы (Release), можно видеть, что быстродействие алгоритмах сопоставимо, но в алгоритме (8), (10), (12), (14), (15), (16), (18) норма невязки меньше в 107 раз!

Определение 1. Векторной нормой Чебышева ||x|| -I|x|| , x 6 R"1 -1 называют величину ||x|| = max IxJ.

Определение 2. Бесконечной нормой

n Л

= max a. .

матрицы a j 6 R

(n2-1)x(",-1)

называется число

j=1

Определение 3. Говорят, что матричная норма ||Л|| согласована с векторной нормой ЦхЦ^, если Ух е ^,\\АхЦ ^ЛЦхЦх .

Т

= u

n -1

Лемма. Матричная норма согласована с векторной нормой ||х|| .

n—1 n—

Доказательство. ||Ax|| = max V a,-.-*.- ^ max V a max Ixj = ||A|| ||x|| . Лемма доказана.

^ i=1,«2-1 j-\ ' i = ,n — j-1 ' = 'П — ^ ^

Определение 4. Говорят, что матрица A монотонна, если обратная к ней матрица имеет неположительные (неотрицательные) элементы [8].

Например, матрица А, определенная формулой (4), является монотонной, а матрица B нет.

Численный эксперимент показывает, что матрица, обратная к A, определенная формулой (4), является монотонной для различных рангов, то есть все элементы матрицы A~l неположительные. Оказывается, что все матрицы с диагональным преобладанием элементов, у которых диагональные и недиагональные элементы имеют противоположные знаки [8], являются монотонными.

Теорема 1. Невырожденная матрица со строгим диагональным преобладанием элементов, у которой диагональные элементы отрицательны (положительны), а недиагональные элементы положительны (отрицательны), является монотонной.

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

База индукции (n = 2).

А =

a 1 < 0, Kl > a 2 I' a2 2 < 0, Ы < 2 A = «2 2«n — a 2я2 1 > An = <

A2 = a 2< 0, A12 = -a2i < 0, Ai = -a 2 < 0, sign (ai + a 2) = sign (a 2) = -1, sign (°n+a 2) = sign (an) = -1,

sign (a" ) = (-1)n = (-1)2 = sign (a22a11 - a1,2 a21) = 1, am,n = —^ < 0, Vm, n = 1,2 .

Видим, что A монотонная, так как все элементы обратной матрицы отрицательны. Кроме того, диагональные алгебраические дополнения прямой матрицы - равносильно диагональные элементы обратной матрицы, больше по модулю недиагональных. База индукции проверена.

Индуктивный переход. Пусть A_j - монотонная, знак диагонального минора sign (Ал1) = (-1) , все диагональные миноры по модулю больше недиагональных миноров. Если amn, A п соответственно элемент и его алгебраическое дополнение, то [11, с. 32]

am,1 Am,1 + am,2 Am,2 + -«m,mAm,m + - + am,nAm,n = An

am,1 As,1 + am,2As,2 + -am,mAS,m + - + am,nAs,n = 0, Vs * m .

(29)

Так как sign (An ) = sign (am1 Am, 1 + am, 2 Am,2 + -am,mAm,m + ■■■ + am, nAm, n ) = Sign mAm, m ) = (—1) ' (—1O^ = (—1)П .

Выражая из формулы (29) диагональное алгебраическое дополнение, получим

sign (a ,A ,+ a nA n + ...a A +... + a A ) = sign (a A ) = sign (a )• sign (A ) = sign (A )

о у m,1 m,1 m,2 m,2 m,m m,m m,n m,n / ö V m,m m,m / ö V m,m / о у m,m / c> \ n /

ö

ö

sign (Am,m ) = sign (An ) sign (am,m ) = (—!)" (—10 = (—!)" + , am,m < 0Vm = 1

1 Am,m

= 1, n, a =—— ö

A

( am!m )=•

''А ^

ö sg ( am,m )= sign

V An У

= sign (Am,m ) sign (An ) = (—if+'(—!)" = —1 .

№ второго уравнения (29) имеем V «m,kAs,k =amA,1 + am,2 As,2 + -^mAm + ... + «m,„As,„ = 0 (m * s) ö

ö a_

n

~m,„As,„ = — V am,iAs,i ö sign (am,mAs,m ) = sign

k=1 k * m

= —sign (am,sAs,s ) = —sign (As,s ) = (—!)"

ö

ö siSn (am,mAs,m ) = ( —1)" ö siSn (am,m ) sg (As,m ) = (—1)" , sg (As,m ) = (—1)" sg (am,m ) = (—1)" • (—1) = (—1)"

k=1

Г A. ^

A-симметрическая

(а-т )= = ^8" () ^8" (Ап ) = ^8" (Кт ) ^8" (Ап ) = (-1)И+1 (-1)" =-1 •

V Ап у

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

Теорема 2. Невырожденная матрица со строгим диагональным преобладанием, у которой диагональные элементы отрицательны (положительны), а недиагональные элементы положительны (отрицательны), имеет верхнюю оценку бесконечной нормы для обратной матрицы

(A) = W-I К» R -

n=1

n* m

l-II

n=1 n* m

1

= min r (A) > 0, A"i < —

m=1,N-1 m V ' 11 II» R

Доказательство. Элемент единичной матрицы ранга N - 1

em,n -8m,n =I am]kak,n = öm11ö1,n + О,^,« + ••• + О^Оу-М = 8m,n . Vm> n = 1 N-1 •

Просуммируем каждое слагаемое в последней сумме по индексу столбца п = 1, ..., N - 1:

N-1 N-1 N-1 N-1 N-1 N-1 N-1 _

Ев = У 8 = 1 = У а\ у а, = а \ У а + а~\ У а +..

т п / , т,п / ' т к / , к,п т 1 / , 1," т 2 / , 2,"

-1n-11 Оу -1, n = 1, Vm = 1, N -1 •

По условию теоремы 2 диагональные и недиагональные элементы матрицы А имеют разные знаки,

f \

N-1

имеем I«m,„ =-

n=1

n=1

n*m J

-rm (A), -alr (A)-am]2r2 (A)-••• -<у-Л-1 (A) = 1,am]n < ö,Vm,n e 1,N -1 о

о a- r (A) + am,2 r2 (A) +... + a-N-1 rN- 1 (A) = 1 , обозначим индекс строки матрицы, обратной к A, сумма

модулей элементов которой максимальна, г0: ||4 || = а- + а- 2 +... + а- > 0 , но и для строки с номером г0 верно

(|a,-! + + ••• + |a-1N-1|rm (A) < Иr (A) + |a^|r (A) +... + |a- - 11^ (A) = rm (A) > ö °

о(|a-д|+|a-,2| +... + |a-У-]) mLrm (A) = ||A" mLrm (A)<1«||A" < . 1 ,,, = ^ •

\l 0' I I 0' I I \f m=1 ,N- 1 " "» m=1,N- 1 4 7 11 "» min r (A) R

m=1N-1 m V ' *

Теорема 2 доказана.

Теорема 3 (достаточные условия устойчивости алгоритма (8)-(18). Пусть выполнены условия: 1) трехдиагональные матрицы Теплица A, B симметрические ранга n -1 ann+1 = ял+1п,Ълл+1 =

= ЪП+1,П . Vn = 2 «1 - 2 :

2) матрица А монотонная с отрицательными диагональными элементами и положительными недиагональными элементами, В - положительная матрица ап-1,п ■ ап п < 0, Ъп т > 0, V" = 2,п1 - 2,У|п - т| < 1;

3) R = min

* m=1,n1 -1

n1 -1

l-I |a

m m | / , | m,n

n=1

n* m J

3 max Fi = F Т

ot=1,"2-1 CO

iii i i i 10 i i 2

= min r (A) = a -2 a , > 2\\B\\ , a > 1-IIB , a , <-IIB ;

и=1,Л1-1 3 3

"1 -1.

= max_ I\fmn\ <+» •

» m=1,n1 -1 n=1

Тогда:

1) IML< 1, Vm = 1, И2 -1;

N-1

k=1

n=1

n=1

n=1

2) формулы прогонки (14), (15), (16) корректны, то есть

v I < m

FT

да

5

<да,m - 1,n2 -1, <(2n2 -3)

F Т

да

5

< +да.

Замечание 6. Нормы матриц A, B, определенных формулой (4), удовлетворяют условиям 2), 3) теоремы 3:

a 1> 1015 «Г 1° > 1°\u |< *« 2 < 2Ы - 2 + 21 - !

Km| >J FIL «| Т>Т I ' rm,m-^ < ~~ « ^ < ^ 'FIL =Т + ^ = 1 ,

3 3

\а I-2\а ,|> 5 «

| m,m | m,m-11 ||

3 3 3

3 6

-1°

-2--2 > 2 . 3

Доказательство теоремы 3 проведем по индукции. Если матрица А симметрическая, то обратная матрица А- также симметрическая. Матрицы А + В и А - В удовлетворяют условиям теоремы 1 (и к ним применима теорема 2), так как имеют строгое диагональное преобладание элементов, а знаки диагональных и недиагональных элементов противоположны. Имеем оценку согласно [4, с. 57]

|У| =11-а~1 и <1 А-1 II и .

II 1 IIда ^ II Ида

1) Базу индукции проверим при к = 1, используя теорему 2. Имеем

Т 2 1

- -a_15 <ai 115 <-115 -

II Ида II Ида Н R И

5

а - 2 а ,

m,m m,m-1

3) 5 1 .<J_!!^-_< 1

215 2 ;

A-F1

<1 A_1 II

F

F Т

да

F1

2 5

5

База индукции проверена.

Индуктивный переход. Пусть верно ||Xt < 1, Vk -1, m -1. Тогда имеем R (A - 5 ) > R (5Х ,+ A) > R (В + A) « |а -b |-2а ,-b ,|> R (5Х ,+A)> |а + b |-2а ,+ b ,|«

* V / * V m-1 ; * у ; | n,n n,n | | n,n-1 n,n-1| * \ m-1 ; | n,n n,n | | n,n-1 n,n-11

|а„,„ | + |b„,„ | - 2 (|а„, ^ | - |b„, ^ |) > R (5X^1 + A) > |а„,„ | - |b„,„ | - 2 (|а„, ^ | + |b„,„-^) . l|Xm| |да-|-( 5Xm_, + A)-1 à <-

«

lB

lB

' R (5Xm-1 + A) |аи^| - |bnn| - 2 (|аи,и_! | +1Ьи,и_1 |)

lB

3) 5

_<-L!—LL

а - 2 а , - 5 2 5 - 5

n,n\ n,n-1 II ||да II II II II

- 1 .

11да II Ида

Первая часть теоремы 1 доказана У к -1, m с 1, n2 - 2.

Обозначим max

m-1,n -1

Fi

F1

n -1.

- max

да m-

^ Çl fmn\ .

2) Так как XI < 1, Vm -1, n2 - 2 , база индукции vJ -

F Т

A-1 F,Т <1 Ai F!' < да

да II Ида да 2 5

F1

5

проверена.

Индуктивный переход. Пусть v J < да, Vk -1, m -1, тогда

(5Хи-1 + A)-1 (Fi - 5vm-11 < ||(5Xm-1 + A)-1[ ||(Fi - 5vm-11 <

( Fi - 5v m-1 ) FТ m

да <.

-1И да

R (5Xm-1 + A)

l5ll да

*

V1 да-

v

m

да

Fl

„

я

+IK-J „ <

F T

„

я

+K-1II „< 2

F T

„

я

+ HV»-2|L< ••• <(m -1)

F T

„

я

+ v, < m

и 11i„

F T

„

(яХ ,+ a)-(Fr ,-Bv ,) < (яХ ,+ A)-1

V n2-2 / \ n2-1 n2-2 / V n2-2 !

F' , - Bv_

1 n-1

я

II 1ют

+1 я „||v

)

< +„ ; (30)

<

Т

<

<

u

= v

n -1

я

-+ v

П -2

<( П2 -1)

я

■ < +».

(31)

По формуле (18) имеем

»uT Л = X мт , + v J < X J u^ + v J < u^ + v J

П2-2| L II "-2 "-1 П2-2|\„ II "2-2\ UN "2-1U „ II "2-2 I1„ II П2 "И l„ II П2-2 II

F Т

„

= ("2 ^^¡Г+ ("2 - = (2" - 3)

я

F Т

„

я rj

F Т

„

я

<+» .

(32)

Теорема 3 доказана.

В программе, написанной на FORTRAN, все массивы и функции заданы с двойной точностью, использован пример (19) с известным аналитическим решением (20). Для вычисления обратных матриц с двойной точностью в алгоритме (8), (10), (12), (14), (15), (16), (18) использована библиотека msImsl [18] c вызовом подпрограммы call dLinrg(n,bb,n,ainv,n) [18].

program puasson;

integer(8),parameter::n1=250,n2=(n1*3)/2,n11=n1-1,n22=n2-1;integer(8)::i,j,k,kk; real(8)::a,b,c,d,pi,h,h1,h2,x,y;real(8)::delta(0:n2,0:n1),max,res(0:n2,0:n1);real(8)::c1,c2,c3,c4,c5,c6,c7,c8; real(8)::nu(n22,n11),lamda(n22,n11,n11),u(0:n2,0:n1),ff(n22,n11); real(8)::aa1(n11,n11),bb(n11,n11),aa(n11,n11),cc(n11,n11),u1(n22,n11);

real(8)::cch,ssh,a1,a2,b1,b2,f,t1,t2,f1,f2,f3,hh2,hh3,hh4,bb1(n11,n11);cch(x)=(dexp(x)+dexp(-x))/2d0;

ssh(x)=(dexp(x)-dexp(-x))/2d0;a1(x,y)=dsin(y);a2(x,y)=dsin(y);b1(x,y)=dsin(x);b2(x,y)=dsin(x);

f(x,y)=dsin(2d0*x)*dsin(3d0*y);f1(x,y)=-13d0*dsin(2d0*x)*dsin(3d0*y)/12d0;

f2(x,y)=97d0*dsin(2d0*x)*dsin(3d0*y)/360d0;f3(x,y)=4d-1*dsin(2d0*x)*dsin(3d0*y);max=-

1000d0;pi=2d0*dasin(1d0);

a=0d0;b=2d0*pi;c=0d0;d=3d0*pi;h1=(b-a)/dfloat(n1);h2=(d-c)/dfloat(n2);call cpu_time(t1); print*,"h1=",h1,"h2=",h2;do i=1,n11;do j=1,n11;if(j==i)then;aa(i,j)=-10d0/3d0;bb(i,j)=2d0/3d0; elseif(j==i+1.or. j==i-1)then;aa(i,j)=2d0/3d0;bb(i,j)=1d0/6d0;else;aa(i,j)=0d0;bb(i,j)=0d0; endif;enddo;enddo;do j=0,n1,1;x=a+h1*dfloat(j);u(0,j)=b1(x,c);u(n2,j)=b2(x,d);enddo; do i=0,n2,1;y=c+h2*dfloat(i);u(i,0)=a1(a,y);u(i,n1)=a2(b,y);enddo;hh2=h1*h1;hh3=hh2*hh2;hh4=hh3*hh2;do j=2,n11-1,1;

x=a+h1*dfloat(j);ff(1,j)=(f(x,c+h2)*hh2+f1(x,c+h2)*hh3+(f2(x,c+h2)+f3(x,c+h2))*hh4)-(1d0/6d0)*(u(0,j-1)+u(0,j+1))-(2d0/3d0)*u(0,j);

ff(n22,j)=(f(x,d-h2)*hh2+f1(x,d-h2)*hh3+(f2(x,d-h2)+f3(x,d-h2))*hh4)-(1d0/6d0)*(u(n2,j-1)+u(n2,j+1))-(2d0/3d0)*u(n2,j);

enddo;do i=2,n22-1,1;y=c+h2*dfloat(i);

ff(i,1)=(f(a+h1,y)*hh2+f1(a+h1,y)*hh3+(f2(a+h1,y)+f3(a+h1,y))*hh4)-(1d0/6d0)*(u(i+1,0)+u(i-1,0))-(2d0/3d0)*u(i,0);

ff(i,n11)=(f(b-h1,y)*hh2+f1(b-h1,y)*hh3+(f2(b-h1,y)+f3(b-h1,y))*hh4)-(1d0/6d0)*(u(i+1,n1)+u(i-1,n1))-

(2d0/3d0)*u(i,n1);

enddo;

ff(1,1)=(f(a+h1,c+h1)*hh2+f1(a+h1,c+h1)*hh3+(f2(a+h1,c+h1)+f3(a+h1,c+h1))*hh4)-(1d0/6d0)*(u(0,0)+u(2,0)+u(0,2))-(2d0/3d0)*(u(1,0)+u(0,1));

ff(n22,1)=(f(a+h1,d-h1)*hh2+f1(a+h1,d-h1)*hh3+(f2(a+h1,d-h1)+f3(a+h1,d-h1))*hh4)-(1d0/6d0)*(u(n2,0)+u(n2,2)+u(n2-2,0))-(2d0/3d0)*(u(n2-1,0)+u(n2,1));

ff(1,n11)=(f(b-h1,c+h2)*hh2+f1(b-h1,c+h2)*hh3+(f2(b-h1,c+h2)+f3(b-h1,c+h2))*hh4)-(1d0/6d0)*(u(0,n1)+u(2,n1)+u(0,n1-2))-(2d0/3d0)*(u(0,n1-1)+u(1,n1)); ff(n22,n11)=(f(b-h1,d-h2)*hh2+f1(b-h1,d-h2)*hh3+(f2(b-h1,d-h2)+f3(b-h1,d-h2))*hh4)-(1d0/6d0)*(u(n2,n1)+u(n2-2,n1)+u(n2,n1-2))-(2d0/3d0)*(u(n2,n1-1)+u(n2-1,n1));

do i=2,n2-2;do j=2,n1-

2;x=a+h1*dfloat(j);y=c+h2*dfloat(i);ff(i,j)=(f(x,y)*hh2+f1(x,y)*hh3+(f2(x,y)+f3(x,y))*hh4); enddo;enddo;call obr(n11,aa,aa1);lamda(1,1:n11,1:n11)=-matmul(aa1(1:n11,1:n11),bb(1:n11,1:n11)); nu(1,1:n11)=matmul(aa1(1:n11,1:n11),ff(1,1:n11));do i=2,n22-1;

bb1(1:n11,1:n11)=matmul(bb(1:n11,1:n11),lamda(i-1,1:n11,1:n11))+aa(1:n11,1:n11);call obr(n11,bb1,aa1); Iamda(i,1:n11,1:n11)=-matmul(aa1(1:n11,1:n11),bb(1:n11,1:n11));

nu(i,1:n11)=matmul(aa1(1:n11,1:n11),(ff(i,1:n11)-matmul(bb(1:n11,1:n11),nu(i-1,1:n11))));enddo; bb1(1:n11,1:n11)=matmul(bb(1:n11,1:n11),lamda(n22-1,1:n11,1:n11))+aa(1:n11,1:n11);call obr(n11,bb1,aa1); u1(n22,1:n11)=matmul(aa1(1:n11,1:n11),(ff(n22,1:n11)-matmul(bb(1:n11,1:n11),nu(n22-1,1:n11))));do i=n22-

1,1,-1;

u1(i, 1:n11)=matmul(lamda(i, 1:n11,1:n11),u1(i+1,1:n11))+nu(i, 1:n11);enddo;do i=1,n22;do

j=1,n11;u(i,j)=u1(i,j);

enddo;enddo;max=-1d3;open(2, file='2.txt');do j=0,n1;do i=0,n2;x=a+h1*dfloat(j);y=c+h2*dfloat(i);

c1= dsin(y)*(cch(x)+ssh(x)*(1d0-cch(2d0*pi))/ssh(2d0*pi));c2= dsin(x)*(cch(y)+ssh(y)*(1d0-

cch(3 d0*pi))/ssh(3d0*pi));

c3= -dsin(2d0*x)*dsin(3d0*y)/13d0;res(i,j)=c1+c2+c3;delta(i,j)= u(i,j)- res(i,j);if( delta(i,j)<=0d0)then; delta(i,j)=-delta(i,j);endif;if( delta(i,j)>max)then;max=delta(i,j);endif;if(mod(i,50)==0.and.mod(j,50)==0)then; 2 write(2,*),y,x;3 write(2,*),res(i,j),u(i,j);endif;enddo;enddo;print*,"norma C =",max,"num="; call cpu_time(t2);print*,"t2-t1=",t2-t1;pause;end program puasson;

subroutine obr(n11,bb,aa1);use msimsl;integer(8)::n11,Lda,Ldainv;real(8)::ainv(n11,n11)

,bb(n11,n11),aa1(n11,n11);

Lda=n11;Ldainv=n11;call dLinrg(n11,bb,n11,ainv,n11);aa1(:,:)=ainv(:,:); end subroutine;

Случай матричной прогонки (h * h2), но с равным числом интервалов n = П для уравнения Пуассона на квадрате (прямоугольнике) с 4 порядком погрешности описан нами в работах [6; 21].

В работе получены результаты:

1) предложен алгоритм прогонки в матричной форме для решения уравнения Пуассона на прямоугольнике с шестым порядком погрешности (формулы (8), (10), (12), (14), (15), (16), (18) за конечное число арифметических операций;

2) показано, что предложенный алгоритм переходит как частный случай в известный алгоритм прогонки в матричной форме Бахвалова Н.С., Жидкова Н.П., Кобелькова Г.М. [3] со вторым порядком погрешности (формулы (22)-(29);

3) на примере (19) и программой, написанной на FORTRAN, показано, что алгоритм (8), (10), (12), (14), (15), (16), (18) имеет шестой порядок погрешности для решения уравнения Пуассона на прямоугольнике и выполняется за конечное число элементарных операций;

4) в теореме 1 доказана монотонность матриц с диагональным преобладанием, у которых элементы главной диагонали отрицательны (положительны), а недиагональные положительны (отрицательны);

5) в теореме 2 получена верхняя оценка бесконечной нормы обратной к монотонной матрице;

6) в теореме 3 получены достаточные условия корректности алгоритма (8), (10), (12), (14), (15), (16), (18);

7) быстродействие приведенного алгоритма в десятки раз превышает быстродействие алгоритма для решения уравнения Пуассона на прямоугольнике методом простой итерации с той же формулой аппроксимации с шестым порядком погрешности и относительной погрешностью.

ЛИТЕРАТУРА

1. Козлов, А.А. Преобразование подобия на множестве полукватернионов / А.А. Козлов, К.С. Сурав-нева, И.Л. Жалейко // Вестник Полоцкого государственного университета. Серия С, Фундаментальные науки. - 2019. - № 4. - С. 115-123.

2. Козлов, А.А. Множество полуоктав / А.А. Козлов // Вестник Полоцкого государственного университета. Серия С, Фундаментальные науки. - 2016. - № 12. - С. 75-85.

3. Бахвалов, Н.С. Численные методы / Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. - 7-е изд. - М. : БИНОМ. Лаборатория знаний, 2011. - 636 с.

4. Бахвалов, Н.С. Численные методы в задачах и упражнениях / Н.С. Бахвалов, А.В. Лапин, Е.В. Чижон-ков. - М. : БИНОМ, 2010. - 240 с.

5. Пастухов, Д.Ф. Аппроксимация уравнения Пуассона на прямоугольнике повышенной точности / Д.Ф. Пастухов, Ю.Ф. Пастухов // Вестник Полоцкого государственного университета. Серия С, Фундаментальные науки. - 2017. - № 12. - С. 62-77.

6. Волосова, Н.К. Модифицированное разностное уравнение К.Н. Волкова для уравнения Пуассона на прямоугольнике с четвертым порядком погрешности // Евразийское Научное Объединение. - 2019. № 6-1 (52). С. 4-11.

7. Александров, П.С. Введение в теорию размерностей / П.С. Александров, Б.А Пасынков. - М. : Наука, 1973. - 577 с.

8. Волков, Ю.С. Оценки норм матриц, обратных к матрицам монотонного вида и вполне неотрицательным матрицам / Ю.С. Волков, В.Л. Мирошниченко // Сиб. мат. журн. - 2009. - Т. 50, №6. -С. 1249 - 1254.

9. Пикулин, В.П. Практический курс по уравнениям математической физики : учеб. пособие / В.П. Пи-кулин, С.И. Похожаев. - М. : Наука, 1995. - 224 с.

10. Волосова, Н.К. Векторный аналог метода прогонки для решения трех- и пятидиагональных матричных уравнений / Волосова [и др.]. // Вестник Полоцкого университета. Серия С. Фундаментальные науки. - 2019. - № 12. - С. 101-115.

11. Ильин, В.А. Линейная алгебра / В.А. Ильин, Э.Г. Поздняк. - М. : Наука : Физматлит - 1978. - 304 с.

12. Пастухов, Д.Ф. Оптимальный порядок аппроксимации разностной схемы волнового уравнения на отрезке / Д.Ф. Пастухов, Ю.Ф. Пастухов, Н.К. Волосова // Вестник Полоцкого государственного университета. Серия С, Фундаментальные науки. - 2018. - № 12. - С. 60-74.

13. Пастухов, Д.Ф. К вопросу о редукции неоднородной краевой задачи Дирихле для волнового уравнения на отрезке / Д.Ф. Пастухов, Ю.Ф. Пастухов, Н.К. Волосова // Вестник Полоцкого государственного университета. Серия С, Фундаментальные науки. - 2018. - № 4. - С. 167-186.

14. Годунов, С.К. Современные аспекты линейной алгебры / С.К. Годунов. - Новосибирск : Научная книга, 1997. - 407 с.

15. Вакуленко, С.П. Способы передачи QR-кода в стеганографии / С.П. Вакуленко, Н.К. Волосова, Д.Ф. Пастухов // Мир транспорта. - 2018. - Т. 16, № 5 (78). - С. 14-25.

16. Пастухов, Д.Ф. Некоторые методы передачи QR-кода в стеганографии / Д.Ф. Пастухов, Н.К. Волосова, А.К. Волосова // Мир транспорта. - 2019. - Т. 17, № 3 (82). - С. 16-39.

17. Волосова, Н.К. Применение преобразования Радона в стеганографии // Некоторые актуальные проблемы современной математики и математического образования : сб. материалов науч. конф., Герценовские чтения - 2018, СПб., 9-13 апр. 2018 г. / Рос. гос. пед. ун-т им. А.И. Герцена. - СПб., 2018. - С. 234-238.

18. Бартеньев, О.В. Фортран для профессионалов. Математическая библиотека IMSL: Ч. 1. - М. : ДИАЛОГ : МИФИ, 2001. - 437 с.

19. Пастухов, Д.Ф. Минимальная разностная схема для уравнения Пуассона в параллелепипеде с шестым порядком погрешности / Д.Ф. Пастухов, Ю.Ф. Пастухов, Н.К. Волосова // Вестник Полоцкого государственного университета. Серия С, Фундаментальные науки. - 2019. - № 4. - С. 154-174.

20. Salih, A. Streamfunction-vorticity formulation // Indian Institute of Space Science and Technology, Department of Aerospace Engineering, Thiruvananthapuram. - 2013.

21. Волосова, Н.К. О решении уравнения Пуассона на прямоугольнике с четвертым порядком погрешности за конечное число элементарных операций // Евразийское Научное Объединение. - 2020. № 2-1 (60). С. 11-17.

Поступила 12.02.2020

ON FINITE METHODS FOR SOLVING THE POISSON EQUATION ON A RECTANGLE WITH THE DIRIHLET BOUNDARY CONDITIO

N. VOLOSOVA, K. VOLOSOV, A. VOLOSOVA, D. PASTUHOV, Y. PASTUHOV

An algorithm for sweeping in matrix form with a sixth order of error for solving the Poisson equation on a rectangle in a finite number of arithmetic operations is proposed. An analytical example and a program using this algorithm confirmed the sixth order of error. Theorem 1 proves the monotonicity of matrices with diagonal dominance, for which the elements of the main diagonal are negative (positive), and the off-diagonal are positive (negative). In Theorem 2, an upper bound is obtained for the infinite norm inverse to a monotonic matrix. In Theorem 3, sufficient conditions for the correctness of the proposed algorithm are obtained. It is shown that the speed of this algorithm is ten times higher than the speed of the algorithm for solving the Poisson equation on a rectangle using the simple iteration method with the same approximation formula with sixth error order.

The Keywords: sweep method in block form, diagonal matrices, monotone matrices, inverse problems mathematical physicists, numerical methods, Poisson equation, translation of an analytical solution into a numerical array.