Zagadnienia

12.1 Metoda elementu skończonego na kwadracie jednostkowym
12.2 Niejednorodny warunek brzegowy
12.3 Zadania

12. Metoda elementu skończonego - wprowadzenia cd. Przypadek dwuwymiarowy.

W tym rozdziale przedstawimy kilka możliwie prostych dyskretyzacji skonstruowanych za pomocą metody elementu skończonego na przykładzie modelowego zadania eliptycznego na kwadracie.

12.1. Metoda elementu skończonego na kwadracie jednostkowym

Rozpatrzmy modelowe zadanie (7.8) z rozdziału 7.2, którego rozwiązanie oznaczymy przez $u_{*}$ .

Słabe sformułowanie zadania (7.8): chcemy znaleźć $u_{*}\in V$ takie, że

$\int _{\Omega}\nabla u_{*}\nabla v+cu\, v\, dx=\int _{\Omega}fv\, dx\qquad\forall v\in V.$

(12.1)

Tutaj $V$ jest odpowiednio dobraną przestrzenią funkcji ciągłych zerujących się na brzegu, dla których obie strony słabego sformułowania mają sens. Np. możemy wziąć domknięcie w odpowiedniej normie przestrzeni funkcji ciągłych zerujących się na brzegu, których słaba pochodna (pochodna w sensie dystrybucyjnym) jest w $L^{2}(\Omega)$ . Później precyzyjniej ustalimy o jaką przestrzeń chodzi.

Jeśli problem (12.1) ma rozwiązanie dostatecznie gładkie, tzn. $u_{*}$ posiada ciągłe pierwsze i drugie pochodne cząstkowe, to $u_{*}$ jest rozwiązaniem wyjściowego zadania. Może się zdarzyć, że istnieje rozwiązanie (12.1), które nie jest nawet ciągłe.

$\$

Rys. 12.1. Modelowa triangulacja złożona z trójkątów $\tau$ na kwadracie $\Omega$ .

12.1.1. Triangulacja obszaru

