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

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

Том 12, № 1, 2007

ПРИМЕНЕНИЕ ИТЕРАЦИОННО-ИНТЕРПОЛЯЦИОННОГО МЕТОДА ДЛЯ РЕШЕНИЯ ТРЕХМЕРНОГО ВОЛНОВОГО УРАВНЕНИЯ*

А. М. Гришин, А. С. Якимов Томский государственный университет, Россия e-mail: [email protected], [email protected]

A numerical algorithm for solving the 3D wave equation based on an iteration-interpolation method is developed. The obtained scheme is absolutely stable. The scheme provides the second order approximation in time and space. Convergence of the scheme is examined and the results of the test calculation is presented.

Введение

В работах [1-7] предложен итерационно-интерполяционный метод (ИИМ) для решения различных уравнений математической физики. Работа [2] переведена на английский язык [3], здесь для одного частного случая утверждалось, что погрешность аппроксимации — O(h2k), где h — шаг разностной сетки, а k = 1, 2 — номер итерации. В [4] даны алгоритмы ИИМ для решения уравнений параболического, гиперболического и эллиптического типа, допускающих точное решение, и приведена оценка точности полученных численных решений, а в [5] — развитие ИИМ для решения трехмерного уравнения теплопроводности. Обобщение ИИМ на трехмерный случай для уравнения параболического типа сделано в [6]. Современное состояние ИИМ изложено в [7].

Суть алгоритма ИИМ заключается в использовании метода последовательных приближений на каждом элементарном отрезке разностной сетки, соответствующей области определения исследуемой краевой задачи, в результате чего точность приближенного решения можно повысить путем как уменьшения шага разностной сетки, так и увеличения числа итераций [1]. В дальнейшем установлена связь ИИМ с теорией сплайнов [3] и даны примеры применения ИИМ для решения некоторых нелинейных краевых задач [4].

В монографии [8] приведен обзор методов решения уравнения с частными производными второго порядка гиперболического типа [9-12]. В [10] решена задача Коши. В работе [11] для двумерного уравнения колебаний имеет место трехслойная аддитивная разностная схема, а в трехмерном пространстве — четырехслойная с погрешностью аппроксимации

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00714).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

0(т + |Л-|2). В [11] отмечено, что вид локально-одномерных схем для гиперболического уравнения второго порядка зависит от числа измерений по пространству.

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

1. Алгоритм обобщения

итерационно-интерполяционного метода для решения трехмерного волнового уравнения

Чтобы упростить дальнейшие выкладки, логическую схему ИИМ для волнового уравнения с постоянными коэффициентами переноса внутри параллелепипеда Я: х = (х1,х2,х3), 0 < хт < ¿т, т =1, 2, 3, при 0 < Ь < Ьк приведем к виду

Для построения разностного аналога уравнения (1.1) введем сетку с координатами узлов хм = гк1 (г = 0,1, 2,...,Ы), х2,з = ]Н2, 3 = 0,1, 2,...,Ы, х3к = кк3, к = 0,1, 2,...,Ы, Ьп = пт, п = 0,1, 2,... , Нт = ¿т/Ы, т = 1, 2, 3, где Нт — шаги по пространству, т — шаг по времени. При написании операторов разностных схем вводятся верхние промежуточные (формальные) слои п, п + 1/3, п + 2/3, а также п — 1, п + 1 — начальный и конечный (фактические) слои по времени Ь.

Используя результаты статьи [6], введем следующие обозначения (одна и две точки "сверху" приписываются частной производной по времени):

Разностные операторы А1аг, Л1пг, А2Ьг,з, Л2пг,з,А3пгз,к, Л3Пгззк внутри области определения Я взяты из статьи [6] и при постоянных шагах по пространству имеют вид

(1.1)

а(Ь,х) = П — у2 д2и/дх\ — ь3д2п/дх\ Ь(Ь,х) = П — у3д 2и/дх\.

.2.

3;

(1.2) (1.3)

