Zagadnienia

11.1 Metoda elementu skończonego dla modelowego zadania eliptycznego w jednym wymiarze
11.2 Zadania

11. Metoda elementu skończonego - wprowadzenie

W tym rozdziale przedstawimy główne idee metody elementu skończonego na przykładzie modelowego zadania eliptycznego rzędu dwa na obszarze jednowymiarowym. Metoda elementu skończonego jest bardziej ogólna od metody różnic skończonych nawet dla zadań różniczkowych zadanych na obszarze w jednym wymiarze. Np. konstrukcje zadań przybliżonych dla warunków brzegowych różnego typu są dużo prostsze niż w przypadku metody różnic skończonych.

11.1. Metoda elementu skończonego dla modelowego zadania eliptycznego w jednym wymiarze

11.1.1. Słabe sformułowanie

Rozpatrzmy modelowe zadanie jednowymiarowe (7.1) z zerowymi warunkami brzegowymi i $c=0$ , którego rozwiązanie oznaczymy $u_{*}$ .

Następnie weźmy dowolną funkcję ciągłą $\phi$ , która jest kawałkami $C^{1}$ na odcinku, tzn. która ma ciągłą pochodną poza skończoną ilością punktów taką, że $\phi(a)=\phi(b)=0$ . Przemnóżmy równanie $Lu=f$ przez tę funkcję. Ze wzoru na całkowanie przez części otrzymujemy:

$\int _{a}^{b}-\frac{d^{2}u_{*}}{dx^{2}}\phi\; dx=\int _{a}^{b}\frac{du_{*}}{dx}\frac{d\phi}{dx}\; dx=\int _{a}^{b}f\phi\; dx.$

Oczywiście tutaj $\frac{d\phi}{dx}$ jest zdefiniowana poza skończoną ilością punktów nieciągłości.

Zamiast zadania (7.1) możemy rozpatrzyć zadanie znalezienia funkcji $u_{*}$ w odpowiedniej przestrzeni $V$ funkcji określonych na odcinku $[a,b]$ zawierających funkcje kawałkami $C^{1}$ i zerujące się w końcach odcinka (na razie nie ustalajmy precyzyjnie o jaką przestrzeń chodzi) takiej, żeby

$\int _{a}^{b}\frac{du_{*}}{dx}\frac{d\phi}{dx}\; dx=\int _{a}^{b}f\phi\; dx\quad\forall\phi\in V.$

(11.1)

Oczywiście rozwiązanie (7.1) spełnia (11.1), a przy odpowiednim doborze $V$ (oraz dostatecznej gładkości $f$ ) można pokazać, że rozwiązanie (11.1) również spełnia (7.1). Zadanie (11.1) nazywamy sformułowaniem uogólnionym (słabym, wariacyjnym) zadania (7.1). Zaś (7.1) nazywamy sformułowaniem klasycznym. Metoda elementu skończonego, która jest szczególnym przypadkiem metody Galerkina polega na tym, że wprowadzamy w specjalny sposób skończenie wymiarową podprzestrzeń przestrzeni $V$ . Następnie szukamy rozwiązania w tej przestrzeni dyskretnej, które spełnia zadanie wariacyjne (11.1) z tym, że przestrzeń $V$ jest zastąpiona przez naszą dyskretną podprzestrzeń.

Proszę zauważyć, że podejście wariacyjne jest inne od metody różnic skończonych, w której konstrukcja rozwiązania określonego na zbiorze dyskretnym (siatce) polega na zastąpieniu odpowiednich pochodnych w równaniu różniczkowym odpowiednimi ilorazami różnicowymi na tej siatce.

11.1.2. Element liniowy

Wprowadźmy podział (triangulacje) odcinka $[a,b]$ na pododcinki (elementy) $T_{h}([a,b])=\{\tau _{k}\}$ , gdzie $\tau _{k}=(x_{k},x_{{k+1}})$ dla $a=x_{0}<\ldots<x_{{N-1}}<x_{N}=b$ . Za parametr tego podziału przyjmijmy $h=\max _{k}|x_{k}-x_{{k-1}}|$ , a punkty $x_{k}$ nazwiemy punktami nodalnymi (węzłami) tego podziału.