Wprowadźmy podział (triangulacje) kwadratu $\overline{\Omega}=[0,1]^{2}$ na trójkąty $T_{h}([0,1]^{2})=\{\tau _{k}\} _{k}$ o jednakowym kształcie i wielkości. Najprościej jest podzielić kwadrat na równe kwadraty $([x^{{(k)}},x^{{(k+1)}}]\times[x^{{(l)}},x^{{(l+1)}}]$ dla $k,l=0,\ldots,N-1$ i $x^{{(k)}}=k*h$ dla $h=1/N$ . Następnie każdy kwadrat dzielimy na dwa trójkąty prowadząc przekątną np. z lewego górnego rogu do dolnego prawego, por. rysunek 12.1.

Zauważmy, że $\bigcup _{{\tau\in T_{h}(\Omega)}}\overline{\tau}=\overline{\Omega}$ i $\partial\tau _{k}\cap\partial\tau _{j}$ jest zbiorem pustym, krawędzią lub wspólnym wierzchołkiem dla dowolnych różnych elementów tej triangulacji.

Za parametr tej dyskretyzacji przyjmujemy $h$ , a punkty $x^{{(kl)}}=(x^{{(k)}},x^{{(l)}})$ nazwiemy punktami nodalnymi tej triangulacji. Proszę zauważyć, że to nie jest tylko jedna możliwa triangulacja kwadratu. Możemy wybierać trójkąty na wiele sposobów tak, aby tylko zachowane zostały warunki, że trójkąty są rozłączne, suma ich domknięć tworzy cały kwadrat, ich wspólne części brzegów to wierzchołek, wspólna krawędź lub zbiór pusty.

Rozpatrujemy w ogólności rodzinę takich triangulacji z $h\rightarrow 0$ .

12.1.2. Element liniowy

Możemy wprowadzić przestrzeń dyskretną:

$V^{h}=\{ u\in C(\overline{\Omega}):u_{{|\tau}}\in P_{1}(\tau)\quad\forall\tau\in T_{h}(\Omega),\quad u(s)=0\quad\forall s\in\partial\Omega\},$

dla $P_{1}(\tau)=\mathrm{span}(1,x_{1},x_{2})$ przestrzeni wielomianów liniowych na $\tau$ . Wielomian liniowy na trójkącie jest zdefiniowany poprzez wartości w wierzchołkach tego trójkąta, które nazywamy punktami swobody tego elementu skończonego, por. rysunek 12.2.

$\$

Rys. 12.2. Wierzchołki trójkąta $\tau$ - punkty nodalne wielomianu liniowego na tym trójkącie.

Analogicznie jak w przypadku jednowymiarowym, możemy wprowadzić funkcje nodalne: definiujemy $\phi _{{k,l}}\in V^{h}$ dla $k,l\in\{ 1,\ldots,N-1\}$ jako funkcję, która spełnia $\phi _{{k,l}}(x^{{(kl)}})=1$ i $\phi _{{k,l}}(x^{{(ij)}})=0$ dla $x^{{(kl)}}\not=x^{{(ij)}}$ dla $x^{{(kl)}}=(x^{{(k)}},x^{{(l)}})$ .

$\$

Rys. 12.3. Wykres funkcji daszkowej $\phi _{{kl}}$ .

$\$

Rys. 12.4. Wykresy dwóch funkcji daszkowych.

$\$

Rys. 12.5. Wykres przykładowej funkcji z przestrzeni elementu skończonego.

W przypadku naszej prostej regularnej triangulacji możemy wyznaczyć wzory na te funkcje na trójkącie $\tau$ . Powiedzmy, że wzory wyznaczymy na tym trójkącie, który jest zaznaczony na rysunku 12.1 - przyjmujemy, że $x^{{(kl)}}=(k*h,l*h)$ jest wierzchołkiem przy kącie prostym:

$\displaystyle\phi _{{k,l}}(x)$	$\displaystyle=$	$\displaystyle 1-\frac{x_{1}-kh}{h}-\frac{x_{2}-lh}{h},$
$\displaystyle\phi _{{k+1,l}}(x)$	$\displaystyle=$	$\displaystyle\frac{x_{1}-k*h}{h},$	(12.2)
$\displaystyle\phi _{{k,l+1}}(x)$	$\displaystyle=$	$\displaystyle\frac{x_{2}-l*h}{h}.$

dla $x=(x_{1},x_{2})$ . Wykresy takiej funkcji dla tej regularnej triangulacji kwadratu lub kilku funkcji możemy obejrzeć na rysunkach 12.3 i 12.4. Na rysunku 12.5 widzimy przykładową funkcję z przestrzeni $V^{h}$ . Również jako zadanie pozostawimy wykazanie, że

$\{\phi _{{k,l}}\} _{{k,l=1,\ldots,N-1}}$

tworzą bazę $V^{h}$ , i że

$u=\sum _{{k,l}}u(x^{{(kl)}})\phi _{{k,l}}.$

Wprowadzamy zadanie dyskretne: chcemy znaleźć $u_{h}\in V^{h}$ takie, że

$\int _{\Omega}\nabla u_{h}\nabla v_{h}+cu_{h}\, v_{h}\, dx=\int _{\Omega}fv_{h}\, dx\qquad\quad\forall v_{h}\in V^{h}.$

(12.3)

Można pokazać, że to zadanie ma jednoznaczne rozwiązanie, i że jeśli $u_{*}\in C^{2}([0,1]^{2})$ to:

	$\displaystyle\\|\nabla(u_{h}-u_{*})\\| _{{L^{2}(\Omega)}}$	$\displaystyle=$	$\displaystyle O(h)$		(12.4)
	$\displaystyle\\| u_{h}-u_{*}\\| _{{L^{2}\Omega)}}$	$\displaystyle=$	$\displaystyle O(h^{2}).$

Zatem w normie $L^{2}$ widzimy oszacowanie zbieżności rzędu dwa analogiczne jak dla dyskretyzacji tego samego zadania przy pomocy metody różnic dzielonych i dyskretnej normy $L^{2}_{h}$ , ale przy dużo słabszych założeniach. Tam musieliśmy założyć, że funkcja jest klasy $C^{4}$ .

Można też zauważyć, że w przypadku bardziej skomplikowanego geometrycznie obszaru (wielokąta) możemy skonstruować analogiczną dyskretyzację wprowadzając triangulację złożoną z trójkątów, a w przypadku różnic dzielonych, jeśli obszar nie jest prostokątem, pojawiają się kłopoty z postawieniem warunku brzegowego.

