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

Аннотация

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

Полный текст

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

$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	неэквивалентны