А^г = (а— + 4аг + а+)/6, Л1 п = ь^п— — 2щ + щ+^/Ь^, г = 1,...,Ы — 1, 0 <хт < ¿т, т = 2, 3, А2Ы,3 = (Ьгз-1 + 4Ьг,з + Ьг,3+1)/6, Л2Щ,з = У2 (щ,— — 2щ,3 + Щ,з+1 )/Н2, г,3 = 1,...,Ы — 1 , 0 <х3 < ¿3, Азпг,з,к = (пг,з,к-1 + 4пг,з,к + щ^к+^/б, Лзпг,з,к = уз(пг,з,к-1 — 2пг,з,к + пг,з,к+1)/кго г,3,к = 1,...,М — 1.

(1.4)

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

Aiu = (w¿-i jk + 4ui,j,k + u¿+i,j,fc)/6, Л1М = vi(ui-i)j)fc - 2uij£ + u¿+i,j,fc)/h?,

A2M = (ui,j-i,k + 4ui,j,k + Ui,j+i,k)/6, Л2и = V2(ui)j-i)fc - 2ui)j)fc + Ui,j+i,k)/h|,

A3U = (ui,j,k-i + 4ui,j,k + ui)j,k+i)/6, Лзu = V3(ui)j)k-i - 2ui,j,k + ui^k+O/hJ,

i,j,k =1,...,N - 1. (1.5)

Следуя алгоритму из [1, 5-7] и используя обозначения (1.2), (1.4) для уравнения (1.1) последовательно запишем следующие соотношения. По направлению координаты xi

a = vid2u/dxi + f

разностные схемы, полученные в результате поочередного интегрирования по пространственным переменным, имеют вид [6]

A i ai ^iun+i/3 + Af-\

i = 1,...,N - 1, 0 <xm < dm, m = 2, 3. (1.6)

По направлению координаты x2

bi = v2d2ui/d^2 + ai, i = 1,..., N - 1

разностные уравнения записываются как

A2bi,j = Л2Un+2/3 + A2ai,j,

i, j = 1,... ,N - 1, 0 < x3 < d3. (1.7)

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

ui,j = v3d2ui;j/dx3 + bi,j, i, j = 1,... , N - 1,

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

A3uiJ;k = Лзun+k + A3biJ;k, i, j, k = 1,..., N - 1. (1.8)

Отметим, что разностные уравнения вида, например, (1.6) получаются непосредственно из логической схемы ИИМ [1, 7], а для нахождения разностных уравнений (1.6)-(1.8) используются операторы из (1.4). Для получения конкретного вида соотношения (1.6) подставим в него значение a из (1.2):

Ai(u - V2d2ui/dx2 - V3d3ui/dx3) = Л^3 + Aifn-\

i = 1,...,N - 1, 0 < xm < dm, m =2, 3.

(1.9)

Из выражений (1.5), (1.9) следует, что оператор Al может действовать на любую сеточную функцию, осредняя ее в направлении xl, но не может действовать на частные производные функций, которые берутся по пространственным переменным с других координатных направлений. Используя (1.5), вместо (1.9) по координатному направлению xl имеем дифференциально-разностную схему

Aiu = Amn+l/3 + A2un-1 + A3un-1 + Aifn-1,

i,j,k =1,...,N - 1. (1.10)

Далее найдем ai из уравнения (1.6), умножая обе части его на A-1 слева. Это можно сделать, так как матрица, соответствующая этому оператору, трехдиагональна (1/6, 4/6, 1/6) при hm = const, m = 1, 2, 3, и имеет диагональное преобладание. Тогда, замечая, что A-lAl = E, получим

ai = A-lAiun+l/3 + fin-1. (1.11)

В левую часть операторного уравнения (1.7) подставим значение b из (1.3), а в правую — значение ai из (1.11). В результате (1.7) примет вид

