Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Numeryczne równania różniczkowe – 8. Teoria zbieżności schematów różnicowych – MIM UW

Zagadnienia

8. Teoria zbieżności schematów różnicowych

W tym rozdziale przedstawimy ogólną teorię zbieżności schematów różnicowych, a następnie pokażemy m.in. zastosowanie tej teorii do przykładów z poprzedniego wykładu.

Osoby zainteresowane obszerniejszym przedstawieniem teorii różnic dzielonych odsyłamy do monografii [27].

8.1. Ogólna teoria zbieżności schematów różnicowych

W tym podrozdziale opiszemy ogólną teorię zbieżności schematów różnicowych. Ograniczymy się do szczegółowego omówienia przypadku schematów liniowych, tzn. aproksymacji równań różniczkowych liniowych.

Teoria ta potrzebna jest zarówno do badania zbieżności schematów różnicowych dla równań eliptycznych, jak i dla schematów dla innych typów równań, np. równań parabolicznych.

Załóżmy, że rozpatrujemy następujące zadanie różniczkowe: chcemy znaleźć u funkcję określoną na obszarze \overline{\Omega} taką, że spełnia równanie różniczkowe z warunkami brzegowymi:

\displaystyle Lu(x) \displaystyle= \displaystyle f(x)\qquad x\in\Omega (8.1)
\displaystyle l_{k}u(x) \displaystyle= \displaystyle g_{k}(x)\qquad x\in\Gamma _{k},\qquad k=1,\ldots,s, (8.2)

gdzie f,g_{k} - to dane funkcje, L - to operator różniczkowy liniowy, l_{k} - to odpowiedni operator różniczkowy brzegowy liniowy określony na \Gamma _{k} dla \Gamma _{k}\subset\partial\Omega.

Będziemy zakładać, że powyższe zadanie jest poprawnie postawione, tzn. że ma jednoznaczne rozwiązanie u\in U dla U przestrzeni liniowej funkcji określonych na \overline{\Omega} z normą \|\cdot\| _{U}. Zakładamy też, że

L:U\rightarrow F

dla F przestrzeni funkcji określonych na \Omega, a

l_{k}:U\rightarrow\Phi _{k}\qquad k=1,\ldots,s

dla \Phi _{k} przestrzeni funkcji określonych na \Gamma _{k}. Wyjściowe zadanie różniczkowe możemy zapisać w postaci operatorowej jako: znaleźć u\in U takie, że

\displaystyle Lu \displaystyle= \displaystyle f (8.3)
\displaystyle l_{k}u \displaystyle= \displaystyle g_{k}\qquad k=1,\ldots,s. (8.4)

Zdefiniujmy \overline{\Omega}_{h} jako siatkę, tzn. zbiór punktów izolowanych, węzłów należących do \overline{\Omega} z parametrem h.

Zakładamy, że istnieje rodzina siatek \{\overline{\Omega}_{h}\} _{h}, czyli rodzina zbiorów punktów izolowanych należących do \overline{\Omega} indeksowanych parametrem h, należącym do pewnego zbioru \omega\subset(0,h_{0}]\subset\mathbb{R}_{+} takim, że 0\in\overline{\omega} (tzn. że istnieje podciąg siatek \Omega _{{h_{k}}} taki, że h_{k}\rightarrow 0).

W praktyce najczęściej stosuje się siatki równomierne, tzn. podzbiory a+h*{\mathbb{Z}}^{d} dla ustalonego punktu a\in{\mathbb{R}}^{d}. Ewentualnie stosuje się siatki o jednolitych krokach w danym kierunku w {\mathbb{R}}^{d}.

Siatkę \overline{\Omega}_{h} przedstawiamy w postaci \overline{\Omega}_{h}=\Omega _{h}\cup\partial\Omega _{h}, gdzie \Omega _{h} będziemy nazywać zbiorem punktów siatkowych wewnętrznych (zazwyczaj zawartych w \Omega), a \partial\Omega _{h} - zbiorem punktów siatkowych brzegowych (zawartych albo leżących w pobliżu \partial\Omega). W zbiorze \partial\Omega _{h} punktów brzegowych możemy dalej wyróżniać podzbiory \Gamma _{{k,h}}. Zakładamy, że rodzina siatek \{\overline{\Omega}_{h}\} _{h} jest gęsta w sensie następującej definicji:

Definicja 8.1

Rodzina siatek \{\overline{\Omega}_{h}\} _{h} jest gęsta (ang. dense) w \overline{\Omega}, gdy dla dowolnego \epsilon>0 istnieje h_{1}\in\omega takie, że dla h<h_{1} i dowolnego x\in\overline{\Omega} kula K(x,\epsilon) o środku w x i promieniu \epsilon zawiera co najmniej jeden punkt y\in\overline{\Omega}_{h}.

