CC BY
8УДК 517.6: 517.958
МИНИМАЛЬНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ ПУАССОНА НА ПАРАЛЛЕЛЕПИПЕДЕ С ШЕСТЫМ ПОРЯДКОМ ПОГРЕШНОСТИ
канд. физ-мат. наук, доц. Д.Ф. ПАСТУХОВ, канд. физ-мат. наук, доц. Ю.Ф. ПАСТУХОВ (Полоцкий государственный университет);
Н.К. ВОЛОСОВА
(Московский государственный технический университет им. Н.Э. Баумана)
Предложен алгоритм решения общей неоднородной краевой задачи Дирихле для трехмерного уравнения Пуассона на параллелепипеде с шестым порядком погрешности и с минимальным 27-точечным шаблоном. Численно проверена устойчивость алгоритма к разрывам первого рода граничных условий на гранях параллелепипеда. Получено разложение невязки задачи в общем виде для неравномерной сетки через производные четного порядка от решения и правой части уравнения и четного порядка по каждой из трех переменных. Написана программа на основе построенного алгоритма и принципа сжатых отображений для явной формулы простой итерации. Решен точно тестовый пример и проведено сравнение его решения с численным решением, которое подтверждает шестой порядок погрешности для формул полученного алгоритма.
Ключевые слова: трехмерное уравнение Пуассона на параллелепипеде, неоднородно-краевая задача Дирихле, принцип сжатых отображений, уравнения математической физики.
Введение. Рассмотрим неоднородную краевую задачу для трехмерного уравнения Пуассона в параллелепипеде с краевыми условиями Дирихле. Известны разностные схемы для уравнения Пуассона с погрешностью второго порядка в двухмерной области А. А. Самарского [1] с погрешностью четвертого порядка на равномерной сетке К.Н. Волкова [2], также схемы большего порядка на прямоугольнике и на равномерной сетке [4]. Во всех случаях для тестирования алгоритма необходимы примеры с точными аналитическими решениями [3]. Однако, как показано в работе [1, с. 57], для схем второго порядка в задаче уравнения Пуассона на прямоугольнике удается достичь равномерно непрерывной нормы погрешности 3 10 при делении его сторон на 160 частей. Отметим работу [11] с формулами для трехмерного лапласиана.
В данной работе предложен простой алгоритм решения трехмерного уравнения Пуассона на параллелепипеде с использованием минимального симметричного шаблона для аппроксимации (27 узлов)
и достижении равномерно непрерывной нормы погрешности 10_9 при делении ребра параллелепипеда на 60 частей и с аппроксимацией шестым порядком погрешности. Сложность трехмерной задачи, по сравнению с двухмерной задачей, заключается в быстром росте количества вычислений и в увеличении общего числа узлов. Например, 10 000 узлов на прямоугольнике и 1 000 000 узлов в параллелепипеде. По сравнению с работой [4], в которой использовались два шаблона, граничный и внутренний на прямоугольнике с восьмым порядком погрешности, в данной работе использование одного универсального
3 3 3 3 3
шаблона позволяет экономить число и время вычислений в п (3 + 5 )/3 = 5,6« число раз за счет выбора одного шаблона, но с шестым порядком погрешности. Кроме того, нужно учитывать, что использование двойного шаблона увеличивает число вычислений в 2-3 раза в каждом узле и уменьшает быстродействие. В работе использована модифицированная формула бинома Ньютона - тринома Ньютона (по аналогии с биномом - прим. авт.).
Полученный алгоритм численного решения трехмерной задачи для уравнения Пуассона может быть применен в различных областях механики [5], кристаллографии, стеганографии [6, 7], для численных задач математической физики, содержащих трехмерный оператора Лапласа, например, волновое уравнение [11, 14-17].
Постановка задачи. Рассмотрим неоднородную краевую задачу для трехмерного уравнения Пуассона в параллелепипеде для достаточно гладкого решения и (х, у, г) как функции трех переменных:
(и^ + иуу + игг) = /(х, у, г), хе (а,Ь), у е (с,й), г е (е, г)
и(а, у, г) = цх(у, г), и(Ь, у, г) = ^(У, г), (у, г) е [с, й]х[е, г] ^
и(х, с, г) = т?(х, г), и(х, й, г) = Ц4(х, г), (х, г) е [а, Ь]х [е, г] и(х, у, е) = ^5 (х, у), и(х, у, г) = Цб (х, у), (х, у) е [а, Ь] х [с, й].
В задаче (1) /(х, у, г) - неоднородная правая часть уравнения Пуассона в параллелепипеде :[а, Ь] х [с, й] х [е, г], у, г), У,г), Цэ(х,г), Ц4(х, г), Ц-5(х, у), Цб(х, у) - неоднородные краевые условия, (х, у, г) — координаты точки. Обозначим внутреннее и граничное множества задачи (1)
О = ( а, Ь)х( с, й )х( е, г), ЭО = [с, й]х [е, г] (х = а, Ь) и [а, Ь] х [е, г] (у = с, й) и [а, Ь]х [с, й] (г = е, г).
В классической постановке задач математической физики необходимо требовать непрерывности граничных условий на пересечении их областей определения, т.е. непрерывность краевых условий на 12 ребрах параллелепипеда, например,
Ц1(0, г) = цз(0, г), "г е [е, г].
Рассмотрим аппроксимацию задачи (1) на шаблоне, содержащем минимальное количество узлов и покрывающим всю область параллелепипеда. Ясно, что на 27-точечном шаблоне, который содержит коробка со сторонами 2/^,2^2,2^3, и решается поставленная задача с минимальным числом узлов в нем. Такой шаблон является универсальным, как на границе рассматриваемого параллелепипеда, так и внутри данной области, и покрывает параллелепипед объема (Ь - а)(й - с)(г - е) кирпичиками с объемом 8/1/2/3 (рисунок).
о * о о * о
О * 0 0 * 0
Рисунок. - Классификация 27 узлов куба (параллелепипеда) по расстоянию от центра и свойствам симметрии для 1, 6, 12, 8 узлов
Докажем несколько вспомогательных утверждений.
Лемма 1. Пусть решение задачи (1) принадлежит классу функций и(х,у,г)е С2Р (О), тогда для суммы узловых значений решения в шести центральных узлах граней куба (параллелепипеда) справедлива формула разложения в ряд Тейлора с центром в узле (т, п, к) :
ит-1,п,к + ит+1,п,к + ит,п-1,к + ит,п+1,к + ит,п,к-1 + ит,п,к+1 6ит,п,к + 2 2
Р 1
=1 (21)!
21 Э и 21 Э и ,21 Э и
/ -тГ + / -гг + / -ТГ-
1 2 ^,21 3 ^21
21 Л
Эх
Эу
Эг2
+о (/2 Р+2+/2 Р+2+/2 Р+2)
(2)
Доказательство. Объединяя попарно противоположные узлы на каждой координатной оси относительно центра (т, п, к) , учитывая, что в силу симметрии сохранятся только частные производные четной степени по соответствующей координате, получим
ит-1,п,к + ит+1,п,к + ит,п-1,к + ит,п+1,к + ит,п,к-1 + ит+1,п,к+1 = = (ит-1,п,к + ит+1,п,к ) + (ит,п-1,к + ит,п+1,к ) + (ит,п,к-1 + ит,п,к+1) =
= 2ит,п,к +
Р 1
( . ъ21 Л
1 (2/)!
1.21
Э и
1 Эх21
+ о (/2 р+2) + 2итлк + 2Х
Р 1
( .->2/. Л
=1 ( 2/)!
2/
2/ Эи
2 Эу2/
+ о
(П2
2 р+2\
+2ит,п,к + 2Х
Р 1
( ^2/ Л
ы (2/)!
2/
2/ Э и
3 Эг21
+ О (/2 р+2 ) =
= 6ит,п,к + 2Х
Р 1
=(2/)!
2/
2/ 2/ Э и 2/ Э и - + п
2/.
+ /
2/
Э 21 Л Эи
1 Эх2/ 2 Э ,2/ 3 Э„ 2/
Эу 2
Эг 2
О (/2 Р+2 + /2 Р+2 + Аз2 Р+2).
+
Что и завершает доказательство формулы (2) Леммы 1. В частности, для равномерной сетки получим из выражения (2) формулу (3)
ит_1,п,к + ит+1,п,к + ит,п_1,к + ит,п+1,к + ит,п,к—1 + ит,п,к+1
= 6ит,п,к + 2Е
Р И
ы (21)!
21 (-ч 21 л21 л21 А
Э и Э и Э и
21 +ЧЛ + л 21
+ О
Эх2' Эу2' Эг2'
(И2Р+2).
(3)
Лемма 2. Пусть решение задачи (1) принадлежит классу функций и (х,у, г)е С2Р (О), тогда для суммы узловых значений решения в 12 узлах - серединах ребер куба (параллелепипеда) справедлива формула разложения в ряд Тейлора с центром в узле с индексами (т, п, к)
ит_1,п+1,к + ит_1,п,к+1 + ит_1,п_1,к + ит_1,п,к—1 + ит+1,п+1,к + ит+1,п,к+1 + +ит+1,п_1,к + ит+1,п,к—1 + ит,п_1,к_1 + ит,п_1,к+1 + ит,п+1,к_1 + ит,п+1,к+1 = Р ' л ( . .. . 32',. _ _ ^21
= 12ит,п,к + 4ЕЕ
1
¡=1 ;=0 (2')!(2' _2')!
32!, ~\21г -\2', А
,2',2'_2' Э и ,1,2'.2'_2' Э и . ,2',21 _2s Э и Й1 ¿2 _ о. о, + Й Й3 _ о. о, + ¿2 Й3
Эх24 Эу2' _2'
Эх2'Эг2' _2'
2'-, 2'_2'
Эу24 Эг
у
+О(1Р=+01(И?Р+2_2гЙ22г' + Й,2Р+2_2гЙ32г' + Аз2Р+2_2гЙ12г')).
(4)
Доказательство. Сгруппируем в формуле (4) суммы 3 четверок узлов, расположенных в плоскостях параллельных граням куба (параллелепипеда) и проходящих через центральный узел с индексами (т, п, к). Для каждой четверки узлов разложим в ряд Тейлора сумму четырёх узловых значений, например, в узлах плоскости 0ху получим разложение по формуле бинома Ньютона:
А = ит_1,п_1,к + ит_1,п+1,к + ит+1,п_1,к + ит+1,п+1,к =
2р ' С' ( Э1и I д'и I
= 4ит,п,к + Е Е -Т л ^ ' _' (_И1)' (_И2 ) _' + Э [_, (_И1 )' ( Й2 ) _' +
'=14=0 '! IЭх4Эу'"'
'.. , а'
Эх4Эу'
( Й1)' (_Й2)' _' ( Й1)' ( Й2 )' ] + О ( Е2=0+1( 12 ^ )) =
2р
'!
= 4ит,п,к + ЕЕ -- ( А )' ( Й2 )' _' ((_1)'+' _' + (_1)' +(_1)'_' +1) .
'=1 '=0
Эх'Эу
Имеем далее
(_ 1)'+'_' +(_1)' +(_1)' _' + 1 = (_1)' (1 +(_1)_' ) + (_1)' + 1 = ((_1)' + 1)(1 +(_1)' ) :
Г 4,' = 2'1 л' = 2'1
[0,' = 2'1 +1V' = 2'1 +1
ЕЕ -21и2',„ 2'_2' Э2'и + О(Ер=+о(й2Р+2_2гй22г)).
(5)
: ('1 ® ', '1 ® ') получим А = 4ит п к + 4Е Е Т^Й2'й2
Поскольку
Ы '=0 (2')! (2')! 1
Эх2' Эу2' _2'
^2'
—21 =_=_
(2')! (2')!(2')!(2' _ 2')! (2')!(2' _ 2')!
ит_1,п_1,к + ит_1,п+1,к + ит+1,п_1,к + ит+1,п+1,к =
, то последняя формула примет вид
= 4ит,п,к + 4ЕЕ
р ' й2'й221_2' э
2'.
£ ;=)(2')!(2' _ 2')! Эх2'Эу2'_2'
+О (ЕР=01( й2 Р+2_2гЙ22г')) .
(6)
+
Две четверки узлов и суммы узловых значений решения в соответствующих узлах в плоскостях 0хг, 0уг имеют ту же симметрию и формулы разложения аналогично формуле (6), тогда, меняя циклически индексы и складывая все 12 узловых значений, получим
Р1 Н12'Н221 -2* Э2'и . Р' й^2' -2* Э2'и
4ит,пЛ + 4УУ/0 шл^мо^ + т,п,к + 4УУ" 1 3 '
/=1 ;=о(2^)!(2/-25)! Эх25Эу21-25 /=1 ;=0(2^)!(2/-2*)! Эх2*Эг21-25
р ' ь 2* г 2/-2* -ч21, , . ,
+4ит п к + 4]Г У ^^--+ о(У Р+01(/2Р+2-2гА22г' + Р+2-2гА32г' + /3Р+2-2г//12г')) =
т'п'к /=1 *=0(2*)!(2/-2*)! Эу2*Эг2'-2* 2 *2 3 13 1 П
р± 1 Э2/и . Э2/и . Э2/и Л
= 12ит,п,к + 4У У "
,2*, 21 -2* Э и . ,2*, 2/-2* Э и . ,2*,2/-2* Э и
+
+О(У Р=+о (/2Р+2"2iА22i' + л|Р+2-2гА32г' + /2Р+2~2\2)).
Лемма 2 и формула (4) доказаны. В частности, для равномерной сетки получим формулу
ит-1,п+1,к + ит-1,п,к+1 + ит-1,п-1,к + ит-1,п,к -1 + ит+1,п+1,к + ит+1,п,к+1 +
+ит+1,п-1,к + ит+1,п,к -1 + ит,п-1,к-1 + ит,п-1,к+1 + ит,п+1,к-1 + ит,п+1,к+1 = Р I 1 ( Э2/ Э2/ Э2/ Л
= 12итпк + 4УА2' У--Т1-— 9Э и о + / и о + / и о + О(/2Р+2). (7)
т,п,к ^ ^ (9„)!(9/-9„)! -, 251 21-2* -.2^21-2* 2/-2* I I у '
Рассмотрим модифицированную формулу Ньютона
(а + Ь + с )п = У ск1'к2'к3 ак1Ьк2 ск3, ск^^3 =---. (8)
1 } о±к£: п п к1!к2!к3! ()
0<к2 <п 0<к3<п кх+к^ +к3 =п
Доказательство формулы (8). Введем переменные
к1 = 7'Ю е [0, г], к2 = г - '0 е [г,0], к3 = п - г, к1 + к2 + к3 = ' + г - ' + п - г = п, поскольку г е [0, п], то области изменения переменных к1 е [0, п], к2 е [п,0], к3 = п - г е [п,0], тогда имеем
п . п г
(а + Ь + с)п = У Сп (а + Ь)гсп-г' = У У С1пС'а'Ь1~'сп~1 =
г=0 г=0 у=0
= У --ак1Ьк2 ск3 = У -п--ак1Ьк2 ск3 =
0<к1<п г!(п - г)! ' !(г - ])! 0<к1<п ' !(г' - № - г')!
0<к2 <п 0<к2 <п
0<к3<п 0<к3<п
к1+к2 +к3 =п к!+к2 +к3 =п
У Ск1,к2,к3 ак1 Ьк2 ск3, =.
п! п!
0<к1<п ' ' '!('" - № -')! к1!к2!к3!
0<к2 <п 0<к3<п к1+к2+к3=п
Формула (8) доказана.
Лемма 3. Пусть решение задачи (1) принадлежит классу функций и (х,у, г)е С2р (О), тогда для суммы узловых значений решения в 8 вершинах куба (параллелепипеда) справедлива формула разложения в ряд Тейлора с центром в узле с индексами (т, п, к)
ит-1,п-1,к-1 + ит-1,п-1,к+1 + ит-1,п+1,к-1 + ит-1,п+1,к+1 + ит+1,п-1,к-1 + ит+1,п-1,к+1 + ит+1,п+1,к-1 + ит+1,п+1,к+1 =
= 8ит,п,к + 8Е Е
а2'1 А2'2 А2'3
Э и
М0<'^£Г (2'1)!(2'2)!(2'3)! Эх2'1 Эу2'2 Эг2'3 0<'2 <' 0<'3<'
'1+'2+'3='
+ О
Е а2'1 И2'2 А2'3
0<'1<'+1 0<'2<'+1 0<'32<'+1 '1+'2+'3='+1
(9)
Доказательство. Разложим в формуле (9) сумму 8 узловых значений в ряд Тейлора в точках, совпадающих с вершинами куба (параллелепипеда) относительно центрального узла с индексами (т, п, к) . Воспользуемся модифицированной формулой бинома Ньютона (8)
(а + Ь + с)п = Е -п1''2''3а'1 Ь'2с'3,пе Сп1''2''3 =
п!
'1!'2!'3!
, '1 + '2 + '3 = п ,
0<'1<п 0<'2 <п 0<'3<п '1+ '2 +'3=п
Подставляя восемь раз формулу (8) в левую часть (9), учитывая С!1-'2-'3 п! 1
и! п!^!^!^! '1!'2!'3! (т, п, к)
, '1 + '2 + '3 = п , получим относительно центрального узла с индексами
um—1,n—1,k _1 + um—1,n—1,k+1 + ит_1,п+1,к _1 + ит_ 1,п+1,к+1 + ит+1,п_1,к _1 + ит+1,п_1,к+1 + ит+1,п+1,к _1 + ит+1,п+1,к+1
= 8ит,п,к + Е Е
2 р ' С'1.'2.'3
2Р ^ С' А'1 Й22 Й'3 Эи '
'=10<'1 <' 0<'2 <' 0<'з<' 1+ '2 +'3 ='
'!
Эх^ау'2 Эг'3
( (_ 1)'1+'2 +'3 + (_1)'1+'2 + (_ 1)'1+ '3 + (_1)1 +
+ (_1) '2 + '3 + (_1)'2 + (_1)'3 + 1) +
О
Е А'1 и22 Й3'3
0<'1 <2 Р+1 0<'2 <2 р+1 0<'32<2Р+1 '1 +'2 +'3 =2 Р+1
Кроме того, преобразуем сумму
(_1)'1+'2 +'3 + (_1)'1+'2 + (_1)'1+'3 + (_1)'1 + (_1)'2 +'3 + (_1)'2 + (_1)'3 +1 = = (_1)'1 ((_1)'2 +1) + (_1)'3 ((_1)'2 +1) + (_1)'1+'3 ((_1)'2 +1) + (_1)'2 +1 = ((_1)'2 +1)((_1)'1 + (_1)'3 + (_1)'1+'3 +1) =
= ((_1)'2 +1)((_1)' (1 + (_1)'3 ) +1 + (_1)'3 ) = ((_1)'2 + 1)((_1)'3 + 1)((_ 1)'1 +1) = 8, '1 = 2'1 л '2 = 2'2 л '3 = 2'3
= • , (Ю)
0, '1 = 2'1 +1V '2 = 2'2 +1V '3 = 2'3 +1.
то, упрощая множитель в двойной сумме используя формулу (10) и возвращаясь к переменным, ' = '1 + '2 + '3 = 2('1 + '2 + '3 ) = 2' , '1 ® '1, '2 ® '2, '3 ® '3, ' ® ' , получим
2 Р I С'1,'2,'3 д'
Е е С_А'1 Й'2 А'3_ _
'=10<'1<' '! Эх'1Эу'2 Эг'3
0<'2 <' 0<'32<' '1+'2 +'3='
( (_1)'1+'2 + '3 + (_1)'1+ '2 + (_1)'1+'3 + (_1)'1 +
+ (-1)*2 + *3 + (-1)*2 + (-1)*3 +1) + О
У
0<*1<2 р+1 0<*2 <2 р+1 0<*3<2 р+1 + +*3 =2 Р+1
И*1 А23 /3
53
р / = 8^ У
/2^ h2sз n2sз
Э и
/ =10< *!</
(2*1)!(2*2)!(2*3)! Эх2*1 Эу2^Эг2*3
0< </ 0< ^ +*2 + *3
+ О
У А22sз 0< р+1 0< < р+1 0< ^ р+1 + *3 = Р+1
(11)
Добавляя к выражению (11) 8ит п к, получим формулу (9). Лемма 3 доказана. В частности, на равномерной сетке формула (9) переходит в формулу (12)
ит-1,п-1,к-1 + ит-1,п-1,к+1 + ит-1,п+1,к-1 + ит-1,п+1,к+1 + ит+1,п-1,к-1 + ит+1,п-1,к+1 + ит+1,п+1,к-1 + ит+1,п+1,к+1 =
Р
1
Эи
= ^ + 8Ул ^ (2*1)!(2*2)!(2*3)! Эх2*Эу2^Эг2*3
0<sз</ +*3
+ О
(/2 р+2).
(12)
Замечание 1. В формуле (12) суммирование происходит по 4 индексам ^52,.3,/, однако в силу уравнения связи . + + *3 = / независимыми остаются только 3 переменных, то есть данная формула представляет тройную сумму по узлам параллелепипеда.
Лемма 4. Пусть решение уравнения Пуассона (1) принадлежит классу функций и(х, у, г) е С4 (О)
и неоднородная правая часть уравнения /(х, у, г) е С2 (О), и пусть наложена связь на коэффициенты
в линейной комбинации производных четвертого порядка от решения уравнения Пуассона: Сумма коэффициентов для несмешанных частных производных равна половине суммы коэффициентов от смешанных производных в линейной комбинации, тогда этого достаточно для преобразования линейной комбинации в сумму частных производных второго порядка от правой части уравнения Пуассона.
Доказательство. Используя уравнение Пуассона, найдём его вторые частные производные от правой части
Г /хх = ^ + ^х + ^х
ихх + иуу + игг
/ (х, у, г) ^
(2) у (2)г
г = и(4) + и(4) + и(4) ¿уу М(2) хиу и(2) у (2)у (2)г
/г = „(2)х + „(2)у + «г4), (2)г (2)г
(13)
где символ и((24))х обозначает частную производную 4 порядка: 2 порядка по переменной х и 2 порядка по у.
(2) у
Тогда для линейной комбинации вторых производных получим с учетом формулы (13)
( Л ( Л (
А/хх + В/уу + С/гг = А
+ и<4> + и(4)
х (2)х (2)х
(2) у (2) г
+ В
„(4) + „(4) + „(4)
(2)х у (2)у
(2) у
+С
(2) г ;
м(4) + м(4) + и4
(2)х (2)у г (2)г (2)г
= А^4 + Ви~у) + Си(4 + м(2) х
(4)
(4)
( а + В) + „(2)х ( а + С) + М(2)у (В + С) . (2) у (2) г (2) г
(14)
Из последней формулы следует, что (А) + (В) + (С) = 1 ((А + В) + (А + С) + (В + С)), то есть сумма
коэффициентов при несмешанных частных производных равна половине суммы коэффициентов при смешанных производных в линейной комбинации производных 4 порядка от решения уравнения Пуассона. Лемма 4 доказана.
Как было показано в Леммах 1-3 аппроксимация трех различных сумм узловых значений решения на 27-точечном шаблоне содержит частные производные от решения только четного порядка, а следовательно, и производные четного порядка от правой части. Следовательно, линейная комбинация производных четвертого порядка и четного порядка по каждой из координат является общим видом разложения в ряд Тейлора для слагаемых 4-й степени по шагам сетки, т.е. утверждение Леммы 4 учитывает общий случай разложения в ряд Тейлора суммы узловых значений для 27-точечного шаблона.
Лемма 5. Пусть решение задачи (1) принадлежит классу функций и(х,у, г)е С6 (О) и неоднородная правая часть уравнения /(х, у, г) е С4 (О), и пусть наложена связь на коэффициенты в линейной комбинации производных шестого порядка от решения уравнения Пуассона: Сумма коэффициентов для несмешанных частных производных и коэффициента при симметричной смешанной производной и((26))х
(2)у (2)г
равна половине суммы коэффициентов от несимметричных смешанных производных шестого порядка. Тогда этого достаточно для преобразования линейной комбинации в сумму частных производных четвертого порядка от правой части уравнения Пуассона.
Доказательство. Используя уравнение Пуассона, найдём его четвертые частные производные от правой части
/4) = 46> + и(61 + и(6)
ихх + иуу + игг = / (х, у, г) ^
/хх = их ) + и(2)х + и(2)х (2) у (2) г
Г = и(4) + и(4) + и(4) ^
/уу и(2)х иу и(2)у (2)у (2)г
/г = и|2)х + и(2)у + иг4)
(2)г (2)г
(4)х (2)у
(4)х (2)г
/ (4) = и(6) + и(б) + и(б) /(2)х и(4)х и(2)х и(2)х
(2) у (2) у (4) у
(2)у (2)г
/ (4) = и(6) + и(6) + и(6) /(2)х и(4)х и(2)х и(2)х
(2)г (2) г
(2)у (2)г
г (4) = и(6) + и(6) + и /у и(2)х иу и(4)у
(4)г (6)
(4)у
(2)г
/(2)у и(2)х и(4)у и(2)у
(2)г
(2) у (2) г
(2)г (4)г
/г<4> = и(61х + и(6) у + игб>.
(2)х (4)г
(4)г
(15)
Тогда для линейной комбинации четвертых производных от правой части уравнения Пуассона получим
( А
(
А/х(4) + В/$Х + + О/^ + + = А
(2) х (2) у
(2)х (2)г
(2) у (2) г
л
и(6 + иг(6) + и%
х (4)х (4)х (2)у (2)г
+В
(6) + (6) + (6) и(4) х + и(2) х + и(2) х
(2)у (4)у
(2) у (2) г
+
+С
(6) + (6) + (6) и(4) х + и(2) х + и(2) х
(2)г
(2) у (2)г
(4)г
+О
и(2)х + иу6 + и(4)у
(4) у
(2)г
+ Е
и(6) + и(6) + и(6) и(2) х + и(4) у + и(2) у
(2) у V (2)г
(2)г (4)г
+ ^
и(2)х + и(2) у + иг } (4)г (4)г
= Аи(6 + Ои(6 + ^гб) + и((46))х (А + В)+ и® (О + В) + и® (С + ^) + и® (С + А) +
(4)х (2) у
(2)х (4) у
(2)х (4) г
(4)х (2)г
+и(26)у (Е + ) + и(46)у (О + Е) + и® (В + С + Е).
(2) у ' (4)г
(4) у 1 (2)г
(2)х (2)у (2)г
(16)
Из последней формулы следует, что
(А) + (О) + (^) + (В + С + Е) = 1 ((А + В) + (О + В) + (С + ^) + (С + А) + (Е + ^) + (О + Е)),
то есть сумма коэффициентов для несмешанных частных производных и коэффициента при симметричной смешанной производной и(2)х равна половине суммы коэффициентов от несимметричных
(2)у (2)г
смешанных производных шестого порядка линейной комбинации от решения уравнения Пуассона. Лемма 5 доказана.
Как было показано в Леммах 1-3 аппроксимация трех различных сумм узловых значений решения на 27-точечном шаблоне содержит частные производные от решения только четного порядка. Поэтому линейная комбинация производных шестого порядка и четного порядка по каждой из координат является общим видом разложения в ряд Тейлора для слагаемых при шестой степени шага сетки. Следовательно, утверждение Леммы 5 учитывает общий случай разложения суммы узловых значений функций для 27-точечного шаблона.
Теорема 1. Пусть решение задачи (1) принадлежит классу функций и(х,у, г)е С6 (О) и неоднородная правая часть уравнения /(х, у, г) е С4 (О), тогда симметричный 27-точечный шаблон на равномерной сетке обеспечивает аппроксимацию 3 мерного уравнения Пуассона на параллелепипеде с шестым
порядком погрешности О (И ). Более того, невозможно получить порядок невязки 27-точечным шаблоном выше шестого. Имеет место разностное уравнение для трёхмерного уравнения Пуассона:
(_ 64 ( ) (
^2 ^ 15 итпк + 15 (ит_1,п,к + ит+1,п,к + ит,п_1,к + ит,п+1,к + ит,п,к_1 + ит,п,к+1) + ю (ит_1,п_1,к + ит_1,п,к_1 + +ит_1,п+1,к + ит_1, п,к+1 + um,n—1,k—1 + ит,п+1,к—1 + um,n—1,k+1 + ит,п+1,к+1 + ит+1,п—1,к + ит+1,п,к—1 + ит+1,п+1,к + ит+1,п,к+1) +
30 (ит_1,п_1,к _
т_1,п_1,к—1 + ит_1,п_1,к+1 + ит_1,п+1,к—1 + ит_1,п+1,к+1 + ит+1,п_1,к—1 + ит+1,п_1,к+1 + ит+1,п+1,к—1 + ит+1,п+1,к+1
+ и
+и
+ ит
+ ит
+ ит
+ ит
) 1 =
/тпк + й2 (-Тхх + /уу + /гг ) + 360
4 (-л4
4 . А
Э/+ЭИ+ЭИ
Эх4 Эу4 Эг4
(4
90
э4 / + э4 / + Э4 /
Эх2Эу2 Эх2Эг2 Эу 2Эг
22
22
+ О (Й
У
(Й6). (17)
Замечание 2. Правая часть формулы (17) за вычетом /тпк равна невязке уравнения Пуассона и в двухмерном случае для прямоугольника переходит в формулу
И
Й2 (/хх + /уу)+ 360
4 ( 4
360
4А
Э/+Э4/
Эх4 Эу4
йй_ 90
4 ( Э4 / А Эх2Эу2
+ О (И6) [4, с. 73].
Доказательство. Заметим, что в силу уравнения Пуассона
Аы(х У, г) ° ихх + иуу + игг = / (х У, г )
операции в правой и левой частях этого уравнения принадлежат одинаковому классу гладкости, и поскольку / (х, у, г)еС4 (О) , то решение имеет класс гладкости на две единицы больше, чем правая часть
уравнения и(х,у, г)е С6 (О).
+
4
А
Оператор Лапласа, как и уравнение Пуассона, является линейным оператором относительно неизвестной функции, поэтому возможно аппроксимировать лапласиан квадратурной формулой линейной относительно узлового значения, используя все 27 узлов шаблона итпк :
Аи = ихх + иуу + игг = 2 (С0итпк + С1 („т-^пД + ит+1,п,к + ит,п-1,к + ит,п+1,к + и\
п,к + ит,п-Хк + ит,п+Хк + „т^к -1 + um'n'k+1
1)-
+С2 („т-^п+^к + ит-1,п,к+1 + ит-Хп-Хк + ит-1,п,к-1 + „т,п-Хк-1 + ит,п~Хк+1 + ит,п+Хк-1 + ^^+1^+1 + +um+1'n-1'k + ит+1,п,к-1 +^+1^+1^ + ит+1,п,к+1) + С3 (ит-1,п-1,к-1 + ит-1,п-1,к+1 + ит-1,п+Хк-1 + ит-1,п+1,к+1 + +^+1^-1^ -1 + „т+^п-^к+1 + ит+1,п+1,к -1 + ^+1^+1^+1)). (18)
Подставим в формулу (18) разложения из Лемм 1-3, т.е. формулы (3), (7), (12) с разложением по степени шага /1 = /2 = /3 = /, 2 Р = 6 включительно, и сгруппируем скобками слагаемые с равной степенью по
/ („(2)х - для краткости будем обозначать, например, частную производную функции и (хт, уп, гк) шесто-
(2) у (2) г
го порядка и по всем координатам х, у, г второго порядка в узле с координатами хт, уп, гк ):
( 1.4 . . 1.2
Аи ( С0итпк + С1 6итпк + /2 (ихх + иуу + игг ) + ( ^ + и<у] + ) + 3г20 („х^ + „у6 + „^ ) + О (^ )
/
/1 I (4) (4) (4) \ / I
их ' + иу/ + и; ; 1+--(+ + иг
12 V х у г ) 36^ х у г
+ С2
12итпк + 4
/2 (ихх+иуу+и77)+/4 — („х4)+„у4))+1 и(4) +—(„х4)+„;4))+1 и(4) +
V хх уу гг) х у ) 4 (2)х х г / 4 (2)х
+ „и 24 х у
и™ +--их + и; +— и,.
4 (2)х 24^ х г > 4 '
(2)у (2)г
+24 (и у4>+-г4')+4
I 48
24 V '"г
+—(иуб) + иг6 | +
и(4)
4 „(2) у 4 (2) г
(
+ kЙ
1
1
,(46)+и(6)+ 720 ^ х у > 48
(6) + (6) „(2)х + „(4) х (4)у (2)у
1
+—и6+„г6)+
720 х г 48
1
(6) + (6) „(2) х + „(4)х (4)г (2)г
- (иуб> + и(6 720^ у г ' 48
(
\\
и(6) + и (6)
„(2)у + и(4)у (4)г (2)г
+ О (/8)
(
+ С3
8итпк + 8
( /2
" I
-2 (их
2у хх + иуу + игг ) + "24 ( ^ + иу4) + иТ' 1 +
4
-44>)-
4
(
(4) + (4) + (4) „(2) х + „(2) х + „(2) у ^ (2)у (2)г (2)г )
Л Л
66 +(„х6)+„уб)+и(6 )+
720^ х у х > 48
(6) + (6) + (6) + (6) + (6) + (6) „(2) х + „(4) х + „(2) х + „(4)х + „(2) у + „(4) у ^ (4)у (2)у (4)г (2)г (4)г (2)г )
6
+ (6) + 8 „(2)х 8 (2) у
(2)г
+ О (К
(/8)
С0 тпк + С1
6итпк + ^ ( „хх + „уу + игг ) + ( 44) + „у^ + ) +
/6 Л ( /4
+ 360 („х6 + "у6) + „г6) ) + О ( /8 ) + С2 12итпк + 4к2 (ихх + иуу + „гг ) + у (+ + ^ ) +
) V
+/4
Л (
и (4) + И(4) + И(4)
(2)х (2)х (2)у _ (2)у (2)г (2)г ^
+
-(„х6) + и(6 + и(6) ) + -9М х у г ) 12
(
\
„(2) + и(6) + „(б) + „(б) + И(6) + „(б)
(2)х (4)х (2)х (4)х (2)у (4)у (4)у (2)у (4)г (2)г (4)г (2)г )
+ О
(//8)
/4
+С3 (8итпк + 4/2 (ихх + иуу + игг ) + Кr (^х^ + и<у) + „г4) ) + 2/4
А
(6) + (6) + (6) + (6) + (6) + (6) „(2) х + „(4) х + „(2) х + „(4) х + „(2) у + „(4) у (4)у (2)у (4)г (2)г (4)г (2)г
Л
(4) + (4) + (4)
„(2) х + „(2) х + „(2) у
V (2)у (2)г (2)г Л
+ К- (и хб) + иУ> + и®) +
+ к6u(:-l + О (/8
(2) х (2) у (2) г
(/8)
итпк
(С0 + 6С1 + 12С2 + 8С3) +
+
+
+
+
4
1
2
/
+
6
/
+
2
6
/
( (А
+ (ихх + иуу + игг )(С + 4С2 + 4С3) + Й2 (-х4) + и!у4 + и^ ) + С2 + С^ +
(
\
+Й4
и(б) + и(б) + и®) ^+^ + С31-х у г 360 90 90 )
С1 +—+ , +
(6) + (6) + (6) + (6) + (6) + (6) и(2) х + и(4) х + и(2)х + и(4)х + и(2) у + и(4) у (4) у (2) у (4)г (2)г (4)г (2)г
(4) + (4) + (4)
и(2)х + и(2)х + и(2) у V (2) у (2) г (2) г
А
12 6
(С2 + 2С3 )
2 + 1 + и® С3 + О(Й6)
2 6 ) (2)х 3 V I
)х (2) у (2) г
(19)
Из формулы (19) следует, что
™|С0 + 6С1 + 12С2 + 8С3| < С0 + 6С1 + 12С2 + 8С3 = 0,
Й Й®0
итпк
то есть получаем первое уравнение связи для слагаемых при Й_2
В задаче аппроксимации трехмерного лапласиана ихх + иуу + игг = / ^ С1 + 4С2 + 4С3 = 1 - полу
чаем второе уравнение связи для слагаемых при Й0 .
2
Согласно Лемме 4 запишем третье уравнение связи для слагаемых при И
3( С-+С32+у У = 23 (С2+2С3)« С1 = 2С2+8С3.
Согласно Лемме 5
четвертое (последнее) уравнение связи при Й4
3!^ + С2 + О. | + С3 = 161 + 1 ^ С1 = 26С2 — 64С3 . V360 90 90| 3 2 \ 12 6 | 1 2 3
Полученная система линейных уравнений имеет единственное решение С0 + 6С1 + 12С2 + 8С3 = 0 С1 + 4С2 + 4С3 = 1 С1 = 2С2 + 8С3 С1 = 26С2 — 64 С3
^^ С0 =--, С1 =-, С2 =-, С3 =-,
0 15 1 15 2 10 3 30
(20)
поэтому первое слагаемое в формуле (19) (С0 + 6С1 + 12С2 + 8С3 ) = 0, в силу уравнения Пуассона
И
второе слагаемое в формуле (19) с учетом системы (20) (ихх + иуу + игг )(С1 + 4С2 + 4С3) = / , используя
Лемму 4, получим третье слагаемое в (19):
( , . (
И
(44)+и(4+44)
) (С
С1 С2 С3 1 ++ — I +
3 3,
2
V
(4) + (4) + (4) и(2) х + и(2)х + и(2) у (2) у (2)г (2)г
Л Л
( С2 + 2С3 )
= Й2 (А/хх + В/уу + С/гг ) =
= 12 ( /хх + /уу + /гг ) , так
как
. с „ С1 С2 С3 7 1 1 15 1
А = В = С = — + — + — =-+ — + — =-= — .
12 3 3 180 30 90 180 12
(21)
Используя Лемму 5, получим четвертое слагаемое в (19):
>) Г с
^ 36(
(-х6) + и(6 + и{6)+ ^ + | +
360 90 90
А/]4) + В/(24))х + С/(24))х + О/(4^ + ЕД4)), + /(4)
'(2) л •'СТ х (2)у (2)г
(6) + (6) + (6) + (6) + (6) + (6) и(2) х + и(4) х + и(2) х + и(4) х + и(2) у + и(4) у V (4)у (2)у (4)г (2)г (4)г (2)г )
44
= (/х(4) + /(4) + /г(4))+ Й_ 360\ х у г ) 90
С2 + С3 ] + и(6) С
^+ У I + и(2)хС3
12 6 У (2)у
(2)г У
'(2)у ' •'у
(2) г
/(4) + /(4) + /(4)
/(3)х /(3)х /(3) у (2) у (2) г (2) г
(22)
+
4
Й
4
Й
дглг-С1 С2 С3 7 1 1 7 + 6 + 2 1 г.^г-г.^г-^ так как А = ^ = F = + —+ — = — + — + =-— = ^7 и В = С = Е, В + С + Е = С3,
360 90 90 5400 900 2700 5400 360
С3 1
В = С = Е = = — .
3 90
Подставляя в левую часть формулы (19) найденные коэффициенты из выражения (20) -
64 ^ 7
11
С0 =-^5", С1 = ^, С2 = —, С3 = —, а в правую часть (19) (невязка аппроксимации трехмерного уравнения Пуассона представляет третье и четвертое слагаемые в правой части (19)) - преобразованные слагаемые через функцию /(х, у, г) - (21), (22), получим следующую формулу:
(- 64 _7. (
2 I г итпк „т-^пк + „т+^пД +
k2V 15 154
ит-Хп,к + ит+1,п,к + ит,п-Хк + ит,п+Хк + ит,п,к-1 + ит,п,к+1) + ю („т-Цп-Цк + ит-Хп,к-1 +
+uш-1'n+1'k + Щт-ХпЛ+Х + ит,п-1,к-1 +uш,n+1'k-1 + ит,п-1,к+1 + „т^!^! + uш^1м-1'k + итЛ,п,к-1 + um+1'n+1'k + um+1'n'k+1) +
+ 30 ( ит-1,п-1,к-1 + ит-1,п-1,к+1 + ит-1,п+1,к-1 + ит-1,п+1,к+1 + ит+1,п-1,к-1 + ит+1,п-1,к+1 + ит+1,п+1,к-1 + ит+1,п+1,к+1 =
/тпк + /2 (/хх + /уу + ./гг ) + 360
4 (Э4 / +Э4/+э4 / Л
Эх4 Эу4 Эг4
(4
90
4Л
Э4/ + Э4/ + Э4/
Эх2Эу2 Эх2Эг2 Эу2Эг2
+ О (/
(к-), (23)
которая совпадает с формулой (17). Из формулы (23) следует, что аппроксимация трехмерного уравнения Пуассона на параллелепипеде 27-точечным симметричным шаблоном имеет шестой порядок погрешности. Всего неизвестных коэффициентов в шаблоне (18) четыре. Тогда добавление условий в систему уравнений (20), аналогично условиям связи в Леммах 4, 5, приведет к несовместности линейной системы уравнений, в которой число неизвестных 4, а число уравнений не менее 5. Другими словами, 27-точечным шаблоном аппроксимировать трехмерное уравнение Пуассона с 8 порядком погрешности невозможно.
Теорема 1 доказана.
Из формулы (23) выразим центральное узловое значение с индексами (т, п, к) и получим формулу простой итерации:
5+1 (15Л( 7
итпк =
64
V15 (
и т-Хп,к + и т+^пк + и т,п-Хк + и m'n+1'k + и т,п,к-1 + и m'n'k+1 ) +
+ 10 I„ т-1,п-1,к + и т-1,п,к-1 + и т-Хп+Хк + и т-ХпЛ+1 + и т,п-Хк-1 + и т,п+Хк-1 + и m'n-1'k+1 + и m'n+1'k+1 + +и т+^п-^к + и ш+Хп,к-1 + и ш+1'n+1'k + и ш+1'n'k+1) + 30 (и ш-1'n-1'k-1 + и т-Хп-Хк+1 + и т-Хп+Хк-1 + и ш-1'n+1'k+1 +
+и ш+1,и-1,к-1 + и т+1,и-1,к+1 + и т+1,п+1,к-1 + и т+1,п+1,к+1) / /тпк ( /хх + /уу + /гг ) 360
1.6 ( ^4
6 (Э^/ Э4/ Э4/Л
.4 + Э .4 + Э_4 )
Н_ 90
3 128 *
Э4/ + Э4/ + Э4/ л Эх2Эу2 Эх2Эг2 Эу2Эг2
Эх Эу Эг
+ О ( к6 )= £ („',„.,„< + ✓ + и"...,, + и-т^ + ✓ + и'„к+1) +
+128(„ m-1'n-1'k + и т-1,п,к-1 + и m-1'n+1'k + и ш-^.^! + и т,п-1,к-1 + и т,п+1,к-1 + и m'n-1'k+1 + и m'n+1'k+1 + ш+1,п-1,к + и ш+1,п,к-1 + и ш+1,п+1,к + и ш+1,п,к+1) + (и т-1,п-1,к-1 + и ш-1,п-1,к+1 + и ш-1,n+1,k-1 + и ш-1,n+1,k+1 +
т+1,п-1,к-1 и т+1,п-1,к+1 и т+1,п+1,к-1
1 + и m+1'n+1'k+1) 64 / /тпк 256 (/хх + /уу + /гг )
5/ 256
(.6 ( э4/ Э4/ Э4/ Л ^ ( Э4
1536
„■ +-1Т+ „
Эх4 Эу4 Эг4
Н_ 384
Э4 / Э4 / Э4 / Л
Эх2Эу2 +Эх2Эг2 +Эу2Эг2
+ О (/
(/8).
(24)
4
/
+
1
5
В формуле простой итерации (24) верхний индекс я обозначает номер итерации, в левой части (24) индекс итерации я = 1 на единицу больше, чем индексы итерации я у всех узловых значений в правой части формулы (24). В формуле (24) значения функции /(х,у,г) и ее частные производные вычисляются в узле с индексами (т, п, к). Формулу простой итерации (24) можно записать в виде
ит+пк = С (итпк, /тпк ) + О ( й8 ), я = ^и— т = 1 п1 —1 п = 1 п2 —1, к = 1, п3 — 1 (25)
где функция С (ит«к, /тпк) полностью совпадает с правой частью формулы (24).
По аналогии с оператором С(и^тк, /тпк) в формуле (25), рассмотрим итерационную последовательность (26) со сжимающим коэффициентом 0 < q < 1, индуцируем следующий вспомогательный оператор С(дитпк, /тпк), если в ее правой части провести замену итпк ® уитпк:
ит+к = С ( ^пк, /тпк ) + О ( й8 ), я = ^и— т = 1 п1 — 1 п = 1 п2 —1, к = 1, п3 — 1 (26)
Введем функцию нормы в пространстве сеточных функций итпк е Р(п1+1) х Р( +1) х Р(п3 +1)
5я =
я пит ехас1
и — и
п1 ,п3 С
= тах
0<т<пу 0<п<«2 0<к <п3
ехас1
т,п,к (и )т,п,к
я пит
ит„,г —(и)
, я = 0,1,2,..., совпадающей с метрической функци-
ей
Л (, пит , ехасЛ \
р (и , и )-
пространстве
, пит , ехас1
и — и
итпк е Р(п1+1) хР(1 +1) хР(«3 +1), определяемой формулой (А.Н. Колмогоров, С.В. Фомин) [8, с. 139], где и'Пт, {иехас*) обозначает соответственно численное решение задачи (27) в узле (т, п, к) на итерационном шаге (я) и след [9, с. 102] точного решения задачи (1) на узел сетки (т, п, к).
Описание численного алгоритма
Введем на параллелепипеде для задачи (1) равномерную сетку с равным шагом
Ь _ а й _ с г _ е Й =-= Й2 =-= Й3 =-= И:
п1 п2 п3
О = [а,Ь]х[с,й]х[е,г], ю«^« ={хт = а + Йт, т = 0,«1, уп = с + Йп, п = 0,п2, гк = е + Йк, к = 0,п3} .
Для сеточной функции ит+Пк запишем систему разностных уравнений, соответствующих задаче в частных производных (1):
<°.к = С ( 9ит,п,к, /т,п,к ) + О ( й8 ) (24Х я = 0,1,2,..., т = 1 п1 —1, п =1, «2 _1, к = 1, п3 _1
лт,п,к, /т,п,к , я+1
Ч+а =т1,0,п,к, С,к =т2,п1,п,к,« = 0, n2, к = 0, т^я = 0,1,2.
ит+0,к =т3,т,0,к,um+1lз,k = , т = 0,n1, к = 0n3, я = 0,1,2.<
(27)
ит+°,0 = т5,т,п,0, и1+° т = т6,т,п,п3 , т = 0, n1, п = 0 n2, я = 0,1,2
т,«2 я+1
т,п,«3
п3
Замечание 3. В разностной задаче (27) # ® 1, иначе изменяются коэффициенты точной аппроксимации (20). Достаточно, чтобы # = 1 _10—15, т.е. # может отличаться от 1 на число, которое в 10 раз большее машинной ошибки округления.
В системе разностных уравнений (27) первое уравнение совпадает с записью формулы (26). За начальное решение итерационного процесса можно выбрать функцию, определяемую граничными условиями [2, с. 428]:
0
ит,п,к
= 1 (т1,0,п,к + т2,п1,п,к + т3,т,0,к + т4,т,п2,к + т5,т,п,0 + т6,т,п,п3 ),т = 1 п1 _ 1,п = 1,«2 _ 1,к = 1,«3 _ 1 (28)
Из формулы (28) видно, что начальное численное решение непрерывно зависит от всех 6 краевых функций как их линейная комбинация, и в классической постановке задачи (1) не может быть разрывов первого рода для начального приближения решения и всех значений индексов (m, n, k).
Теорема 2 (об устойчивости итерационных формул (23), (26). Итерационная последовательность (26) ufmUk = G (qusmnk, fmnk) + O (h8) является сжимающим отображением в метрическом пространстве
usmnk £ R(П1+1) х R(n2 +1 xR(П3 +1, то есть для любой итерационной последовательности численных решений (26) существует единственное предельное значение lim usmnk, равное точному решению задачи (1),
q®1-0 h®0
удовлетворяющее разложению (23): lim usmnk = lim (uexact) , "m = Ö^, n = Ön", k = 0^, lim (uexact) = G | lim (uexact) , lim fmnk
s®+¥ h®0^ 'mnk h®0^ 'mnk Vh®0^ 'mnk h®0
q®1-0 h®0
Доказательство. Покажем, что итерационная последовательность (26) является сжимающим отображением.
Обозначим ?>sm,n,k = y"m,n,k - xSm,n,k, "m = 0, n1, n = 0 n2, k = 0 n3 . Рассмотрим 2 произвольные сеточные функции и получим:
"xm+k, ym+1 £ R(ni+1)x R(n2 +1)x R(n3+1) : p( x^, y^ ) = p( G ( a^nk, fmnk ) , G ( qyLk, fmnk ))« |xm+i,k - yJS+t,
G ( qxmnk, fmnk ) - G ( aymnk, fmnk )
I (26)
(24) (la\
— 7 a °s + °s + ds + ds + ds + ds
— I - ■ 1° m-1,n,k +0 m+1,n,k +0 m,n-1,k + 0 m,n+1,k +0 m,n,k-1 +0 m,n,k+11
3q 128
641
°s + °s + °s + °s + °s + °s + °s + °s
° m-1,n-1,k +° m-1,n,k-1 +° m-1,n+1,k +° m-1,n,k+1 +° m,n-1,k-1 +° m,n+1,k-1 +° m,n-1,k+1 +° m,n+1,k+1 '
s
s
s
m+1,n-1,k ^ u m+1,n,k -1 ^ u m+1,n+1,k ^ u m+1,n,k+1
s
128
s
s
m-1,n-1,k-1 ^ u m-1,n-1,k+1 ^ u m-1,n+1,k -1
+ °Sm-U+1,k+1 + °Sm+1,n-1,k-1 + °Sm+1,n+1,k-1 + °Sm+1,n+1,k-1 + °Sm+1,n+1,k+1 ) - max
' 0-n-n1
0—m—n2 0—k—n3
°S
s
m,n,k
(7q • 6 3a-12 a • 8 ^ (84 + 36 + 8
xl —-— + —-+ -— I = q I-I • max
I 64 128 128 J V 128 J 0—n—n1
0—m—n2 0—k—n3
°s+1
° m,n,k
= a max
0—n—n1
m,n,k
= ap(xmnk, ySmnk ), 0 < a < ^
что равносильно
max
0—n—n 0—m—n2 0—k—n3
°s+1
° m,n,k
— a max
0—n—n1 0—m—n2 0—k—n32
0—m—n2 0—k—n32
m,n,k
(29)
при достаточно малом шаге й, то есть р (41°, у^) < (хтпк, ут«к).
Применяя определение 1 [8, с. 74] по А.Н. Колмогорову, получим, что последовательность (26) является сжимающей и сходится к единственному пределу, который обозначим и*« (теорема - принцип сжимающих отображений [8, с. 75])
s * * / * \/8\
lim umnk ° umnk : umnk = G ( aumnk, fmnk ) + O ( h )
+
a
+
x
Но точное решение удовлетворяет формулам (23), (25), из которых получим предельный переход lim (uexact) = GI lim (uexact) , lim fmnk I, где оператор G(qusmnk, fmnk ) является линейным и непре-
h®o\ 'mnk Vh®0^ 'mnk h®0 ) V '
рывным в формуле (23) по umnk и по q. Далее перейдем к следующему пределу в (30) по переменным q и h:
( \
* . / * \ lim Umnk = G (qumnk, fmnk )
q®1-0 q®1-0 v '
h®0 h®0
=G
lim „ umnk, fmnk q®1-0 h®0
^ lim umnk = lim (uexact)
q®1-0 mnk h®0^ )
q®1-0 h®0
mnk
Теорема 2 доказана.
Рассмотрим тестовый пример
uxx + uyy + uzz = sin(x)sin( y)sin(z), x e (0, я), y e (0, я), z e (0, я) u(0, y, z) = u(p, y, z) = sin(y)sin(z),( y, z) e [0, я]х[0, я] u(x,0, z) = u(x, я, z) = sin(x)sin(2z),(x, z) e [0, я]х [0, я] u(x, y,0) = u(x, y, я) = sin(3x) sin( y),(x, y) e [0, я] х [0, я].
(31)
Постановка задачи в примере (31) является классической, так как все шесть краевых функций непрерывны не только на своей грани куба, но и на всех 12 ребрах куба равны друг другу (нулю).
Задача (31) является линейной относительно неизвестной функции и(х, у, г), поэтому воспользуемся методом редукции линейной задачи и( х, у, г) = „х( х, у, г) + „2( х, у, г) + „3( х, у, г) + „4( х, у, г), где функции „1 (х, у, г), „2 (х, у, г), „3 (х, у, г), „4 (х, у, г) являются решениями соответственно следующих частных задач:
1)
2)
3)
u1xx + u1yy + u1zz = sin(x)sin( y)sin(z), x e (0, я), y e (0, я), z e (0, я)
u1 (0, y, z) = u1 (я, y, z) = 0,(y, z) e [0, я] х [0, я]
%(x,0, z) = %(x, я, z) = 0,(x, z) e [0, я] х [0, я]
u1 (x, y, 0) = u1 (x, y, я) = 0, (x, y) e [0, я] х [0, я]
u2xx + u2yy + u2zz = 0, x e (0, я), y e (0, я), z e (0, я)
u2(0, y, z) = u2 (я, y, z) = sin( y)sin(z),( y, z) e [0, я]х[0, я]
u2 (x,0, z) = u2 (x, я, z) = 0,(x, z) e [0, я] х [0, я]
u2 (x, y, 0) = u2 (x, y, я) = 0,(x, y) e [0, я] х [0, я]
u3xx + u3yy + u3zz = 0, xe (0,я), y e (0,я), z e (0,я)
u3 (0, y, z) = ^(я, y, z) = 0,( y, z) e [0, я] х [0, я]
u3(x,0, z) = u3(x, я, z) = sin(x)sin(2z),(x, z) e [0, я] х[0, я]
u3 (x, y, 0) = u3 (x, y, я) = 0,( x, y) e [0, я] х [0, я]
u4xc + u4yy + u4zz = 0, x e (0, я), y e (0, я), z e (0, я)
u4 (0, y, z) = ^(я, y, z) = 0,( y, z) e [0, я]х[0, я]
u4 (x,0, z) = u4 (x, я, z) = 0,(x, z) e [0, я] х[0, я]
u4(x, y,0) = u4(x, y,я) = sin(3x)sin(y),(x, y)e [0,я]х[0,я].
Ищем решение задачи 1) в виде их (х, у, г) = С sin(х^ш( у^ш(г), которое удовлетворяет однородным краевым условиям системы уравнений 1):
их(0, у, г) = и^я, у, г) = их( х,0, г) = х, я, г) = „х( х, у,0) = их( х, у, я) = 0.
Подставим х, у, г) в первое уравнение системы 1):
1
1 .
„1хх + „хуу + „хгг = -3Сап(х)ап(у)ап(г) = sin(x)sin(y)sin(г) ^ С = --, „х(х, у, г) = --sin(x)sin(y)sin(г).
Ищем решение задачи 2) в виде „2( х, у, г) = f (x)sin( y)sin( г), которое удовлетворяет однородным краевым условиям системы уравнений 2): „2(х,0, г) = „2(х, я, г) = „2(х, у,0) = „2(х, у, я) = 0. Подставим „2( х, у, г) = f (х^т( y)sin( г) в первое и второе уравнения системы 2):
„2хх + „2уу + „2гг = f (x)sin(у^ш(г) - 3f(x)sin(у^ш(г) = 0 ^ f (х) - 2f (х) = 0 и2(0, у, г) = f (0)sin(у) sin(г) = sin(уг),"(у, г) е [0, я] х[0, я] ^ f (0) = 1 „2 (я, у, г) = f (я^т(y)sin(г) = sin(y)sin(г),"(у, г) е [0, я]х [0, я] ^ f (я) = 1.
Общее решение краевой задачи с обыкновенным дифференциальным уравнением (ОДУ) второго порядка есть
f" (х) - 2 f (х) = 0, f (х) = АзЪ (ч/2х) + ВеЬ (л/2х), f (0) = 1 ^ В = 1, f (я) = 1 ^ 1 = АэЪ (Т2я) + еЬ (Т2я),
\ - еЬ (л/2я)4
1 - еЬ (Т2я)
А =-, V. , , f(х) = еЬ(л/2х)
sЬ (\/2я)
„2( х, у, г) =
еЬ (л/2х) +
1 - еЬ (72я) sЬ (л/2я)
sh (>/2я) sЬ (72х)
sh (-Лх),
sin( y)sin( г).
Ищем решение задачи 3) в виде „3(х,у,г) = f(у^т(x)sin(3г), которое удовлетворяет однородным краевым условиям системы уравнений 3):
„3 (0, у, г) = „3 (я, у, г) = „3 (х, у, 0) = „3 (х, у, я) = 0. Подставим „3(х, у, г) = f (у^ш( х^т(2г) в первое и третье уравнения системы 3):
и3хх + „3 уу + „3 гг
= f (у)яп(х)яп(2г)-5f(у^ш(х)яп(2г) = 0f (у)-5f(у) = 0 • и3 (х, 0, г) = f (0) sin(х) sin(3г) = sin(х) sin(3г),"(х, г) е [0, я] х[0, я] ^ f (0) = 1 „3(х, я, г) = f (я^ш(x)sin(3г) = sin(х^т(2г),"(х, г) е [0, я]х [0, я] ^ f (я) = 1.
Общее решение краевой задачи с ОДУ второго порядка есть
f" (у) - 5 f (у) = 0, f (у) = А sh (45х) + В еЬ (45х), f (0) = 1 ^ В = 1, f (я) = 1 ^ 1 = А sh (л/5я) + еЬ (л/5я)
1-еЬ(Т5я)
А =-г^-1, f(у) = еЬ(л/5у)-
sЬ (л/5я)
„3 (х, у, г) =
еЬ (75 у) +
1-еЬ(Т5я)
sh (Т5я)
- еЬ (75я) sь (Т5я)
sь (75 у),
sь (75 у)
sin( х) sin(3 г).
Ищем решение задачи 4) в виде „4(х,у,г) = f(г^т(3х^т(у), которое удовлетворяет однородным краевым условиям системы уравнений 3):
„4(0, у, г) = „4 (я, у, г) = „4( х,0, г) = „4( х, я, г) = 0.
Подставим и4 (х, у, г) = / (г)8т(3х)8т( у) в первое и четвертое уравнения системы 4):
-4хх + и4уу + и4гг = / (г)8ш(3х)8ш(у) — 10/(г)8ш(3х)8т(у) = 0 ^ / (г) —10/(г) = 0 • и4(х, у,0) = /(0) 8ш(3х) 8Ш( у) = 8ш(3х)8ш( у),"(х, у) е [0, я] х [0, я] ^ /(0) = 1 и4 (х, у, я) = /(я) 8ш(3х) 8ш( у) = 8т(3х) 8ш( у),"(х, у) е [0, я] х [0, я] ^ /(я) = 1.
Общее решение краевой задачи с ОДУ второго порядка есть
/" (г) —10/(г), /(г) = А8Ь (<Л0г) + ВеЬ (>/10г), /(0) = 1« В = 1, /(я) = 1«1 = А8Ь (-Л0я) + еЬ (>/10я)
(1 — еЬ (у/Щ
1 — еЬ Шя) . _ .
А =-, _ , ;, /(х) = еЬ (у/10г)-
8Ь (л/Шя)
и4 ( х, у, г) =
еЬ ^л/10г) +
1 — еЬ(75я) 8Ь (л/5я)
8Ь (л/10я) 8Ь (у/10 г)
8Ь (>Я0г),
8Ш(3х)8Ш( у) .
Решение примера (31) равно (в силу линейной редукции) сумме решений систем уравнений 1) - 4): и(х, у, г) = -1(х, у, г) + и2(х, у, г) + и3(х, у, г) + и4(х, у, г) = — х)8ш(у)8ш(г) +
еЬ (л/2х) +
1 — еЬ (л/2я) 8Ь (>/2я)
8Ь (л/2х)
еЬ (л/10г) +
8ш( у)8т( г) +
1—еЬ(75я)
8Ь (Т5я)
еЬ(75 у) + 8Ь (710 г)
1 — еЬ (л/5я) 8Ь (Т5я)
8Ь (л/5 у)
8т( х)8т(2 г) +
8Ш(3х)8Ш( у).
(32)
В таблице показаны результаты численного решения примера 1 с параметрами программы «1 = 60, «2 = 60, «3 = 60, т = 13500 (первые три столбца координаты - х, у, г; четвертый и пятый столбцы -численное и точное решение соответственно).
Таблица
х У 2 питепеа1 ехае!
1 2 3 4 5
0.261799387799149 0.261799387799149 0.261799387799149 0.19326654098 0.19326654078
1.30899693899575 0.261799387799149 0.261799387799149 0.18328544131 0.18328544145
2.35619449019234 0.261799387799149 0.261799387799149 0.28564460479 0.2856446045
0.261799387799149 1.30899693899575 0.261799387799149 0.4608673730 0.46086737233
1.30899693899575 1.30899693899575 0.261799387799149 -0.28784462945 -0.28784462878
2.35619449019234 1.30899693899575 0.261799387799149 0.35457814658 0.35457814586
0.261799387799149 2.35619449019234 0.261799387799149 0.35373212267 0.35373212214
1.30899693899575 2.35619449019234 0.261799387799149 -0.14968371807 -0.14968371760
2.35619449019234 2.35619449019234 0.261799387799149 0.30424860387 0.30424860333
0.261799387799149 0.261799387799149 1.30899693899575 0.22895617234 0.22895617228
1.30899693899575 0.261799387799149 1.30899693899575 0.242830568967 0.242830568968
2.35619449019234 0.261799387799149 1.30899693899575 0.23202215878 0.23202215871
0.261799387799149 1.30899693899575 1.30899693899575 0.59406289056 0.59406289039
1.30899693899575 1.30899693899575 1.30899693899575 -6.56191877Е-002 -6.56191875Е-002
2.35619449019234 1.30899693899575 1.30899693899575 0.15445645365 0.15445645348
0.261799387799149 2.35619449019234 1.30899693899575 0.45123800883 0.45123800870
1.30899693899575 2.35619449019234 1.30899693899575 1.299659600Е-002 1.299659609Е-002
2.35619449019234 2.35619449019234 1.30899693899575 0.15774935695 0.15774935681
0.261799387799149 0.261799387799149 2.35619449019234 -1.683548459Е-002 -1.683548468Е-002
1.30899693899575 0.261799387799149 2.35619449019234 -0.57131703876 -0.57131703855
2.35619449019234 0.261799387799149 2.35619449019234 -0.35628466793 -0.35628466797
0.261799387799149 1.30899693899575 2.35619449019234 0.45799547075 0.45799547037
+
+
Окончание таблицы
1 2 3 4 5
1.30899693899575 1.30899693899575 2.35619449019234 -0.18843379843 -0.18843379794
2.35619449019234 1.30899693899575 2.35619449019234 9.32515 8014E-002 9.325157979E-002
0.261799387799149 2.35619449019234 2.35619449019234 0.30256837412 0.30256837385
1.30899693899575 2.35619449019234 2.35619449019234 -0.26000946507 -0.26000946466
2.35619449019234 2.35619449019234 2.35619449019234 -2.109387609E-002 -2.109387631E-002
С нормой невязки \\unum -uexact\\C = 1.0854397936022005e-009 (« = «2 = «3 = 60), где равномерно непрерывная норма разности двух сеточных функций численного и аналитического решений имеет порядок 10-9 и определяется по формуле \\unum -uexact\\c = max ыП^к -uTrnk , т.е. программа работает
c 0<n<n1 ' ' ' '
0<m<«2 0<k <n3
с использованием алгоритма согласно формулам (19)-(28) с относительной точностью REAL (4), лучше чем 10-8 [10].
При «1 = 10, «2 = 10, «3 = 10, m = 6000 \\unum -uexact\\c = 5.019113151139010E-005, а при «1 = 20, «2 = 20, «3 = 20, m = 6000 (уменьшение всех шагов сетки в 2 раза) \\u«um -uexact\\c = 7.891130154868975E-007 , то алгебраический порядок погрешности p равен 6 при использовании формул (25), (27) - q = 1, так как норма невязки уменьшается в 63,4 раз: II и 1
||u«um-uexact\\c 5.019113151139010E-005 _ . 06 ОР „ - - 63.4 » 2 = 2F, p = 6.
u - u II 1/2 7.891130154868975E-007
\\il«um 11 exact ||c
При использовании алгоритма (26), (27) q = 1 - eps = 1 - 5 -10-8 = 0.99999995
1
uum-uexact\\c 5.013128145980872E-005 пл n o6 оР „ c - -71.7 » 2 = 2p, p = 6.
II,. - u II 1/2 6.993329699778172E-007
11 "Mum Mexact | lc
Более высокое значение для порядка погрешности во втором случае объясняется экспоненциальным затуханием ошибки округления в правой части формулы (27) в слагаемых вида usm_1«_1k-1,— c параметром q = 0.99999995.
Лемма 6. Пусть начальная норма разности численного и точного решения конечна, тогда после s шагов итерационной формулы (26) конечная норма затухает по экспоненциальному закону от s и справедлива оценка
d 0 exact _„ 0, «um exact D , • "0 „„/'„\
R0 = um,«,k - um,«,k = max Vm m - u«,m,k , Rs = lim-exP (-es), q = 1 (33)
c 0<«<«1 e®0 e
0<m<«2 0<k <«3
Доказательство. Согласно теореме о неподвижной точке [8, с. 75] итерационная последовательность (26) имеет предел, и после s итераций получим оценку погрешности [8, с. 75]:
-es
R qs R qs R (1 -epe R0 , N
Rs < , q ® 1 - 0 <^e = 1 - q ®+0, Rs < R0q- = ---> -°exp (-es ) =
1 - q 1 - q e e®o e
s= -1 ,,d>0 „ / -d \ z^-1
Rq , , e1+d R0exP (-e d) R0 ed
= —exp (-es) = ---- =-0 -
e e
eexp, , ' ed
R z1/d R z[1/d]+1 R
Ro = z Ro <_Z_—® 0. (34)
1 exP (z)e®oexp (z)- z[1/d]+1 z[1/d]+2 z®¥
1/d expV Z)Z®¥--1--
Z1 d ([1/d ]+1)! ([1/d ]+2)!
Лемма 6 и формула (33) доказаны.
Замечание 4. Из формулы (33) следует, если s e = 1, Rs ® sR^e 1, то ошибка округления в алго-
e®0
ритме (26) как и в (25) прямо пропорциональна числу итераций s. При выполнении условий в алгоритме (26): e = 1 -q ® 0, se = e-d s(e) = —+d-> d >0, Rs ® 0, Rs мала.
e1+d e®0
Быстродействие алгоритма в секундах составляет для ni = 30, П2 = 30, пз = 30, m = 6000
II II —8
(t2 - ti = 12.85448 с, с нормой погрешности \\unum — uemct\\c = 6-10 ). В то время как авторы работы [11]
получили время счета 15 минут (п1 = 30, п2 = 30, п3 = 30) с алгоритмами вида O (й4) и погрешность
в процентах, т.е. с точностью порядка 10—2.
Программа написана на языке FORTRAN и представлена ниже с использованием всех функций и переменных двойной точности, также занесены данные из тестового примера 1.
module stolb;use dfimsl;integer(8),parameter:: n1=30,n2=30,n3=30,m=6000;end module program one;use dfimsl;use stolb;real(8)::a,b,c,d,e,r,h1,h2,h3,ss integer(8)::i,j,k;real(8)::res(0:n1+2,0:n2+2,0:n3+2)
real(8)::res0(0:n1+2,0:n2+2,0:m+2),l(n1+2),nu(n1+2),res1(0:n1+2,0:n2+2,0:n3+2,0:m+1),max,delta(0:n1+2,0:n2+2,0:n3+1)
real(8)::x,y,z,t,ch,sh,f,u1x,u2x,u1y,u2y,u1z,u2z,du,d4u,d4xyzu,pi,h22,h44,h66
ch(x)=(dexp(x)+dexp(-x))/2d0;sh(x)=(dexp(x)-dexp(-x))/2d0;f(x,y,z)=dsin(x)*dsin(y)*dsin(z);
u1x(y,z)= dsin(y)*dsin(z);u2x(y,z)= dsin(y)*dsin(z);u1y(x,z)=dsin(x)*dsin(2d0*z)
u2y(x,z)=dsin(x)*dsin(2d0*z);u1z(x,y)=dsin(3d0*x)*dsin(y)
u2z(x,y)=dsin(3d0*x)*dsin(y);du(x,y,z)=-3d0*dsin(x)*dsin(y)*dsin(z)
d4u(x,y,z)=3d0*dsin(x)*dsin(y)*dsin(z);d4xyzu(x,y,z)=3d0*dsin(x)*dsin(y)*dsin(z)
a=0d0;c=0d0;e=0d0;pi=2d0*dasin(1d0);b=pi;d=pi;r=pi
h1=(b-a)/dfloat(n1);h2=(d-c)/dfloat(n2);h3=(r-e)/dfloat(n3)
max=-1d3;do k=0,n3;z=e+h3*dfloat(k);do j=0,n2;y=c+h2*dfloat(j)
do i=0,n1;x=a+h1*dfloat(i);call uu(x,y,z,t);res(i,j,k)=t;if(mod(i,5)==0.and.mod(j,5)==0.and.mod(k,5)==0)then;endif if(max<dabs(res(i,j,k)))then;max=dabs(res(i,j,k));endif;
enddo;enddo;enddo;do i=1,n1-1;x=a+h1*dfloat(i);do j=1,n2-1;y=c+h2*dfloat(j);do k=1,n3-1;z=r+h3*dfloat(k)
res0(i,j,k)=(u1x(y,z)+u2x(y,z)+u1y(x,z)+u2y(x,z)+u1z(x,y)+u2z(x,y))/6d0
enddo;enddo;enddo;do s=0,m;do j=0,n2;do k=0,n3;y=c+h2*dfloat(j)
z=e+h3*dfloat(k);res1(0,j,k,s)=u1x(y,z);res1(n1,j,k,s)=u2x(y,z)
enddo;enddo;do i=0,n1;do k=0,n3;x=a+h1*dfloat(i);z=e+h3*dfloat(k)
res1(i,0,k,s)=u1y(x,z);res1(i,n2,k,s)=u2y(x,z);enddo;enddo;do i=0,n1;do j=0,n2;x=a+h1*dfloat(i);y=c+h2*dfloat(j); res1(i,j,0,s)=u1z(x,y);res1(i,j,n3,s)=u2z(x,y);enddo;enddo;enddo h22=h1*h1;h44=h22*h22;h66=h44*h22;do s=0,m-1;do k=1,n3-1;z=e+h3*dfloat(k); do j=1,n2-1;y=c+h2*dfloat(j);do i=1,n1-1;ss=0d0;x=a+h1*dfloat(i)
ss=(7d0/64d0)*(res1(i-1,j,k,s)+res1(i+1,j,k,s)+res1(i,j-1,k,s)+res1(i,j+1,k,s)+res1(i,j,k-1,s)+res1(i,j,k+1,s)) ss=ss+(3d0/128d0)*(res1(i-1,j+1,k,s)+res1(i-1,j,k+1,s)+res1(i-1,j-1,k,s)+res1(i-1,j,k-1,s)+res1(i,j-1,k-1,s)+res1(i,j-1,k+1,s)+res1(i,j+1,k-1,s)+res1(i,j+1,k+1,s)+res1(i+1,j-1,k,s)+res1(i+1,j,k-1,s)+res1(i+1,j+1,k,s)+res1(i+1,j,k+1,s))
ss=ss+(1d0/128d0)*(res1(i-1,j-1,k-1,s)+res1(i-1,j-1,k+1,s)+res1(i-1,j+1,k-1,s)+res1(i-1,j+1,k+1,s)+res1(i+1,j-1,k-1,s)+res1(i+1,j-1,k+1,s)+res1(i+1,j+1,k-1,s)+res1(i+1,j+1,k+1,s)) ss=ss-(15d0/64d0)*f(x,y,z)*h22-(5d0*h44/256d0)*du(x,y,z)-(h66/1536d0)*(d4u(x,y,z))-(h66/384d0)*(d4xyzu(x,y,z))
res1(i,j,k,s+1)=ss;enddo;end do;enddo;enddo;max=-1d3;do k=0,n3;do j=0,n2;do i=0,n1 delta(i,j,k)=dabs(res(i,j,k)-res1(i,j,k,m));if(mod(i,2)==0.and.mod(j,2)==0.and.mod(k,2)==0)then endif;if(max<=delta(i,j,k))then;max=delta(i,j,k);endif;enddo;enddo;enddo; print*,"norma C",max; end program one;
subroutine uu(x,y,z,t);real(8): :x,y,z,t,pi,c1,c2,c3,c4,ch,sh
ch(x)=(dexp(x)+dexp(-x))/2d0;sh(x)=(dexp(x)-dexp(-x))/2d0
pi=2d0*dasin(1d0);c1=-dsin(x)*dsin(y)*dsin(z)/3d0
c2=dsin(y)*dsin(z)*(ch(x*dsqrt(2d0))+sh(x*dsqrt(2d0))*(1d0-ch(pi*dsqrt(2d0)))/sh(pi*dsqrt(2d0))) c3=dsin(x)*dsin(2d0*z)*(ch(y*dsqrt(5d0))+sh(y*dsqrt(5d0))*(1d0-ch(pi*dsqrt(5d0)))/sh(pi*dsqrt(5d0))) c4=dsin(3d0*x)*dsin(y)*(ch(z*dsqrt(10d0))+sh(z*dsqrt(10d0))*(1d0-ch(pi*dsqrt(10d0)))/sh(pi*dsqrt(10d0))) t=c1+c2+c3+c4;end subroutine
Особую сложность в методах решения имеют дифференциальные уравнения с разрывами первого рода в краевых и начальных условиях [12, 13, 18], а также с разрывами первого рода в правой части
уравнения Пуасона [19, 20]. Предложенный в работе алгоритм (26) с формулой (24) устойчив по отношению к разрывам первого рода в краевых условиях, например, к прибавлению постоянной к краевой функции одной грани (точки разрыва первого рода краевых условий появляются на все 4 ребрах данной грани).
В работе получены результаты:
1. Предложен минимальный симметричный 27-точечный шаблон для аппроксимации трёхмерного уравнения Пуассона на параллелепипеде.
2. Доказано, что на 27-точечном шаблоне и равномерной сетке невозможна аппроксимация уравнения Пуассона выше шестого порядка.
3. В леммах 1-3 доказано, что коэффициенты разложения в ряд Тейлора суммы узловых значений в группах узлов, расположенных симметрично и на одинаковом расстоянии от центра куба, содержат только частные производные четного порядка и четного порядка по каждой из трех переменных.
4. В леммах 4, 5 получено достаточное условие на коэффициенты аппроксимации уравнения Пуассона линейной квадратурной формулы, как следствие уравнения Пуассона, по одному уравнению на каждую четную степень шага сетки.
5. Найдена явная разностная схема аппроксимирующая трехмерное уравнение Пуассона на параллелепипеде с шестым порядком погрешности, доказана ее устойчивость.
6. Численно показано с использованием точного аналитического решения тестирующего примера (28) и программы, что явная формула простой итерации (24) с алгоритмом (26) имеют шестой порядок погрешности.
7. Итерационная формула (25) устойчива относительно ошибки округления, исключая ее экспоненциальный рост, но допускает линейный рост ошибки округления. В то время как сжимающий алгоритм (27) и формула (26) подавляют даже линейный рост ошибки округления. Разность между начальной итерацией задачи (1) и ее точным решением падает экспоненциально. Алгоритм (27) обеспечивает единственность и существование численного решения задачи (1), совпадающего с точным при любой начальной допустимой итерации в пространстве сеточных функций. Он минимизирует ошибку до невязки аппроксимации О (й8).
ЛИТЕРАТУРА
1. Самарский, А. А. Численные методы решения обратных задач математической физики : учеб. пособие / А.А. Самарский, П.Н. Вабищевич. - Изд. 3-е. - М. :Изд-во ЛКИ, 2009. - 480 с.
2. Методы ускорения газодинамических расчетов на неструктурированных сетках / К.Н. Волков [и др.]. - М. : Физматлит, 2013. - 709 с.
3. Пикулин, В.П. Практический курс по уравнениям математической физики : учеб. пособие / В .П. Пикулин, С .И. Похожаев. - М. : Наука, 1995. - 224 с.
4. Пастухов, Д.Ф. Аппроксимация уравнения Пуассона на прямоугольнике повышенной точности / Д.Ф. Пастухов, Ю.Ф. Пастухов // Вестник Полоцкого университета. Серия С, Фундаментальные науки. - 2017. - № 12. - С. 62-77.
5. Кирхгоф, Г.Р. Механика : лекции по математической физике / Г.Р. Кирхгоф ; пер. с 4 нем. изд. - М. : ДомКнига, 2014. - 392 с.
6. Вакуленко, С.П. Способы передачи QR-кода в стеганографии / С.П. Вакуленко, Н.К. Волосова, Д.Ф. Пастухов // Мир транспорта. - 2018. - Т.16, № 5(78). - С. 14-25.
7. Волосова, Н.К. Преобразование Радона и краевой задачи для уравнения Пуассона в стеганографии / Н.К. Волосова // Тезисы докладов Международной конференции по дифференциальным уравнениям и динамическим системам, Суздаль, 6-11 июля 2018 г. / МФТИ. - Суздаль, 2018. - С. 61.
8. Колмогоров, А.Н. Элементы теории функций и функционального анализа / А.Н. Колмогоров, С.В. Фомин. - М. : Наука, 1976. - 543 с.
9. Бахвалов, Н.С. Численные методы в задачах и упражнениях / Н.С. Бахвалов, А.В. Лапин, Е.В. Чи-жонков. - М. : БИНОМ, 2010. - 240 с.
10. Бартеньев, О.В. Современный Фортран / О.В. Бартеньев. - М. : ДИАЛОГ - МИФИ, 2000. - 450 с.
11. Гришин, А.М. Об одном методе решения трехмерного эллиптического уравнения общего вида / А.М. Гришин, А.С. Якимов // Вычислительные технологии. - 2000. - Т. 5, № 5. - С. 38-52.
12. Козлов, А.А. Об управлении показателями Ляпунова двумерных линейных систем с локально интегрируемыми коэффициентами / А.А. Козлов // Дифференциальные уравнения. - 2008. - Т. 44, № 10. -С. 1319-1335.
13. Козлов, А.А. Об управлении показателями Ляпунова линейных систем в невырожденном случае / А.А. Козлов // Дифференциальные уравнения. - 2007. - Т. 43, № 5. - С. 621-627.
14. Вакуленко, С.П. К методу оценки состояния железнодорожного полотна / С.П. Вакуленко, К.А. Во-лосов, Н.К. Волосова // Мир транспорта. - 2016. - Т. 14, № 3 (64). - С. 20-35.
15. Вакуленко С.П. К вопросу о нелинейных волнах в стержнях / С.П. Вакуленко, А.К. Волосова, Н.К. Волосова // Мир транспорта. - 2018. - Т. 16, № 3 (76). - С. 6-17.
16. Пастухов, Д.Ф. Оптимальный порядок аппроксимации разностной схемы волнового уравнения на отрезке / Д.Ф. Пастухов, Ю.Ф. Пастухов, Н.К. Волосова // Вестник Полоцкого университета. Серия С, Фундаментальные науки. - 2018. - № 4. - С. 167-186.
17. Пастухов, Ю.Ф. Необходимые условия в обратной вариационной задаче / Ю.Ф. Пастухов // Фундаментальная и прикладная математика. - 2001. 7:1. С. 285-288.
18. Козлов, А.А. О равномерной глобальной достижимости двумерных линейных систем с локально интегрируемыми коэффициентами / А.А. Козлов, И.В. Инц // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. - 2017. - Т. 27, № 2. - С. 178-192.
19. Решение уравнения Пуассона в целых числах по модулю Р с кусочно-разрывной правой частью стеганографии / Н.К. Волосова, К.А. Волосов, Д.Ф. Пастухов, Ю.Ф. Пастухов // Евразийское научное объединение. - 2019. - Т. 1, № 1 (47). - С. 4-9.
20. Эффективная итерационная формула для краевой задачи уравнения Пуассона со сложно распределенными источниками / Н.К. Волосова [и др.] // Герценовские чтения : сб. ЬХХП Междунар. конф. по дифференциальным уравнениям и динамическим системам, СПб., 9-13 апр. 2019 г. / Рос. пед. ун-т им. А.И. Герцена. - СПб., 2019. - С. 234-238.
Поступила 01.03.2019
MINIMUM SCHEME OF THE DIFFERENCES FOR EQUATION OF THE POISSON ON BOX WITH SIXTH RATHER INACCURACY
D. PASTUHOV, Y. PASTUHOV, N. VOLOSOVA
The Offered algorithm of the decision of the general lumpy marginal problem Dirihle for three-dimensional equation of the Poisson on parallelepiped with sixth rather inaccuracy with minimum 27 point patterns. Numerically stability of the algorithm is checked to breakup of the first sort of the border conditions on side parallelepiped. It Is Received decomposition to inaccuracy of the problem in general type for uneven net through derived even order from decision and right part of equation and even order on each of three variable. The Writtenned program on base of the built algorithm and principle of the compressed images for evident formula iteration idle time. It Is Solved exactly and numerically test example, which comparison confirms the sixth order to inaccuracy for molded the numerical algorithm.
Keywords: three-dimensional equation of the Poisson on parallelepiped, lumpy-marginal task Dirihle.
