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

Вычислительные технологии

Том 6, № 2, 2001

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

А. М. Гришин, А. С. Якимов Томский государственный университет, Россия e-mail: rector@tsu.ru

Stable difference schemes have been obtained by the generalization of iterative-interpolation method with the conventional approximation error for the solution of a three-dimensional equation of a general type. The convergence of the iterative process is studied. Calculation results are presented for the Poisson equation and at the constant transport coefficients the convergence rate has been theoretically estimated.

Известны итерационные методы [1-4] решения многомерного уравнения Пуассона. Эти методы иногда используют параметр релаксации ш (0 < ш < 2), который необходимо подбирать при решении эллиптического уравнения общего вида. Их не удается применить при численном решении трехмерного уравнения в частных производных, за исключением, может быть, универсального локально-одномерного метода расщепления (разностные схемы суммарной погрешности аппроксимации) [5, 6] или метода Дугласа [7, 8] (схема стабилизирующей поправки с погрешностью аппроксимации в обычном смысле).

В работах [2, 4-6, 8] отмечено, что итерационные схемы можно трактовать как методы установления (при t ^ ж) для соответствующего нестационарного уравнения со стационарными (не зависящими от времени) граничными условиями. Однако имеется одна особенность [2]: в нестационарных задачах для обеспечения точности решения шаги по времени т должны быть достаточно малы, в стационарных же задачах оптимальные итерационные параметры то выбираются из условия, при котором разностное решение выйдет на стационарное за наименьшее число шагов K и параметр т0 может принимать относительно большие значения.

Следуя [6], можно показать, что симметричный вариант локально-одномерной схемы расщепления для уравнения Пуассона в трехмерном случае имеет сходимость O(N) при hm = h = 1/N (h — шаг разностной сетки по пространству, m = 1, 2, 3). Но, как справедливо отмечено в [4], аддитивные схемы обладают свойством "близости" к решаемому уравнению только при ограничении на итерационный параметр т (например, т ^ 0 при n ^ ж, n — номер итерации). В общем случае свойство "близости" итерационной разностной схемы к решаемому уравнению сформулировано в [8] как условие полной аппроксимации.

* Работа выполнена при финансовой поддержке Федеральных целевых программ "Университеты России — фундаментальные исследования", "Интеграция" (проект "Академический университет") и Российского фонда фундаментальных исследований, гранты № 99-01-00352, № 99-01-00363.

© А. М. Гришин, А. С. Якимов, 2001.

Цель данной работы — показать, что итерационно-интерполяционный метод (ИИМ) [9] свободен от отмеченных выше недостатков, а также исследовать аппроксимацию, устойчивость и найти теоретическую оценку скорости сходимости для трехмерного уравнения Пуассона при постоянных коэффициентах переноса.

1. Алгоритм метода

Для простоты дальнейших выкладок логическую схему метода приведем с постоянными коэффициентами переноса внутри параллелепипеда R {x = (x^x2,x3), 0 < xm < dm, m = 1, 2, 3}, С4 > 0:

3 д 2T

£ CmTTY - с,Т(х) + f(x) = 0 . (1.1)

m=1 dXm

Вместо (1.1) рассмотрим нестационарное уравнение

дТ 3 д 2T с4Т + — + /(*)• (1-2)

Ui m=1 dxm

Для построении разностного аналога уравнения (1.2) введем неравномерную сетку x1,i (i = 0,1, 2, ..., N1), X2j (j = 0,1, 2, ..., N2), X3,k (k = 0,1, 2, ..., N3). Введем верхние индексы n, n + 1/3, n + 2/3, n +1, соответствующие начальному и последующим этапам итераций.

Используя результаты работ [9, 10], введем обозначения:

a(t, x) = Т + С4Т - С2д2Т/дх2 - с^д2Т/дх2 , (1.3)

b(t,x) = Т - С3д2Т/дх2ъ + С4Т , (1.4)