12.1.3. Element kwadratowy i kubiczny

W tym rozdziale przedstawimy dwa typy przestrzeni elementu skończonego wyższych rzędów - element kwadratowy i kubiczny. Rozpatrzmy triangulację kwadratu jednostkowego $T_{h}([0,1]^{2})$ jak w Rozdziale 12.1.1. Wtedy definiujemy:

$V^{h}_{p}=\{ u\in C(\overline{\Omega}):u_{{|\tau}}\in P_{p}(\tau)\quad\forall\tau\in T_{h}(\Omega),\quad u(s)=0\quad\forall s\in\partial\Omega\}.$

dla $p=1,2,3,\ldots$ . Tutaj $P_{p}(\tau)=\mathrm{span}(x_{1}^{k}x_{2}^{l})_{{0\leq k+l\leq p}}$ to przestrzeń wielomianów na $\tau$ stopnia nie większego od $p$ .

Przypadek $V^{h}_{2}$ określamy jako przestrzeń elementu skończonego kwadratowego, a $V^{h}_{3}$ kubicznego.

Dla $p=2$ widzimy, że wielomian kwadratowy na trójkącie jest określony jednoznacznie poprzez swoje wartości w trzech wierzchołkach i trzech środkach krawędzi. Wszystkie te punkty wewnątrz $\Omega$ określamy jako punkty nodalne $V^{h}_{2}$ . Z każdym takim punktem $x$ wiążemy funkcję, która jest równa jeden w tym punkcie, a zero we wszystkich pozostałych wierzchołkach i punktach środkowych krawędzi, por. rysunek 12.6.

$\$

Rys. 12.6. Wierzchołki i punkty środkowe krawędzi trójkąta $\tau$ -punkty swobody wielomianu kwadratowego na tym trójkącie.

Z kolei dla $p=3$ wielomian kubiczny na trójkącie jest jednoznacznie określony poprzez swoje wartości w wierzchołkach, w środku ciężkości oraz w dwu punktach wewnątrz każdej krawędzi dzielących ją na trzy równe odcinki, por. rysunek 12.7.

$\$

Rys. 12.7. Punkty swobody wielomianu kubicznego na trójkącie $\tau$ .

W przypadku przestrzeni $V^{h}_{p}$ możemy wprowadzić analogiczne bazy nodalne złożone z funkcji, które przyjmują wartość jeden w ustalonym punkcie nodalnym, a zero - w pozostałych. Wzory takich funkcji są coraz bardziej skomplikowane wraz ze wzrostem $p$ . Dla $p=2$ wypiszemy wzory na obcięcie funkcji nodalnej na ustalonym trójkącie $\tau$ . Niech punkty $x^{{(k)}}$ $k=0,1,2$ będą trzema wierzchołkami tego trójkąta $\tau$ i niech punkt $m^{{(k)}}$ $k=0,1,2$ będzie środkiem odcinka między $x^{{(k)}}$ a $x^{{(k+1)}}$ . Przy czym utożsamiamy $0$ z $3$ . Niech $q_{k}$ będzie funkcją liniową taką, że $q_{k}(x^{{(j)}})=0$ dla $j\not=k$ i $q_{k}(x^{{(k)}})=1$ . Wtedy funkcja nodalna związana z wierzchołkiem elementu $x^{{(k)}}$ dla elementu kwadratowego wynosi na trójkącie $\tau$ :

$\phi _{{x^{{(k)}}}}=q_{k}(q_{k}-q_{{k-1}}-q_{{k+1}})\qquad k=0,1,2\; mod\; 3,$

(12.5)

a funkcja nodalna związana z punktem środkowym krawędzi elementu na $\tau$ :

$\phi _{{m^{{(k)}}}}=4*q_{k}*q_{{k+1}}\qquad k=0,1,2\; mod\; 3.$

(12.6)

Jako zadanie pozostawiamy sprawdzenie tych wzorów i znalezienie analogicznych dla funkcji bazowych elementu kubicznego. Zastosowanie elementów: kwadratowego ( $p=2$ ), kubicznego ( $p=3$ ) czy nawet dla większych $p$ jest korzystne, o ile $u_{*}$ - czyli rozwiązanie zadania różniczkowego (12.1) jest bardziej regularne, tzn. należy do $C^{{s+1}}(\overline{\Omega})$ dla $s$ większych od jeden do $s=p$ . Tzn. wtedy można wykazać, że