A2(uj - ьзд2и%]/dx\) = А2пП+ 2/3 + A2(A-lAiun+1/3 + fin,'1) =

= Aj + (A2A-1 )(AU+l/3) + A2%r\

i,j = 1,...,N - 1, 0 < x3 < d3. (1.12)

Используя обозначения (1.5) и соотношение A2A-1 = E (при hm = const Al = A2 = A3), аналогично предыдущему (1.9), (1.10) из (1.12) получим

A2U = Aiun+l/3 + A2Un+2/3 + A3un-1 + A2fn-1,

i,j,k =1,...,N - 1. (1.13)

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

A3uhJ,k = A3 + A33(A-2lA2unjl/3 + A-Aiun+l3 + fjl) =

= ^un+l + (A3A-l)(A2un+2k/3) + ^A-XAiun+f3) + A3 ft-. (1.14)

Замечая, что A3A-1 = E, A3A-1 = E, из операторного уравнения (1.14) окончательно получим следующее уравнение:

A311 = Aiun+l/3 + A2un+2/3 + A3un+l + A3P-1,

i,j,k =1,...,N - 1. (1.15)

Таким образом, для решения трехмерного волнового уравнения (1.1) внутри области определения R получены разностные уравнения (1.10), (1.13), (1.15), аппроксимация и устойчивость которых будет рассмотрена ниже.

2. Постановка задачи и метод построения разностной схемы

Пусть требуется решить трехмерное уравнение в частных производных второго порядка гиперболического типа [9-12]

, . d2u , ч du v^ d

g(t'x) at* + s(t'x) dt = ml dX

m=1

, . du Vm(x,t)-

+

3 r\2 3 r\

d2u , s du

+wj]d-d— + Wm(x,t) d— + / (x,t) (2.1)

dxmdxk dxm

m=1 m=1

k=m

в параллелепипеде R: x = (x1,x2,x3), 0 < xm < dm, m =1, 2, 3, для w = const, |w| < 1 [11] с начальными условиями

du

u|t=0 = Pi(x), — = P2 (x) (2.2)

dt t=o

и для простоты анализа при граничных условиях первого рода

u|Xl=dl = s1(t,x2,x3), u|x2=d2 = s2(t,x1,x3), u|x3=d3 = s3 (t,x1,x2), u|Xl=0 = q1(t,x2,x3), u|x2=0 = q2(t,x1,x3), u|x3=0 = q3(t,x1,x2). (2.3)

В общем случае при наличии смешанных производных имеем 27-точечный шаблон в пространстве, а при их отсутствии — семиточечный крест в параллелепипеде R. Используя результаты, полученные в разд. 1 и в [5], запишем разностные операторы внутри области определения R. Для простоты дальнейшего анализа положим g = const, s = const (g = s = 1), vm = const, wm = const, m =1, 2, 3 в (2.1). При переменных коэффициентах переноса в (2.1) разностные операторы подробно выписаны в [6, 7]. Тогда имеем

u j = т-2 «+ - 2u j + tj = (unn^+k - un:-;k)/(2r),

A1ui;j,fc = W1(ui+1)j)k - ui-1)j)k)/(2h1),

Лl)2U = 2w(ui+1)j+1)k - ui+1)j-1)k - ui-1)j+1)k + ui-1)j-1)k)/(4h1^2),

Лl)зu = 2w(ui+1ijik+1 - ui+1)j-fc-1 - uj-1J;fc+1 + ui-1)j-fc-1)/(4h1h3),

Л2)зu = 2w(ujJ+1;fc+1 - ui)j+1)k-1 - ujj-1)fc+1 + ui)j-1)k-1)/(4h2h3),

A/ = (/¿-1,j,k + 4/i)j-fc + /¿+1,j,k )/6, Лl-зu = (Л12 + Л1,3 + Л2,3 )u,

=1,...,N - 1. (2.4)

Операторы Am, Am, m = 2, 3, по другим координатным направлениям x2,x3 получаются циклической заменой индексов [5]. На n временном слое решение можно найти из формулы Тэйлора и уравнения (2.1), если воспользоваться разложением из [13], которое в общем случае перепишется так:

dun-k

uij,k = ui,j,fc + т dt +

П—1

о / г\ \ г\ 3 с\2 3 с\

d dui j fc\ dui j k d2Ui j k du

dxm V m dxm ) dt ^ dxmdxs ^ m dxr.