h1,i-1 = x1,i — x1,i-1 , h2,j-1 = x2,j - x2,j-1 , h 1,i = 0.5(h1,i-1 + h1,i),

h1,i = x1,i+1 - x1,i , h3,k-1 = x3,k - x3,k-1 ,

где точка над символом обозначает частную производную по времени.

Разностные операторы внутри области определения R из работы [10] имеют вид

Aa = [(h1a)i-1 + 4(h 1a)i + h^a^]/6h 1,i ,

Л1Т = c1[(Ti+1 - Ti)/h1,i + (Ti-1 - Ti)/h1,i-1]/h 1,i,

i = 1,...,N1 - 1 , 0 < xm < dm , m =2, 3 . (1.5)

Разностные операторы внутри области определения R для нового подхода записываются следующим образом:

A1T = [hM

1,j,k + 4h 1,iTi,j,k + h1,iTi+1,j,k] /6h 1,i ,

Л1T = c1[(Ti+1,j,k - Ti,j,k)/h1,i + (Ti-1,j,k - Ti,j,k)/h1,i-1]/h 1,i,

i,j,k =1,...,Nm - 1, m =1, 2, 3. (1.6)

Отметим, что операторы A2T, A3T, A2T, A3T в уравнениях (1.13) - (1.16), а также AT2, AT3 в приведенных в п. 2 уравнениях подробно выписаны в [9].

Следуя алгоритму из [10] и используя обозначения (1.3), (1.5), для уравнения (1.2) последовательно получим по направлению координаты х1

а = ад 2Т/дх\ + /,

разностные схемы поочередного интегрирования по пространственным переменным на основе ИИМ имеют вид

Аа = А1Т™+1/3 + А/? , г = 1,...,^ - 1, 0 <хт <dm, т =2, 3; (1.7)

по направлению оси х2

Ы = С2д2Тг/дх22 + аг, г = 1, ..., N - 1, разностные уравнения записываются так:

А2Ы,, = Л2Т?+2/3 + А2аг,, ,

г = 1,...,^ - 1 , 3 = 1, ...,N2 - 1, 0 <хз <dз; (1.8)

окончательно по направлению координаты х3

Т+ С4Т= сзд2Тг^/дх23 + Ы^, г = 1, ..., N1 - 1 , j = 1, ..., N2 - 1

разностные схемы ИИМ принимают вид

Аз (Тг ,3,к + с4Тг,], к) = ЛзТ?+ + АзЬ

г, 2, к = 1,...,^ - 1 , т =1, 2, 3. (1.9)

Отметим, что разностные уравнения вида, например (1.7), получаются непосредственно из логической схемы ИИМ [11]. Для определения явной зависимости разностного соотношения (1.7) от неизвестного Т подставим в него значение а из (1.3) и получим

А1(Тг + С4Т - С2д2Тг/дх22 - С2д3Тг/дх2) = Л^1^ + АЛ? ,

г = 1,...,^ - 1 , 0 <хт <dm , т = 2, 3. (1.10)

Воспользовавшись обозначениями (1.6), имеем вместо (1.10) по координатному направлению х1 дифференциально-разностную схему (А1Т = А1Тг, Л1Т = Л1Тг, так как х2, хз в интервале 0 < хт < dm, т = 2, 3 изменяются непрерывным образом в (1.5)):

А1(Т + сАТп+1/з) = Л1Тп+1/з + Л2Тп + ЛзТп + А1/п ,

г, 2,к =1,...^т - 1 , т =1, 2, 3. (1.11)

Из уравнений (1.7), умножая обе его части на А-1 слева, найдем аг. Это можно сделать, так как матрица, соответствующая этому оператору, трехдиагональная (1/6, 4/6, 1/6 ) при hm = Л, т = 1, 2, 3 ис диагональным преобладанием. Тогда, отмечая, что А-1А1 = Е, получим

аг = А-1Л1Т™+1/з + /п . (1.12)