Proszę zauważyć, że rodzina siatek zdefiniowana w Rozdziale 7.1 jest w sposób oczywisty gęsta dla [a,b].

Wprowadzamy teraz rodzinę zadań przybliżonych (schematów różnicowych), które dają się zapisać w następujący sposób: chcemy znaleźć funkcję u_{h} określoną na \overline{\Omega}_{h} taką, że

\displaystyle L_{h}u_{h}(x) \displaystyle= \displaystyle f_{h}(x)\qquad x\in\Omega _{h}
\displaystyle l_{{k,h}}u(x) \displaystyle= \displaystyle g_{{k,h}}(x)\qquad x\in\Gamma _{{k,h}},\qquad k=1,\ldots,s,

czy inaczej - operatorowo

\displaystyle L_{h}u_{h} \displaystyle= \displaystyle f_{h} (8.5)
\displaystyle l_{{k,h}}u \displaystyle= \displaystyle g_{{k,h}}\qquad k=1,\ldots,s. (8.6)

Zakładamy, że L_{h}:U_{h}\rightarrow F_{h} i l_{{k,h}}:U_{h}\rightarrow\Phi _{{k,h}} dla k=1,\ldots,s, gdzie:

  1. U_{h} jest przestrzenią liniową unormowaną funkcji określonych na \overline{\Omega}_{h} z normą \|\cdot\| _{{U_{h}}},

  2. F_{h} jest przestrzenią liniową unormowaną funkcji określonych na \Omega _{h} z normą \|\cdot\| _{{F_{h}}},

  3. \Phi _{{k,h}} jest przestrzenią liniową unormowaną funkcji określonych na \Gamma _{{k,h}} z normą \|\cdot\| _{{\Phi _{{k,h}}}}.

Zazwyczaj wszystkie rozpatrywane przestrzenie są zupełne, tzn. są przestrzeniami Banacha. W przypadku gdy \Omega jest ograniczony, są one też przestrzeniami skończenie wymiarowymi. Jeśli L_{h} i wszystkie l_{{k,h}} są operatorami liniowymi, to mówimy, że rozpatrujemy zadanie przybliżone (dyskretne) liniowe, czy schemat różnicowy liniowy. W przeciwnym razie - gdy choć jeden z operatorów jest nieliniowy, to mamy do czynienia z zadaniem przybliżonym nieliniowym, czy schematem różnicowym nieliniowym.

Proszę zauważyć, że rozpatrujemy rodzinę zadań przybliżonych, parametryzowanych przez h. Tak, jak w przykładzie w Rozdziale 7.1, aby mówić o zbieżności rozwiązania zadania dyskretnego u_{h}\in U_{h} do u\in U musimy mieć możliwość porównania obu funkcji. Dlatego zakładamy, że istnieje rodzina operatorów obcięcia (ang. restriction) r_{h}^{U}:U\rightarrow U_{h}, które są liniowe i ograniczone jednostajnie (ang.uniformly bounded) względem h, tzn. \exists K>0\quad\forall h\in\omega\quad\forall u\in U

\| r_{h}^{U}u\| _{{U_{h}}}\leq K\| u\| _{U}.

Operator obcięcia pozwala porównywać rozwiązania w normie przestrzeni dyskretnej, ale możemy porównywać je również w normie przestrzeni U. W tym celu musimy wprowadzić rodzinę operatorów liniowych przedłużenia p_{h}^{U}:U_{h}\rightarrow U. Najczęściej za operatory przedłużenia bierze się odpowiednie operatory interpolacji.

Uwaga 8.1

Można wprowadzić pojęcie zbieżności aproksymacji przestrzeni. Tzn. rodzinę trójek (u_{h},r_{h}^{U},p_{h}^{U}) nazywamy aproksymacją przestrzeni U i mówimy, że ta aproksymacja jest zbieżna, jeśli dla dowolnego u\in U zachodzi zbieżność

\| p_{h}^{U}r_{h}^{U}u-u\| _{U}\rightarrow 0\qquad h\rightarrow 0.

W teorii zbieżności metod różnicowych najczęściej nie stosuje się operatorów przedłużenia, a za to wprowadza się warunek zgodności norm:

Definicja 8.2

Jeżeli dla danej przestrzeni unormowanej (U,\|\cdot\| _{U}) i rodziny przestrzeni unormowanych z odpowiednimi operatorami obcięcia (U_{h},\|\cdot\| _{{U_{h}}},r_{h}^{U}) zachodzi zbieżność

