Asymptotic equivalence of two-stage difference schemes

Abstract

Value estimations for fluctuations of difference schemes’ matrix operators, which provide asymptotic equivalence of these schemes, are obtained. Numeric experiments for the confirmation of theoretic results are described.

Full Text

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

$Lu+f\left(x\right)=0$

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

$Au=f.$

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

$\frac{\partial u}{\partial t}=Lu+f\left(x,t\right),$

исследованию асимптотического поведения разностных схем для которых и посвящена эта работа.

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

$B_0\left(t_n\right)u\left(t_{n+1}\right)=\sum _{s=1}^{r-1}{C}_s\left(t_n\right)u\left(t_{n+1-s}\right)+f\left(t_n\right),n\ge r-\mathrm{1,}$

где $B_0\left(t\right),C_s\left(t\right)$ – линейные операторы, действующие из некоторого линейного нормированного пространства ℋ в ℋ, которое в свою очередь зависит от некоторого параметра ℎ. Заметим, что эти линейные операторы также зависят от τ и ℎ.

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

$B_0{u}_{n+1}+B_1{u}_n=\tau {\tilde{f}}_n,n\ge 0.$

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

$B\left(t_n\right)\frac{u_{n+1}-u_n}{\tau }+A\left(t_n\right)u_n={\tilde{f}}_n,n\ge 0;$

$B{u}_t+Au=\tilde{f}\left(t\right).\text{ (1)}$

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

${\text{u}}_t-u_{xx}=\mathrm{0,}$

$u\left(x\mathrm{,0}\right)=\varphi \left(x\right),0\le x\le L,$

$u\left(\mathrm{0,}t\right)={\psi }_1\left(t\right),0\le t\le T,\text{ (2)}$

$u\left(\mathrm{0,}t\right)={\psi }_2\left(t\right),0\le t\le T.$

Пусть 𝜏 – шаг сетки по 𝑡, ℎ – шаг сетки по 𝑥. Тогда для этой задачи (1) примет вид

$B{u}_t+Au=\mathrm{0,}\text{ B=E+}\tau A,A=-\Lambda ,\text{(3)}$

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

$\Lambda \text{=}\left(\begin{array}{ccccc}-2& 1& 0& \dots & 0\\ 1& -2& 1& \dots & 0\\ 0& 1& -2& \dots & 0\\ & & \dots & & \\ 0& \dots & 0& 1& -2\end{array}\right)\frac{1}{h^2}.$

Разрешая (3) через $u_{n+1}$ – значение на верхнем слое, – получим

$u_{n+1}=R{u}_n,R=E-\tau B^{-1}A,$

где оператор 𝑅 называется оператором перехода, а 𝐵⁻¹ – обратный к 𝐵 оператор.

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

$u_{n+1}=R_1{u}_n,u_{n+1}=R_2{u}_n,R_1=R_h^{\left(1\right)}\ne R_h^{\left(2\right)}=R_2.$

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

$\begin{array}{l}||R_h^{{\left(1\right)}^n}-R_h^{{\left(2\right)}^n}{P}_2||\le C{h}^l,\\ ||R_h^{{\left(1\right)}^n}{P}_1-R_h^{{\left(2\right)}^n}||\le C{h}^l,\end{array}$

где 𝑃₁,𝑃₂ – некоторые операторы, воздействующие на начальные условия, соответствующие первой и второй разностной схеме, 𝑙 – порядок аппроксимации.

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

$P_1\equiv P_2\equiv E.$

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

$||R_1^n-R_2^n||\le C{h}^l.\text{ (4)}$

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

Для явной схемы 𝐵≡𝐸, поэтому оператор перехода имеет вид

$R=E-\tau A.$

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

$A_2=A_1+\text{Q, (5)}$

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

$\text{Q=}\left(\begin{array}{ccccc}\epsilon & \epsilon & 0& \dots & 0\\ \epsilon & \epsilon & \epsilon & \dots & 0\\ 0& \epsilon & \epsilon & \dots & 0\\ & & \dots & & \\ 0& \dots & 0& \epsilon & \epsilon \end{array}\right).\text{ (6)}$

Тогда имеем два оператора перехода 𝑅₁ и 𝑅₂:

$R_1=E-\tau A_1,R_2=E-\tau A_2.$

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

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

$||R_1^n-R_2^n||=||\left(R_1-R_2\right)\left(R_1^{n-1}+R_1^{n-2}{R}_2+\dots +R_1{R}_2^{n-2}+R_2^{n-1}\right)||.\text{ (7)}$

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

$||A||={||A||}_1=\underset{i}{\max }\sum _j\left|a_{ij}\right|.$

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

$||R_1^n-R_2^n||=||D_1{D}_2||\text{, (8)}$

$D_1=R_1-R_2,D_2=R_1^{n-1}+R_1^{n-2}{R}_2+\dots +R_1{R}_2^{n-2}+R_2^{n-1}.$

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

$||D_1{D}_2||\le ||D_1||\cdot ||D_2||.\text{ (9)}$

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

$||D_2||=||R_1^{n-1}+R_1^{n-2}{R}_2+\dots +R_1{R}_2^{n-2}+R_2^{n-1}||.$

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

$||R_1^{n-1}||\le M_1,||R_2^{n-1}||\le M_2,$

$||R_1^m{R}_2^{n-1-m}||\le \max \left(M_1,M_2,M_3,M_4\right)=M,\text{ 0}\le \text{m}\le \text{n-1.}$

Получим

$||D_2||=||R_1^{n-1}+R_1^{n-2}{R}_2+\dots +R_1{R}_2^{n-2}+R_2^{n-1}||\le nM.\text{ (10)}$

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

$||D_1||=||E-\tau A_1-E+\tau A_2||=||\tau \left(A_2-A_1\right)||\le \left|\tau \right|\cdot ||A_2-A_1||$

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

$||D_1||\le \left|\tau \right|\cdot ||A_2-A_1||=\tau ||A_1+\text{Q}-A_1||=\tau ||\text{Q}||.\text{ (11)}$

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

$||\text{Q}||=\underset{i}{\max }{\sum }_j\left|q_{ij}\right|\underset{i}{\max }\left(2\epsilon \mathrm{,3}\epsilon \right)=3\epsilon .\text{ (12)}$

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

$||D_1||\le \tau ||\text{Q}||=3\tau \epsilon .\text{ (13)}$

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

$||R_1^n-R_2^n||\le 3\tau \epsilon nM.$

Выполняя замену 𝜏=𝑟ℎ²,𝑛=ℎ⁻¹, имеем:

$||R_1^n-R_2^n||\le 3Mr\epsilon h.$

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

$3Mr\epsilon h\le C{h}^l,$

Или

$\epsilon \le \tilde{C}{h}^{l-1},\tilde{C}=\frac{c}{3Mr}.\text{ (14)}$

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

Рассмотрим два различных оператора 𝐴₁ и 𝐴₂ аналогично (5). Получим два оператора перехода 𝑅₁ и 𝑅₂:

$R_1=E-\tau B^{-1}{A}_1,R_2=E-\tau B^{-1}{A}_2.$

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

$\begin{array}{l}||D_1||=||E-\tau B^{-1}{A}_1-E+\tau B^{-1}{A}_2||=\\ =||\tau B^{-1}\left(A_2-A_1\right)||\le \tau ||B^{-1}||\cdot ||\text{Q}||=3\tau \epsilon ||B^{-1}||.\text{ (15)}\end{array}$

Для оценки ‖𝐵⁻¹‖ воспользуемся тем, что 𝐵 – монотонная матрица. Действительно, используя (3):

$\begin{array}{l}B=E-\tau A=\left(\begin{array}{ccccc}1& 0& 0& \dots & 0\\ 0& 1& 0& \dots & 0\\ 0& 0& 1& \dots & 0\\ & & \dots & & \\ 0& \dots & 0& 0& 1\end{array}\right)-\left(\begin{array}{ccccc}2& -1& 0& \dots & 0\\ -1& 2& -1& \dots & 0\\ 0& -1& 2& \dots & 0\\ & & \dots & & \\ 0& \dots & 0& -1& 2\end{array}\right)\frac{\tau }{h^2}=\\ =\left(\begin{array}{ccccc}1-2r& r& 0& \dots & 0\\ r& 1-2r& r& \dots & 0\\ 0& r& 1-2r& \dots & 0\\ & & \dots & & \\ 0& \dots & 0& r& 1-2r\end{array}\right).\end{array}$

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

