Асимптотическая эквивалентность двухслойных разностных схем

Краевые задачи для дифференциальных уравнений типа

Lu + f(x) = 0

можно трактовать как операторные уравнения:

Au = f.

Это справедливо в частности для уравнений параболического типа:

du

— = Lu + f(x, t), исследованию асимптотического поведения разностных схем для которых и посвящена эта работа.

Рассмотрим искомую функцию u = u(x, t) как некоторую абстрактную функцию u(t) со значениями в некотором банаховом пространстве. Если ввести для t равномерную сетку с шагом т, то можно получить выражение для r-слойной операторно-разностной схемы:

г-1

B₀ ⁽ t_n ⁾ u(t_n +i) = ^ C_s ⁽ t_n ⁾ u ⁽ t_n +i- _s' ⁾ + f(t_n),    n>r-1,

S = 1

где B ₀(t'),C_s(t') - линейные операторы, действующие из некоторого линейного нормированного пространства К в К, которое в свою очередь зависит от некоторого параметра h. Заметим, что эти линейные операторы также зависят от г и h.

Тогда для двухслойной схемы имеем:

Boun+! + B1un = rfn,n>0.

Выполняя замену, предложенную в [1], получим каноническую форму двухслойных схем:

B(tn) Un+\ Un + A(tn)Un =fn,n>

Bu_t + Au = f(t).

Рассмотрим первую краевую задачу для одномерного однородного уравнения теплопроводности:

Ut   UXX

u(x, 0) = ф(х),    0 < x < L,

u(0,t) = ^(t),   0 < t

u(0, t) = ^2(t),    0 < t

Пусть г - шаг сетки по t, h - шаг сетки по х. Тогда для этой задачи (1) примет вид

But + Au = 0, B = Е + тА, А = -Л,(3)

где Л - оператор вторых конечных разностей:

(-2   1    0    ...0\

1   -2   1        0\ 1

0   1-2  ...  0т^.

h2

0.01-2/

Разрешая (3) через un+1 - значение на верхнем слое, - получим un+1 = Run,    R = Е- tB-1A, где оператор R называется оператором перехода, а B-1 - обратный к B оператор.

Рассмотрим два разностные схемы:

^un + 1 = ^R1^un,    ^un + 1 = ^R2^un,    ^R1 = ^^ ^ ^Rh^ = ^R2.

Согласно [2], условие асимптотической эквивалентности разностных схем при условии согласования норм можно записать как

\\Rh}" -^"Ц^'.

|R^^,l'f₁-R^^!>"|,, где P1, P₂ - некоторые операторы, воздействующие на начальные условия, соответствующие первой и второй разностной схеме, I - порядок аппроксимации.

Не теряя общности, можно положить

P1=P₂= E.

Тогда получим иную форму условия асимптотической эквивалентности:

\\R?-R?\\l.                                 (4)

Получим некоторые оценки для нормы выше, при которых (4) справедливо.

Для явной схемы В = E, поэтому оператор перехода имеет вид

R = E -тА.

Рассмотрим два различных оператора А1 и А₂ такие, что:

А₂=А1 + Q,                                 (5)

где Q - оператор возмущения. В этой работе будем полагать, что возмущение постоянно и одинаково для каждого узла сетки:

	(£ £	£ £	0 £	...    °\ ...     0
Q =	0	£	£		0
	\0		0	£	г)

Тогда имеем два оператора перехода R1 и R₂: R₁=E-tA₁,   R₂=E-tA₂.

Получим оценку для возмущения £, при котором схемы, порождаемые операторами R1 и R₂, являются асимптотически эквивалентными.

Подставляя операторы в левую часть условия (4) и используя формулу разности степеней:

I\R? - R^\ = \(R - R2XR?-1 + R^2R2 + - + RiR^2 + Rr^\.(7)

Здесь и далее будем рассматривать m-норму:

||A|| = ||A||i = max^ la^.

Обозначим множители в правой части (7) как D1 и D₂, то есть:

\R? - R?W = WD1D2W,(8)

D1 = R1- R₂,    D₂ = Ri^- + R1~²R₂ + - + R1R2^- + R?^-1.

По свойству нормы:

WD1D2W < WDiW • WD2W.(9)

Сначала оценим вторую норму, потому что далее эта оценка будет справедлива для всех рассматриваемых случаев. Имеем:

WD2W = WR?-¹ + RV²R2 + - + RiR^² + R1^-1W.

По теореме о достаточных условиях устойчивости двухслойных схем в линейных нормированных пространствах [1] можно положить:

HR?^-1H 1,    HR?^-1H < М2,

Н R”R?^-1-”| | < max(M1,M₂,M1M₂) = М,   0<т<п-1.

Получим

НО2Н = ||R?-1 + R?-2R2 + -" + R1R?-2 + R?-1W < пМ.

Теперь оценим первую норму

^|D1\ = ^Н^ ^{— т4}1 -^{Е + т4}2^Н = ^042 ^— ^1^)Н< 1^Т1 • ^Н^2 ^— ^1\

Учитывая, что т > 0:

НМ < |т| • Н42 -41Н = тН41 + Q -41Н = тН^Н.

Легко вычислить норму Q:

HQH = maxS7- Iqijl = тах(2г, 3г) = 3г.

Подставляя (12) в (11):

НМК THQH = 3тг.

Объединяя (8), (9), (10) и (13), получим:

HR? -R?H < 3тгпМ.

Выполняя замену т = rh², п = h^-1, имеем:

HR? -R?H < 3Mrгh.

Используя правую часть условия (4):

3Mrгh < Ch1, или

г < Ch^1-1, C = —.

,       3Mr

Теперь получим аналогичную оценку для неявной схемы, когда возмущениям подвержен только оператор 4.

Рассмотрим два различных оператора 41 и 4₂ аналогично (5). Получим два оператора перехода R1 и R₂:

R1 = E-tB-141,   R2=E-tB-142.

Как было сказано ранее, оценка для D₂ из (8) остаётся прежней, поэтому нас интересует оценка D1:

HD1H = НЕ — тВ^-141 - Е + тВ^-142Н =

= НтВ^-1(42 -41)Н < тНВ^-1Н • HQH = 3тгНВ^-1Н.                (15)

Для оценки ||В-1Н воспользуемся тем, что В - монотонная матрица. Действительно, используя (3):

10 1 0    1 В = E - гА = 1 0  0	0   .    0        2   -1   0    .     0 0  ...  0|  1-1   2   -1   ...   0) T 1  ...  0-0  -1  2   ...  0 77 = 1                                         h2
\0	0  0  1/   \0   ...  0-1  2/

1 - 2r     г        0      ...      0

I    r = 1    0	1 - 2r      r      ...      0    I r     1 - 2r  .     0   I.
\0	...        0      r   1 - 2r/

В работе [3] представлена теорема, согласно которой для монотонной матрицы G

1 HG-1Hm<-—, R* (G ) где R»(G) - минимальная величина диагонального преобладания матрицы G. Заметим, что здесь используется -норма, но для симметричной матрицы они взаимозаменяемы:

Ик

= max^ laiy I =

max^ |ау| = H^H„.

i