Oczywiście najprostszym podziałem jest podział równomierny, jeśli bierzemy $x_{k}=k*h$ dla $h=(b-a)/N$ .

Zakładamy, że rozpatrujemy rodzinę podziałów z $h$ dążącym do zera.

Teraz na bazie danego podziału $T_{h}([a,b])$ możemy wprowadzić przestrzeń dyskretną:

$V^{h}=\{ u\in C([a,b]):u_{{|\tau}}\in P_{1}\quad\forall\tau\in T_{h}([a,b]),\quad u(a)=u(b)=0\},$

gdzie $P_{1}=span(1,x)$ jest przestrzenią wielomianów liniowych, tzn. stopnia nie większego niż jeden. Funkcje z tej przestrzeni to funkcje kawałkami liniowe (czyli klasyczne splajny liniowe) z zerowymi warunkami brzegowymi na brzegu, czyli w szczególności są to funkcje ciągłe kawałkami klasy $C^{1}$ , zatem $V^{h}\subset V$ .

Zadanie dyskretne polega na znalezieniu $u_{h}\in V^{h}$ takiego, że

$\int _{a}^{b}\frac{du_{h}}{dx}\frac{d\phi}{dx}\; dx=\int _{a}^{b}f\phi\; dx\quad\forall\phi\in V^{h}.$

(11.2)

Pozostawimy jako zadanie wykazanie istnienia jednoznacznego rozwiązania (11.2).

$\$

Rys. 11.1. Baza daszkowa dla $N=5$ .

Na razie wprowadźmy tzw. funkcję nodalną: niech $\phi _{k}\in V^{h}$ dla $k\in\{ 1,\ldots,N-1\}$ będzie taką funkcją, że $\phi _{k}(x_{k})=1$ i $\phi _{k}(x_{j})=0$ dla $j\not=k$ . Możemy podać wzór na taką funkcję:

$\phi(x)=\left\{\begin{array}[]{lcl}0&&x\not\in[x_{{k-1}},x_{{k+1}}],\\ \frac{x-x_{{k-1}}}{x_{k}-x_{{k-1}}}&&x\in[x_{{k-1}},x_{k}],\\ 1-\frac{x-x_{k}}{x_{{k+1}}-x_{k}}&&x\in[x_{k},x_{{k+1}}].\end{array}\right.$

(11.3)

Widzimy wykres kilku takich funkcji, por. rysunek 11.1.

Nietrudno pokazać, że $\phi _{k}$ jest elementem $V^{h}$ , i że $(\phi _{1},\ldots,\phi _{{N-1}})$ tworzy bazę $V^{h}$ taką, że jeśli $u\in V^{h}$ , to

$u=\sum _{{k=1}}^{{N-1}}u(x_{k})\phi _{k},$

gdzie $u(x_{k})$ wartość funkcji $u$ w punkcie nodalnym $x_{k}$ . Wstawiając $u_{h}=\sum _{{k=1}}^{{N-1}}u_{h}(x_{k})\phi _{k}$ do (11.2) otrzymujemy następujący układ równań liniowych:

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

(11.4)

z $\vec{u}=(u(x_{1}),\ldots,u(x_{{N-1}}))^{T}$ ,

$A_{h}=(a_{{kl}})_{{k,l=1}}^{{N-1}},\qquad a_{{kl}}=\int _{a}^{b}\frac{d\phi _{k}}{dx}\frac{d\phi _{l}}{dx}\; dx,\qquad\vec{f}=(f_{1},\ldots,f_{{N-1}})^{T}$

dla $f_{k}=\int _{a}^{b}f\phi _{k}\, dx$ .

Zauważmy, że macierz $A_{h}$ jest symetryczna i trójdiagonalna. Można wykazać, że jest dodatnio określona, więc można powyższy układ rozwiązać metodą przeganiania lub odpowiednim wariantem metody Choleskiego kosztem rzędu $c*N$ dla stałej $c$ niezależnej od $N$ .

11.1.3. Zbieżność

Zastanówmy się nad zbieżnością $u_{h}$ do rozwiązania $u_{*}$ . Najpierw trzeba ustalić w jakiej normie chcemy wykazać zbieżność.

Naturalną normą jest norma energetyczna związana z formą dwuliniową w słabym sformułowaniu (11.1):

$\| u\| _{V}:=\sqrt{\int _{a}^{b}\left|\frac{du}{dx}\right|^{2}\, dx}.$

Wprowadzając oznaczenie $e_{h}=u_{*}-u_{h}$ i odejmując (11.2) od (11.1) otrzymujemy:

$\displaystyle\int _{a}^{b}\frac{de_{h}}{dx}\frac{d\phi}{dx}\, dx=\int _{a}^{b}\frac{d}{dx}(u_{*}-u_{h})\frac{d\phi}{dx}\, dx=0\qquad\forall\phi\in V^{h}.$

Następnie dla dowolnego $v_{h}\in V^{h}$ widzimy, że

	$\displaystyle\int _{a}^{b}\left\|\frac{de_{h}}{dx}\right\|^{2}\, dx$	$\displaystyle=$	$\displaystyle\int _{a}^{b}\frac{de_{h}}{dx}\frac{du_{}}{dx}\, dx=\int _{a}^{b}\frac{de_{h}}{dx}\frac{du_{}}{dx}\, dx-\int _{a}^{b}\frac{de_{h}}{dx}\frac{dv_{h}}{dx}\, dx$
		$\displaystyle=$	$\displaystyle\int _{a}^{b}\frac{de_{h}}{dx}\frac{d}{dx}(u_{*}-v_{h})\, dx.$

Korzystając z nierówności Schwarza w $L^{2}(a,b)$ otrzymujemy:

$\| e_{h}\| _{V}^{2}\leq\| e_{h}\| _{V}\| u_{*}-v_{h}\| _{V},$

czyli

$\| e_{h}\| _{V}\leq\| u_{*}-v_{h}\| _{V}\qquad\forall v_{h}\in V^{h}.$

Przyjmując za $v_{h}$ liniowy interpolator $u_{*}$ w punktach nodalnych, tzn. $I_{h}u_{*}:=\sum _{k}u_{*}(x_{k})\phi _{k}$ , można wykazać:

$\left\|\frac{d^{k}}{dx^{k}}(u_{*}-I_{h}u_{*})\right\| _{{\infty,[a,b]}}\leq C_{k}h^{{2-k}}\left\|\frac{d^{2}u_{*}}{dx^{2}}\right\| _{{\infty,[a,b]}}\qquad k=0,1$

(11.5)

dla pewnych stałych $C_{k}$ niezależnych od $h$ , $u_{*}$ i $[a,b]$ . Dowód tego oszacowania pozostawiamy jako zadanie, wynika on wprost z oszacowań błędu interpolacji dla splajnów liniowych.

Wtedy od razu otrzymujemy:

$\| e_{h}\| _{V}\leq\| u_{*}-I_{h}u_{*}\| _{V}\leq C_{1}h\left\|\frac{d^{2}u_{*}}{dx^{2}}\right\| _{{\infty,[a,b]}},$

czyli dla funkcji klasy $C^{2}$ błąd w normie energetycznej zachowuje się jak $O(h)$ .

Można też wykazać, że w normie $L^{2}(a,b)$ zachodzi:

$\| e_{h}\| _{{L^{2}(a,b)}}\leq Ch^{2}\left\|\frac{d^{2}u_{*}}{dx^{2}}\right\| _{{\infty,[a,b]}},$

czyli błąd zachowuje się jako $O(h^{2})$ - co nie jest oczywiste.

Porównajmy to oszacowanie z oszacowaniem błędu z metody różnic skończonych (MRS) na siatce równomiernej dla tego samego zadania różniczkowego. Można pokazać wtedy zbieżność dyskretną $O(h^{2})$ w dyskretnej normie $L^{2}_{h}$ , por. (10.4), która w przybliżeniu odpowiada normie $L^{2}$ , czyli możemy powiedzieć, że szybkość zbieżności w tym przypadku metody elementu skończonego i metody różnic skończonych jest tego samego rzędu. Ale w MRS musieliśmy założyć równomierność siatki i wyższą gładkość rozwiązania ( $u_{*}\in C^{4}((a,b))$ ).

11.1.4. Inne przestrzenie elementu skończonego

Przestrzeń dyskretną $V^{h}$ można też zdefiniować inaczej.

Dla danego podziału $T_{h}([a,b])$ zdefiniujemy następujące przestrzenie elementu skończonego dla dowolnego $p=1,2,3,4,\ldots$ :

$V^{h}_{p}=\{ u\in C([a,b]):u_{{|\tau}}\in P_{p}\quad\forall\tau\in T_{h}([a,b]),\quad u(a)=u(b)=0\}$

gdzie $P_{p}$ jest przestrzenią wielomianów stopnia nie przekraczającego $p$ .

Widzimy, że $V^{h}=V^{h}_{1}$ . Przestrzeń $V^{h}_{2}$ nazywamy przestrzenią elementu skończonego funkcji ciągłych kawałkami kwadratowych, a przestrzeń $V^{h}_{3}$ - przestrzenią elementu skończonego funkcji ciągłych kawałkami kubicznych, czy inaczej - metodą elementu skończonego typu Lagrange'a kwadratową lub kubiczną. Możemy teraz postawić zadanie dyskretne, jak poprzednio, tzn. szukamy $u_{{h,p}}\in V^{h}_{p}$ takiego, że

$\int _{a}^{b}\frac{du_{{h,p}}}{dx}\frac{d\phi}{dx}\; dx=\int _{a}^{b}f\phi\; dx\quad\forall\phi\in V^{h}_{p}.$

(11.6)

Zadanie to ma jednoznaczne rozwiązanie. Analogicznie jak dla elementu liniowego możemy wprowadzić tu tzw. bazę nodalną w $V^{h}_{p}$ . Wprowadzamy dodatkowe punkty wewnątrz odcinka $[x_{k},x_{{k+1}})$ dla $k=0,\ldots,N-1$ :

$x_{{k,j}}=x_{k}+\frac{j}{p}(x_{{k+1}}-x_{k}),\quad j=0,\ldots,p-1.$

Oczywiście $x_{{k,0}}=x_{k}$ .

Dla każdego punktu $x_{{k,j}}$ oprócz $x_{{0,0}}=x_{0}=a$ wprowadzamy funkcję bazową $\phi _{{k,j}}\in V^{h}_{p}$ taką, że

$\phi _{{k,j}}(x_{{l,i}})=\left\{\begin{array}[]{lcl}1&&l=j\quad i=j\\ 0&&\mathrm{w\; przeciwnym\; przypadku}\end{array}\right.$

(11.7)

Można pokazać, że $(\phi _{{k,j}})_{{k,j}}$ jest bazą i to taką, że

$\mathrm{supp}(\phi _{{k,j}})=\left\{\begin{array}[]{lcl}{[x_{{k-1}},x_{{k+1}}]}&&j=0,k>0\\ {[x_{k},x_{{k+1}}]}&&j\not=0.\end{array}\right.$

Powstaje pytanie: po co stosować przestrzenie $V^{h}_{p}$ dla $p>1$ ?

Jeśli rozwiązanie $u_{*}$ jest bardziej regularne, tzn. należy do $C^{{p+1}}([a,b]),$ to można wykazać, że dla $I_{{h,p}}=\sum _{{k,j}}u(x_{{k,j}})\phi _{{k,j}}\in V^{h}_{p}$ :

$\left\|\frac{d^{k}}{dx^{k}}(u_{*}-I_{{h,p}}u_{*})\right\| _{{\infty,[a,b]}}\leq C_{k}h^{{p+1-k}}\left\|\frac{d^{{p+1}}u_{*}}{dx^{{p+1}}}\right\| _{{\infty,[a,b]}}\qquad k=0,1$

(11.8)

i stąd, jak poprzednio dla elementu liniowego, otrzymujemy, że

$\| u_{*}-u_{{h,p}}\| _{V}\leq C_{p}h^{p}\left\|\frac{d^{{p+1}}u_{*}}{dx^{{p+1}}}\right\| _{{\infty,[a,b]}},$

(11.9)

czyli błąd zachowuje się jak $O(h^{p})$ co oznacza, że w tej normie zachodzi zbieżność rzędu $p$ .

11.2. Zadania

Ćwiczenie 11.1

Pokaż, że $\phi _{k}$ zdefiniowana w (11.3) jest w $V^{h}$ , i że $(\phi _{1},\ldots,\phi _{{N-1}})$ tworzy bazę $V^{h}$ .

Ćwiczenie 11.2

Wykaż (11.4). Policz wszystkie różne od zera elementy macierzy $A_{h}$ układu (11.4). Pokaż, że dla $c=0$ i równomiernego podziału odcinka tzn. $x_{k}=a+k*h$ macierz ta jest równa macierzy dyskretyzacji metodą różnic skończonych dla tego samego zadania, pomnożonej przez parametr $h$ , tzn. jest macierzą układu równań liniowych (7.7). Czy oba układy równań liniowych po przeskalowaniu przez $h$ (7.7) są wtedy identyczne?

Rozwiązanie:

Układy nie są identyczne, ponieważ prawe strony mogą być różne, jakkolwiek wtedy prawą stronę (7.7) możemy uznać za prostą aproksymację całek z prawej strony (11.4).

Ćwiczenie 11.3

Pokaż, że macierz $A_{h}$ w (11.4) jest zawsze trójdiagonalna, symetryczna i dodatnio określona.

Wskazówka:

Pokaż, że dla dowolnych funkcji $u=\sum _{k}u_{k}\phi _{k},v=\sum _{k}v_{k}\phi _{k}\in V^{h}$ zachodzi $\vec{u}^{T}A_{h}\vec{v}=\int _{a}^{b}\frac{du}{dx}\frac{dv}{dx}\, dx$ , dla $\vec{u}=(u_{1},\ldots,u_{{N-1}})^{T},\vec{v}=(v_{1},\ldots,v_{{N-1}})^{T}$ .

Ćwiczenie 11.4

Zaproponuj metodę rozwiązywania układu równań (11.4) będącą odpowiednią wersją metody eliminacji Gaussa dla macierzy symetrycznej trójdiagonalnej dodatnio określonej, której koszt wynosi $O(N)$ .

Ćwiczenie 11.5 (laboratoryjne)

Dla podziału równomiernego na odcinku $[-1,1]$ rozwiąż w octavie (11.4) dla znanego rozwiązania $u(x)=\sin(\pi*x)$ , czyli dla $f=-\pi^{2}*\sin(x)$ . Prawą stronę możemy policzyć odpowiednią funkcją octave'a. Policz rozwiązania dyskretne dla $2h$ i $h$ . Następnie policz normy dyskretne maksimum dla błędów $u_{*}-u_{h}$ i $u_{*}-u_{{2h}}$ (czyli maksima błędów w punktach nodalnych) i ich stosunek. (Zakładając, że błąd dla $h$ wynosi $O(h^{p})$ , stosunek ten powinien w przybliżeniu wynosić $2^{p}$ ).

Ćwiczenie 11.6 (częściowo laboratoryjne)

Udowodnij, że $\phi _{{k,j}}$ z (11.7) są w $V^{h}_{p}$ , i że stanowią bazę tej przestrzeni. Wyprowadź bezpośrednie wzory na $\phi _{{k,j}}$ . Narysuj w octave wykresy wszystkich $\phi _{{k,j}}$ na odcinku $[x_{k},x_{{k+1}}]$ dla $[a,b]=[0,2]$ , $h=1$ i $p=2,3$ przy użyciu funkcji octave'a plot().

Ćwiczenie 11.7

Korzystając z (11.8) wykaż (11.9) tzn., że

$\| u_{*}-u_{{h,p}}\| _{V}=O(h^{p})$

o ile $u_{*}$ jest dostatecznie gładka.

Wskazówka:

Dowód przebiega identycznie jak w przypadku $p=1$ , tzn. liniowego elementu skończonego.

Ćwiczenie 11.8

Udowodnij jednowymiarową nierówność Friedrichsa, a mianowicie, że jeśli $f$ jest funkcją ciągłą kawałkami klasy $C^{1}$ na $[a,b]$ i $u(a)=0$ , to

$\int _{a}^{b}|f|^{2}\, dx\leq(b-a)^{2}\int _{a}^{b}\left|\frac{df}{dx}\right|^{2}\, dx.$

Ćwiczenie 11.9

Pokaż, że $(u,v)_{V}=\int _{a}^{b}\frac{du}{dx}\frac{dv}{dx}\, dx$ jest iloczynem skalarnym na $V^{h}_{p}$ i - ogólniej na dowolnej przestrzeni funkcyjnej zawartej w przestrzeni funkcji ciągłych kawałkami $C^{1}$ zerujących się w końcach odcinka $[a,b]$ . W szczególności $\| u\| _{V}$ jest normą na $V^{h}_{p}$ .

Ćwiczenie 11.10

Wyprowadź układ równań liniowych

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

którego rozwiązaniem jest wektor współczynników w bazie $\{\phi _{{k,j}}\}$ , por. (11.7), rozwiązania $u_{{h,p}}$ zadania (11.6). Określ ilość elementów różnych od zera w macierzy, czy przy odpowiednim porządku indeksów funkcji bazy ta macierz może być jest pasmowa? Jeśli tak, to znajdź wielkość pasma. Czy jest symetryczna i dodatnio określona? Zaproponuj algorytm bezpośredni rozwiązywania tego układu kosztem $O(n)$ , dla $n$ wymiaru $V^{h}_{p}$ .

Ćwiczenie 11.11

Udowodnij (11.5).

Wskazówka:

Oszacowanie dla $k=0$ wprost wynika z oszacowania błędu interpolacji Lagrange'a, por. np. [14], [18] lub w języku angielskim [17] lub [25]. Natomiast w przypadku $k=1$ wystarczy zauważyć, że w każdym przedziale $(x_{k},x_{{k+1}})$ dla $0\leq k<N$ istnieje punkt $\xi _{k}$ taki, że $u_{*}(\xi _{k})=I_{h}u_{*}(\xi _{k})$ , a następnie skorzystać z twierdzenia o wartości średniej.

Ćwiczenie 11.12

Udowodnij (11.8).

Wskazówka:

Ponownie jak w poprzednim zadaniu oszacowanie dla $k=0$ wprost wynika z oszacowania błędu interpolacji Lagrange'a, por. np. [14], [18] lub w języku angielskim [17] lub [25].

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

Numeryczne rozwiązywanie równań różniczkowych

Zagadnienia

11. Metoda elementu skończonego - wprowadzenie

11.1. Metoda elementu skończonego dla modelowego zadania eliptycznego w jednym wymiarze

11.1.1. Słabe sformułowanie

11.1.2. Element liniowy

11.1.3. Zbieżność

11.1.4. Inne przestrzenie elementu skończonego

11.2. Zadania

Ćwiczenie 11.1

Ćwiczenie 11.2

Ćwiczenie 11.3

Ćwiczenie 11.4

Ćwiczenie 11.5 (laboratoryjne)

Ćwiczenie 11.6 (częściowo laboratoryjne)

Ćwiczenie 11.7

Ćwiczenie 11.8

Ćwiczenie 11.9

Ćwiczenie 11.10

Ćwiczenie 11.11

Ćwiczenie 11.12