$\|\nabla(u_{h}-u_{*})\| _{{L^{2}(\Omega)}}=O(h^{s}).$

12.1.4. Metoda elementu skończonego z podziałem obszaru na prostokąty

$\$

Rys. 12.8. Przykład obszarów będących sumą prostokątów.

W przypadku, gdy nasz obszar jest prostokątem lub kwadratem, możemy wprowadzić przestrzeń elementu skończonego: tzw. elementów skończonych prostokątnych. Taki element możemy w ogólności zastosować, jeśli obszar jest sumą prostokątów o brzegach będących sumą odcinków równoległych do osi współrzędnych; np. tzw. $L$ -obszarem, por. rysunek 12.8.

Opiszemy tę metodę dla naszego modelowego zadania postawionego na kwadracie. Dla $\Omega=(0,1)^{2}$ wprowadźmy triangulacje złożoną z kwadratów, por. rysunek 12.9:

$T_{h}(\Omega)=\{\tau _{{kl}}\} _{{k,l=0,\ldots,N-1}}$

dla $h=1/N$ i dla kwadratów:

$\tau _{{kl}}=(k*h,(k+1)*h)\times(l*h,(l+1)*h),$

$\$

Rys. 12.9. Przykład triangulacji kwadratu jednostkowego złożonej z kwadratów.

$\$

Rys. 12.10. Wierzchołki kwadratu $\tau$ -punkty swobody wielomianu biliniowego na tym kwadracie.

Wtedy przestrzeń dyskretną definiujemy jako przestrzeń funkcji ciągłych biliniowych na kwadratach:

$V^{h}_{{bi}}=\{ u\in C(\overline{\Omega}):u_{{|\tau}}\in Q_{1}(\tau)\quad\forall\tau\in T_{h}(\Omega),\quad u(s)=0\quad\forall s\in\partial\Omega\}.$

gdzie

$Q_{1}(\tau _{{kl}})=P_{1}((k*h,(k+1)*h))\otimes P_{1}((l*h,(l+1)*h))=\mathrm{Span}(1,x_{1},x_{2},x_{1}\, x_{2})$

przestrzeń wielomianów biliniowych określonych na $\tau _{{kl}}$ , czyli liniowych ze względu na każdą zmienną z osobna.

Wielomian biliniowy na prostokącie jest zdefiniowany poprzez wartości w wierzchołkach tego prostokąta. Nazywamy je punktami swobody tego elementu skończonego, por. rysunek 12.10. Analogicznie do przypadku elementów trójkątnych możemy tu wprowadzić funkcje nodalne: definiując $\phi _{{k,l}}\in V^{h}_{{bi}}$ dla $k,l\in\{ 1,\ldots,N-1\}$ jako funkcję, która spełnia $\phi _{k}(x^{{(kl)}})=1$ i $\phi _{k}(x^{{(ij)}})=0$ dla $x^{{(kl)}}\not=x^{{(ij)}}$ dla $x^{{(kl)}}=(x^{{(k)}},x^{{(l)}})$ . Podanie wzoru na te funkcje jest nietrudnym zadaniem z uwagi na prostą geometrię elementu. Pozostawiamy to jako zadanie, por. rysunek 12.11

$\$

Rys. 12.11. Wykres funkcji bazowej przestrzeni MES biliniowej na kwadracie.

Następnie definiujemy zadanie dyskretne: chcemy znaleźć $u_{{hb}}\in V^{h}_{{bi}}$ takie, że

$\int _{\Omega}\nabla u_{{hb}}\nabla v_{h}+cu_{h}\, v_{h}\, dx=\int _{\Omega}fv_{h}\, dx\qquad\quad\forall v_{h}\in V^{h}_{{bi}}.$

(12.7)

W tym przypadku możemy analogicznie wykazać, że zachodzi oszacowanie błędu takie, jak w przypadku elementu liniowego na trójkącie. Tzn. można wykazać, że to zadanie ma jednoznaczne rozwiązanie, i że jeśli $u_{*}\in C^{2}([0,1]^{2})$ , to

	$\displaystyle\\|\nabla(u_{{hb}}-u_{*})\\| _{{L^{2}(\Omega)}}$	$\displaystyle=$	$\displaystyle O(h)$		(12.8)
	$\displaystyle\\| u_{{hb}}-u_{*}\\| _{{L^{2}\Omega)}}$	$\displaystyle=$	$\displaystyle O(h^{2}).$