Л + &i)un+1/3 + Л + &2)un+2/3 + (Лз + &3)un+1 + Ai-3Un-1 + Ä3fn-1,

m=1 4 7 m=1 m=1

s=m

i,j,k =1,..., N- 1. (2.5)

Используя обозначения (1.5), (2.4), получим следующую систему разностных уравнений внутри параллелепипеда R (при написании разностных схем предполагается, что (un-1)1/3 = (un-1)2/3 = (un-1)1 = un-\ (un)1/3 = (un)2/3 = (un)1 = un)):

T-2Ä1(un+1/3 - 2un + un-1) + (2t)-1 Ä1(un+1/3 - un-1) =

= (Л1 + A1)un+1/3 + (Л2 + A2)un-1 + (A3 + A3)un-1 + A1-3un-1 + Äfn-1; (2.6)

T-2Ä2(un+2/3 - 2un + un-1) + (2t)-1 Ä2(un+2/3 - un-1) = (Л1 + A1)un+1/3 + (Л2 + A2)un+2/3 + (Л33 + A3)un-1 + A—un-1 + Ä2fn-1; (2.7) t-2Ä3(un+1 - 2un + un-1) + (2t)-1 Ä3(un+1 - un-1)

i,j,k =1,...,N - 1. (2.8)

Здесь un-1 задано из первого начального (2.2), а un находится по формуле (2.5) при использовании разностных операторов Лm, Am, m = 1, 2, 3, Л1>2u, Л1^, Л^^ из (1.5), (2.4).

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

T-2(Ä2un+2/3 - Älun+1/3) + (2t)-1(Ä2un+2/3 - Älun+1/3) = = (Л2 + A2)(un+2/3 - un-1) + T-2Ä2(2un - un-1) + (2t)-1Ä2un-1-

-[t-2Ä1(2un - un-1) + (2t)-1Ä1un-1] + Ä2fn-1 - Ä1fn-1, (2.9)

T-2(Ä3un+1 - Ä2un+2/3) + (2t)-1(Ä3un+1 - Ä2un+2/3) = = (Л3 + A3)(un+1 - un-1) + T-2Ä3(2un - un-1) + (2t)-1Ä3un-1--[t-2Ä2(2un - un-1) + (2t)-1Ä2un-1] + Ä3fn-1 - Ä2fn-1,

i,j,k =1,...,N - 1. (2.10)

На гранях x1 = 0 , x1 = d1 имеем

uN+j = sSL> j = 0,...,N, k = 0,...,N,

u

П,5 = СЙ> 3 = к = 0,...,М. (2.11)

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

3. Устойчивость разностной схемы ИИМ по начальным данным

В работе [5] в приближении замороженных коэффициентов [9, 13] методом Фурье найдено условие безусловной устойчивости явно-неявных разностных уравнений (2.6)-(2.8) или (2.6), (2.9), (2.10) по начальным данным при g = f = 0, s = 1, vm = wm = const, m =1, 2, 3. Представляет интерес оценка величины n — множителя роста гармоники при переходе со слоя на слой [13] при g = 0. Сделаем стандартную подстановку [13]

Un+1 = ПИ™-1, «П-fc = еХр[1 (glXl,i + g2X2j + g3X3,fc)], (3.1)

где / = \J — 1, gm = 0, ±1, ±2,..., m = 1, 2, 3 — гармоники. Подставим (3.1) в (2.6), (2.9), (2.10) при g = 1, s = f = p2 = 0 и для сокращения дальнейших выкладок введем обозначения

0m = 4т2Vmhm2 sin2(gmhm/2), Xm = sin(gmfrm), 6m = T^W^Xm, rm = [1 + 2cos2(gmhm/2)]/3, m = 1, 2, 3,

Cl = T2(hih2)-1XiX2, C2 = T2 (h2h3)-1X2X3, C3 = T2(h hi)-1 X1X3,

333

c = 2wJ^cm, a =^2 Om, e = bm, Y = в — a — c. (3.2)

m= 1 m= 1 m= 1

Тогда по каждому из направлений получим

u

