Zagadnienia

9. Metoda różnic skończonych - stabilność schematów dla zadań eliptycznych w normie maksimum

W tym rozdziale zajmiemy się przedstawieniem metod badania stabilności schematów różnicowych dla zadań liniowych w dyskretnej normie maksimum.

Będziemy badali stabilność schematu zapisanego w formie (8.5)-(8.6). Dla x\in\Omega _{h} możemy zapisać (8.5) jako:

L_{h}u(x)\equiv=\sum _{{y\in N_{h}(x)}}A(x,y)u_{h}(y)=f_{h}(x), (9.1)

gdzie N_{h}(x) jest podzbiorem \overline{\Omega}_{h} punktów, dla których A(x,y)\not=0, czyli uwzględnionych w równaniu dla tego x.

Jeśli x\in\Gamma _{{k,h}}, to dla (8.6) zachodzi:

l_{{k,h}}u(x)\equiv\sum _{{y\in N_{h}(x)}}A(x,y)u_{h}(y)=g_{{k,h}}(x),

gdzie N_{h}(x) jest zdefiniowane analogicznie jak poprzednio. N_{h}(x) jest zdefiniowane jednoznacznie.

N_{h}(x) nazywamy otoczeniem siatkowym punktu x. Wprowadzimy również otoczenie siatkowe nakłute: N_{h}^{{\prime}}(x)=N_{h}(x)\setminus\{ x\}. Oczywiście N_{h}(x) może być jednopunktowe, wtedy N_{h}^{{\prime}}(x) jest zbiorem pustym.

Zapiszmy schemat (8.5)-(8.6) jako:

\mathcal{L}_{h}u_{h}(x)\equiv\sum _{{y\in N_{h}(x)}}A(x,y)u_{h}(y)=\psi _{h}(x)\qquad x\in\overline{\Omega}_{h}, (9.2)

gdzie