W przypadku zadania na kwadracie jednostkowym element biliniowy wydaje się bardziej naturalny, ale przestrzeń elementu skończonego jest bardziej skomplikowana, ponieważ na każdym kwadracie triangulacji lokalna przestrzeń MESu zawiera wielomiany wyższego stopnia. Czy możemy uzyskać lepsze oszacowania błędu? Nie - jeśli chodzi o rząd zbieżności, tzn. dla elementu biliniowego otrzymamy co najwyżej $O(h)$ w normie $\|\nabla\cdot\| _{{L^{2}(\Omega)}}$ , ale oszacowanie błędu jest ostrzejsze, por. rozdział 4.6 w [4].

12.2. Niejednorodny warunek brzegowy

Można się zastanowić co się dzieje, jeśli w (7.8) w warunku brzegowym Dirichleta wartość prawej strony jest różna od zera, tzn. $u_{*}=g$ na brzegu. Tak jak poprzednio, otrzymujemy nowe słabe sformułowanie: chcemy znaleźć $u_{*}$ takie, że $u_{*}=g$ na brzegu $\Omega$ spełniające (12.1). Jeśli założymy, że znamy funkcję $\tilde{g}$ określoną na $\Omega$ taką, że $\tilde{g}=g$ na brzegu, to definiując $\hat{u}_{*}=u_{*}-\tilde{g}$ otrzymujemy:

$a(\hat{u}_{*},v)=(f,v)-a(\tilde{g},v)=F(v)\qquad\forall v\in V,$

dla $V$ jak w (12.1) i $\hat{u}^{*}\in V$ , czyli (12.1) ale z prawą stroną zależną dodatkowo od $\tilde{g}$ .

Następnie możemy wprowadzić zadanie dyskretne tak, jak dla zerowych warunków brzegowych.

12.3. Zadania

Ćwiczenie 12.1

Udowodnij, że funkcja bazowa przestrzeni elementu liniowego na trójkącie spełnia wzory (12.1.2).

Ćwiczenie 12.2

Definiujemy rozwiązanie (12.1) jako

$u_{h}=\sum _{{k,l=0}}^{{N-2}}u_{{k+(N-1)*l}}\;\phi _{{k+1,l+1}}$

z $u(x^{{(k+1,l+1)}})=u_{{k+(N-1)*l}}$ . Wstawiając $u_{h}$ w tej postaci do (12.1) otrzymujemy układ równań liniowych na wektor współrzędnych $\vec{u}=(u_{{k+(N-1)*l}})_{{k,l=0}}^{{N-2}}$ . Policz różne od zera elementy macierzy tego układu $A_{h}$ i elementy wektora prawej strony $f_{h}$ . Czy elementy na diagonali tej macierzy zależą od $h$ ? Pokaż, że dla równomiernego podziału, tzn. $x^{{(k)}}=k*h$ , macierz ta jest równa macierzy dyskretyzacji metodą różnic skończonych dla tego samego zadania, pomnożonej przez parametr $h^{2}$ , tzn. macierzą układu równań liniowych dla ćwiczenia 7.5 (dla $c=0$ ). Czy oba układy po przeskalowaniu przez $h$ są wtedy identyczne?

Rozwiązanie:

Układy równań liniowych nie są identyczne, ponieważ prawe strony są różne, tak jak w przypadku jednowymiarowym.

Ćwiczenie 12.3

Pokaż, że macierz $A_{h}$ z ćwiczenia 12.2 jest zawsze pasmowa (wyznacz wielkość pasma, wyznacz czy zależy ono od $N$ , czyli odpowiednio od $h$ ), symetryczna i dodatnio określona.

Ćwiczenie 12.4

Zakładając, że posiadamy procedurę z metodą iteracyjną, która rozwiązuje układ równań z ćwiczenia 12.2 w $I(N)$ iteracji, i dla której w każdej iteracji wykonujemy jedno mnożenie przez macierz $A_{h}$ oraz $10*N^{2}$ operacji algebraicznych określ, ile musi wynosić $I(N)$ , aby metoda ta była tańsza od odpowiedniej taśmowej wersji eliminacji Gaussa dla macierzy symetrycznej dodatnio określonej.

