Zagadnienia

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*V takie, że

Ωu*v+cuvdx=ΩfvdxvV. (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 L2Ω. 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 τ na kwadracie Ω.

12.1.1. Triangulacja obszaru

Wprowadźmy podział (triangulacje) kwadratu Ω¯=0,12 na trójkąty Th0,12=τkk o jednakowym kształcie i wielkości. Najprościej jest podzielić kwadrat na równe kwadraty ([xk,xk+1]×[xl,xl+1] dla k,l=0,,N-1 i xk=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 τThΩτ¯=Ω¯ i τkτ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 xkl=xk,xl 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 h0.

12.1.2. Element liniowy

Możemy wprowadzić przestrzeń dyskretną:

Vh=uCΩ¯:u|τP1ττThΩ,us=0sΩ,

dla P1τ=span1,x1,x2 przestrzeni wielomianów liniowych na τ. 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 τ - punkty nodalne wielomianu liniowego na tym trójkącie.

Analogicznie jak w przypadku jednowymiarowym, możemy wprowadzić funkcje nodalne: definiujemy ϕk,lVh dla k,l1,,N-1 jako funkcję, która spełnia ϕk,lxkl=1 i ϕk,lxij=0 dla xklxij dla xkl=xk,xl.

\
Rys. 12.3. Wykres funkcji daszkowej ϕ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 τ. Powiedzmy, że wzory wyznaczymy na tym trójkącie, który jest zaznaczony na rysunku 12.1 - przyjmujemy, że xkl=k*h,l*h jest wierzchołkiem przy kącie prostym:

ϕk,lx=1-x1-k*hh-x2-l*hh,
ϕk+1,lx=x1-k*hh, (12.2)
ϕk,l+1x=x2-l*hh.

dla x=x1,x2. Wykresy takiej funkcji dla tej regularnej triangulacji kwadratu lub kilku funkcji możemy obejrzeć na rysunkach 12.312.4. Na rysunku 12.5 widzimy przykładową funkcję z przestrzeni Vh. Również jako zadanie pozostawimy wykazanie, że

ϕk,lk,l=1,,N-1

tworzą bazę Vh, i że

u=k,luxklϕk,l.

Wprowadzamy zadanie dyskretne: chcemy znaleźć uhVh takie, że

Ωuhvh+cuhvhdx=ΩfvhdxvhVh. (12.3)

Można pokazać, że to zadanie ma jednoznaczne rozwiązanie, i że jeśli u*C20,12 to:

uh-u*L2Ω=Oh (12.4)
uh-u*L2Ω)=Oh2.

Zatem w normie L2 widzimy oszacowanie zbieżności rzędu dwa analogiczne jak dla dyskretyzacji tego samego zadania przy pomocy metody różnic dzielonych i dyskretnej normy Lh2, ale przy dużo słabszych założeniach. Tam musieliśmy założyć, że funkcja jest klasy C4.

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 Th0,12 jak w Rozdziale 12.1.1. Wtedy definiujemy:

Vph=uCΩ¯:u|τPpττThΩ,us=0sΩ.

dla p=1,2,3,. Tutaj Ppτ=spanx1kx2l0k+lp to przestrzeń wielomianów na τ stopnia nie większego od p.

Przypadek V2h określamy jako przestrzeń elementu skończonego kwadratowego, a V3h 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 Ω określamy jako punkty nodalne V2h. 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 τ -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 τ.

W przypadku przestrzeni Vph 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 τ. Niech punkty xk k=0,1,2 będą trzema wierzchołkami tego trójkąta τ i niech punkt mk k=0,1,2 będzie środkiem odcinka między xk a xk+1. Przy czym utożsamiamy 0 z 3. Niech qk będzie funkcją liniową taką, że qkxj=0 dla jk i qkxk=1. Wtedy funkcja nodalna związana z wierzchołkiem elementu xk dla elementu kwadratowego wynosi na trójkącie τ:

ϕxk=qkqk-qk-1-qk+1k=0,1,2mod 3, (12.5)

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

ϕmk=4qkqk+1k=0,1,2mod 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 Cs+1Ω¯ dla s większych od jeden do s=p. Tzn. wtedy można wykazać, że

uh-u*L2Ω=Ohs.

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 Ω=0,12 wprowadźmy triangulacje złożoną z kwadratów, por. rysunek 12.9:

ThΩ=τklk,l=0,,N-1

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

τkl=k*h,k+1*h×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 τ -punkty swobody wielomianu biliniowego na tym kwadracie.

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

Vbih=uCΩ¯:u|τQ1ττThΩ,us=0sΩ.

gdzie

Q1τkl=P1k*h,k+1*hP1l*h,l+1*h=Span1,x1,x2,x1x2

przestrzeń wielomianów biliniowych określonych na τ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 ϕk,lVbih dla k,l1,,N-1 jako funkcję, która spełnia ϕkxkl=1 i ϕkxij=0 dla xklxij dla xkl=xk,xl. 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źć uhbVbih takie, że

Ωuhbvh+cuhvhdx=ΩfvhdxvhVbih. (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*C20,12, to

uhb-u*L2Ω=Oh (12.8)
uhb-u*L2Ω)=Oh2.

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 Oh w normie L2Ω, 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 Ω spełniające (12.1). Jeśli założymy, że znamy funkcję g~ określoną na Ω taką, że g~=g na brzegu, to definiując u*=u*-g~ otrzymujemy:

au*,v=f,v-ag~,v=FvvV,

dla V jak w (12.1) i u*V, czyli (12.1) ale z prawą stroną zależną dodatkowo od 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

uh=k,l=0N-2uk+N-1*lϕk+1,l+1

z uxk+1,l+1=uk+N-1*l. Wstawiając uh w tej postaci do (12.1) otrzymujemy układ równań liniowych na wektor współrzędnych u=uk+N-1*lk,l=0N-2. Policz różne od zera elementy macierzy tego układu Ah i elementy wektora prawej strony fh. Czy elementy na diagonali tej macierzy zależą od h? Pokaż, że dla równomiernego podziału, tzn. xk=k*h, macierz ta jest równa macierzy dyskretyzacji metodą różnic skończonych dla tego samego zadania, pomnożonej przez parametr h2, 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 Ah 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 IN iteracji, i dla której w każdej iteracji wykonujemy jedno mnożenie przez macierz Ah oraz 10*N2 operacji algebraicznych określ, ile musi wynosić IN, 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 uhVh, por. rozdział 12.1.2, zachodzi tzw. nierówność odwrotna:

uhL2ΩCh-1uhL2Ω.

Tutaj stała C jest niezależna od parametru h i uh.

Wskazówka: 

Wystarczy pokazać to oszacowanie na dowolnym elemencie triangulacji τTh i wykorzystać fakt, że uhL2τ=uh-cL2τ 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 P2 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 P3 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 V2h wyprowadź układ równań liniowych:

Ah,2u=f

taki, że jego rozwiązaniem są współczynniki rozwiązania przybliżonego uh,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 qk.

Ćwiczenie 12.12

Wyznacz wzory na współczynniki bazowych funkcji nodalnych dla elementu biliniowego na kwadracie xkl+0,h2 w bazie 1,x1-xk,x2-xl,x1-xk*x2-xl. 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 Vbih, wyprowadź układ równań liniowych:

Ah,biu=f

taki, że jego rozwiązaniem są współczynniki rozwiązania przybliżonego uhb 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.

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.