${||G^{-1}||}_{\infty }\le \frac{1}{R_{\ast }\left(G\right)},$

где 𝑅_∗(𝐺) – минимальная величина диагонального преобладания матрицы 𝐺. Заметим, что здесь используется -норма, но для симметричной матрицы они взаимозаменяемы:

${||A||}_1=\underset{i}{\max }\sum _j\left|a_{ij}\right|\underset{j}{\max }\sum _i\left|a_{ij}\right|={||A||}_{\infty }.$

Вычислим 𝑅_∗(𝐵):

𝑅_∗(𝐵)=1+2𝑟−𝑟−𝑟=1.

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

$||B^{-1}||\le \frac{1}{R_{\ast }\left(B\right)}=\frac{1}{1}=1.\text{ (16)}$

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

$||D_1||\le 3\tau \epsilon ||B^{-1}||\le 3\tau \epsilon .\text{ (17)}$

Так как оценка (17) совпадает с оценкой (13), а оценка (10) одинакова для обоих случаев, то делаем вывод, что для неявной схемы в случае 𝐴₁≠𝐴₂ оценка возмущения одинакова и выражается (14).

Наконец рассмотри случай, когда для неявной схемы 𝐴₁≡𝐴₂,𝐵₁≠𝐵₂. Имеем:

$R_1=E-\tau B_1^{-1}A,R_2=E-\tau B_2^{-1}A,$

$B_2=B_1+\text{Q},$

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

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

$\begin{array}{l}||D_1||=||E-\tau B_1^{-1}A-E+\tau B_2^{-1}A||\le \\ \le ||\tau A\left(B_2^{-1}-B_1^{-1}\right)||\le \tau ||A||\cdot ||B_2^{-1}-B_1^{-1}||.\text{ (18)}\end{array}$

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

$||A||=\underset{i}{\max }\sum _j\left|q_{ij}\right|=\underset{j}{\max }(2h^{-2}+h^{-2}\mathrm{,2}{h}^{-2}+h^{-2}+h^{-2})=4h^{-2}.$

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

$||D_1||\le \tau \cdot 4h^{-2}||B_2^{-1}-B_1^{-1}||.$

Для оценки $||B_2^{-1}-B_1^{-1}||$ используем теорему о матрице, обратной сумме матриц [4], в частности её формой как тождества Хуа [5]:

${(A+B)}^{-1}=A^{-1}-{(A+A{B}^{-1}A)}^{-1}.$

Оценим

$\begin{array}{l}||B_2^{-1}-B_1^{-1}||=||{(B_1+\text{Q)}}^{-1}-B_1^{-1}||=||B_1^{-1}-{(B_1+B_1{\text{Q}}^{-1}{\text{B}}_1\text{)}}^{-1}-B_1^{-1}||=\\ =||{(B_1+B_1{\text{Q}}^{-1}{\text{B}}_1\text{)}}^{-1}||=||{\left(B_1\right(E+{\text{Q}}^{-1}{B}_1\left)\right)}^{-1}||.\end{array}$

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

$||{\left(B_1\right(E+{\text{Q}}^{-1}{B}_1\left)\right)}^{-1}||=||{(E+{\text{Q}}^{-1}{B}_1)}^{-1}{B}_1^{-1}||\le ||{(E+{\text{Q}}^{-1}{B}_1)}^{-1}||||B_1^{-1}||$

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

$\begin{array}{l}||{(E+{\text{Q}}^{-1}{B}_1)}^{-1}||\cdot ||B_1^{-1}||\le ||{(E+{\text{Q}}^{-1}{B}_1)}^{-1}||=||{\left({\text{Q}}^{-1}\right(B_1+\text{Q}\left)\right)}^{-1}||=\\ =||{(B_1+\text{Q)}}^{-1}\text{Q}||\le ||{(B_1+\text{Q)}}^{-1}||\cdot ||\text{Q}||.\end{array}$

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

$||{(B_1+\text{Q)}}^{-1}||\cdot ||\text{Q}||\le 3\epsilon ||{(B_1+\text{Q)}}^{-1}||.\text{ (19)}$