Ćwiczenie 12.5

Pokaż, że dla dowolnej funkcji przestrzeni liniowej elementu skończonego $u_{h}\in V^{h}$ , por. rozdział 12.1.2, zachodzi tzw. nierówność odwrotna:

$\|\nabla u_{h}\| _{{L^{2}(\Omega)}}\leq C\, h^{{-1}}\| u_{h}\| _{{L^{2}(\Omega)}}.$

Tutaj stała $C$ jest niezależna od parametru $h$ i $u_{h}$ .

Wskazówka:

Wystarczy pokazać to oszacowanie na dowolnym elemencie triangulacji $\tau\in T_{h}$ i wykorzystać fakt, że $\|\nabla u_{h}\| _{{L^{2}(\tau)}}=\|\nabla(u_{h}-c)\| _{{L^{2}(\tau)}}$ dla dowolnej stałej $c$ .

Ćwiczenie 12.6

(laboratoryjne) Napisz w octave odpowiednią wersję eliminacji Gaussa dla macierzy pasmowej symetrycznej dodatnio określonej. Zastosuj ją do układu równań liniowych z ćwiczenia 12.2 dla $N$ różnej wielkości. Porównaj czas w octavie z czasem dla standardowej metody rozwiązywania równań liniowych octave'a zarówno, gdy macierz układu jest w formacie pełnym, jak i formacie rzadkim (można użyć funkcje octave: sparse(), tic() i toc() i operator octave: \).

Ćwiczenie 12.7

Udowodnij, że funkcja z $P_{2}$ jest jednoznacznie wyznaczona przez określenie jej wartości w wierzchołkach i środkach krawędzi trójkąta (jak opisano w rozdziale 12.1.3, por. rysunek 12.6).

Ćwiczenie 12.8

Udowodnij, że funkcja z $P_{3}$ jest jednoznacznie wyznaczona przez określenie jej wartości jak opisano w rozdziale 12.1.3, por. rysunek 12.7.

Ćwiczenie 12.9

Dla przestrzeni $V^{h}_{2}$ wyprowadź układ równań liniowych:

$A_{{h,2}}\vec{u}=\vec{f}$

taki, że jego rozwiązaniem są współczynniki rozwiązania przybliżonego $u_{{h,p}}$ w bazie nodalnej tej przestrzeni elementu skończonego, por. rozdział 12.1.3. Czy macierz tego układu jest pasmowa? Jeśli tak, to znajdź wielkość pasma. Czy jest symetryczna i dodatnio określona?

Ćwiczenie 12.10

Udowodnij wzory na bazowe funkcje elementu kwadratowego na trójkącie (12.5) i (12.6).

Ćwiczenie 12.11

Wyznacz wzory na bazowe funkcje nodalne dla elementu kubicznego, analogiczne do wzorów (12.5) i (12.6) dla funkcji bazowych elementu kwadratowego, tzn. wzory na te funkcje w zależności od funkcji $q_{k}$ .

Ćwiczenie 12.12

Wyznacz wzory na współczynniki bazowych funkcji nodalnych dla elementu biliniowego na kwadracie $x^{{(kl)}}+(0,h)^{2}$ w bazie $(1,(x_{1}-x^{{(k)}}),(x_{2}-x^{{(l)}}),(x_{1}-x^{{(k)}})*(x_{2}-x^{{(l)}}))$ . Oblicz elementy macierzy układu równań liniowych dla dyskretyzacji metodą elementu skończonego biliniowego dla zadania dyskretnego (12.7) w tej bazie nodalnej.

Ćwiczenie 12.13

Dla przestrzeni $V^{h}_{{bi}}$ , wyprowadź układ równań liniowych:

$A_{{h,bi}}\vec{u}=\vec{f}$

taki, że jego rozwiązaniem są współczynniki rozwiązania przybliżonego $u_{{hb}}$ w bazie nodalnej związanej z wierzchołkami elementów kwadratowych. Czy macierz tego układu jest pasmowa? Jeśli tak, to znajdź wielkość pasma. Czy jest symetryczna i dodatnio określona?

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