\lim _{{h\rightarrow 0}}\| r_{h}^{U}u\| _{{U_{h}}}=\| u\| _{U}\qquad\forall u\in U,

to mówimy, że normy dyskretne \|\cdot\| _{{U_{h}}} są zgodne (ang. consistent) z normą \|\cdot\| _{U},

Od tej pory będziemy zakładali zgodność norm dyskretnych z normą w U, według powyższej definicji.

Definicja 8.3

Zadanie przybliżone (zadanie dyskretne, schemat różnicowy) (8.5)-(8.6) jest zbieżne (czasami używa się terminu zbieżne dyskretnie) jeśli

\| r_{h}^{U}u-u_{h}\| _{{U_{h}}}\rightarrow 0\qquad h\rightarrow 0,

gdzie u - to rozwiązanie zadania (8.3)-(8.4), a u_{h}\in U_{h} - to rozwiązanie dyskretne zadania przybliżonego (8.5)-(8.6).

Jeśli dodatkowo zachodzi

\| r_{h}^{U}u-u_{h}\| _{{U_{h}}}\leq O(h^{p})

to mówimy o zbieżności (dyskretnej) rzędu p.

Wielkość \| r_{h}^{U}u-u_{h}\| _{{U_{h}}} będziemy nazywać błędem dyskretnym dla zadania przybliżonego (ang. dicrete error).

Kolejnym krokiem jest wprowadzenie pojęcia aproksymacji zadania ciągłego (wyjściowego zadania różniczkowego) przez zadanie dyskretne.

Definicja 8.4 (aproksymacja; rząd schematu; (ang. consistency))

Mówimy, że zadanie przybliżone (8.5)-(8.6) aproksymuje zadanie (8.3)-(8.4), jeśli lokalne błędy aproksymacji zdefiniowane jako

e_{{0,h}}:=\| L_{h}r_{h}^{U}u-f_{h}\| _{{F_{h}}}\qquad e_{{k,h}}=\| l_{{k,h}}r_{h}^{U}u-g_{{k,h}}\| _{{\Phi _{{k,h}}}}\qquad k=1,\ldots,s,

dążą do zera dla h\rightarrow 0. Tutaj u jest rozwiązaniem zadania (8.3)-(8.4), a f_{h},g_{{k,h}} są z zadania dyskretnego (8.5)-(8.6). Jeśli dodatkowo zachodzi:

e_{{k,h}}=O(h^{p})\qquad k=0,1,\ldots,s,

to mówimy, że schemat aproksymuje (8.3)-(8.4) z rzędem p (ang. local truncation error is of order p), (inaczej, że lokalne błędy aproksymacji są rzędu p, dane zadanie przybliżone lub schemat różnicowy ma rząd p, rząd aproksymacji schematu wynosi p).

Drugim ważnym pojęciem jest stabilność zadania dyskretnego. Tu podamy definicje stabilności dla schematu liniowego:

Definicja 8.5 (stabilność; (ang. stability))

