1 laboratorium rozwiązanie (bardzo) dokładne MRS: gęsta siatka u
2 u rozwiązanie metodą elementów skończonych z liniowymi funkcjami kształtu na dziewięciu węzłach rozwiązanie MES w tej wersji (liniowe fcje kształtu 1D) jest dokładne w węzłach odchylenia: tylko między węzłami
3 dla porównania: MRS na 9 – ciu węzłach MES uwaga – na rysunku - dla MRS punkty łączymy liniami tylko dla ilustracji dla MES z liniowymi funkcjami kształtu połączenie punktów: ma znaczenie dosłowne
4 u odstępstwo między wynikiem MES a dokładnym
5 u rozwiązanie MES w naszej bazie jest odcinkami liniowe a rozwiązanie dokładne jest liniowe tam gdzie =0 wniosek: tam gdzie =0 wystarczy jeden element!
6 u rozwiązanie MES w naszej bazie jest odcinkami liniowe a rozwiązanie dokładne jest liniowe tam gdzie =0 wniosek: tam gdzie =0 wystarczy jeden element! pomysł: przesunąć wszystkie węzły poza brzegowymi do obszaru gdzie nie znika gęstość ładunku – tam gdzie u zaokrąglone. wiemy już, że przesuwanie czerwonych punktów pójdzie po krzywej dokładnej.
7 Kryterium wyboru węzłów? (bx) i=2,8 x1= -x9 = -1zacieśniamy węzły wokół x=0 W przypadku ogólnym zawsze można policzyć pozostałość Lu-f=r i badać np. (r,r).
8 Kryterium wyboru węzłów? (bx) przy okazji dyskusji metod relaksacyjnych dowiedzieliśmy się, że najbliższe dokładnego jest rozwiązanie, które minimalizuje funkcjonał całki działania wykorzystajmy działanie jako kryterium jakości rozwiązania w metodzie elementów skończonych i=2,8 x1= -x9 = -1zacieśniamy węzły wokół x=0 W przypadku ogólnym zawsze można policzyć pozostałość Lu-f=r i badać np. (r,r). dla równania Poissona mamy znacznie lepsze kryterium:
9 do oceny jakości wyboru węzłów użyjemy macierzy A i wektora F, które uprzednio i tak musimy wyznaczyć aby wyliczyć rozwiązanie c. jednorodny warunek brzegowy Dirichleta: c 1 =c 9 =0 rozpoznajemy:
10 0.000.250.500.751.00 bx -0.0120 -0.0118 -0.0116 -0.0114 -0.0112 a 100001500020000250003000035000 równoodleg³e 201 pkt działanie w relaksacji MRS działanie w MES w funkcji parametru bx numer iteracji MRS funkcjonał działania a wybór położeń węzłów:
11 rozwiązanie dla optymalnego bx: czerwone: MES niebieskie: wynik dokładny czarne: niejednorodność równania
12 Wybór węzłów: przez optymalizację funkcjonału... metoda Galerkina, w tym w wersji elementów skończonych ma charakter wariacyjny metoda Galerkina dla dowolnej bazy jest równoważna metodzie Reyleigha-Ritza gdy ta stosowalna metoda Reyleigha-Ritza: rozwiązanie w bazie funkcyjnej dla ustalonych funkcji bazowych (w naszym przykładzie: dla ustalonych węzłów) c wyznaczone przez warunek minimum a. (pokazać, że warunek min a produkuje Ac=F)
13 Wybór węzłów: przez optymalizację funkcjonału... metoda elementów skończonych ma charakter wariacyjny w ogóle: metoda Galerkina dla dowolnej bazy jest równoważna metodzie Reyleigha-Ritza gdy ta stosowalna metoda Reyleigha-Ritza: rozwiązanie w bazie funkcyjnej dla ustalonych funkcji bazowych (w naszym przykładzie: dla ustalonych węzłów) c wyznaczone przez warunek minimum a. (pokazać, że warunek min a produkuje Ac=F) uwaga: w naszym przykładzie : dodatkowo optymalizowaliśmy funkcje bazowe (położenie węzłów). Zasada wariacyjna (najmniejszego działania) wykorzystana została więc dwukrotnie. jeśli tylko znamy funkcjonał dla równania różniczkowego: przyda się do optymalizacji kształtu elementów (2D)
14 Równanie Poissona, funkcje kształtu liniowe wynik MES dokładny w węzłach w MRS wszystko co możemy otrzymać, to wartości w węzłach, które są dokładne TYLKO w granicy x 0 !!! Wniosek: cały rachunek na 2 elementach: +1q f(q) niezależnie od wyboru q: f(q) da dokładne rozwiązanie równania wracamy do zagadnienia ścisłości rozwiązania MES w węzłach
15 q = 0.1 q=-0.6 q możemy ustawić gdziekolwiek, zawsze dostaniemy rozwiązanie dokładne wniosek: wystarczy przeskanować q przez pudło aby uzyskać dokładny wynik
16 +1q f(q) jeden węzeł wewnątrz pudłą funkcja kształtu v 1 = (x+1)/(q+1) dla xq rozwiązanie przybliżone: v=f(q)v 1 (x) Lu= F(v)=½(v,Lv)-( ,v) funkcjonał: wyliczymy F jako funkcję f(q) z warunku min F(v) wyznaczymy f(q) wiemy że MES wyprodukuje takie f(q) aby F było minimalne
17 v 1 = (x+1)/(q+1) dla xq v=f(q)v 1 (x) F(v)=½(v,Lv)-( ,v) v znika na brzegach:
18 v 1 = (x+1)/(q+1) dla xq v=f(q)v 1 (x) F(v)=½(v,Lv)-( ,v) v znika na brzegach:
19 =0
20 policzmy pochodne f(q) po q
21 =0 policzmy pochodne f(q) po q
22 =0 policzmy pochodne f(q) po q spełnia silną formę równania ! u’’= - stąd wynik dokładny dla u w x=q
23 silna a słaba forma równania: różnica u(x)’’= - x +1q v(q) v’(q) v’’(q) funkcja v nie spełnia silnej formy równania różniczkowego, tylko słabą: (v,v’’)f(q)=-(v, ) druga pochodna potencjału = ładunek druga pochodna potencjału: delta Diraca, niezależnie od tego jak wygląda niejednorodność równania Poissona
24 Podobny zabieg dla MRS: 3 węzły. Metoda różnic skończonych, siatka nierównomierna Iloraz różnicowy drugiej pochodnej dla nierównej siatki: ll pp + tracimy jeden rząd dokładności w porównaniu z siatką równomierną Problem rozwiązany w metodzie elementów skończonych. Wzór trójpunktowy W MES: nie ma problemu bo pochodne i całki liczymy dokładnie!
25 p=1-x l=x+1 +1x v(x) Co się stanie jeśli taki zabieg powtórzymy w MRS?
26 p=1-q l=q+1 +1q v(q) u(q) rozwiązanie dokładne Co się stanie jeśli taki zabieg powtórzymy w MRS?
27 liniowe fcje kształtu a warunki Neumanna tylko v 1 oraz v 2 wnoszą przyczynek do pochodnej na lewym końcu: pierwszy wiersz macierzy S (-1/h 2 1/h 2 0 0 0.... ) prawa strona pierwszy wiersz F : C warunki mieszane [Robina] u’(x 1 )+Du(x 1 )=E: (-1/h 2 +D 1/h 2 0 0 0... )
28 wybrane narzędzia MES umożliwiające jej automatyzację w więcej niż 1D: 1) macierze sztywności pojedynczych elementów oraz ich 2) składanie do globalnej macierzy sztywności 3) przestrzeń odniesienia i jej mapowanie do przestrzeni fizycznej
29 Przestrzeń referencyjna [odniesienia] Problem fizycznie zadany jest na siatce [x 1,x 2,x 3,...x N ] Rachunki (całkowanie elementów macierzowych) dla każdego elementu chcemy przenieść do przedziału (-1,1) w 1D x= -1 1 y= -1 y=1 element w przestrzeni fizycznej element w przestrzeni odniesienia x 1 x 2 x 3 x 4 x 5 -1 1 itd przestrzeń fizyczna przestrzeń odniesienia
30 Przestrzeń referencyjna [odniesienia] Problem fizycznie zadany jest na siatce [x 1,x 2,x 3,...x N ] Rachunki (całkowanie elementów macierzowych) dla każdego elementu chcemy przenieść do przedziału (-1,1) Element K m =(x m-1, x m ) (-1,1) mapowanie z (-1,1) do K m : x (x m +x m-1 )/2+(x m -x m-1 )/2 gdzie z przedziału (-1,1) Modelowy operator w 1D x= -1 1 y= -1 y=1 będziemy całkować jego elementy macierzowe w przestrzeni odniesienia element w przestrzeni fizycznej element w przestrzeni odniesienia
31 x( ) (x m-1 +x m )/2+(x m -x m-1 )/2 skala transformacji m-tego elementu: (czynnik skali, jakobian) przy transformacji: granice całki zmieniają się na –1,1, poza tym dx=J m d transformacja pochodnych: Całkowanie macierzy sztywności w przestrzeni referencyjnej 1D: J nie zależy od w 2D: zobaczymy, że nie zawsze tak jest [gdy element zmienia swój kształt w mapowaniu. w 1D: odcinek -> odcinek] całkę i pochodne przenosimy do przestrzeni odniesienia: element macierzowy całkowany w elemencie [fizycznym] J m =(x m -x m-1 )/2 pole elementu fizycznego / pole elementu odniesienia
32 x( ) (x m-1 +x m )/2+(x m -x m-1 )/2 całkowanie wektora sztywności: całka (f,v j ) transformuje się jak wyraz z a 0. Całkowanie macierzy sztywności w przestrzeni referencyjnej
33 odcinkowo liniowe funkcje kształtu w przestrzeni odniesienia x( ) (x i +x i+1 )/2+(x i+1 - x i )/2 v i (x( )) =1/2 –1/2 v i+1 (x( )) =1/2 W elemencie i+1 dwie funkcje kształtu x v i (x) 1 fcja kształtu xixi x i+1 x i-1 K i+1 KiKi fizyczna x i+2 v i+1 (x) K i+1 -1 1 v i v i+1 odniesienia
34 Przykład całkowanie w przestrzeni odniesienia dla bazy odcinkami liniowej (całka po elemencie K i+1 ) v i (x( )) =1/2 –1/2 v i+1 (x( )) =1/2 J i+1 =(x i+1 -x i )/2 pole elementu fizycznego / pole elementu odniesienia tu prim to pochodna po x a tu po ten wynik już znamy z całkowania w przestrzeni fizycznej
35 Macierz sztywności pojedynczego elementu składanie macierzy globalnej Zmieniamy punkt widzenia: (z funkcji kształtu na elementy) x m-1 x m -1 1 J m =h m /2 u m ( )=u 1 m 1 ( ) +u 2 m 2 ( ) u 1 u 2 (parametry węzłowe niewiadome) v i (x( )) =1/2 –1/2 v i+1 (x( )) =1/2 1 =1/2 –1/2 2 =1/2 +1/2 element [ funkcje bazowe : ważone parametrami węzłowymi ]
36 Zmieniamy punkt widzenia: (z funkcji kształtu na elementy) x m-1 x m -1 1 J m =h m /2 macierz sztywności elementu m [wymiar taki jak liczba funkcji kształtu na element] u m ( )=u 1 m 1 ( ) +u 2 m 2 ( ) u 1 u 2 (parametry węzłowe niewiadome) 1 =1/2 –1/2 2 =1/2 +1/2 element Em=Em= E m 11 E m 12 E m 21 E m 22 zależność od m w J m : m x( ) (x m-1 +x m )/2+(x m -x m-1 )/2
37 na przekrywających się węzłach:suma Em=Em= E m 11 E m 12 E m 21 E m 22 węzły na granicy elementów obsługują więcej niż jeden element Składanie (assembly) globalnej macierzy sztywności macierz globalna S (rozmiar = liczbie węzłów) U 1 U 2 U 3 U 4 1 2 3 u m ( )=u 1 m 1 ( ) +u 2 m 2 ( ) u 2 1 = u 1 2 na U 2 opiera się rozwiązanie w dwóch sąsiednich elementach globalna [U] i lokalna [u] numeracja węzłów do macierzy S wchodzą całki po sąsiednich elementach z funkcją kształtu wspólną dla sąsiadów
38 u(x=-1)=0 u(x=1)=0 Przedział (-1,1) Podzielony na 7 elementów (8 węzłów) x m-1 x m -1 1 J m =h m /2 u( )=u 1 1 ( ) +u 2 2 ( ) u 1 u 2 1 =1/2 –1/2 2 =1/2 +1/2 case study
39 1234 h2h2 h3h3 h4h4 Składanie (assembly) macierzy sztywności z całek po elementach Element 2 Element 3 Forma już znana dodajemy elementy z różnych macierzy lokalnych które odpowiadają temu samemu węzłowi
40 Wektor obciążeń pojedynczego elementu/składanie globalnego po elemencie K i po K i+1 xixi x i-1 x i+1 U i-1 U i U i+1 druga funkcja elementu i oraz pierwsza elementu i+1 = to ta sama funkcja v i (x)
41 Składanie globalnej macierzy sztywności – buchalteria węzłów m=1 m=2m=3 m-numeruje elementy g=1 g=2 g=3 g=4 g – globalna numeracja węzłów 1 2 1 2 1 2 lokalne numery węzłów nr (k,m) – przyporządkowanie numeru globalnego węzłowi o lokalnym numerze k w elemencie m 1 = nr (1,1) 2 = nr (2,1) = nr (1,2) 3 = nr (2,2) = nr (1,3) 4 = nr (2,3)
42 pętla po elementach m=1,M pętla po węzłach lokalnych k=1,N pętla po węzłach lokalnych l=1,N identyfikacja numeru globalnego węzła i=nr(k,m) j=nr(l,m) S(i,j)=S(i,j)+E(m,k,l) Składanie globalnej macierzy sztywności – buchalteria węzłów m=1 m=2m=3 m-numeruje elementy g=1 g=2 g=3 g=4 g – globalna numeracja węzłów 1 2 1 2 1 2 lokalne numery węzłów nr (k,m) – przyporządkowanie numeru globalnego węzłowi o lokalnym numerze k w elemencie m 1 = nr (1,1) 2 = nr (2,1) = nr (1,2) 3 = nr (2,2) = nr (1,3) 4 = nr (2,3) identycznie składa się macierze dla wyższych funkcji kształtu i w więcej niż 1D
43 pętla po elementach m=1,M pętla po węzłach lokalnych k=1,N pętla po węzłach lokalnych l=1,N identyfikacja numeru globalnego węzła i=nr(k,m) j=nr(l,m) S(i,j)=S(i,j)+E(m,k,l) Składanie globalnej macierzy sztywności – buchalteria węzłów m=1 m=2m=3 g=1 g=2 g=3 g=4 1 2 1 2 1 2 1 = nr (1,1) 2 = nr (2,1) = nr (1,2) 3 = nr (2,2) = nr (1,3) 4 = nr (2,3)
44 pętla po elementach m=1,M pętla po węzłach lokalnych k=1,N identyfikacja numeru globalnego węzła i=nr(k,m) F(i)=F(i)+P(k,m) składanie globalnego wektora obciążeń z lokalnych:
45 o potrzebie używania wyższych funkcji kształtu (i o laboratorium): u(1)=u(-1)=0 rozwiązanie dokładne: równanie Poissona z odcinkowo stałą prawą stroną równania
46 z liniowymi funkcjami kształtu: poza węzłami nie uzyskamy dokładnego rozwiązania tego równania (nigdy nie uzyskamy rozwiązania silnej postaci równania, druga pochodna wewnątrz elementów jest zawsze równa zeru a ma być równa niejednorodności dla równań Poissona niejednorodność to źródło potencjału, dla równania przew. ciepl. – źródło ciepła) o potrzebie używania wyższych funkcji kształtu (i o laboratorium): u(1)=u(-1)=0 rozwiązanie dokładne: równanie Poissona z odcinkowo stałą prawą stroną równania rozwiązanie dokładne:
47 całka działania a rozkład elementów dla funkcji odcinkowo linowych: dokładne działanie optymalne rozwiązanie odcinkami liniowe rozwiązanie dokładne i=2,8
48 Funkcje kształtu Lagrange’a wyższych rzędów: u1u1 u2u2 u4u4 funkcje kształtu jeden element, cztery (n) węzły u3u3 wielomian stopnia n-1, taki, że
49 Funkcje kształtu Lagrange’a wyższych rzędów: u1u1 u2u2 u4u4 funkcje kształtu jeden element, cztery (n) węzły u3u3 wielomian stopnia n-1, taki, że wiemy jak go wskazać: wielomian węzłowy Lagrange’a
50 Funkcje kształtu wyższych rzędów: u1u1 u2u2 u4u4 funkcje kształtu jeden element, cztery (n) węzły u3u3 wielomian stopnia n-1, taki, że wiemy jak go wskazać: wielomian węzłowy Lagrange’a funkcje kształtu Lagrange’a: rozwiązanie interpolowane wielomianowo w każdym z elementów. jedynie ciągłość rozwiązania między elementami. w przeciwieństwie do problemów z KSN: wartości funkcji w węzłach nie są znane. należy je wyliczyć. istota FEM.
51 Elementy wyższych rzędów: Jeden element, trzy funkcje bazowe, 3 parametry węzłowe Funkcje bazowe : w danym węźle tylko jedna z nich niezerowa (co min. gwarantuje liniową niezależność funkcji bazowych) 1 = ( -1)/2 2 =-( -1)( +1) 3 = ( +1)/2 funkcje wierzchołkowe (vertex functions) 1 na krawędziach elementu funkcja bąbelkowa (bubble function) węzeł wewnątrz elementu -1 0 1 u1u1 u2u2 u3u3
52 węzły= granice elementów liniowa baza Lagrange’a kwadratowa baza Lagrange’a Funkcje kształtu Lagrange’a: odcinkowo liniowe i kwadratowe każda strzałka to węzeł granice elementów czerwone
53 1 = ( -1)/2 2 =-( -1)( +1) 3 = ( +1)/2 liczone numerycznie całki wyliczone analitycznie: ilu punktowego Gaussa należałoby użyć aby dokładnie scałkować m.sztywności numerycznie? przy równym podziale przedziału E takie samo dla każdego elementu lecz P nie! [inny zakres x( )] x( ) (x m +x m+1 )/2+(x m+1 - x m )/2 Macierz sztywności dla kwadratowych f. Lagrange’a
54 Liczba wierszy: 2n+1 (n-liczba elementów) Składanie globalnej macierzy sztywności i wektora obciążeń dla kwadratowych funkcji Lagrange’a lokalne globalne
55 Laboratorium: 4 elementy: x x x x 2 centralne elementy o długości bx krzyżyki: węzły bąbelkowe dokładne MES dla bx=0.5 (równoodległe węzły) widzimy: dokładne dla węzłów granicznych
56 Laboratorium: 4 elementy: dokładne 2 centralne elementy o długości bx krzyżyki: węzły bąbelkowe x x x x a na laboratorium zobaczymy, że potencjał dokładny odtworzony
57 Wyniki dla problemu modelowego funkcje odcinkami liniowe 9 węzłów (8 elementów) funkcje kwadratowe 9 węzłów (4 elementy) rozwiązanie dokładne
58 Funkcje liniowe i kwadratowe funkcje odcinkami liniowe 9 węzłów (8 elementów) funkcje kwadratowe 9 węzłów (4 elementy) [ błąd =0 tylko na węzłach wierzchołkowych ] jeden rząd funkcji kształtu więcej: maksymalny błąd zmniejszony trzykrotnie rozmiar problemu liniowego bez zmian, ale S ma więcej niezerowych elementów
59 dla kwadratowej bazy Lagrange’a pochodna jest nieciągła na granicy elementów nieciągła pochodna
60
61 -1 0 1 u1u1 u2u2 u3u3 Co zrobiliśmy : poprowadziliśmy przez każdy element wielomian interpolacyjny Lagrange’a. Zabieg zakończył się sukcesem. Lepsza dokładność prawie tej samej złożoności obliczeniowej w porównaniu z liniowymi funkcjami kształtu. Rozmiar URL bez zmian, ale macierz układu – więcej niezerowych elementów. Chcemy podnieść rząd wielomianu interpolacyjnego. Czy równomiernie rozłożenie większej ilości węzłów na elemencie jest dobrym pomysłem ? NIE
62 Błąd interpolacji Lagrange’a (przypomnienie): x 0,x 1,..., x n – n+1 różnych węzłów f(x) –gładka funkcja interpolowana (klasy co najmniej n+1) x w przedziale interpolacji należy do (najmniejszego) przedziału, w którym mieszczą się punkty x i norma nieskończoność: ||g(x)|| = max |g(x)| w przedziale (a,b) odchylenie funkcji interpolowanej od wielomianu Lagrange’a
63 -1 1 4 punkty 0 -a a punkty rozłożone równomiernie norma nieskończoność wielomianu węzłowego równomierne rozłożenie nie jest optymalne dla celów aproksymacyjnych
64 Efekt Rungego nieoptymalność interpolacji na równoodległych węzłach robi się drastyczna dla wysokiego rzędu wielomianu interpolacyjnego 3 punkty 5 punktów 11 punktów im wyższy stopień wielomianu interpolacyjnego tym gorsze przybliżenie [większa norma nieskończoność błędu] - szczególniej przy brzegach przedziału
65 Węzły Czebyszewa: bliskie optymalnym 3 punkty 5 punktów 11 punktów 21 punktów więcej węzłów przy brzegach [w.czebyszewa: węzły wielomianów ortogonalnych na przedziale (–1,1) z wagą 1/sqrt(1-x**2) [! waga gęsto punkty przy brzegu, błąd nie urośnie ]
66 -1 -0.5 0.5 1 u1u1 u2u2 u4u4 u3u3 Kubiczne funkcje kształtu Lagrange’a z węzłami Czebyszewa cos( j/3), j=0,1,2,3,: -1,-0.5,0.5,1 Jeden element, 4 funkcje bąbelki wierzchołek
67 błąd: funkcje odcinkami liniowe 9 węzłów (8 elementów) funkcje kwadratowe 9 węzłów (4 elementy) funkcje kubiczne 10 węzłów (3 elementy) zwiększenie stopnia wielomianów kształtu o jeden: max. odchylenie wyniku od dokładnego zmniejsza się 3 krotnie Wyniki dla problemu modelowego
68 MES używa jako funkcji bazowych określonych na elemencie wielomianów potrafimy je numerycznie różniczkować i całkować dokładnie x 1 1 1 1 x 2 2x 2x+ x 2x 2x x 3 3x 2 3x 2 +3x x+ x 2 3x 2 + x 2 3x 2 x 4 4x 3 4x 3 +6x 2 x+4x x 2 + x 3 4x 3 +4x x 2 4x 3 x 5 5x 4 5x 4 +10x 3 x+10x 2 x 2 +5x x 3 + x 4 5x 4 +10x 2 x 2 + x 4 5x 4 -4 x 4 u’(x) (błąd na czerwono) u(x) C=1/2C=-1/6 różniczkowanie: a całkowanie... Gaussa
69 kwadratury Gaussa-Legendra do całkowania elementów macierzowych Gauss= najbardziej efektywna metoda dla MES funkcje kształtu są wielomianami(!), a Gauss całkuje je dokładnie ważona suma funkcji podcałkowej w wybranych punktach x i Chcemy wybrać tak wagi i punkty aby kwadratura była dokładna dla wielomianu jak najwyższego stopnia (funkcje kształtu będą wielomianami) Na pewno uda nam się skonstruować kwadraturę dokładną dla wielomianu stopnia n-1
70 Wybieramy wagi i punkty Gaussa, tak aby dokładnie scałkować wielomian stopnia 2n-1 [2n współczynników, 2n wag i punktów] Przykład: n=2 – dokładnie scałkujemy wielomian stopnia 3 f(x)=a+bx+cx 2 +dx 3 kwadratury Gaussa-Legendra
71 f(x)=a+bx+cx 2 +dx 3 a,b,c,d – dowolne. Każda z powyższych całek musi zostać policzona dokładnie. wstawiamy po kolei 1 za jeden z a,b,c,d=reszta 0. [kwadratura ma działać również dla f(-x)] x 1 oraz x 2 będą rozłożone symetrycznie względem 0 ( x 1 =-x 2 ) wtedy z (2) w 1 =w 2 =1 (z 1) (4) - zawsze spełnione 2/3=x 2 2 +x 2 2 z (3) x 2 = (1/3) 1/2 x 1 =-x 2
72 kwadratura Gaussa dokładna dla wielomianów stopnia 3: w 1 =w 2 =1, x 1 =(1/3) 1/2 x 2 =- (1/3) 1/2 wystarczy dodać wartości funkcji w dwóch punktach aby uzyskać dokładną całkę dla wszystkich wielomianów stopnia 3
74 Na przedziale –1,1 wybieramy (dowolnie) n – punktów i prowadzimy przez nie wielomian interpolacyjny Lagrange’a funkcji f(x) Próbkując funkcję w n dowolnych punktach: na pewno uda się skonstruować kwadraturę dokładną dla wielomianu stopnia n-1 Jeśli f(x) – wielomian stopnia nie większego niż n-1 f(x)=y(x) (interpolując wielomian dostaniemy ten sam wielomian) na wyborze punktów x i można zyskać dokładność dla n stopni więcej
75 Dalej o wyborze punktów Gaussa: Tw. Jakobiego: kwadratura oparta na wielomianie interpolacyjnym Lagrange’a jest dokładna dla wielomianów stopnia 2n-1, jeśli punkty x i wybrane tak, że wielomian stopnia n jest ortogonalny do wszystkich wielomianów stopnia (n-1) zobaczmy, że tak jest: kwadratury Gaussa-Legendra
76 dla dowolnego wielomianu stopnia n i dowolnej liczby r istnieje taki wielomian o stopniu o jeden niższym i taka liczba R, że: P n (x)=(x-r) P n-1 (x)+R przykład: 1+x+x 2 =(x-2) (ax+b)+c=c-2b+(b-2a)x+ax 2 – wyliczymy sobie a,b, oraz c f 2n-1 (x)=z n (x)f n-1 (x)+q n-1 (x) kwadratury Gaussa-Legendra f 2n-1 (x)=(x-x 1 ) f 2n-2 (x)+r 0 f 2n-1 (x)=(x-x 1 ) [(x-x 2 ) f 2n-3 (x)+r 1 ]+r 0 =(x-x 1 )(x-x 2 ) f 2n-3 (x)+r 0 +r 1 (x-x 1 ) q 1 (x)
77 całka oparta o przepis interpolacyjny na n punktach będzie dokładna dla każdego wielomianu stopnia n-1 Problem: jak wybrać wielomian stopnia n z(x) tak aby ortogonalny dla każdego wiel. stopnia n-1 f 2n-1 (x)=z n (x)f n-1 (x)+q n-1 (x)
78 Problem: jak wybrać z(x) aby ortogonalny dla każdego wiel. stopnia n-1 każdy wielomian można zapisać w postaci sfaktoryzowanej wybrać zera znaczy wybrać wielomian (co do stałej multiplikatywnej) wielomian Legendre’a stopnia n -ortogonalny na przedziale [–1,1] do wszystkich wielomianów stopnia n-1. -zera tego wielomianu wyznaczą optymalne punkty Gaussa
79 Przedział [-1,1]. Mamy zbiór niezaleznych liniowo funkcji h 0 =1, h 1 =x, h 2 =x 2, h 3 =x 3,... które nie są ortogonalne [iloczyn skalarny określony z funkcją wagową w(x)]. Chcemy skonstruować bazę wielomianów ortogonalnych. funkcje bazowe dla tego przedziału, z wagą w(x)=1 są to wielomiany Legendre’a. u 0 = 1 u 1 =a+x Jakie a aby (u 0,u 1 )=0 ?: odp.: a=0 u 1 =x u 2 =x 2 +bx+c (u 2,u 0 )= 2/3+2c=0 (u 2,u 1 )=0 b=0 u 2 =(x 2 -1/3) W literaturze wielomiany Legendre’a normalizowane tak aby P k (1)=1 : 1,x,3/2 (x 2 -1/3) itd. Ortogonalizacja Grama-Schmidta
80 Punkty Gaussa zapewniające maksymalną dokładność (do wielomianu stopnia 2n-1): zera n-tego wielomianu Legendra Dla 2n-1=3 [ punkty Gaussa tam gdzie wcześniej wyliczyliśmy ] l 1 =(x+1/sqrt(3))/(2/sqrt(3)). całka z niego od od –1 do 1 = 1 kwadratury Gaussa-Legendra W bazie P 0,P 1,..., P n-1 można opisać wszystkie wielomiany stopnia n-1, P n ortogonalny do wszystkich wektorów bazy, więc i do wszystkich wielomianów stopnia n-1
81 Dokładne do wielomianów stopnia 3 stopnia 5 stopnia 11 Wagi i punkty Gaussa
82 MES z liniowymi funkcjami kształtu: przykład zastosowania nr. 2: równanie oscylatora harmonicznego niewiadome: funkcja własna u n (x) oraz wartość własna E n E n =1/2+n rysunek z Wikipedii
83 x xixi x i+1 x i-1 v i (x) 1 element K i długości i =x i -x i-1 element K i+1 długości i+1 = x i+1 -x i węzły baza funkcyjna
84 Hc=EOc tzw. uogólnione macierzowe równanie własne [„zwykłe” równanie własne gdy : O=1] macierze H,O – trójprzekątniowa, symetryczna nasz operator=samosprzężony (hermitowski) całki się liczy analitycznie [wielomiany stopnia najwyżej 4]
85 równomierny rozkład 21 węzłów (elementy równej długości) siatka od –6.2 do 6.2 E dok E num błąd 5.56.048 0.548 4.54.878 0.378 3.53.736 0.235 2.52.625 0.125 1.51.549 0.049 0.50.510 0.01 MRS: przy tej samej liczbie węzłów E num błąd 4.772 -0.728 4.027 -0.473 3.220 -0.280 2.358 -0.172 1.446 -0.054 0.489 -0.011 względna przewaga MES rośnie dla stanów o wyższej energii
86 kropy: wartości z MES linie: dokładne funkcje własne MES nie jest ściśle dokładna w węzłach jak dla Poissona, ale niezła różnica dokładny=MES dla wektora własnego o najniższej wartości własnej
87 przydałaby się siatka nierównomierna: więcej węzłów tam gdzie wektory własne przyjmują duże wartości, mniej na ogonach n=21 punktów =1 =1.5
88 oscylator harmoniczny: optymalizacja rozkładu elementów E n =n+1/2 n=0 n=1 2 3 4 dla wszystkich stanów optymalny jest wykładnik 1.5 =1 =1.5 E
89 Metoda różnic skończonych, siatka nierównomierna Iloraz różnicowy drugiej pochodnej dla nierównej siatki: ll pp + tracimy jeden rząd dokładności w porównaniu z siatką równomierną Problem rozwiązany w metodzie elementów skończonych. Wzór trójpunktowy W MES: nie ma problemu bo pochodne i całki liczymy dokładnie!
90 k-1 k k+1 [różnica z MES: potencjał tylko do diagonali, O=1] Uwaga: metoda produkuje niesymetryczną macierz H: komplikacja, np. dla niesymetrycznych macierzy nie ma gwarancji, że wartości własne będą rzeczywiste, a wektory własne wzajemnie ortogonalne
91 Metoda różnic skończonych, siatka nierównomierna n=0 n=1 2 wartości własne mniejsze niż dokładne. dla n=0, optymalna jest siatka równomierna. Pewna poprawa (wzrost) jest dla wzbudzonych. Ale (!!) wiemy, że wzrost to poprawa tylko dlatego, że rozwiązanie dokładne jest znane. W typowej sytuacji - nie znamy rozwiązania dokładnego. Nie będziemy wiedzieli czy mamy do czynienia z poprawą czy pogorszeniem. w MES: przeciwnie, metoda ma charakter wariacyjny! Lepsza siatka znaczy mniejsza wartość własna! Parametr dla optymalizacji siatki! E n =n+1/2
92 Lu=f eliptyczne rozwiązanie poszukiwane w bazie funkcyjnej Hu=Eu własne (wstawić rozwiązanie próbne, wyrzutować na wektory bazy) (układ równań algebraicznych) Lc=fHc=EOc (uogólniony macierzowy problem własny) funkcja próbna (trial) równanie eliptyczne a różniczkowe równanie własne w metodzie elementów skończonych 1D: np. L=-d 2 /dx 2 np. H=-d 2 /dx 2 +x 2
93 wariacyjny charakter metody dla równania eliptycznego i własnego MES (i Galerkin w ogóle) równoważna metodzie Reyleigha-Ritza, gdy ta stosowalna Lu=fHu=Eu wartościom c i, które spełniają równania algebraiczne odpowiada minimum funkcjonału:
94 dokładne rozwiązania równania oscylatora są nie tylko ciągłe, ale również ciągłe z pochodną baza Lagrange’a: gwarantuje tylko ciągłość funkcji na granicy elementów ciągłość pochodnej nie jest gwarantowana jeśli zainwestujemy w funkcje kształtu tak aby zapewnić ciągłość z pochodną - niższe wartości własne wyjść powinny kiedy ciągłość pochodnej ważna 1)dla równania Poissona: nieciągłe pole elektryczne (pochodna potencjału) gdy rozkład gęstości ładunku o formie delty Diraca 2) równanie przewodnictwa cieplnego: strumień ciepła proporcjonalny do gradientu temperatury = gdy ten nieciągły na granicy elementów stan nie może być ustalony
95 MES: baza funkcji ciągłych z pochodnymi na granicy elementów k-1 k k+1 kk k+1 f g f przenosi ciągłą wartość funkcji przez granice elementów: f: 1 w węźle k, 0 w sąsiednich f’ zero we wszystkich węzłach g przenosi ciągłą pochodną przez dwa elementy: g’: 1 w węźle k, 0 w sąsiednich g zero we wszystkich węzłach parametry węzłowe: 4 na element – wartości i pochodne na obydwu końcach przedziału zapewniały ciągłość funkcji: 2 nie znikające f na element, jedna na węzeł 4 funkcje kształtu na element po 2 funkcje na węzeł z węzłem k wiążemy dwie funkcje:
96 parametry węzłowe: tutaj u i u’ współczynniki rozwinięcia (niewiadome) : wartości funkcji i pochodnej w węzłach Metoda ES da nam ekstra oszacowanie pochodnych rozwiązania
97 konstruujemy funkcję f k : stopnia musi trzeciego wielomian k-1 k k+1 kk k+1 f współczynniki a k, b k,c k, d k : f(x k )=1, f(x k-1 )=0, f’(x k )=0, f’(x k-1 )=0 podobnie
98 f f’ przykładowy przebieg funkcji f niosącej ciągłość rozwiązania przez granicę elementów
99 konstruujemy funkcję g k : k-1 k k+1 kk k+1 g g(x k-1 )=g(x k )=g’(x k-1 )=0 g’(x k )=1
100 g g’ przykładowy przebieg funkcji g niosącej ciągłość pochodnej rozwiązania przez granicę elementów
101 funkcje f k g k są ortogonalne jeśli węzły równoodległe symetria, ale (f k,g k+1 ) już nie 0 w teorii interpolacji: wielomiany 3-go rzędu interpolujące wartości i pochodne - sklejki Hermite’a
102 Hc=EOc (f,Hf) (g,Hg) (g,Hf) (f,Hg) macierze 2n na 2n każda z ćwiartek trójprzekątniowa sklejka Hermita dla równania własnego: struktura macierzy
103 4.502 3.510 2.500 1.500 0.500 równoodległe węzły zaznaczone kropami rewelacyjny wynik 13 węzłów/26 równań liniowe 23 węzły / (23 równania) 2.625 1.549 0.51
104 optymalizacja elementów (więcej węzłów na środku ) wynik dla sklejek Hermite’a -funkcje kształtu są tak elastyczne, że optymalizacja elementów staje się zbędna n=0 wynik dla funkcji liniowych =1 =1.5
105 funkcje kształtu Hermite’a w przestrzeni referencyjnej zastosowanie dla problemów eliptycznych 1D -1 1 u 1 =u(x 1 )u 2 =u(x 2 ) Pochodne funkcji (po x nie po ) = nowe elementy węzłowe = oprócz ciągłości przez granicę elementu dostaniemy oszacowanie pochodnych na granicy
106 4 równania: 2 wartości, 2 pochodne Szukamy funkcji kształtu dla rozwiązań ciągłych z pochodną (funkcji kształtu Hermite’a) x( ) (x m-1 +x m )/2+(x m -x m-1 )/2
107 Funkcje kształtu Hermita na rysunku funkcje związane z pochodnymi / Jakobian każda z funkcji kształtu Hermite’a: odpowiedzialna za wartość lub pochodną, żadna nie przeszkadza pozostałym (pełen podział kompetencji)
108 Funkcje kształtu Hermite’a: interpolacja jednomianów, przedział [–1:1]: dokładna funkcja v( )= u 1 0 = -1, u 2 0 =1, u 1 1 =1, u 2 1 =1 u= wartości pochodne dokładna funkcja v( )= u 1 0 = 1, u 2 0 =1, u 1 1 =-2, u 2 1 =2 u=1 dokładna funkcja v( )= u 1 0 = -1, u 2 0 =1, u 1 1 =3, u 2 1 =3 u= dokładna funkcja v( )= u 1 0 = 1, u 2 0 =1, u 1 1 =-4, u 2 1 =4 u=-1+2 J m =1
109 x4x4 (-1+2x 2 ) kubiczna sklejka Hermita a wielomian 4-tego rzędu wartość funkcji i pochodnej OK, ale sukces interpolacji umiarkowany
110 Trzeba zawęzić przedział x4x4 (-1+2x 2 )
111 zawęzić przedział x 4 od [–1,1] do [0,1] aby nie przejmować się czynnikiem skali J m : trzymamy [-1,1] a przekształcamy funkcję [(x+1)/2] 4 dokładna funkcja v( )=[( u 1 0 = 0, u 2 0 =1, u 1 1 =0, u 2 1 =2 u= itd.
112 u(x=-1)=0 u(x=1)=0 Przedział (-1,1) Podzielony na 8 elementów (9 węzłów, 18 parametrów węzłowych) J m =h m /2 funkcje kształtu Hermita: zastosowanie -1 1 x( ) (x m-1 +x m )/2+(x m -x m-1 )/2 1 2 34 macierz sztywności dla pojedynczego elementu
113 składanie macierzy sztywności dla f.k. Hermita dwa elementy wspólne parametry węzłowe: wartość pochodna wspólne również elementy pozadiagonalne dla tego węzła x=-1x=+1 u(lewy) u’(lewy) u(prawy) u’(prawy)
114 Prawa strona/jeden element
115 Warunki brzegowe: u(x=-1)=0 u(x=1)=0 u(-1) u’(-1) u(0) u’(0) u(1) u’(1)
116 Warunki brzegowe: u(x=-1)=0 u(x=1)=0 u(-1) u’(-1) u(0) u’(0) u(1) u’(1) tu damy zera tu damy jedynki na diagonali Su=F
117 dwa elementy / 3 węzły / 6 parametrów węzłowych dokładne FEM z FKH pochodne uzyskane z rozwiązania układu równań w węzłach: -0.387, 0.387, -0.387 zamiast dokładnej -1/ , 1/ , -1/ -0.318, 0.318, -0.318 Wyniki dla problemu modelowego
118 3 elementy / 4 węzły / 8 parametrów węzłowych dokładne FKH+ 2E FKH+ 3E Wyniki dla problemu modelowego
119 Przy tej samej złożoności: w porównaniu z bazą Lagrange’a w bazie Hermita błąd mniejszy (i ciągły, bo rozwiązanie przybliżone ciągłe) Baza kubiczna Lagrange’a i baza Hermita porównanie dokładności
120 funkcje kształtu Hermita a warunki brzegowe Neumanna u(x=-1)=0 u’(x=1)=-1/ Przykład dla dwóch elementów 1 1 0 co tutaj? odpowiedzialna za pochodną na końcu przedziału mamy ustawić tak, żeby pochodne pozostałych funkcji znikają w 1
121 1 1 0 -1 / warunek Neumanna dla funkcji kształtu Hilberta
122 -1.00-0.500.000.501.00 -0.20 0.00 0.20 0.40 0.60 0.80 Wynik umiarkowany sukces pochodna jak trzeba ale wartość całkiem nie taka co jest nie tak? warunek Neumanna dla funkcji kształtu Hermite’a
123 macierz sztywności pojedynczego elementu Przypomnijmy sobie skąd się wzięła w metodzie Galerkina v(x) jest funkcją bazową. dla warunków Dirichleta u(1)=u(-1)=0 usuwaliśmy z bazy funkcje, które nie znikały na prawym brzegu. v(x)=0 Teraz mamy parę funkcji dla których wyraz nie znika. Jaką? warunek Neumana dla funkcji kształtu Hilberta
124 Teraz mamy parę funkcji, dla których wyraz nie znika. Jaką?
125 no i co?
126 1 1 to już mamy to trzeba dodać S 56 =S 56 +1 warunek Neumana dla funkcji kształtu Hermita
127 Wynik: prawy warunek brzegowy Neumanna: 2 elementy Znacznie lepiej prawie dobrze, pochodna idealnie, wartość nie (poprzednio – przy warunkach Dirichleta – wartość była super pochodna taka sobie) 2 elementy gdy więcej?
128 Wyniki : 4 elementy Wynik: prawy warunek brzegowy Neumanna