В левую часть операторного уравнения (1.8) подставим значение Ь из (1.4), а в правую — значение аг из (1.12), тогда она примет вид

А2(С4Т,, + Т, - Сзд2Тг/дх2) = Л2Т?+2/з + А2А-1 ЛТ?^/* + /?) =

= Л2ТП+2/з + (А2АГ1)(Л1ТП+1/з) + А2/п , г = 1,...,т - 1, 2 = 1, ...,N2 - 1, 0 <хз <dз. (1.13)

Используя обозначения (1.6) и соотношение А2А-1 = Е, аналогично (1.10), (1.11) из

(1.13) получим

А2(сТп+2/з + Т) = Л1Тп+1/з + Л2Тп+2/з + ЛзТп + А/,

г,2,к =1, ...,Nm - 1, т = 1, 2, 3. (1.14)

Наконец, в правую часть уравнения (1.9) подставим значение Ьг^, полученное из разностного соотношения (1.8) и предварительно умноженное слева на А-1. Тогда разностная схема (1.9) перепишется так:

Аз(с4Тг ¿к + Т ¿к) = ЛзТП+51 + Аз(А2-1Л2ТП+2/з + А^ЛТ?^ + /п) =

= (АзА-1)^^) + (АзА-1)(Л1Т;П+1/з) + Аз/^к. (1.15)

Учитывая, что АзА-1 = Е, АзА-1 = Е, из операторного уравнения (1.15) окончательно следует уравнение

Аз(с4Тп+1 + Т)= ЛТ п+1/з + Л2Т п+2/з + ЛзТп+1 + Аз/п,

г,2,к =1,...,^ - 1, т = 1, 2, 3. (1.16)

Таким образом, для решения трехмерного эллиптического (1.1) или нестационарного (1.2) уравнений внутри области определения R получены разностные уравнения (1.11),

(1.14), (1.16), аппроксимация и устойчивость которых будет исследована ниже. Отметим, что в логической схеме ИИМ (1.7)-(1.16) неявно присутствует подход Дугласа для решения трехмерного уравнения теплопроводности [7, 8]. Поэтому полученные разностные схемы (1.11), (1.14), (1.16) имеют погрешность аппроксимации в обычном смысле. Так как в итоге получается схема стабилизирующей поправки [8], то двухслойная разностная схема по времени для отмеченных уравнений абсолютно устойчива для любых начальных данных.

2. Постановка задачи и метод построения разностного решения

Пусть требуется решить трехмерное уравнение в частных производных [12, 13] з д дТ з д2Т з дТ

^ жг + ш ^ + Е «-аг - *т + / = о (2.1)

m=1 ^т \ ^^т/ т=1 т=1 и^т

к=т

в параллелепипеде R при w = const и граничных условиях первого рода

Т\Xl=dl = Pl(d1,x2,x3), Т\Х2 =d2 = P2(x1,d2,x3),

Т\x3=d3 = Р3 (x1, x2, d3) , (2.2)

Т \x1=0 = s1(x2,x3), Т\x2=0 = s2(x1,x3),

Т\x3=0 = s3(x1,x2). (2.3)

Наряду с основной краевой задачей (2.1)-(2.3) рассмотрим вспомогательную нестационарную задачу

— - У — (v —1 + wy Q2T +У и—-с T+f-0 (2 4)

dt ¿1V mdxm) ¿1 дхтдхк mdxm 4

k=m

при произвольном начальном условии Т\t=0 = p(x) и с теми же стационарными граничными значениями.

Уравнения типа (2.1) применяются в механике сплошных сред [13]. В общем случае при наличии смешанных производных имеем 27-точечный шаблон в пространстве, а при их отсутствии — 7-точечный крест.

Разностные операторы внутри области определения R записываются в виде

A1T = т 1[h1,i-1Ti-1,j,k + 4h 1,iTi,j,k + h1,iTi+1,j,k]/6h 1,h Л1Т = [V ,i+1,j,k + v1,i,j,k)(Ti+1,j,k - Ti,j,k)/h1,i + (v1,i-1,j,k + v1,i,j,k)x

x(Ti-1,j,k - Ti,j,k)/h1,i-1]/2h 1,u

Л12Т = 2w[(Ti+1 ,j+1,k - Ti+1,j 1,j+1,k + Ti-1,j,k) /h2,j +

+ (Ti+1,j,k - Ti+1,j + Ti -1,j-1,k)/h2,j- 1]/4h 1i,

Л1,3Т = 2w[(Ti+1,j,k+1 - Ti+1,j,k - Ti-1,j,k+1 + Ti-1,j,k)/h3,k + + (Ti+1,j,k - Ti+1,

+ Ti

)/h33,k-1]/4h Ц,

Л233Т = 2w[(Ti ,j+1,k+1 - Ti,j+1,k - Ti,j-1,k+1 + Ti,j-1,k)/h3,k +

+ (Ti j-1,k + Ti,j-1,k-1)/h33,k-1 ]/4h2j,

Л1-3Т = (Л1,2 + Л2,3 + Л1,3)Т,

A1T = [(П1 ,i+1,j,k + 2u1,i,j,k)(Ti+1,j,k - Ti,j,k) + (u1,i-1,j,k + 2u1,i,j,k)x

x (Ti,j,k - Ti-1,j,k)]/6hц. (2.5)

В случае переменных коэффициентов vm = vm(x,t), um = um(x,t) и hm = const (m = 1, 2, 3), используя операторные обозначения (2.5), получим следующую систему разностных уравнений внутри параллелепипеда R:

A1(Tn+1/3 - Тп) = (Л1 + A1)Tn+1/3 + (Л2 + Д2)Тга+

+(Л3 + Д3)Тп + (Л12 + Л2,3 + Л133)Тп + rA1fn , (2.6)

A2(Tn+2/3 - Тп) = (Л1 + Д1)Тп+1/3 + (Л2 + Д2)Тп+2/3+

+(Л3 + Д3)Тп + Л-Тп + rA2fn , (2.7)

Аэ(Тга+1 - Тп) = (Л! + Д1)Тга+1/3 + (Л2 + А2)Тп+2/3+ +(Лз + Дз)Тп+1 + Л1-зТп + тАзР ,

г, к = 1,...,^ - 1, т =1, 2, 3. (2.8)

Практически при исследовании аппроксимации разностных схем (2.6)-(2.8) и программировании алгоритма расчета удобно пользоваться вместо двух последних уравнений (2.7), (2.8) более простыми разностными уравнениями [8]. Вычтем из второго уравнения первое, а из третьего второе, тогда вместо двух последних уравнений системы (2.6) - (2.8) получим

А2тп+2/3 - а1Тп+1/3 = (Л2 + Д2)(Тп+2/3 - Тп)+

+А2(Т + т/)п - А1(Т + т/)п, (2.9)

АзТп+1 - А2Тп+2/3 = (Лз + Дз)(Тп+1 - Тп) + Аз(Т + т/)п-

-А2(Т + т/)п, г, 2,к = 1,...,^ - 1, т = 1, 2, 3. (2.10)

На гранях параллелепипеда хт = 0 , хт = dm (т = 1, 2, 3) имеем

Т^Хи = РЛ^г., 2 = 0,...,N2, к = 0,...,N3,

ТспЙ = 8п+к, 2 = 0,...,N2, к = 0,...,N3. (2.11)

На остальных гранях формулы вида (2.11) получаются одна из другой круговой заменой индексов.

3. Об экономичности метода

Представляют интерес оценка арифметических действий для нахождения численного решения и аппроксимация (абсолютная устойчивость по начальным данным доказана в [9]). Внутри области R получается = Пт=1 - 1) разностных уравнений (2.6)-(2.8) или (2.6), (2.9), (2.10) для определения Б1 неизвестных Т.

На гранях хт = 0, хт = dm (т = 1, 2, 3) имеем также Б2 = + + ЩЩ)

конечных алгебраических выражений (2.11) для определения Б2 заданных функций из (2.2),

(2.3).

С учетом абсолютной устойчивости из [9] расчет по однородным разностным схемам (2.6) - (2.8) экономичен, так как для получения численного решения в узлах области определения понадобится 0(Б1) арифметических действий, пропорциональное числу узлов R.

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

проксимации разностных схем (1.11), (1.14), (1.16) для уравнения (1.2) при c4 = f = 0 и hm = const (m = 1, 2, 3), тогда имеем

т-1 Ai (Tn+1/3 - Tn) = Л1Тп+1/3 + Л2Гп + Л3Тn, (3.1)

т-1(A2Tn+2/3 - AiTn+1/3) = Л2(Тп+2/3 - Tn) + т-1(A2Tn - AiTn), (3.2)

т-1(A3Tn+1 - A2Tn+2/3) = Л3(Тп+1 - Tn) + т-1(A3Tn - A2Tn), (3.3)

i,j,k =1,...,Nm - 1, m = 1, 2, 3.

Из выражений (3.2), (3.3) последовательно находим

г-1А2Тп+2/з = г-1Л3Тп+1 - Лз(Тп+1 - Тп) - т-1(АзТп - А2Тп), (3.4)

т-1А1Тп+1/3 = т-1А2Тп+2/з - Л2(Тп+2/3 - Тп) - т-1(Л2Тп - АгТп) = = т-1АзТп+1 - Лз(Тп+1 - Тп) - т-1(АзТп - А2Тп) - Л2[А-1А3Тп+1 --тА-1Л3(Тп+1 - Тп) - А-1(АзТп - А2Тп) - Тп] - т-1(А2Тп - АхТп). (3.5)

Умножим слева равенства (3.4), (3.5) на А-1, А-1 соответственно, тогда они перепишутся так:

т-1Тп+2/3 = т-1А-1АзТп+1 - А-1Лз(Тп+1 - Тп) - т-1А-1(АзТп - А2Тп), (3.6)

т-1Тп+1/з = А-1[т-1АзТп+1 - Лз(Тп+1 - Тп) - т-1(АзТп - А2Тп)]--А-1Л2[А-1АзТп+1 - тА-1Лз(Тп+1 - Тп) - А-1(АзТп - А2Тп) - Тп]-

-т-1А-1(А2Тп - А1Тп). (3.7)

Найдем разностную схему в целых шагах и покажем, что разностные схемы ИИМ удовлетворяют условию полной аппроксимации [8]. Для этого сложим почленно уравнения (3.1)-(3.3) и получим

т-1(АзТп+1 - АТп) = ЛТп+1/з + Л2Тп+2/з + ЛзТп+1 + т-1(АзТп - А{Тп). (3.8)

Для нахождения погрешности аппроксимации разностной схемы 1 обобщенного ИИМ (3.8) исключим из нее величины Тп+1/з, Тп+2/з при помощи формул (3.6), (3.7), тогда она перепишется так:

т-1Аз(Тп+1 - Тп) = (Л1 + Л2 + Лз)Тп+1 - т[Л1А-1Лз(Тп+1 - Тп)+

+Л1А-1Л2[Тп+1 - Тп - тА-1Лз(Тп+1 - Тп)] + Л2А-1Лз(Тп+1 - Тп)}.

В результате для уравнения (1.2) при f = 0 и с4 = 0 имеем окончательно

дTг,3,k/дt + 0.5тТ^к = £ стд2Тг,3,к/дх2т + о(т + ^ .

Таким образом, погрешность аппроксимации разностной схемы (3.8) на решениях уравнения (1.2) имеет первый порядок по времени и второй — по пространству.

Для уравнения (2.4) локальная погрешность аппроксимации разностной схемы (2.6)-(2.8) есть O[Y,a=i(ha,v — ha,v) + т] (v = i при a =1, v = j при a = 2, v = k при a = 3), для hm = cm (cm = const, m = 1, 2, 3) — 0(т + hf + hf + h|), а при w = um = 0

и Vm = Cm - 0[т + £m=1 (hm,q — hm,qhm,q-1 + ^q-i)] для hi,i — hi— = h2,j — hf,j-i =

h3,k — h3kk-i. Граничные условия типа (2.11) точно аппроксимируют соответствующие граничные условия из (2.2), (2.3), так как они заданы как функции пространственных координат xm (m = 1, 2, 3) явно на целом слое.

В работе [9] в приближении замороженных коэффициентов [6, 8] найдено условие безусловной устойчивости явно-неявных разностных уравнений (2.6)-(2.8) или (2.6), (2.9), (2.10) по начальным данным. Для n(T,hm,qm) — коэффициента возрастания гармоники exp[I(qiX\ti + q2X2,j + дзХз,к)], I = л/—Т, qm = 0, ±1, ±2, ..., т = 1, 2, 3 — спектральное уело-

вие Неймана при w = с4 = 0 имеет вид п = (A - IB)/(C - ID), \п\ = [(A2 + B 2)/(C2 + D2)]0-5, где

A = F + s; C = F + a; D = E + b; B = E + 12 = -1;

F = Г3Г2Г1 + ri(a2a3 - b2h) + r2(aia3 - Ь^з) + r3(a2ai-

-b2bi) + aia2a3 - aib2b3 - a2bib3 - a3b2bi;

s = air2(r3 - ri) + ria2(r3 - Г2); a = r3r2ai + r3a2ri + a3^ri;

E = ri(b2a3 + a2b3) + r2(bia3 + a^) + r3(b2ai + a2bi) +

+a2a3bi + aia3b2 + aia2b3 - b^b3;

# = bir2(r3 - ri) + rib2(r3 - r2); b = r3r2bi + r3b2ri + b3r2ri;

am = irvmhm2 sin2 qmhm/2; bm = TUmh^ sin qmhm; rm =(1 + 2 cos2 qmhm/2)/3, m = 1, 2, 3.

Как было отмечено во введении, при получении стационарного решения методом установления представляет интерес нахождение оптимального шага по времени, при котором численное решение выйдет на стационарное за минимальное число итераций K. Следуя известному алгоритму из монографии [6, с. 405-407], получим оценку K в кубе Q: (0<xm< 1, m =1, 2, 3) для (2.6), (2.9), (2.10) при vm = hm = const, w = c4 = 0, m =1, 2, 3.

Собственные функции разностного оператора из (2.6), (2.9), (2.10) в кубе Q на равномерной сетке равны

3

Wq^qs (x) = П sin(^qmXm), 1 < qm < Nm - 1, m =1, 2, 3. (3.9)

m=i

Подставим (3.9) в разностную схему (2.6), (2.9), (2.10) и, полагая w(П+^З (x) = nqw^^w (x), получим множитель роста гармоник nq = nq (т, hm, qm) аналогично [9], но для простоты дальнейших выкладок будем считать, что vm = v, um = u, hm = h, qm = q, m = 1, 2, 3. Так как нас интересуют асимптотические оценки, то, считая Nm = N достаточно большими, можно положить AmT ^ ЛmT (h-2 ^ h-i , v ~ u) в (2.6), (2.9), (2.10). Тогда

nq = (r3 + 3ra2 + a3)(r3 + 3r2a + 3ra2 + a3)-i , (3.10)

a = 4Tvh-2 sin2(nqh/2), r = [1 + 2 cos2(nqh/2)]/3, m = 1, 2, 3.

В [6] показано, какие гармоники затухают наиболее медленно и тем самым больше всего препятствуют выходу на стационарный режим. Значение nq меняется в пределах 0 < nq < 1. Сначала nq монотонно убывает при увеличении номера q, а затем возрастает при N ^ ж. При этом, как следует из дальнейших выкладок, наибольшим значение nq может быть либо при q =1, либо при q = N - 1. Сделаем оценки для a и r в (3.10):

при q =1 sin(nh/2) = sin(n/2N) & n/2N, ai = Tvn2, r = 1;

при q = N - 1 sin[n(N - 1)h/2] = sin[n(N - 1)/2N] « 1,

aN-i = 4tvN2, r =1/3.

В результате экстремальные множители можно представить в виде

ni = (1 + 3ai)(1 + 3ai + 3ai)-i « (1 + 3a2)(1 - 3ai) « 1 - 3ai. (3.11)

щ-1 = (1/27 + _ + а%-1)(1/27 + _ + а%_ + aN-1/З)-1 =

= (1 + 1/aN _1)(1 + 1/aN-1 + 1/3aN _1)_\

щ-1 ~ (1 + _1)[1 - (1/aN-1 + 1/3а%-1)] « 1 - 4/(3aN-1). (3.12)

(считается, что (ту)2 ^ (ту) при ту < 1).

Тогда то выбирается так, чтобы п1 (то) = nN_1 (то). Из (3.11), (3.12) получим то: 3а1 = 4/(3aN_1) или

то = [у(6п)2/3]-1 N _4/3 = у_^_4/3/7.08. (3.13)

Следовательно, минимальное необходимое число итераций, при котором разностное решение выйдет на стационарное с заданной точностью е [6], есть

[ехр(-4/^ _1)]К = е

или

К (то) = 3N4/31п(1/е)/4 п. (3.14)

4. Результаты численных расчетов

Оценку скорости сходимости (3.14) проведем на решении модельных (тестовых) задач. Сначала найдем численное решение эллиптического уравнения общего вида

3 3 3

ci £ дТ/дхт + с £ д2T/3x2m + сз ]Т д2T/dxm3xk - f = 0, (4.1)

m=1 m=1 m=1,k=m

с постоянными коэффициентами переноса и граничными условиями первого рода, у которых известно точное аналитическое решение, например, T =1 + 2m=i xzm:

33

f = x— + c2z(z - xш2,

m=1 m=1

T\x 1=0 = (1 + x3 + x2), T |x2=0 = (1+ x1 + x3), T 1x3=0 = (1 + x1 + x2), T\x 1=1 = (2 + x3 + x2), T\x2=1 = (2 + x1 + x3), T\хз=1 = (2 + x1 + x2). (4.2)

Для решения задачи (4.1), (4.2) воспользуемся разностными уравнениями (2.6), (2.9), (2.10). Были взяты следующие значения входных данных для опорного варианта c1 = c2 = 1, c3 = 2, z = 2 и ^ = 1. Программа составлена на языке Фортран-77 (Power Station-4), расчет проводился на ПЭВМ Pentium-2 (133 МГц) с двойной точностью. Число итераций K(т0) отслеживалось по относительному изменению нормы вектора погрешности

\\yn\\ = max \(T - Tn)/T\,

где T — заданное точное решение дифференциальной задачи (4.1), (4.2); Tn — численное решение. В табл. 1 приведены результаты расчетов K(т0) с \\yn\\ < £, £ = 0, 1 для двух вариантов: Tн = 1, Tн = 0 при c2 = 1 и различных c1, c3.

Видно, что выбор начального приближения может влиять на скорость сходимости итерационного процесса. Отметим, что время расчета варианта при N = 81 составило 30 минут. Однако расчет при c1 = 1, c3 = 2 и прочих одинаковых входных данных общего

уравнения (4.1) по количеству итераций изменяется незначительно. При увеличении в (4.1) первого слагаемого (с1 = 500) для N = 81, Тн = 1, с3 = 2 имеем К = 32, что также не ухудшает предложенный алгоритм.

Таблица 1

Ci = с3 = 0 С1 = 1 сз = 2

Тн 0 1 1 1

N 21 41 81 21 41 81 21 41 81

К 24 55 133 17 40 96 18 43 104

Для проверки эффективности метода и оценки теоретической формулы (3.14) рассматривалась разностная задача Дирихле в кубе для уравнения Пуассона с постоянными и переменными коэффициентами:

дТ '

3 д m=1 Xm

'дхг.

TI:

f, 0,

xm

е (0,1),

(4.3)

(4.4)

v1(x) = 1 + C ^ (xm — 0.5)2, C = const,

m=1

^(х) = 1 + С[0.5 - (х1 - 0.5)2 - (х2 - 0.5)2 + (хз - 0.5)2], 1>з(х) = 1 + С[0.5 + (Х1 - 0.5)2 + (Х2 - 0.5)2 - (хз - 0.5)2].

Источник f в правой части уравнения (4.3), как и в (4.1), подбирался, чтобы пробная функция Т = Пт=1 хт(1 - хт) была точным решением. Отметим, что в двумерном случае постановка задачи (4.3), (4.4) совпадает с тестовой задачей, приведенной в [1].

В табл. 2 К отвечает численному решению краевой задачи (4.3), (4.4) при \\уп\\ < е, а К1, К2 соответствуют теоретической оценке скорости сходимости (3.14) для е1 = 0.1, е2 = 0.01 и С = 0. Видно, что ^ = (К - ^)/К{, ^ = 1, 2) улучшается с ростом N и е. Это говорит о сходимости алгоритма и правильности полученных формул (3.13), (3.14).

Таблица 2

v

£i = 0.1 £i = 0.01

N = 21 N = 41 N = 81 N = 21 N = 41 N = 81

К 36 83 203 К 73 165 399

Кг 32 78 193 k2 64 156 386

h 0.13 0.064 0.051 ¿2 0.14 0.061 0.035

В табл. 3 приведено численное решение этой задачи при vm = const, е = 0.1 для C =2, C = 62 и прочих одинаковых входных данных; параметр т0 из (3.13) находился при v = (1 + 0.5C)/2. Таким образом, и при переменных коэффициентах переноса сходимость итерационного процесса при помощи ИИМ не ухудшается.

Таблица 3

C = 2 C= 62

N 21 41 81 21 41 81

К 22 51 124 29 67 161

Выводы

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

2. На тестовых примерах показана сходимость итераций, и при постоянных коэффициентах переноса для уравнения Пуассона найдена теоретическая оценка скорости сходимости итерационного метода.

Список литературы

[1] Самарский А. АНиколаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.

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

[3] Ильин В. П. Разностные методы решения эллиптических уравнений. Новосибирск: Наука, 1970.

[4] Коновалов А. Н. Численное решение задач теории упругости. Новосибирск: Наука, 1968.

[5] Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971.

[6] Калиткин Н. Н. Численные методы. М.: Наука, 1978.

[7] Douglas J. Alternating direction methods for three space variables // Numerische Math. 1962. Vol. 4, No. 6. P. 41-63.

[8] Яненко Н. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

[9] Гришин А. М, Якимов А.С. Обобщение итерационно-интерполяционного метода для решения трехмерного параболического уравнения общего вида // Вычисл. технологии. 1999. Т. 4, № 2. С. 26-41.

[10] Якимов А. С. Об одном методе расщепления // Числ. методы механики сплошной среды. 1985. Т. 16, № 2. С. 144-161.

[11] Гришин А. М., Берцун В. Н., Зинченко В. И. Итерационно-интерполяционный метод и его приложения. Томск: Изд-во ТГУ, 1981.

[12] Ладыженская О. А, Уральцева Н. Н. Линейные и квазилинейные уравнения эллиптического типа. М.: Наука, 1973.

[13] Лойцянский Л. Г. Механика жидкости и газа. М.: Наука, 1973.

Поступила в редакцию 23 ноября 1999 г. в переработанном виде — 11 ноября 2000 г.