\psi _{h}(x)=\left\{\begin{array}[]{lcl}f_{h}(x)&&x\in\Omega _{h}\\
g_{{k,h}}(x)&&x\in\Gamma _{{k,h}}\quad s=1,\ldots,s.\end{array}\right.

Wtedy zachodzi następujące twierdzenie, pozwalające na wykazanie stabilności niektórych schematów w normie dyskretnej maksimum:

Twierdzenie 9.1

Niech \mathcal{L}_{h} dla (8.5)-(8.6) będzie w formie (9.2). Załóżmy, że dla pewnej stałej \alpha>0 i dla h\leq h_{0}:

|A(x,x)|-\sum _{{y\in N_{h}^{{\prime}}(x)}}|A(x,y)|\geq\alpha\qquad\forall x\in\overline{\Omega}_{h}.

Wtedy

\| u_{h}\| _{{\infty,h}}\leq\frac{1}{\alpha}\max(\| f_{h}\| _{{\infty,h,\Omega _{h}}}+\sum _{{k=1}}^{s}\| g_{{k,h}}\| _{{\infty,h,\Gamma _{{k,h}}}}).

Widzimy, że \| u_{h}\| _{{\infty,h}}=\max _{{x\in\overline{\Omega}_{h}}}|u_{h}(x)|=|u_{h}(x_{0})| dla pewnego x_{0}\in\overline{\Omega}_{h}.

Rozpatrzmy równanie ze schematu dla tego punktu:

\displaystyle|\psi _{h}(x_{0})| \displaystyle= \displaystyle|\sum _{{y\in N_{h}(x_{0})}}A(x_{0},y)u_{h}(y)|
\displaystyle\geq \displaystyle|A(x_{0},x_{0})||u_{h}(x_{0})|-\sum _{{y\in N_{h}^{{\prime}}(x_{0})}}|A(x_{0},y)||u_{h}(y)|
\displaystyle\geq \displaystyle|A(x_{0},x_{0})||u_{h}(x_{0})|-\sum _{{y\in N_{h}^{{\prime}}(x_{0})}}|A(x_{0},y)||u_{h}(x_{0})|
\displaystyle= \displaystyle\left(|A(x_{0},x_{0})|-\sum _{{y\in N_{h}^{{\prime}}(x_{0})}}|A(x_{0},y)|\right)|u_{h}(x_{0})|\geq\alpha|u_{h}(x_{0})|,

czyli

\| u_{h}\| _{{\infty,h}}=|u_{h}(x_{0})|\leq\frac{1}{\alpha}\max _{{x\in\overline{\Omega}_{h}}}|\psi _{h}(x)|.

Jak widzimy, jest to proste kryterium. Sprawdźmy je na naszym modelowym zadaniu (7.5)-(7.6):

Przykład 9.1

Dla zadania (7.5)-(7.6) otrzymujemy następujący układ:

\displaystyle-\frac{1}{h^{2}}u(x_{{k-1}})+(\frac{2}{h^{2}}+c)u(x_{k})-\frac{1}{h^{2}}u(x_{{k+1}}) \displaystyle= \displaystyle f(x_{k})\qquad k=1,\ldots,N-1
\displaystyle u(x_{k}) \displaystyle= \displaystyle g(x_{k})\qquad k=0,N.

dla x_{k}=a+k*h.

Zatem N_{h}(x_{k})=\{ x_{{k-1}},x_{k},x_{{k+1}}\} dla k=1,\ldots,N-1 i N_{h}(x)=\{ x_{k}\} dla x_{k}\in\{ a,b\}, tzn. dla k=0,N. Sprawdzamy założenie twierdzenia:

|A(x_{k},x_{k})|-\sum _{{y\in N_{h}(x_{k})\setminus\{ x_{k}\}}}|A(x_{k},y)|=\left\{\begin{array}[]{lcr}c&&k=1,\ldots,N-1\\
1&&k=0,N.\end{array}\right.

Zatem \alpha=\min\{ 1,c\} i z naszego kryterium, tzn. z twierdzenia 9.1, otrzymujemy stabilność zadania przybliżonego w normie dyskretnej supremum tylko w przypadku c>0 ze stałą \frac{1}{\alpha}=\max\{ 1,\frac{1}{c}\}.

Powyższe oszacowanie sugeruje, że jeśli c=0, to schemat może nie być stabilny w dyskretnej normie maksimum. Okaże się, że istnieją jednak inne kryteria badania stabilności, które są bardziej precyzyjne. Przedstawimy je poniżej.

9.1. Różnicowa zasada maksimum

Jak wiadomo, por. np. rozdział 6.4 w [11], dla równania eliptycznego spełnionych jest szereg zasad maksimum. Okaże się, że odpowiednio skonstruowane schematy różnicowe, czyli problemy przybliżone (różnicowe), spełniają analogiczne różnicowe zasady maksimum. Korzystając z tych zasad będziemy mogli wykazać stabilność tychże schematów.

Załóżmy, że operator \mathcal{L}_{h} określony na \overline{\Omega}_{h} jest w formie (9.2).

Definicja 9.1

Operator \mathcal{L}_{h} w postaci (9.2) będziemy nazywać operatorem dodatniego typu (ang. positive operator) w \overline{\Omega}_{h}, jeśli dla dowolnego x\in\overline{\Omega}_{h}

  1. A(x,x)>0,

  2. A(x,y)<0\qquad\forall y\in N_{h}^{{\prime}}(x) ,

  3. \sum _{{y\in N_{h}(x)}}A(x,y)\geq 0.

Dodatkowo dla operatora typu dodatniego przedstawiamy siatkę \Omega _{h} jako dwa rozłączne zbiory \overline{\Omega}_{h}=\sum _{{k=1,2}}\Omega _{h}^{{(1)}}\cup\Omega _{h}^{{(2)}} zdefiniowane jako:

\Omega _{h}^{{(1)}}=\{ x\in\overline{\Omega}_{h}:\sum _{{y\in N_{h}(x)}}A(x,y)=0\}

i

\Omega _{h}^{{(2)}}=\overline{\Omega}_{h}\setminus\Omega _{h}^{{(1)}}=\{ x\in\overline{\Omega}_{h}:\sum _{{y\in N_{h}(x)}}A(x,y)>0\}.

Wprowadzamy jeszcze jedną definicję:

Definicja 9.2

Załóżmy, że \mathcal{L}_{h} w postaci (9.2) jest operatorem dodatniego typu w \overline{\Omega}_{h}, dla którego zachodzi warunek: \Omega _{h}^{{(2)}}\not=\emptyset i \overline{\Omega}_{h} jest zbiorem skończonym. Wtedy powiemy, że \overline{\Omega}_{h} spełnia warunek spójności siatki (ang. mesh connectivity condition, mesh is connected), jeśli dla dowolnego x\in\Omega _{h}^{{(1)}} istnieje ciąg elementów siatki \{ x_{i}\} _{{i=1}}^{N}\subset\Omega _{h}^{{(1)}} i y\in\Omega _{h}^{{(2)}} taki, że x_{1}=x, i x_{{i+1}}\in N_{h}(x_{i}) dla i=1,\ldots,N-1 i y\in N_{h}(x_{N}).

Wtedy zachodzi następująca różnicowa zasada maksimum:

Twierdzenie 9.2 (Różnicowa zasada maksimum - ang. finite difference maximum principle)

Załóżmy, że \mathcal{L}_{h} w postaci (9.2) jest operatorem dodatniego typu w \overline{\Omega}_{h}, i że \overline{\Omega}_{h} jest zbiorem skończonym spełniającym warunek spójności siatki. Wtedy, jeśli

\mathcal{L}_{h}u_{h}(x)\geq 0\qquad\forall x\in\overline{\Omega}_{h},

to

u_{h}(x)\geq 0\qquad\forall x\in\overline{\Omega}_{h}.

Dowód można znaleźć w Rozdziale 10 w [10].

Wniosek 9.1

Załóżmy, że spełnione są założenia twierdzenia 9.2. Wtedy zadanie (9.2) ma jednoznaczne rozwiązanie.

Załóżmy, że zadanie (9.2) ma dwa różne rozwiązania u_{k} dla k=1,2. Wtedy z twierdzenia 9.2 wynika, że \mathcal{L}_{h}(u_{1}-u_{2})=0 zatem (u_{1}-u_{2})\geq 0 ale i (u_{2}-u_{1})\geq 0, czyli u_{1}=u_{2}. Z kolei zauważmy, że (9.2) jest układem równań liniowych, więc jednoznaczność rozwiązania z prawą stroną równą zero jest równoważna istnieniu rozwiązania dla dowolnego \psi _{h}.

Jako kolejny wniosek z różnicowej zasady maksimum otrzymujemy następujące kryterium porównawcze:

Twierdzenie 9.3

Załóżmy, że spełnione są założenia twierdzenia 9.2 oraz niech

\mathcal{L}_{h}u_{h}(x)=f_{h}(x)\qquad\mathcal{L}_{h}v_{h}(x)=g_{h}(x)\qquad x\in\overline{\Omega}_{h}.

Wtedy, jeśli

|f_{h}(x)|\leq g_{h}(x)\qquad x\in\overline{\Omega}_{h},

to

|u_{h}(x)|\leq v_{h}(x)\qquad x\in\overline{\Omega}_{h}.

Niech z_{h}=u_{h}-v_{h}, a w_{h}=u_{h}+v_{h}. Stąd

\mathcal{L}_{h}(-z_{h}(x))=-f_{h}(x)+g_{h}(x)\geq 0,\qquad\mathcal{L}_{h}w_{h}(x)=f_{h}(x)+g_{h}(x)\geq 0\qquad x\in\overline{\Omega}_{h}.

Zatem z twierdzenia 9.2 otrzymujemy:

z_{h}(x)\leq 0,\qquad w_{h}(x)\geq 0\qquad x\in\overline{\Omega}_{h},

a stąd otrzymujemy |u_{h}(x)|\leq v_{h}(x) dla x\in\overline{\Omega}_{h}.

Z ostatniego twierdzenia otrzymujemy następujące kryterium badania stabilności w dyskretnej normie maksimum:

Twierdzenie 9.4 (kryterium stabilności z różnicowej zasady maksimum)

Załóżmy, że spełnione są założenia twierdzenia 9.2

oraz, że istnieje nieujemna funkcja v_{h} określona na \overline{\Omega}_{h} taka,że

0\leq v_{h}\leq M,\qquad\mathcal{L}_{h}v_{h}\geq 1.

Wtedy u_{h} - rozwiązanie (9.2) z prawą stroną \psi _{h}, spełnia:

\| u_{h}\| _{{\infty,h}}=\max _{{x\in\overline{\Omega}_{h}}}|u_{h}(x)|\leq M\|\psi _{h}\| _{{\infty,h}}.

Dla prostoty załóżmy, że \|\psi _{h}\| _{{\infty,h}}=1 (\mathcal{L}_{h} jest liniowe, więc zawsze możemy przeskalować u_{h} i \psi _{h} przez stałą różną od zera). Wtedy

|\mathcal{L}_{h}u_{h}(x)|=|\psi _{h}(x)|\leq 1\leq\mathcal{L}_{h}v_{h}\qquad x\in\overline{\Omega}_{h},

zatem z twierdzenia 9.3 otrzymujemy:

|u_{h}(x)|\leq v_{h}(x)\leq M=M\|\psi _{h}\| _{{\infty,h}}\qquad x\in\overline{\Omega}_{h}.
Przykład 9.2

Powróćmy do dyskretyzacji naszego modelowego zadania, tzn. do (7.5)-(7.6). Pozostawiamy jako proste zadanie sprawdzenie, że operator \mathcal{L}_{h} w tym przypadku jest operatorem dodatniego typu, i że siatka spełnia warunek spójności.

Aby pokazać oszacowanie stabilności korzystając z naszego kryterium należy znaleźć funkcję nieujemną \psi określoną na \overline{\Omega}, czyli w szczególności na każdej siatce ograniczonej, taką że

\mathcal{L}_{h}\psi\geq 1.

Na brzegu widzimy, że \mathcal{L}_{h}\psi(x)=\psi(x) dla x\in\{ a,b\}, więc wystarczy przyjąć \psi takie, że \psi\geq 1 na brzegu \Omega.

Najprościej będzie znaleźć funkcję \psi taką, że L\psi=-\frac{d^{2}\psi}{dx^{2}}\geq 2. Następnie, korzystając z tego, że rząd aproksymacji zadania przybliżonego jest dwa w każdym punkcie siatki, tzn.

|(L\psi-\mathcal{L}_{h}r_{h}\psi)(x)|=|(L\psi-L_{h}r_{h}\psi)(x)|=O(h^{2})\qquad x\in\Omega _{h},

możemy wywnioskować, że istnieje stała h_{0} taka, że dla h\leq h_{0} funkcja

\psi _{h}(x)=\psi(x)=r_{h}\psi(x)\qquad x\in\overline{\Omega}_{h}

spełnia \mathcal{L}_{h}\psi _{h}\geq 1.

W naszym przypadku np. dla c\geq 0 wystarczy zdefiniować:

\psi(x)=1+\left(\frac{b-a}{2}\right)+(x-a)*(x-b).

Wtedy 0\leq\psi\leq 1+\frac{b-a}{2}+\frac{(b-a)^{2}}{4}=M i -\frac{d^{2}\psi}{dx^{2}}+c*\psi\geq-\frac{d^{2}\psi}{dx^{2}}=2. Zatem z naszego kryterium otrzymujemy dla h\leq h_{0}, że

\| u_{h}\| _{{\infty,h}}\leq\left(1+\frac{b-a}{2}+\frac{(b-a)^{2}}{4}\right)\max\{|g(a)|,|g(b)|,\| f_{h}\| _{{\infty,\Omega _{h},h}}\},

czyli stabilność w dyskretnej normie maksimum.

Proszę zauważyć, że stała w oszacowaniu nie zależy od stałej c, za to - inaczej niż w przypadku poprzedniego prostszego kryterium, zależy od długości odcinka (a,b).

Jeśli rozwiązanie (7.1) jest w C^{4}([a,b]), to otrzymujemy, że:

\|\mathcal{L}_{h}r_{h}u-Lu\| _{{\infty,h,\Omega _{h}}}=O(h^{2})

dla \mathcal{L}_{h} z (8.7), a warunki brzegowe spełnione są dokładnie. Zatem, korzystając z twierdzenia 8.1, otrzymujemy:

\| r_{h}u-u_{h}\| _{{\infty,h}}=O(h^{2}).

9.2. Zadania

Ćwiczenie 9.1

Zbadaj stabilność schematu (8.10) dyskretyzacji modelowego problemu dwuwymiarowego w dyskretnej normie maksimum dla c>0.

Ćwiczenie 9.2

Sprawdź, czy operator z (8.10) jest dodatniego typu i zbadaj stabilność schematu (8.10) dyskretyzacji modelowego problemu dwuwymiarowego w dyskretnej normie maksimum dla c=0 korzystając z różnicowej zasady maksimum.

Ćwiczenie 9.3

Rozpatrzmy problem -\frac{du}{dt}(t)+c*u(t)=f(t) dla t\in(0,1) i f gładkiej funkcji z warunkiem brzegowym Neumanna \frac{du}{dt}(s)=0 dla s\in\{ 0,1\} oraz schemat różnicowy na siatce jednorodnej \overline{\Omega}_{h}=\{ x_{k}\} _{{k=0}}^{N} dla x_{k}=k*h:

-\overline{\partial}\partial _{h}u_{h}(x_{k})+cu_{h}(x_{k})=f(x_{k})\qquad k=1,\ldots,N-1

z \partial _{h}u_{h}(x_{0})=\overline{\partial}_{h}u_{h}(x_{N})=0. Czy to zadania wyjściowe oraz zadanie dyskretne mają jednoznaczne rozwiązanie dla c>0?

Zbadaj rząd tego schematu oraz stabilność w dyskretnej normie maksimum dla stałej c>0. Podaj oszacowanie błędu dyskretnego w dyskretnej normie maksimum w terminach O(h^{p}).

Ćwiczenie 9.4

Zbadaj rząd i stabilność w normie maksimum schematu skonstruowanego analogicznie jak schemat (8.10) dyskretyzacji modelowego problemu dwuwymiarowego: -\triangle u+c*u=f w \Omega=(0,1)^{2} z zerowym warunkiem Dirichleta na brzegu kwadratu oprócz krawędzi \Gamma _{{1}}={0}\times(0,1), gdzie jest postawiony zerowy warunek brzegowy Neumanna tzn. u(s)=0 dla s\in\partial\Omega\setminus\Gamma _{1} i \frac{\partial u}{\partial n}(0,s)=-\frac{\partial u}{\partial x}(0,s)=0 na \Gamma _{1}.

Warunek brzegowy na \Gamma _{1} przybliżamy w schemacie różnicowym przez odpowiednią różnicę skończoną wprzód, tzn. przez \partial _{1}u_{h}(0,k*h) dla k=1,\ldots,N-1 z h=\frac{1}{N}.

Ćwiczenie 9.5

Zbadaj stabilność w dyskretnej normie maksimum schematu z ćwiczenia 8.6.

Ćwiczenie 9.6

Zbadaj stabilność w dyskretnej normie maksimum schematu z ćwiczenia 8.9.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.