CC BY
106
Вычислительные технологии
Том 5, № 1, 2000
ГРАНИЧНЫЕ УСЛОВИЯ ПОВЫШЕННОЙ ТОЧНОСТИ В ПОЛЮСАХ КООРДИНАТНЫХ СИСТЕМ *
В. И. ПААСОНЕН Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: [email protected]
The method of setting the improved difference boundary conditions in poles (in points of singularities) of polar, cylindrical and spherical coordinates is developed, with the purpose of their application in high-order schemes for symmetric and nonsymmetrical boundary value problems. Difference equations in poles representing special approximations of differential equations in Cartesian coordinates, fitted with an order of accuracy and simplifications under the form are obtained. The problem on the realization of boundary conditions in the implicit schemes of an approximate factorization and in iterative processes is investigated.
В работах [4, 5, 8] построены схемы четвертого порядка точности для уравнений второго порядка в любых ортогональных системах координат. Для того чтобы воспользоваться этими результатами применительно к системам координат с особенностями (линиями или точками неоднозначности преобразования координат, для краткости называемыми в дальнейшем полюсами), необходимо с достаточной точностью поставить разностные граничные условия в узлах, соответствующих этим особенностям. Так, например, в случае краевой задачи для цилиндра на его оси необходимо поставить разностные граничные условия, порядок погрешности которых соответствовал бы порядку аппроксимации используемой схемы. В симметричном случае это линия r = 0 прямоугольной области изменения радиуса и высоты (r, z), а для произвольной задачи — плоскость r = 0 трехмерной области (r, z, ф), где ф — полярный угол. В полярных координатах, в зависимости от того, симметрична задача или нет, полюс r = 0 есть точка или прямая соответственно, а в сферических — точка или плоскость. Общим моментом в этих примерах является неоднозначность преобразования координат в полюсах и отсутствие в исходных краевых задачах для дифференциального уравнения граничных условий в них, тогда как разностная задача непременно требует их постановки.
Для симметричных задач в случае схем второго порядка точности вопросы постановки граничных условий, сходимости и реализации алгоритмов хорошо разработаны (см., например, монографии [2, 6]). В случае же схем повышенной точности построение удовлетворительных соотношений в полюсах даже для симметричных краевых задач является
* Работа выполнена при финансовой поддержке программы Интеграционных фундаментальных исследований СО РАН, проект №43, и Российского фонда фундаментальных исследований, грант №97-01-00819.
© В. И. Паасонен, 2000.
проблемой. Несимметричные разностные краевые задачи в этом отношении еще сложнее даже при традиционных схемах, не говоря уже о высокоточных. Решению этих проблем и посвящена данная статья.
Работа с аналогичным содержанием одновременно направлена в Russian Journal of Numerical Analysis and Mathematical Modelling, для наших англоязычных коллег. Предлагаемая работа является по отношению к ней сокращенным вариантом — в ней опущены некоторые подробности.
1. Схема повышенной точности
Рассмотрим краевую задачу в полярных, цилиндрических или сферических координатах для параболического (или соответствующего ему эллиптического) уравнения
ди
— = grad (Л div и) - /. (1.1)
В любых ортогональных координатах уравнение (1.1) можно записать через коэффициенты Ламе:
ди=£ £ (л*а - /• (1.2)
1=\ 4 7
Выражения коэффициентов Ь, * для конкретных криволинейных координат опускаются, так как они общеизвестны. Заметим, что в виде (1.2) естественным образом можно записать и симметричные случаи, при этом понижается размерность задачи. Каждый одномерный дифференциальный оператор
д ди
Ьи = дхгх) дх
входящий в уравнение (1.2), будем аппроксимировать разностным
Ли(х) = —- ,—--—-и(х), (1.3)
п
где Д± — разделенные разности "вперед" и "назад", а коэффициенты с(х ± П/2) выбраны так, чтобы выполнялось разложение
л = Ь + 12 ь(ф ь) +0(п4). (1.4)
Функция с(х), определяющая коэффициенты оператора (1.3), зависит от переменного коэффициента ф(х) дифференциального оператора Ь и, таким образом, представляет собой зависимость с(ф), называемую шаблонным функционалом. В работе [4] сформулирован необходимый и достаточный критерий1 справедливости разложения вида (1.4). Он представляет собой условие на первые два члена разложения функционала с(ф):
h2
Ф) = ф(ж) +—
^ - ф'(х)(ф(х))-1ф'(х)
+ O(h4). (1.5)
:В работе [8] критерий приведен в иной форме, эквивалентной ввиду тождества ф(ф 1)'ф' + ф'ф ф' = 0.
Условию (1.5) удовлетворяет, например, выбор
к\ 2
с х ± .
2 3
х ± М + ф(х)ф(х ± к)
2) ф(х) + ф(х ± к)_
который, впрочем, не определяется однозначно, так как важны только первые два члена разложения по степеням к2.
Пусть одномерные разностные операторы Лг вида (1.3), аппроксимирующие
г 9 ф 9 ф Х
Г = -К~фг д—, фг = а^Х ,
^Ух г ихг
выбраны так, что для каждого из них выполняется критерий (1.5) по соответствующей переменной хг. Тогда схема приближенной факторизации
т 7/П+1 — ип
1Т Аг--— = Пип - Рп (1.6)
1 т
г=1
аппроксимирует исходное уравнение (1.2) с погрешностью 0(т2 + к4), если сомножители оператора верхнего слоя имеют вид
1 к2<г 1
Аг = Е — тЬгЛг(агЕ), аг =---— , 1г = -—— ,
г г гУ г ^ г 2 12т г ХагЬг'
где Е — тождественный оператор, т — шаг по времени, п — номер временного шага (пт = ¿), а оператор на нижнем слое и правая часть имеют вид
т к2 < П = £ ЬгЛг(Е + Ь] Л] ),
г=1 ]=г
£ п+1 I £ п т к 2
Рп = 1 2+1 + ^72 ЬгЛгШп)
2 12
г=1
(подробности см. в работах [4, 5, 8]).
В стационарном случае схему вида (1.6) формально можно считать итерационным процессом решения разностного уравнения четвертого порядка точности
т Ь2
= р, р = f + ЬлШ),
г=1
представляющего собой обобщение схемы Микеладзе [3]. Но, разумеется, аг и т в этом случае следует считать независимыми итерационными параметрами, отказавшись от представленной выше связи аг с шагами сетки.
Итак, будем считать, что во внутренних узлах сеточной области исходное дифференциальное уравнение аппроксимируется с четвертым порядком относительно шагов пространственной сетки и, для динамических задач, со вторым порядком по времени. Будем предполагать, кроме того, что условия на внешних границах также аппроксимируются с повышенной точностью. Например, в условиях второго и третьего рода — путем записи специальных компактных аналогов уравнения на границах [7] или с помощью многоточечных одномерных аналогов нормальных производных [1]. Поэтому сосредоточимся исключительно на предмете данной работы — построении специальных разностных уравнений в полюсах координатных систем.
2. Одномерные симметричные задачи
Рассмотрим сначала простые примеры — краевые задачи в круге и шаре для эллиптического уравнения с симметричными граничными условиями. Коэффициенты для простоты здесь и в дальнейшем будем считать постоянными. В обоих случаях входные данные и решение зависят только от радиуса, и (1.1) превращается в обыкновенное дифференциальное уравнение
1 д 1 du f (2 1)
r1 dr dr '
где r — радиус, l — индекс задачи, равный 1 для круга и 2 для шара.
Построим граничное условие в центре круга. С этой целью запишем схему Микеладзе на квадратной сетке в начале декартовых координатах (x,y):
h2 h2 Лиоо = Rfoo, Л = Лх + Лу + — ЛхЛу, R = E + —(Лж + Лу ). (2.2)
Здесь Лх, Лу — обычные аппроксимации вторых производных, а двойной индекс (ij) указывает на точку (ih,jh), в данном случае на начало координат (0, 0). Уравнение (2.2) в индексах запишется в виде
1 2 1
[U11 + 4 Uio - 20 uoo] = зfoo + — Fio, (2.3)
где вследствие симметрии
U11 = u11 + u1-1 + u-11 + u-1-1 = 4 u(h л/2),
Uio = U1o + Uo1 + u-1o + Uo-1 = 4 u(h), Fw = f1o + fo1 + f-1o + fo-1 = 4f (h).
Здесь и в дальнейшем за функцией v(r) = vx2 + y2) = u(x,y) условимся сохранять обозначение u. Из симметрии немедленно следует связь между значениями решения в точках 0, h, hv^ оси r:
u(h^2) + 4u(h) - 5u(0) = f[f (h) + 2f(0)]. (2.4)
Сложность заключается в том, что формулой (2.4) невозможно воспользоваться непосредственно как граничным условием, так как равномерная по r сетка не содержит точку h\/2 в качестве своего узла. Тем не менее соотношение (2.4) оказывается полезным для построения граничного условия четвертого порядка аппроксимации. С этой целью сначала запишем схему (2.2) еще раз, но с шагом h\/2, а затем исключим из этих двух уравнений значение функции u в посторонней, не принадлежащей сетке, точке В итоге получим
трехточечное граничное условие, содержащее значения искомой функции исключительно в узлах сетки:
u
(2h) - 16u(h) + 15u(0) = h2[f (hv^2) - 2f (h) - 2f (0)]. (2.5)
Возможны и другие подходы. Например, можно строить граничное условие исходя не из компактной схемы Микеладзе, а из схемы типа "большой крест" на основе пятиточечных аппроксимаций вторых производных. Тогда вместо (2.5) получается граничное условие
Ц2й) - 16и(й) + 15и(0) = -3й2/(0), (2.6)
с точностью до членов шестого порядка малости не противоречащее (2.5). Это следует из справедливости для достаточно гладкой четной функции равенства
/(Ь^2) = 2/(Ь) - /(0) + 0(Ь4)
Иной способ исключения неудобной точки — интерполяция. Для четной функции с погрешностью 0(Ь6) имеем
Г- 1 4 1
м(^2) = -и(2Ь) + -и(Ь) — и(0). (2.7)
6 3 2
Подстановка (2.7) в (2.4) дает еще один эквивалентный по точности вариант граничного условия в центре круга.
В случае центрально-симметричной задачи поступаем аналогично — используем в центре шара схему Микеладзе на кубической сетке, которая в индексах имеет вид
^[ии + 2^10 - 24иооо] = 1 /000 + 112^о, (2.8)
где выражения и11 ,и10 содержат не по 4, как в плоском случае, а 12 или 6 слагаемых, причем в силу симметрии слагаемые в каждой из этих сумм совпадают. Отсюда вытекает равенство, справедливое с погрешностью шестого порядка
«(Ь^2) + и(Ь) - 2и(0) = (0) + /(Ь)]. (2.9)
Здесь та же проблема, и решается она вышеизложенным способом: уравнение (2.9) пишется еще раз с шагом Ь\/2, и после исключения и(Ь\/2) получается граничное условие
и(2Ь) - и(Ь) = (0) + 2/(Ь^2) - /(Ь)],
причем ввиду четности функции / выражение в скобках можно заменить выражением
[3/(Ь) - / (0)].
Заметим, что если исходить из схемы четвертого порядка точности на шаблоне "большой крест", то можно получить иное трехточечное условие, также достаточно точное. Третий вариант граничного условия получится, если к первому слагаемому (2.9) применить интерполяционную формулу (2.7).
3. Симметричная задача в цилиндрических координатах
Рассмотрим теперь осесимметричную задачу в цилиндрических координатах. Для вывода граничного условия при г = 0 запишем схему четвертого порядка точности в декартовых координатах в точках оси х = у = 0, г = причем шаги по переменным х и у примем одинаковыми и равными Ь, а по переменной г — равным д. Эта схема имеет вид
2
Ь2 Ь2 + д
Лх + Лу + лх + —ЛхЛу + 4 (ЛЖЛ^ + ЛуЛх)
6 у 12
где 2
12 д
и00к = 0к, (3.1)
К
00к =
Ь2 ~2
Е + —(Лх + Лу
/00к.
12 12
Сгруппируем в (3.1) отдельно от всех остальных слагаемые, содержащие Лх. В результате получим уравнение
Ли00к + ВУ00к = ^00к, ^00к = Лх и00к, (3.2)
где Л — ранее определенный двумерный оператор,
Ъ2 + п2
В = Е + + Лу),
а индекс (г]к) соответствует точке (х^,у^,гк).
Из условия симметрии вытекает возможность сделать следующие упрощения:
Лпоок = зЪ^ [ик(Ъ^/2) + (Ъ) - Ъпи(0)],
Вп,
ООк
1 + Р2
Е + —ЪА+Г
п
Лг ик (0), р = -, Ъ
Я
00 к
2
Ъп
Е +—А+г + — Л2 3 + 12 2
/к(0), (3.3)
где А+г — оператор разделенной разности "вперед" по радиальной переменной. Шаблон уравнения (3.2) после проведенных упрощений становится семиточечным — он содержит по три точки на оси г = 0 и на прямой г = Ъ, отстоящей на шаг от оси, и седьмую точку (г = Ъл/2, г = ^к), не принадлежащую сеточной области.
Таким образом, подобно рассмотренным выше одномерным случаям, полученное соотношение содержит значение искомой функции в нерегулярной точке (г, г) = (Ъ^,гк). Но в отличие от них повторная запись уравнения (3.2) с шагом Ъ\[2 с целью исключения из пары уравнений неудобного значения ничего не дает. Этот путь приводит лишь к появлению вместо одного двух других, столь же неудобных, значений на соседних к г = гк плоскостях г = £к±ъ Если вместо компактной схемы для построения граничного условия использовать схему на основе пятиточечных аппроксимаций вторых производных, то это приводит к проблеме "законтурных" точек при записи уравнения в точках оси, отстоящих на один шаг от оснований цилиндра.
В силе остается третий путь — замена значения функции и(г, г) в точке г = Ъ\/2 результатом ее интерполяции по значениям в точках г = 0, Ъ, 2Ъ. Заметим, что значение Пк (Ъл/2) = и входит в уравнение (3.2) с множителем 0(Ъ 2). Поэтому для со-
хранения четвертого порядка точности интерполяцию необходимо провести по меньшей мере с погрешностью шестого порядка. А это именно тот порядок, который для четных функций обеспечивает трехточечная интерполяционная формула (2.7).
Подправленный оператор, действующий на регулярной сетке, имеет вид
Лпоок = 9Ъ2 [пк(2Ъ) + 32пк(Ъ) - 33пк(0)].
Теперь подстановка (3.3) в (3.2) с учетом этой поправки даст уравнение, записанное на семиточечном регулярном шаблоне, в котором неудобная точка (Ъ\/2, гк) заменена точкой (2Ъ, гк), принадлежащей сетке.
Реализация граничного условия при использовании итерационной схемы в форме (1.6) может быть осуществлена различными способами. Обычно это предполагает переход к двухслойной форме не только во внутренних узлах сетки, но и в граничном условии, а также какое-нибудь расщепление условия (3.2), которое выделяло бы на верхнем итерационном слое одномерный разностный оператор, действующий ортогонально границе, т. е. по радиусу. При этом обычно точное граничное условие (3.2) выполняется лишь при п ^ то. Эта техника хорошо известна, и улучшенная аппроксимация не вносит в порядок действий принципиальных отличий в сравнении с применяемым при использовании традиционных схем. Различие состоит лишь в степени сложности шаблона на нижнем итерационном слое.
4. Несимметричная краевая задача
Рассмотрим сначала несимметричную задачу в полярных координатах. Будем предполагать, что на промежутке [0, п/2] по угловой переменной ф укладывается целое число шагов. Это ограничение не является обременительным и необходимо для того, чтобы лучи
ф = 0, п/2, п, 3п/2,
соответствующие полуосям декартовой системы координат, совпали бы с некоторыми линиями сетки.
Если в качестве граничного условия записать в центре круга схему четвертого порядка точности вида (1.6) в декартовых координатах ж, у, то, перейдя к индексным выражениям, получим уравнение
^ [ЗДЬ) + 4 5хо(Л) - 20 и(0,0)] = 3/(0, 0) + -2(4.1)
где выражения вида
£10(Л) = и(Л, 0) + и(-Л, 0) + и(0, Л) + и(0, -Л),
$п(Л) = Л) + и(-Л, Л) + -Л) + и(-Л, -Л)
не удается упростить так легко, как в симметричных задачах. Однако следует заметить, что ввиду сделанного выше предположения по крайней мере слагаемые функции $ю(Л) определены на полярной сетке — это значения в узлах (г, ф) при г = к и ф = 0, п/2, п, 3п/2.
Проблему создает функция 511(Л), содержащая в качестве слагаемых значения функции в точках (Л\/2, п/4 + кп/2), не принадлежащих сетке. Однако здесь, как и в симметричном случае, оказывается возможным эквивалентное с точностью до допустимой погрешности переопределение функции 511(Л). Основой для такого преобразования служит свойство четности функции 511(Л), которое очевидным образом вытекает из самой формы ее записи.
Легко получить ее разложение в окрестности нуля, разложив каждое из слагаемых и(±Л, ±Л) в ряд Тейлора по обеим переменным. Для наших целей достаточно удержать лишь три члена разложения. Тогда с погрешностью шестого порядка получим
ЗД) = 4« + + £ ( ^ + 0 + 6 + «(П (4.2)
где А — оператор Лапласа в координатах ж, у.
Задача состоит в переопределении функции 5л (Л) так, чтобы в новом представлении она зависела бы от значений функции и (ж, у) только на осях ж, у в узлах полярной сетки. Это означает необходимость избавиться от смешанной производной, что можно сделать, используя следствия дифференциального уравнения. На гладких решениях уравнения Пуассона имеем
д4и 1 4 „ 1 (д4и д4ич -А/ -- — +
<9ж23у2 2 ' 2\ <9ж4 <9у4
Заменим в (4.2) таким образом вычисленную смешанную производную, воспользуемся равенством Аи = /, а затем аппроксимируем четвертые производные пятиточечными разностными отношениями. Тогда с эквивалентной погрешностью получим
5П(Л) = 4и + 2Л2/ + у А/ - + Л2 )и. (4.3)
Ясно, что в (4.3) вполне можно было заменить производные от / разностными отношениями, и тогда получилась бы несколько иная формула, аналогичная (4.3). Если же не проводить в (4.2) замену Ди через правую часть /, то в итоге получился бы другой вариант приближенного представления функции Бц(П), также с погрешностью шестого порядка:
Бп(П) = 4и + 2П2(Лл + Лу)и - у (ЛХ + Л;)и + у (Лх + Лу)/. (4.4)
Теперь осталось подставить функцию (4.3) или (4.4) в (4.1), чтобы получить регулярное граничное условие на оси г = 0. В осях (х,у) это условие пишется на шаблоне "большой крест" в начале координат, а в осях (г, ф) имеет шаблон с тремя точками г = 0, П, 2П на четырех параллельных прямых ф = 0, п/2, п, 3п/2. Ясно, что при г = 0 точки с разными координатами ф отождествляются.
Рассмотрим теперь задачу в сферических координатах (г, ф, 0), не предполагая симметрии решения. Пусть как и выше, шаги сетки по угловым переменным выбраны соизмеримыми с п/2, что обеспечивает принадлежность всех осей декартовой системы к линиям сетки. Записав схему четвертого порядка точности в начале координат х,у, г, получим уравнение вида
(ип(П) + 2Цю(П) - 24и(0, 0, 0)) = ^ш, (4.5)
где
иш(П) = и(П, 0,0) + и(-П, 0,0) + и(0, п, 0) + и(0, -П, 0) + и(0,0, П) + и(0,0, -П),
ип(П) = Бг (П) + Бу (П) + Бх(П),
/о00 — правая часть схемы, а функции Бх, Бу, Бх аналогичны функции Бп(П), определенной в (4.1). Каждая из них представляет собой сумму четырех значений функции и(х, у, г) при одном нулевом аргументе и двух других аргументах, равных ±П. Индекс означает фиксированную переменную. Так, например, функция
Бг(П) = и(П, П, 0) + и(П, -П, 0) + и(-П, П, 0) + и(-П, -П, 0)
имеет фиксированный третий аргумент г = 0.
Для получения регулярного граничного условия функцию [/ц(П) нужно переопределить так, чтобы она зависела только от значений функции и(х,у,г) на осях. Из результатов предыдущего раздела следует, что Бх (П) разлагается в ряд в полном соответствии с формулой (4.2)
Б,(П)=4и + 2П2 (дХи + 0) + £ (1X4 + 6^ + ) +0(П6).
Разложения двух других функций Бх(П), 5у(П) получаются круговой заменой х,у,г. Суммируя их, получим
П4
Ц7ц(П) = 12и + 4П2Ди + —
3
д4и д4и д4и ( д4и д4и д4и
+ ^ + ^ + 3 п пп п + „ „„ „ +
дх4 ду4 дг4 \дж2 ду2 дх2дг2 ду2дг2
+ 0(П6), (4.6)
где А — оператор Лапласа в пространстве х,у,г. Исключая все смешанные производные порядка 2 + 2 с привлечением следствий уравнения Пуассона по формуле
д4 и д4и д4и 1
+ ТОТ^ + тт-г^ = тгА/ - (ЛХ + ЛУ + Л>] + 0(-2)
дх2ду2 дх2дг2 ду2дг2 2 и заменяя Аи правой частью /, получим с погрешностью 0(-6)
адь) = 12и + 4-2/ + -6 [3А/ - (ЛХ + л; + Л2)и]. (4.7)
Подстановка (4.7) в (4.5) дает окончательное выражение граничного условия в центре сферической системы координат. Здесь также возможны иные эквивалентные по точности варианты формулы (4.7). Например, можно заменить А/ его разностным аналогом или отказаться от замены Аи правой частью /, а непосредственно аппроксимировать с погрешностью 0(-4) функцию Аи на шаблоне "большой крест".
5. Задача для параболического уравнения
Построение граничных условий в полюсах криволинейных систем координат в случае нестационарных краевых задач отличается от рассмотренных выше случаев тем, что в качестве исходного объекта для преобразований необходимо взять схему повышенной точности для параболического уравнения, записанную в декартовых координатах в точках неоднозначности преобразования координат. При этом важно построить именно схему без факторизации по х и у, так как у нее оператор, действующий на верхнем слое, имеет более простой шаблон, чем у факторизованной схемы (простой "крест" вместо "ящика"), в то время как различие в их разложении находится в пределах погрешности аппроксимации.
Рассмотрим сначала краевую задачу в полярных координатах. В начале координат запишем схему, аппроксимирующую с погрешностью 0(т2 + -4) уравнение теплопроводности в осях х, у:
ип+1 — ип /п+1 + /п -2
В-= Лип - = / 0+ / + -(Лх + Л;)/п, (5.1)
т 2 12
1 -2
где
В = Вх; = Е - ат(Лх + Л;), а
2 12т
а двумерный оператор Л определен ранее. В индексах оператор Л дает девятиточечное выражение на шаблоне "ящик", уже обсуждавшееся выше. В симметричном случае оно сворачивается и после применения интерполяции (2.7) переходит в трехточечное выражение на регулярном шаблоне (0, -, 2- ). В произвольном случае оператор Л переопределяется, как описано выше, эквивалентным по точности выражением на шаблоне "большой крест". Отличие состоит лишь в том, что в параболическом случае Аи = / + и^, поэтому в выражениях из разд. 4, полученных с использованием уравнения и его следствий, необходимо заменить / на сумму / + и4. Учитывая это замечание и подставляя (4.4) в выражение
Лип = ^[£п(-) + 4 йо(-) - 20 ип(0,0)], б-2
после замены производной по времени разностным отношением получим с погрешностью 0(т2 + к4)
где
к2
Лип = Л ип + —(ЛХ + Лу) 12
к2
Iп +
ип+1 - ип
т
(5.2)
л = Лх + Лу - —(ЛХ + Л2)
Из (5.2) следует граничное условие вида (5.1), в котором Л заменено на Л, а в В и Рп учтены поправки, связанные со вторым слагаемым (5.2)
т 1п+1 + 1п
В = Е - 2(Лх + ЛУ), Рп = ] ^ .
Аналогично строится граничное условие в центре сферической системы координат. Отличие состоит лишь в размерности шаблона "крест".
В цилиндрических координатах схему на оси цилиндра запишем в виде
В
ип+1 - ип
I п+1 + Iп к2 о2
1 2+ 1 + кк2(Лх +Лу )1п + ^Лг Г, (5.3)
где
к2 + 02
П = Л + Л,, С = Е + _4 (Лх + Лу),
В = Е - втЛг - ат (ЛХ + Лу),
12 1
к2
а
2 12т'
в
1
о2
2 12т
Действия по преобразованию выражения Пип к регулярной форме практически сводятся к описанным выше преобразованиям в Лип, так как операторы В, С и правая часть определены на шаблоне "крест", не содержащем точек (±к, ±к), и поэтому в симметричном случае дают двухточечные выражения, а в несимметричном не нуждаются в преобразованиях. Модификация П влечет, как и выше, модификацию правой части Рп и оператора В. Новые их выражения имеют вид
П = Л + Лг С, Рп
I п+1 + 1п 02
1—+1 + Лг г 2 12 -1
т
В = Е - втЛг - ^(Лх + Лу). Построение завершается неполной факторизацией оператора
В = В г Вху =[Е - вт Лг ] Е - ^(Лх + Лу)
т
При полной факторизации, т. е. при условии
Вху ВхВу
Е - т ЛХ
2х
Е - т Лу 2 у
шаблон пополнился бы еще четырьмя точками (±к, ±к), что привело бы к потере регулярности. Это соображение свидетельствует о правильности способа факторизации в цилиндрических координатах. То же можно утверждать относительно выбора оператора без факторизации в полярных и сферических координатах.
6. Реализация граничных условий
Рассмотрим вопрос о реализации алгоритмов на примере краевой задачи в цилиндре для уравнения теплопроводности. Сначала рассмотрим симметричную задачу. Схема повышенной точности типа (1.3) запишется в виде
„ п+1 _ „ п
Аг Аг- = С(„П,/П)
т
(здесь О есть результат вычислительной работы над правой частью и искомой функцией на предыдущем слое), а соответствующая ей схема в дробных шагах:
Аг ип+ 2 = О(ип,/п),
Граничные условия на оси
„п+1_ „п
Аги-^ = „п+1. (6.1)
т
ип+1 - ип
Бг Бху-= £(ип,/п)
также представим в расщепленном виде:
Бг ип = Д(ип,/п), ип+1 _ ип
Бху и-— = ип (6.2)
т
Последовательность действий при решении разностной краевой задачи (6.1) - (6.2) следующая. Сначала прогонками параллельно оси решается первое уравнение системы (6.1) с граничными условиями на основаниях цилиндра Г:
1, л ип+1|г - ип|г
ип+ 2 |г = А,
т
Затем прогонкой по оси решается первое уравнение системы (6.2) с аналогичными граничными условиями в центрах оснований цилиндра, вытекающими из определения величины ип. Решение задачи на слое завершают радиальные прогонки для второго уравнения системы (6.1). При этом в качестве граничных условий на оси г = 0 выступает второе уравнение системы (6.2) (выше установлено, что оно двухточечное), а на правой границе — заданные условия на внешней цилиндрической поверхности.
Заметим, что для реализации схемы (6.1) важна правильная последовательность операторов в произведении Аг Аг. Перестановка операторов привела бы к необходимости ставить граничные условия для промежуточной величины ип+1/2 не на основаниях, где это удается просто сделать, а на оси цилиндра, где их постановка затруднена. Неверная постановка граничных условий для промежуточных величин приводит к неэквивалентности схемы в дробных шагах исходной схеме и к потере точности.
Перейдем теперь к рассмотрению вопроса о реализации схемы в произвольном случае. Вместо (6.1) имеем схему с записью в дробных шагах:
Афип+3 = О(ип,/п),
2 1 Аг 2 = 3
(6.3)
т
Здесь также важно, чтобы последним был оператор Аг. Иначе на оси было бы необходимо ставить граничное условие для промежуточной величины.
Первое одномерное уравнение системы (6.3) решается прогонкой по окружностям с условиями периодичности на границах ф = 0 и ф = 2п, второе — по всем внутренним линиям сетки, параллельным оси цилиндра, за исключением самой оси, с условиями на основаниях цилиндра, вытекающими из третьего разностного уравнения (6.3). Затем на оси решается первое из уравнений (6.2) способом, описанным выше. Тем самым трехмерное граничное условие на оси сводится к плоскому, записанному на шаблоне "крест", т. е. ко второму из уравнений (6.2).
Наиболее сложен третий этап — решение системы одномерных уравнений по радиусам цилиндра. Для этого в каждой плоскости г = ^, параллельной основанию цилиндра, по четырем направлениям ф = 0, п/2, п, 3п/2, соответствующим полуосям х и у, осуществляются навстречу друг другу четыре прямые прогонки от периферии к оси цилиндра. Затем в точке встречи решается система пяти уравнений с пятью неизвестными. Система включает в себя второе из уравнений (6.2) и четыре связи между значениями на оси и значениями в ближайших точках (х = ±Л,, у = ±Л.), которые определяются благодаря уже вычисленным прогоночным коэффициентам. Эта очень простая система позволяет вычислить значение решения на оси г.
Только теперь, с учетом того обстоятельства, что решение на оси уже вычислено, появляется возможность осуществить прогонки решения абсолютно по всем направлениям угла ф с условиями первого рода на оси цилиндра.
В сферической системе координат реализация аналогична, и даже проще в том отношении, что нет необходимости расщеплять граничное условие. Различие лишь в том, что в центральной точке условие записано не на плоском, а на трехмерном "кресте", и для вычисления решения в центре необходимо провести не четыре, а шесть встречных прогонок по полуосям. В точке пересечения возникает необходимость решить систему из семи уравнений, по числу точек объемного "креста". Несмотря на большую размерность, она довольно простая и значение искомого решения в центре легко выражается через значения в соседних узлах.
В стационарном случае возможны различные варианты реализации граничных условий. Если в записи итерационного граничного условия на верхнем слое шаблон представляет собой простой "крест", то реализация ничем не отличается от случая нестационарной задачи. Если же это "большой крест" (например, при точном выполнении граничного условия на каждой итерации), то для замыкания системы, определяющей решение в полюсе, необходимо к прогоночным соотношениям в ближайших к полюсу точках подключить соотношения в предыдущих точках каждой из четырех (шести) прогонок. Тогда для цилиндра получится система из девяти уравнений, а для шара — из тринадцати, по числу точек в шаблоне.
7. Заключение
Положенные в основу данной работы идеи представляются автору естественными, дополнительные требования к исходным данным практически сводятся только к требованию
повышенной гладкости, а полученные в результате алгоритмы являются легко реализуемыми в рамках традиционных подходов к решению многомерных задач. Очевидно, аналогичные построения могут быть проведены для иных ортогональных систем координат с простыми полюсами, в частности для тороидальной. Более сложными для реализации граничных условий, но отнюдь не безнадежными, выглядят координатные системы с двумя фокусами, например эллиптическая, эллипсоидальная или бифокальная.
Есть все основания полагать, что аналогичным образом можно строить специальные аппроксимации в полюсах более произвольных сеток, включая адаптивные. Это важно в том отношении, что далеко не всегда произвольные области целесообразно отображать на прямоугольник. В частности, если граница физической области гладкая, то ее отображение на прямоугольник не может быть качественным. Ему следует предпочесть гладкое отображение на область с полюсом, а для этого необходимо преодолевать неоднозначность преобразования, используя в численных моделях соотношения, аналогичные рассмотренным здесь.
Выражаю признательность моей дочери Юлии за помощь при оформлении настоящей работы.
Список литературы
[1] ЖДАН С. А., ПААСОНЕН В. И. Схема повышенной точности для задачи о непотенциальном течении идеальной жидкости в плоских каналах. Численные методы механики сплошной среды, Новосибирск, 3, №5, 1972, 35-39.
[2] Калиткин Н. Н. Численные методы. Наука, М., 1978.
[3] МиКЕЛАДЗЕ Ш. Е. О численном интегрировании уравнений эллиптического и параболического типов. Известия АН СССР. Сер. матем., 5, №1, 1941, 57-74.
[4] ПААСОНЕН В. И. Компактные схемы для систем уравнений второго порядка с конвективными членами. Вычислительные технологии, 3, №1, 1998, 55-66.
[5] ПААСОНЕН В. И. Обобщение методов повышенной точности для нелинейных уравнений второго порядка в ортогональных системах координат. Численные методы механики сплошной среды, Новосибирск, 8, №2, 1977, 94-99.
[6] САМАРСКИЙ А. А. Теория разностных схем. Наука, М., 1983.
[7] ФРЯЗИНОВ И. В. О схеме повышенного порядка точности решения третьей краевой задачи для уравнения Au — qu = — f в прямоугольнике. Журн. вычисл. матем. и матем. физики, 11, №2, 1971, 515-517.
[8] PAASONEN V. I. Compact schemes for system of second-order equations without mixed derivatives. Russ. J. Numer. Anal. Math. Modelling, 13, №4, 1998, 335-344.
Поступила в редакцию 24 декабря 1998 г.