Для оценки $||{(B_1+\text{Q)}}^{-1}||$ снова воспользуемся теоремой из [3]. Чтобы однозначно вычислить 𝑅_∗(𝐵₁+𝑄), введём условие

$r>\epsilon \text{ (20)}$.

Тогда

$\begin{array}{l}{R}_{\ast }(B_1+\text{Q)}=1+2r+\epsilon -(r-\epsilon )-(r-\epsilon )=\\ =1+2r+\epsilon -2r+2\epsilon =1+3\epsilon .\end{array}$

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

$||{(B_1+\text{Q)}}^{-1}||\cdot ||\text{Q}||\le \frac{1}{R_{\ast }(B_1+\text{Q)}}=\frac{1}{1+3\epsilon }.$

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

$||{(B_1+\text{Q)}}^{-1}||\cdot ||\text{Q}||\le 3\epsilon ||{(B_1+\text{Q)}}^{-1}||\le 3\epsilon \frac{1}{1+3\epsilon }\le 3\epsilon .$

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

$||D_1||\le \tau \cdot 4h^{-2}||B_2^{-1}-B_1^{-1}||\le \tau \cdot 4h^{-2}\cdot 3\epsilon =12\tau \epsilon h^{-2}.$

Выполняя замену 𝜏=𝑟ℎ² и объединяя с оценкой нормы 𝐷₂:

$||R_1^n-R_2^n||=||D_1{D}_2||\le 12Mr\epsilon h^{-1}\le C{h}^l,$

или

$\epsilon \le \tilde{C}{h}^{l+1},\tilde{C}=\frac{c}{12Mr}.\text{ (21)}$

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

$\varphi \text{(x)=sin}\pi x,{\psi }_1\left(t\right)\equiv \mathrm{0,}{\psi }_2\left(t\right)\equiv \mathrm{0,}\text{ L=1,}$

$h=\mathrm{0.01,}\text{ r=0.4.}$

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

$u\text{(}x,t{\text{)=e}}^{-{\pi }^{2}t}\sin \pi x.$

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

$\epsilon \le \tilde{C}{h}^0=\tilde{C}.$

Неявная схема имеет второй порядок аппроксимации (𝑙=2), поэтому (14) для неявной схемы (𝐴₁≠𝐴₂) имеет вид

$\epsilon \le \tilde{C}h,$

а (21) для неявной схемы (𝐵₁≠𝐵₂) имеет вид

$\epsilon \le \tilde{C}{h}^3.$

В рамках эксперимента будем рассматривать возмущения = [0, ℎ³, ℎ², ℎ, ℎ^0.5, ℎ⁰] = [0, 0.000001, 0.0001, 0.01, 0.1, 1.0], но заметим, что в виду (20) для неявной схемы (B₁ ≠ B₂) мы не можем полагать ε = 1.0. Результаты эксперимента представлены в таблице 1. Как можно видеть, все полученные оценки корректны.

Таблица 1

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

Схема	ε	$||R_1^n-R_2^n||$	Итог
явная	0.000000	0.000004	эквивалентны
явная	0.000001	0.000004	эквивалентны
явная	0.000100	0.000006	эквивалентны
явная	0.010000	0.000120	эквивалентны
явная	0.100000	0.001157	эквивалентны
явная	1.000000	0.011472	эквивалентны
неявная (А1 ≠ А2)	0.000000	0.000011	эквивалентны
неявная (А1 ≠ А2)	0.000001	0.000011	эквивалентны
неявная (А1  ≠ А2)	0.000100	0.000010	эквивалентны
неявная (А1 ≠ А2)	0.010000	0.000028	эквивалентны
неявная (А1 ≠ А2)	0.100000	0.000374	эквивалентны
неявная (А1 ≠ А2)	1.000000	0.003827	неэквивалентны
неявная (B1 ≠ В2)	0.000000	0.000011	эквивалентны
неявная (B1 ≠ В2)	0.000001	0.000181	эквивалентны
неявная (B1 ≠ В2)	0.000100	0.019000	неэквивалентны
неявная (B1 ≠ В2)	0.010000	0.825713	неэквивалентны
неявная (B1 ≠ В2)	0.100000	0.983464	неэквивалентны