n+1/3(r1 + O1 — /61) = 2unr1 + un-1(—Г1 — O2 — O3 + Ib2 + /63 — C); (3.3)

ига+2/3(г2 + «2 - /62) = Г1Мга+1/3 + 2иП(г2 - Г1) + Мга-1(Г1 - Г2 + «2 - /62); (3.4) ига+1(гз + «3 - /63) = Г2Мга+2/3 + 2ип(гз - Г2) + Мга-1(Г2 - Г3 + «3 - /63), (3.5)

= ип-1[1 + 0, 5(/в - а - с)].

Выразим ига+1/3 из (3.3) и подставим в (3.4), затем найдем ига+2/3 из полученного соотношения и внесем его в (3.5). В результате для определения п окончательно имеем уравнение

5 + и - / (Т + V)

n

W-/D

= [(S + и )2 + (Г + V )2Р (36)

|П| = (W2 + D2)0,5 , (36)

где

/2 = — 1,S = F + p — Г1Г2С; W = F + о; D = E + 6; T = E + q;

F = Г3Г2Г1 + Г1 (0203 — 6263) + Г2 (0103 — 6163) + r3(o2O1 — 6261) + O1O2O3 — O16263 — O26163 — O36261; p = О1Г2 (Г3 — Г1) + Г102(Г3 — Г2), о = Г3Г2О1 + Г3О2Г1 + О3Г2Г1;

E = Г1 (62O3 + 0263) + Г2 (61O3 + 0163) + Г3 (62O1 + 0261) + 02 0361 + O1O362 + O1O263 — 616263; q = 61Г2(Г3 — Г1) + Г162(r3 — Г2), 6 = Г3Г261 + Г362Г1 + 63Г2Г1;

Ф = 61 (Г3Г2 — Г1Г2) + 62(Г3Г1 — Г1Г2) + (Г3 — Г2) (0162 — 6102);

ф = Г1Г2Г3 + 01 (Г3Г2 — Г1Г2) + 02(Г3Г1 — Г1Г2) + (Г2 — Г3Х0201 + 6162);

U = y(Ф + Ф); V = 1(Ф - Ф). (3.7)

Необходимое условие устойчивости Неймана разностных схем записывается как Inl < 1 [9, 13]. В силу того что | singmhml < 1, | cos gmhml < 1, найдем оценку сверху для п В частности, знаменатель будет минимален в (3.6) при gmhm/2 = kn, m = 1, 2, 3, k = 1, 2,... Тогда из выражений (3.2), (3.6), (3.7) имеем rm = 1, singmhm = singmhm/2 = 0, m = 1, 2, 3, p = q = c = b = a = 0, E = T = D = 0, U = V = 0, S = W = F, F = 1 и

Inl = F = 1. (3.8)

Из (3.6), (3.8) видно, что разностные схемы (2.6), (2.9), (2.10) безусловно устойчивы для т > 0.

4. Об экономичности метода и погрешности аппроксимации разностной схемы

Представляют интерес число арифметических действий для нахождения численного решения и аппроксимация (абсолютная устойчивость по начальным данным доказана ниже). Внутри R: x Е (0 < xm < dm, m = 1, 2, 3) получается S\ = (N — I)3 разностных уравнений (2.6)-(2.8) или (2.6), (2.9), (2.10) для определения S\ неизвестных u.

На гранях xm = 0, xm = dm, m = 1, 2, 3, имеем также S2 = 6N2 конечных алгебраических выражений (2.11) для определения S2 заданных функций из (2.3).

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

Пусть существуют ограниченные вплоть до четвертого порядка производные по пространству и по времени от искомой функции. Не умаляя общности, найдем погрешность аппроксимации разностных схем (2.6), (2.9), (2.10) для уравнения (2.1) при s = f = w = 0, vm = const, wm = 0, g = vm = 1, m = 1, 2, 3, тогда имеем

т-2Ai(un+1/3 — 2un + un-1) = Лгип+1/3 + A2un-1 + Лзип(4.1)

т-2(A2un+2/3 — Aiun+1/3) = A2(un+2/3—un-1)+T-\2(A2un — Aiun) — (A2un-1 — Aiun-1)]; (4.2)