Вычислим R*(By.

R»(B) = 1 + 2г-г-r = 1.

Таким образом

Нв-1Н <    = 1=1.

R» (В)    1

Подставляя (16) в (15):

НМ < 3^llB-1|| < 3гг.

Так как оценка (17) совпадает с оценкой (13), а оценка (10) одинакова для обоих

случаев, то делаем вывод, что для неявной схемы в случае А1 ^ А2 оценка возмущения

одинакова и выражается (14).

Наконец рассмотри случай, когда для неявной схемы А1 = А2, B1 ^ B2. Имеем:

R1 = Е - тВ-1 А,    R2= Е - tB-1A,

B2 = Bi + Q,

причём Q по-прежнему имеет вид (6).

Рассмотрим оценку D1:

НМ1 = HE-iB—1A-E + iB—1AH <

< HtA(b2-1 - B1-1)H < тНАН • Нв-1 - B1-1H.

Легко вычислить норму А:

НАН = max i

^ RijI = max(2^ 2

+ h-2, 2h-2 + h-2 + h-2) = 4h-2.

Подставляя её в (18):

IIM < т • 4h-2HB-1 — В--1Н.

Для оценки IB-¹ — В—¹1| используем теорему о матрице, обратной сумме матриц [4], в частности её формой как тождества Хуа [5]:

(А + В)-1 = А-1 — (А + АВ-1 А)-1.

Оценим

IB-1 — В-1 II = \(В1 + Q)-1 — В-1 II = ||В1-1 — (В1 + B1Q-1B1)-1 — В1-1|| =

= |(В1 + B1Q-1B1)-1| = |(В1(Е + е-1В1))-1|.

По свойству обратной матрицы:

IKB^F + Мв^Г1!! = IK£ + Q"1B1)"1B"1H < Н(я + Q-1B1)-1H • НВ-1!

Оценка нормы невозмущенной матрицы В представлена выражением (16). Подставим её вместо ||В-1Н:

Н(Е + Q-1B1)-1H • НВ1-1Н < Н(я + Q-1B1)-1H = |(Q-1(B1 + Q))-1| =

= H(B1 + Q)-1QH

Используя оценку нормы оператора Q, представленной (12), имеем

Н(В1 + Q)-1H • HQH <3

Для оценки Н(В1 + Q)-1H снова воспользуемся теоремой из [3]. Чтобы однозначно вычислить RS(B1 + Q), введём условие г > £.                                        (20)

Тогда

R*(B1 + Q) = 1 + 2г + £ — (г — £) — (г — £) =

= 1 + 2г + £ — 2г + 2£ = 1 + 3£.

Таким образом

11 "№ + Q)H

Подставляя в (19), получаем

Н(В1 + Q)-1H • HQH < 3£Н(В1 + Q)-1H < 3£^< 3£.

Объединяя рассуждения выше:

UMI < т • 4h-2HB-1 — В-1! < т • 4h-2 • 3£ = 12T£h-2.

Выполняя замену т = гh2 и объединяя с оценкой нормы D2:

||R” — r*|| = HD1D2H < 12Mг£h-1 < Ch1, или

£ < Chl+1, C = —.                          (21)

12МГ                                     v 7

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

ф(х) = sin?rx,    ^1(t) = 0,    ^2(t) = 0, L = 1,

h = 0.01, r = 0.4.

Эта задача имеет точное решение:

u(x,t) = e-^^ sin ttx.

Для решения с помощью неявной схемы будем использовать обычный метод прогонки. Явная схема имеет первый порядок аппроксимации (I = 1), поэтому (14) для явной схемы имеет вид

£

Неявная схема имеет второй порядок аппроксимации (I = 2), поэтому (14) для неявной схемы (41 ^ Л₂) имеет вид

£ < Ch, а (21) для неявной схемы (В1 Ф В2) имеет вид

£ < Ch3.

Таблица 1

Результаты численного эксперимента

Схема	£	Н^-^Н	Итог
явная	0.000000	0.000004	эквивалентны
явная	0.000001	0.000004	эквивалентны
явная	0.000100	0.000006	эквивалентны
явная	0.010000	0.000120	эквивалентны
явная	0.100000	0.001157	эквивалентны
явная	1.000000	0.011472	эквивалентны
неявная (41 Ф 42)	0.000000	0.000011	эквивалентны
неявная (41 Ф 42)	0.000001	0.000011	эквивалентны
неявная (41 Ф 42)	0.000100	0.000010	эквивалентны
неявная (41 Ф 42)	0.010000	0.000028	эквивалентны
неявная (41 Ф 42)	0.100000	0.000374	эквивалентны
неявная (41 Ф 42)	1.000000	0.003827	неэквивалентны
неявная (В1 Ф В2)	0.000000	0.000011	эквивалентны
неявная (В1 Ф В2)	0.000001	0.000181	эквивалентны
неявная (В1 Ф В2)	0.000100	0.019000	неэквивалентны
неявная (В1 Ф В2)	0.010000	0.825713	неэквивалентны
неявная (В1 Ф В2)	0.100000	0.983464	неэквивалентны

В рамках эксперимента будем рассматривать возмущения

г = [0,h3,h2,h,h°'5,h°] = [0,0.000001,0.0001,0.01,0.1,1.0], но заметим, что в виду (20) для неявной схемы (В1 ^ В2) мы не можем полагать г = 1.0.

Результаты эксперимента представлены в таблице 1. Как можно видеть, все полученные оценки корректны.