Zagadnienia

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.

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.