т-2(A3un+1 — A2un+2/3) = A3(un+1 — un-1) + т-2[2(A3un — A2un) — (A3un-1 — A2un-1)], (4.3)

i,j,k = 1,...,N — 1. Из выражений (4.3), (4.2) последовательно находим [7]

т-2 A2un+2/3 = т-2A3un+1 — A3(un+1 — un-1) + т-2[A3un-1 — A2un-1 — 2(A3un — A2un)],

un+2/3 = (a-1 A3)un+1 — т 2A-lA3(un+l — un-1) + A-1(A3un-1 —

—A2un-1) — 2A-1(A3un — A2un), (4.4)

un+1/3 = (A-1A2)un+2/3 — т2А-1А2 (un+2/3 — un-1) + A-1[A2un-1 — —Axun-1 — 2(A2un — A1un)],

un+i/3 = (аг1А2)[(А-1АЗК+1 - т2А-1Лз(мга+1 - ига-1)+ +А"1(Азига"1 - A2un-1) - 2А"1(Азига - ^u")]--т2А-1Л2[(А2"1Аз)мга+1 - т2А2"1Лз(ига+1 - un-1)+ +А"1 (Азип-1 - Ли""1) - 2А"1 (Азип - Ли") - un-1]+

+A-1[(A2Ura"1 - Aiu""1) - 2(A2un - Aiun)]. (4.5)

Найдем разностную схему в целых шагах и покажем, что разностные схемы ИИМ — схемы с абсолютной погрешностью аппроксимации на уравнении. Для этого сложим почленно уравнения (4.1)-(4.3). Тогда получим

т"2[Азмга+1 + Ai (и""1 - 2un)] = Л1Мга+1/3 + Л2Мга+2/3 + Лзмга+1+

+2т"2(Азмга - Aiun) - т"2(Азмга-1 - Aiun-1). (4.6)

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

т"2Аз(мга+1 - 2un + и""1) = Л1((А"1Аз)мга+1 - т2A"%(ura+1 - un-1)+

+А"1(Азип"1 - A2un-1) - 2А-1(Азмга - A2un) - т2А"%х х [(А-1Аз)мга+1 - т2А"1Лз(мга+1 - un-1)+ +А"1 (Ази""1 - A2un-1) - 2А"1 (Ази" - A2un) - un-1]+ +A"1(A2Un"1 - AiMra-1) - 2A"1 (A2«ra - Ai«ra)}+ +Л2[(А-1АзК+1 - т2А"1Лз(мга+1 - и""1) + А"1(Азмга"1 -

-А2Мга-1) - 2А"1(Азига - А2«га)] + Лзмга+1. (4.7)

Используя свойство переместительности [14] симметричных матриц Am, Лт, а также

-1 А

m Aq

Am1Aq = E, m, q = 1, 2, 3, из (4.7) для уравнения (2.1) окончательно получим

+т2 ^ = ± ^+42+£ *

m=1 \ m=1

Таким образом, погрешность аппроксимации разностных схем (4.1)-(4.3) на уравнении (2.1) при s = f = w = 0, vm = const, wm = 0, g = vm =1, m =1, 2, 3, имеет второй порядок по времени и по пространству. На уравнении (2.1) локальная погрешность аппроксимации

разностной схемы (2.6)-(2.8) — O ( т2 + ) при vm = const, wm = const, m = 1, 2, 3.

m=1

Первое начальное условие из (2.2) удовлетворяется точно, а второе — с погрешностью аппроксимации 0(т2) согласно (2.5) [13]. Граничные уравнения из (2.11) на соответствующей грани из (2.3) аппроксимируются точно, так как они заданы от пространственных координат т = 1, 2, 3, и времени явно на целом слое.

5. Сходимость и примеры применения метода

Сходимость разностных схем ИИМ установим на конкретном примере при решении частного трехмерного уравнения (5.1) с начальными условиями (5.2) и граничными значениями первого рода (5.3) в области R : (0 < xm < b, m = 1, 2, 3), 0 < t < tk, с постоянными коэффициентами переноса Cm = const, m = 1, 2,33, 4, на последовательности сгущающихся сеток по т и hm, m = 1, 2, 3:

_ d2u ^ du ^ ^ d2u ^ d2u ^ ^ du . ,

C W + C at = C £ sXm + w E aXXmaXXk: + C £ dxm + F t); <5Л)

m=l m m=\ m=l

k=m

du

u\t=0 = (xi - b)z + x2 + xz3 + 1 = U, — |t=0 = 0; (5.2)

u\Xl=o = (1 + ty)[(-b)z + x2 + xZ + 1], u\Xl=b = (1 + ty)x + xZ + 1), u\X2=o = (1 + ty)[xi - b)z + xZ + 1], u\x2=b = (1 + ty)[bz + (xi - b)z + xZ + 1], u\x3=o = (1 + ty)[(xi - b)z + x2> + 1], u\x3=b = (1 + ty)[bz + (xi - b)z + x2> + 1]. (5.3) Легко видеть, что задача (5.1)-(5.3) имеет точное решение

u = (1 + ty) [(xi - b)z + x2 + xZ + 1] (5.4)

при условии, что источник F в уравнении (5.1) взят в виде

F = Uy[CA(y - 1)ty-2 + CZty-i]-

-(1 + ty)z{C2 x (z - 1)[xZZ-2 + (xi - b)z-2 + xz2-2] + Ci[xZ-i + (xi - b)z-i + xz2-i]}.

При исследовании сходимости разностных уравнений (2.6), (2.9), (2.10) для решения тестового примера воспользуемся следующими значениями опорных данных: C4 = 10-4, C2 = 1, Ci = w = 0, z = 4, y = 4, b = 1, N = 21, hm = 0, 05, m = 1, 2, 3. Программа составлена на языке Фортран-90, расчет проводился на ПЭВМ Pentium-4 (2.5 ГГц, транслятор Power Station-5) с двойной точностью.

В табл. 1 дается максимальная относительная погрешность е = (u - u) 100 %/u (u — точное решение (5.4), u — приближенное численное решение) в момент времени tk = 3 при различных значениях CZ. Видно, что фактический порядок точности разностных схем по т лучше теоретического, полученного в четвертом разделе статьи на гладких функциях.

Таблица 1.

Cz 0 0.05

т 0.004 0.002 0.001 0.0005 0.004 0.002 0.001 0.0005

59 17.5 1.54 0.1 42.4 11.2 1.1 0.26

Таблица 2.

h 0.2 0.1 0.05 0.025

£, % 2.54 0.77 0.28 0.17

В табл. 2 приведена та же погрешность е, что и в табл. 1, для различных шагов по пространству при С4 = 10-2, С3 = С = т = 0, у = 2, г = 6, т = 2 • 10-3 и прочих одинаковых входных данных в момент времени ¿^ = 1. Наибольшее время расчета варианта из таблиц составляет не более одной минуты. Из табл. 2 видно, что численное решение задачи сходится к точному при измельчении шагов разностной сетки.

Представляет интерес влияние числа Куранта Ки = т(С2/С4)0'5/Л, на точность численного решения задачи [13]. В табл. 3 дана погрешность е для опорных данных из табл. 2, различных шагов по пространству и различных С2: 1; 10; 100, С4 = 0, 01 при п = 2 • 10-3 и т2 = 10-3. Как видно из табл. 3, для улучшения точности расчета целесообразно уменьшать шаг по времени т. Наконец, в табл. 4 рассмотрена погрешность е для опорных данных из табл. 1 при т = 10-3 и различных значениях С1, Сз, т.

Таким образом, наличие смешанных производных в (5.1) (варианты под номерами 2, 3 и 4 в табл. 4) не ухудшает точность расчета на конкретном тестовом примере.

Рассмотрим решение задачи, аналогичной задаче нахождения колебания прямоугольной мембраны 0 < < Ь, т =1, 2, 3, с закрепленным краем

и|г = 0, (5.5)

вызванные начальным отклонением

Ч=о = Х1Х2Х3, (5.6) где Хт = вт(пжт/Ь), т = 1, 2, 3, и начальной скоростью, равной нулю

Ж = 0. (5.7)

ОТ *=0

Если взять точное решение краевой задачи (5.1), (5.5)-(5.7) в виде

и = сс^ап^Х^^^ то источник ^ в уравнении (5.1) перепишется при С3 = С1 = т = 0 как

^ = п2 сс8(ап£) Г 3С - с^2 ) Х1Х2Х3.

Таблица 3.

Ки 0.2 0.64 2 0.4 1.28 4 0.8 2.56 8

Н 0.1 0.1 0.1 0.05 0.05 0.05 0.025 0.025 0.025

£1, % 0.77 0.79 5 0.28 0.61 13.7 0.17 3.8 43.1

Ки 0.1 0.32 1 0.2 0.64 2 0.4 1.28 4

Н 0.1 0.1 0.1 0.05 0.05 0.05 0.025 0.025 0.025

£2, % 0.72 0.73 0.83 0.232 0.229 1.25 0.1 0.28 22.7

Таблица 4.

№ С1 С3 ш £, %

1 1 0 0 1.59

2 0 0 0.5 1.57

3 0 0.05 0.5 1.12

4 1 0.05 0.5 1.17

Таблица 5.

Ku 1 3.2 10 2 6.4 20 4 12.8 40

C2 1 10 100 1 10 100 1 10 100

h 0.1 0.1 0.1 0.05 0.05 0.05 0.025 0.025 0.025

£, % 0.81 0.82 2.65 0.21 0.21 2.01 0.052 0.056 1.85

В табл. 5 дана погрешность е для С4 = 10-4, а = 1, Ь = 1 и различных С2 при т = 10-3 в момент времени = 1. Как видно, при Ки < 40 (к = 0.025, С2 = 100) существует устойчивое численное решение модельной краевой задачи (5.1), (5.5)-(5.7).

Заключение

На основании вышеизложенного можно сделать следующие выводы.

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

( з

ных сеток внутри параллелепипеда 0 1т2 + ^ кт

\ т=1

2. При постоянных коэффициентах переноса разностные уравнения безусловно устойчивы.

3. На тестовых примерах показана практическая сходимость разностной схемы ИИМ.

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

[1] Гришин А.М. Об одном видоизменении метода М.Е. Швеца // Инженерно-физический журн. 1970. Т. 19, № 1. С. 84-93.

[2] Гришин А.М., Берцун В.Н. Итерационно-интерполяционный метод и теория сплайнов // Докл. АН СССР. 1974. Т. 214, № 4. С. 701-704.

[3] Grishin A.M., Beroün V.N. The iterational-interpolation method and spline theory // Sovet. Math. Dokl. 1974. Vol. 15, N 1. P. 222-227.

[4] Гришин А.М. Об одном итерационно-интерполяционном методе // Тр. НИИ прикл. мат. и механики при ТГУ. Томск: Изд-во ТГУ, 1973. С. 45-48.

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

[6] Якимов А.С. Об одном методе расщепления // Числ. методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1985. Т. 16, № 2. С. 144-161.

[7] ИтЕрАционно-интЕрполяционный метод и его приложения / А.М. Гришин, В.И. Зин-ченко, К.Н. Ефимов, А.Н. Субботин, А.С. Якимов. Томск: Изд-во ТГУ, 2004.

[8] Марчук Г.И. Методы расщепления. М.: Наука, 1988.

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

[10] Коновалов А.Н. Метод дробных шагов решения задачи Коши для многомерного уравнения колебаний // Докл. АН СССР. 1962. Т. 147, № 1. С. 25-27.

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

[12] Дьяконов Е.Г. О применении разностных схем с расщепляющимся оператором для гиперболических уравнений с переменными коэффициентами // Докл. АН СССР. 1963. Т. 151, № 4. С. 762-765.

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

[14] ДЕмидович Б.П., Марон И.А. Основы вычислительной математики. М.: Наука, 1966.

Поступила в редакцию 18 апреля 2005 г., в переработанном виде — 22 ноября 2006 г.