Liniowe zadanie przybliżone (8.5)-(8.6) jest stabilne (poprawnie postawione), jeśli istnieje stała h_{0} taka, że dla dowolnego h\in\omega, h\leq h_{0} i dla dowolnych f_{h}\in F_{h} i g_{{k,h}}\in\Phi _{{k,h}} k=1,\ldots,s zachodzą:

  1. istnieje jednoznacznie wyznaczone rozwiązanie u_{h}\in U_{h} spełniające (8.5)-(8.6),

  2. rozwiązanie to spełnia następującą nierówność:

    \| u_{h}\| _{{U_{h}}}\leq C(\| f_{h}\| _{{F_{h}}}+\sum _{{k=1}}^{S}\| g_{{k,h}}\| _{{\Phi _{{k,h}}}},

    gdzie C - to dodatnia stała niezależna od h (ani oczywiście od f_{h},g_{{k,h}}).

Uwaga 8.2

W literaturze czasami za stabilność zadania przybliżonego przyjmuje się tylko warunek (2) z definicji 8.5.

Proszę zauważyć, że stabilność zadania przybliżonego jest samoistną cechą związaną tylko z definicją samego zadania dyskretnego- ona nie zależy w żaden sposób od rozwiązania równania różniczkowego. Dodatkowo warto też zauważyć, że jeśli U_{h} jest przestrzenią skończenie wymiarową, istnieje rozwiązanie (8.5)-(8.6) i spełniony jest warunek (2) z definicji 8.5, to wtedy to rozwiązanie jest jednoznaczne.

Sformułujemy teraz następujące twierdzenie o zbieżności zadania przybliżonego:

Twierdzenie 8.1 (Lax-Filipow)

Jeśli liniowe zadanie przybliżone (8.5)-(8.6) jest stabilne oraz aproksymuje zadanie (8.3)-(8.4), którego rozwiązaniem jest u, wtedy zadanie przybliżone jest zbieżne i

\| r_{h}^{U}u-u_{h}\| _{{U_{h}}}\leq C(\sum _{{k=0}}^{s}e_{{k,h}}).

Tutaj u_{h} - to rozwiązanie zadania przybliżonego (8.5)-(8.6).

Z powyższego twierdzenia otrzymujemy od razu następujący wniosek:

Wniosek 8.1

Jeśli zadanie przybliżone (8.5)-(8.6) jest stabilne oraz aproksymuje zadanie (8.3)-(8.4) z rzędem p, to

\| r_{h}^{U}u-u_{h}\| _{{U_{h}}}=O(h^{p}).

Oznaczmy z_{h}=r_{h}^{U}u-u_{h}. Z liniowości L_{h} i l_{{k,h}} dla k=1,\ldots,s wynika, że z_{h} spełnia zadanie przybliżone z odpowiednimi prawymi stronami:

L_{h}z_{h}=L_{h}r_{h}^{U}-f_{h}=w_{h},\qquad l_{{k,h}}z_{h}=l_{{k,h}}r_{h}^{U}u-g_{{k,h}}=w_{{k,h}}\qquad k=1,\ldots,s,

zatem z definicji stabilności zadania przybliżonego otrzymujemy następujące oszacowanie:

\displaystyle\| z_{h}\| _{{U_{h}}} \displaystyle\leq \displaystyle C(\| g_{h}\| _{{F_{h}}}+\sum _{{k=1}}^{s}\| w_{{k,h}}\| _{{\Phi _{{k,h}}}})=C(\| L_{h}r_{h}^{U}-f_{h}\| _{{F_{h}}}+\sum _{{k=1}}^{s}\| l_{{k,h}}r_{h}^{U}u-g_{{k,h}}\| _{{\Phi _{{k,h}}}}).

Następnie z faktu aproksymacji zadania (8.3)-(8.4) przez zadanie przybliżone (8.5)-(8.6) otrzymujemy ostatecznie oszacowanie:

\| r_{h}^{U}u-u_{h}\| _{{U_{h}}}=\| z_{h}\| _{{U_{h}}}\leq C(\sum _{{k=0}}^{s}e_{{k,h}})=O(h^{p})\rightarrow 0\qquad h\rightarrow 0.

Powyższe twierdzenie można krótko podsumować, że aby otrzymać schemat zbieżny z rzędem p musi być on stabilny i posiadać rząd aproksymacji p.

Proszę zauważyć, że twierdzenie jest bardzo ogólne, a dowód jest prosty. Pojawia się pytanie: jak dobrać odpowiednie przestrzenie i operatory, aby zadania przybliżone (schematy) były stabilne i miały możliwie wysoki rząd aproksymacji.

Uwaga 8.3

Proszę zauważyć, że powyższa teoria zbieżności może zostać zastosowana do zadań różniczkowych rożnego typu - zarówno eliptycznych, jak i parabolicznych, czy hiperbolicznych.

8.2. Zastosowanie teorii zbieżności do prostych schematów jedno- i dwuwymiarowych

8.2.1. Przypadek jednowymiarowy

Wracamy teraz do dyskretyzacji modelowego zadania jednowymiarowego (7.5)-(7.6). Za przestrzeń U weźmy przestrzeń funkcji ciągłych na \overline{\Omega}, czyli C([a,b]) z normą supremum \| u\| _{\infty}=\sup _{{t\in[a,b]}}|u(t)|. Oznaczmy przez C_{h}(K_{h}) zbiór funkcji określonych na dowolnym podzbiorze K_{h} siatki \overline{\Omega}_{h} z normą \| u\| _{{\infty,h,K_{h}}}=\max _{{x\in K_{h}}}|u(x)|. Za przestrzeń dyskretną przyjmijmy U_{h}=C(\overline{\Omega}_{h}). Jeżeli wprowadzimy operatory L_{h}:U_{h}\rightarrow F_{h} i l_{h}:U_{h}\rightarrow\Phi _{h} zdefiniowane jako (por. (7.3)):

\displaystyle L_{h}u_{h}(x) \displaystyle= \displaystyle-\partial\overline{\partial}u_{h}(x)+c*u_{h}(x)\qquad x\in\Omega _{h}
\displaystyle l_{h}u_{h}(x) \displaystyle= \displaystyle u(x)\qquad x\in\partial\Omega _{h}=\overline{\Omega}_{h}\setminus\Omega _{h}

dla F_{h}=C_{h}(\Omega _{h}) i \Phi _{h}=C_{h}(\partial\Omega _{h}), to zadanie (7.5)-(7.6) możemy zapisać w formie operatorowej jako:

\displaystyle L_{h}u_{h}=f_{h}, (8.7)
\displaystyle l_{h}u_{h}=g_{h}

dla f_{h}\in F_{h} zdefiniowanego jako (f_{h})(x)=f(x) dla x\in\Omega _{h} oraz g_{h}\in\Phi _{h} z (g_{h})(x)=g(x) dla x\in\partial\Omega _{h}.

W tym przypadku dla funkcji ciągłej przekształceniem obcięcia (ang. restriction) jest r_{h}^{{C([a,b])}}:U\rightarrow U_{h} zdefiniowany jako

(r_{h}^{{C([a,b])}}u)(x)=r_{h}u(x)=u(x)\qquad x\in\overline{\Omega}_{h}.

Możemy teraz zbadać zbieżność błędu dyskretnego:

\| r_{h}u-u_{h}\| _{{\infty,h}}

dla h\rightarrow 0, co jest równoważne badaniu zbieżności w punktach siatki. Zauważmy, że otrzymujemy

\| r_{h}u\| _{{\infty,h}}\leq\| u\| _{\infty}\qquad\forall u\in U,

co oznacza jednostajną ograniczoności operatorów obcięcia.

Można też w tym przypadku łatwo wprowadzić operator przedłużenia (ang. prolongation) p_{h}:U_{h}\rightarrow U. Np. niech p_{h} będzie funkcją ciągłą liniowo interpolującą wartości u_{h} pomiędzy punktami siatki tj.

p_{h}(x)=u_{h}(x_{k})+\frac{u_{h}(x_{{k+1}})-u_{h}(x_{k})}{h}(x-x_{k})\qquad x_{k}\leq x\leq x_{{k+1}}.

Następnie możemy badać zbieżność błędu \| p_{h}u_{h}-u\| _{\infty} dla h\rightarrow 0. Jeśli błąd zbiega do zera, to mówimy o zbieżności schematu w normie supremum.

Zauważmy, że w naszym przypadku dodatkowo zachodzi

\| p_{h}r_{h}u-u\| _{\infty}\rightarrow 0\qquad h\rightarrow 0\qquad\forall u\in U, (8.8)

czyli zbieżność aproksymacji przestrzeni wyjściowej przez przestrzeń dyskretną oraz

\lim _{{h\rightarrow 0}}\| r_{h}u\| _{{\infty,h}}=\| u\| _{\infty}\qquad\forall u\in U, (8.9)

czyli zachodzi zgodność rodziny norm przestrzeni dyskretnych (U_{h},\|\cdot\| _{{\infty,h}}) z normą przestrzeni wyjściowej (U,\|\cdot\| _{\infty}). Wykazanie tego pozostawiamy jako zadanie, por. ćwiczenie 8.1.

Innym wyborem przestrzeni i norm jest badanie zbieżności i błędu w normie L^{2}(a,b), czy odpowiednio dyskretnej normie L^{2}_{h} definiowanej dla u\in U_{h} jako:

\| u\| _{{0,h}}=\sqrt{h\sum _{{x\in\overline{\Omega}_{h}}}|u(x)|^{2}},

gdzie U_{h}=L^{2}_{h}(\overline{\Omega}_{h}) - to przestrzeń wszystkich funkcji określonych na \overline{\Omega}_{h}. Oczywiście zmieniliśmy oznaczenie przestrzeni funkcji dyskretnych na siatce. Jest to ten sam zbiór funkcji określonych na siatce, ale zmieniła się norma dyskretna.

Aby otrzymać zgodność norm powinniśmy inaczej zdefiniować obcięcia np. poprzez uśrednienia, czyli r_{h}^{{L^{2}}}:L^{2}(\Omega)\rightarrow L^{2}_{h}(\overline{\Omega}_{h}) dla x\in\overline{\Omega}_{h} definiujemy:

r_{h}^{{L^{2}}}u(x)=\frac{1}{K(x,h)\cap\Omega}\int _{{B(x,h)\cap\Omega}}u(x)\, dx.

dla K(x,h) kuli o środku w x i promieniu h.

Inna możliwość to rozważenie zbioru funkcji ciągłych U=C([a,b]) ale z normą L^{2}, oraz normy dyskretnej typu L^{2}_{h} na U_{h}. Następnie możemy przeprowadzić analizę z obcięciem r_{h}:=r_{h}^{{C([a,b])}}u. Zbiór funkcji U=C([a,b]) z normą L^{2} nie jest przestrzenią zupełną, ale jest gęstą podprzestrzenią przestrzeni L^{2}(a,b).

Nietrudno zauważyć, że problem przybliżony aproksymuje problem wyjściowy z rzędem dwa, o ile rozwiązanie należy do C^{4}(\overline{\Omega}), w obu powyżej przedstawionych przestrzeniach dyskretnych, czyli w odpowiednich normach dyskretnych. Wykazanie, że schemat jest stabilny zarówno w C(\overline{\Omega}_{h}) jak i L^{2}(\overline{\Omega}_{h}) jest trudniejsze. Zajmiemy się tym w kolejnych wykładach.

8.2.2. Przypadek dwuwymiarowy

Rozpatrzmy ponownie modelowe zadanie dwuwymiarowe na kwadracie jednostkowym (7.8). Analogicznie, jak w przypadku jednowymiarowym, niech U=C([0,1]^{2}) z normą supremum \| u\| _{\infty}=\sup _{{t\in[0,1]^{2}}}|u(t)| i C_{h}(K_{h}) będzie przestrzenią funkcji określonych na podzbiorze K_{h} siatki \overline{\Omega}_{h} (por. (7.10)) z normą \| u\| _{{\infty,h,K_{h}}}=\max _{{x\in K_{h}}}|u(x)|. Przestrzeń dyskretną definiujemy jako U_{h}=C(\overline{\Omega}_{h}).

Operator siatkowy (ang. mesh operator or discrete operator) L_{h}:U_{h}\rightarrow F_{h} i brzegowy l_{h}:U_{h}\rightarrow\Phi _{h} definiujemy jako:

\displaystyle(L_{h}u_{h})(x) \displaystyle= \displaystyle(-\sum _{{k=1,2}}\partial\overline{\partial}_{k}+c)u_{h}(x),\qquad x\in\Omega _{h}
\displaystyle(l_{h}u_{h})(x) \displaystyle= \displaystyle u(x),\qquad x\in\partial\Omega _{h}=\overline{\Omega}_{h}\setminus\Omega _{h}

dla F_{h}=C_{h}(\Omega _{h}) i \Phi _{h}=C_{h}(\partial\Omega _{h}). Teraz zadanie (7.11) możemy zapisać w formie operatorowej jako

\displaystyle\left\{\begin{array}[]{rcl}L_{h}u_{h}&=&f_{h}\\
l_{h}u_{h}&=&g_{h}\end{array}\right. (8.10)

dla funkcji prawej strony f_{h}\in F_{h} oraz g_{h}\in\Phi _{h} zdefiniowanych jako f_{h}(x)=f(x) dla x\in\Omega _{h} i g_{h}(x)=g(x) dla x\in\partial\Omega _{h}. Operatorem obcięcia (ang. restriction) jest r_{h}:U\rightarrow U_{h}, zdefiniowany jako

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

Tak samo jak w przypadku jednowymiarowym badamy błąd: \| r_{h}u-u_{h}\| _{{\infty,h}}, lub w dyskretnej normie L^{2}, tzn. w

\| u\| _{{0,h}}=\sqrt{h^{2}\sum _{{x\in\overline{\Omega}_{h}}}|u(x)|^{2}}.

Tu U_{h}=L^{2}_{h}(\overline{\Omega}_{h}) jest zdefiniowana jako przestrzeń wszystkich funkcji określonych na \overline{\Omega}_{h}. Oczywiście zbiór funkcji siatkowych jest ten sam, zmieniła się tylko norma.

Można pokazać, że zachodzi zgodność norm dyskretnych z odpowiednimi normami, oraz że schemat (7.11) posiada rząd aproksymacji dwa i jest stabilny w obu normach dyskretnych. Wykazanie rzędu aproksymacji jest prostym zadaniem, natomiast pokazanie stabilności jest trudniejsze, por. rozdziały 9 i rozdział 10.

8.3. Zadania

Ćwiczenie 8.1

Udowodnij (8.8) i (8.9).

Ćwiczenie 8.2

Zbadaj rząd lokalnych błędów aproksymacji schematu (8.7) dyskretyzacji modelowego problemu jednowymiarowego w obu normach dyskretnych.

Ćwiczenie 8.3

Wykaż, że rząd aproksymacji schematu (7.11) w dyskretnych normach maksimum i L^{2}_{h} wynosi dwa, o ile rozwiązania wyjściowego zadania różniczkowego są dostatecznie gładkie.

Ćwiczenie 8.4

Zbadaj rząd lokalnych błędów aproksymacji schematu (8.10) dyskretyzacji modelowego problemu dwuwymiarowego w obu normach dyskretnych.

Ćwiczenie 8.5

(Przybliżony warunek brzegowy) Rozpatrzmy modelowe zadanie jednowymiarowe z warunkiem brzegowym Dirichleta:

-u^{{\prime\prime}}(x)=f(x)\qquad x\in(0,1)\qquad x(0)=a\quad x(1)=b

dla f\in C^{\infty}.

Rozpatrzmy następującą dyskretyzację zbudowaną na siatce \overline{\Omega}_{h}=\{ x_{k}\} _{{k=0,\ldots,N-1}} dla x_{k}=k*h z k=0,\ldots,N-1 z (N-1)*h<1<N*h. Definiujemy \Omega _{k}=\{ k*h\} _{{k=1,\ldots,N-2}} i \partial\Omega _{h}=\{ 0,(N-1)*h\}, oraz operatory:

L_{h}u_{h}(x)=\partial\overline{\partial}u_{h}(x)\qquad x\in\Omega _{h}

i l_{h}u_{h}(0)=a,\; l_{h}u_{h}(x_{{N-1}})=b. Zbadaj rząd lokalnego błędu aproksymacji tej dyskretyzacji w dyskretnej normie maksimum.

Wskazówka: 

Wystarczy zbadać błąd operatora brzegowego w punkcie x_{{N-1}}. W pozostałych punktach błąd jest jak w schemacie (8.7).

Ćwiczenie 8.6

Rozważmy modelowe zadanie, jak i siatkę niezawierającą prawy koniec obszaru, tak jak w poprzednim ćwiczeniu.

Operatory L_{h} i l_{h}(0)=a definiujemy tak samo, natomiast zmodyfikujmy operator brzegowy l_{h} w prawym końcu, tzn. w punkcie x_{r}:=x_{{N-1}}<1.

Rozpatrzmy tzw. aproksymację Collatza, tzn. niech wartość l_{h}(x_{r})=b będzie liniowo interpolowała warunek brzegowy w końcu obszaru:

l_{h}u_{h}(x_{r})=u_{h}(x_{r})+\frac{u_{h}(x_{r})-u_{h}(x_{r}-h)}{h}(1-x_{r})=(1+\frac{\tilde{h}}{h})u_{h}(x_{r})-\frac{\tilde{h}}{h}u_{h}(x_{r}-h)=b,

gdzie \tilde{h}=1-x_{r}<h.

Zbadaj lokalny błąd aproksymacji tego schematu w normach dyskretnych maksimum i L^{2}_{h} i jego rząd, tzn. czy zachowuje się jak O(h^{p}) dla pewnego p naturalnego.

Ćwiczenie 8.7

Rozpatrzmy zadanie z poprzedniego ćwiczenia, ale z siatką nierównomierną: \overline{\Omega}_{h}=\{ x_{k}\} _{{k=0,\ldots,N-1}}\cup\{ 1\} dla x_{k}=k*h z k=0,\ldots,N-1 z (N-1)*h<1<N*h. Operator l_{h} możemy zdefiniować jako

l_{h}u(0)=l_{h}u(\tilde{x})=0,

gdzie \tilde{x}=1, ale musimy zmodyfikować definicję L_{h}, tzn.

L_{h}u(x_{k})=-\partial\overline{\partial}u(x_{k})\qquad k=1,\ldots,N-2

i

L_{h}u(x_{{N-1}})=au(x_{{N-2}})+bu(x_{{N-1}})+cu(\tilde{x}).

Wyznacz a,b,c w zależności od wartości h i \tilde{h}=1-x_{{N-1}}<h tak, aby lokalny błąd aproksymacji schematu był możliwie mały.

Ćwiczenie 8.8

Rozpatrzmy zadanie jednowymiarowe -\frac{d^{2}u}{dt^{2}}+u=f na [(0,1) z warunkiem Neumanna \frac{du}{dt}(0)=\frac{du}{dt}(1)=0. Rozpatrzmy następującą dyskretyzację zbudowaną na siatce \overline{\Omega}_{h}=\{ x_{k}\} _{{k=0,\ldots,N-1}} dla x_{k}=k*h z h=1/N.

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

oraz

l_{1}u(0)=\partial _{h}u(0)=0,\qquad l_{2}u(0)=\overline{\partial}_{h}u(1)=0.

Zbadaj rząd lokalnego błędu aproksymacji tego schematu względem parametru siatki h w dyskretnej normie maksimum i dyskretnej normie L^{2}_{h}.

Ćwiczenie 8.9

Rozpatrzmy zadanie jednowymiarowe -\frac{d^{2}u^{*}}{dt^{2}}+u^{*}=f na (0,1) z warunkiem Neumanna \frac{du^{*}}{dt}(0)=\frac{du^{*}}{dt}(1)=0. Rozpatrzmy następującą dyskretyzację o podwyższonym rzędzie zbudowaną na siatce \overline{\Omega}_{h}=\{ x_{k}\} _{{k=0,\ldots,N-1}} dla x_{k}=k*h z h=1/N. W punktach wewnętrznych siatki stosujemy standardowo aproksymacje na trzech punktach:

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

natomiast na brzegu podnosimy rząd schematu, a dokładniej zakładamy, że równanie jest spełnione w punktach brzegu, tzn. funkcja f jest określona na \overline{\Omega}=[0,1] i -\frac{d^{2}u^{*}(x)}{dt^{2}}+u^{*}(x)=f(x) dla x\in\{ 0,1\}. Rozpatrzmy lewy punkt brzegu x=0. Widzimy, że

\partial _{h}u^{*}(0)=\frac{du^{*}}{dt}(0)+\frac{d^{2}u}{dt^{2}}(0)*\frac{h}{2}+O(h^{2})=\frac{du^{*}}{dt}(0)+(u^{*}(0)-f(0))\frac{h}{2}+O(h^{2}),

o ile u jest dostatecznie gładka. Zatem - korzystając z obu faktów - możemy skonstruować równanie różnicowe:

\partial _{h}u_{h}(0)-u_{h}(0)\frac{h}{2}=-f(0)\frac{h}{2}

przybliżające warunek Neumanna w punkcie x=0 z wyższym rzędem.

Skonstruuj analogiczne równanie różnicowe przybliżające warunek Neumanna w punkcie x=1 z wyższym rzędem. Pokaż, że rząd lokalnego błędu aproksymacji tego schematu względem parametru siatki h wynosi dwa w dyskretnej normie maksimum i dyskretnej normie L^{2}_{h} dla odpowiednio gładkiego rozwiązania. Przetestuj w octavie rząd lokalnego błędu schematu w normie dyskretnej maksimum dla u=\cos(\pi x) metodą połowienia kroków.

Ćwiczenie 8.10

Rozpatrzmy modelowe zadanie dwuwymiarowe na kole o średnicy jeden tzn. (7.8) dla \overline{\Omega}=\overline{K}(0,1). Dobierzmy siatkę na płaszczyźnie o parametrze h równomierną zawierającą punkt (0,0), tzn. \{(k*h,l*h)\} _{{k,l}}.

Za \Omega _{h} uznajmy wszystkie punkty siatki, które należą do \Omega i wszystkie punkty przecięcia prostych zadających siatkę z brzegiem \Omega. Te punkty przecięcia uznajemy za brzegowe punkty siatki. Otrzymujemy oczywiście siatkę nierównomierną, bo odległość między brzegowym punktem siatki, a jego sąsiadem wewnętrznym jest mniejsza od h (poza ewentualnie pojedynczymi punktami).

Tu warunek brzegowy możemy zadać dokładnie. Pojawia się pytanie: jak przybliżyć drugą pochodną w punktach wewnętrznych siatki, których punkty sąsiednie są na brzegu?

Definiujemy w takim punkcie x\in\Omega _{h} (załóżmy, że tylko jego prawy sąsiad x_{r}=x+\tilde{h}\vec{e}_{1} jest na brzegu):

L_{h}u(x)=a*u(x-h\vec{e}_{1})+bu(x)+cu(x+\tilde{h}\vec{e}_{1})+-\partial\overline{\partial}_{2}u(x).

dla pewnych parametrów a,b,c.

Jeśli x ma dwa punkty sąsiednie leżące na brzegu (powiedzmy prawy i dolny punkt sąsiedni), tzn. x_{p}=x+\tilde{h}\vec{e}_{1},x_{d}=x-\hat{h}\vec{e}_{1} są na brzegu, to oczywiście musimy wyznaczyć całe równanie różnicowe:

L_{h}u(x)=a*u(x-h\vec{e}_{1})+bu(x)+cu(x+\tilde{h}\vec{e}_{1})+\alpha*u(x-h\vec{e}_{1})+\beta u(x)+\gamma u(x-\hat{h}\vec{e}_{2}).

Zadanie: Wyznacz odpowiednie parametry a,b,c, czy \alpha,\beta,\gamma tak, aby lokalny błąd schematu |L_{h}u(x)-Lu(x)| był możliwie mały, tzn. żeby schemat posiadał możliwie wysoki rząd lokalnego błędu aproksymacji względem h w dyskretnej normie maksimum.

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.