Debreceni Egyetem Informatika Kar
ELLIPTIKUS EGYENLETEK MEGOLDÁSA VÉGESELEM MÓDSZERREL
Témavezető:
Készítette:
Dr Baran Ágnes
Szuetta Judit
egyetemi adjunktus
informatika tanár
Debrecen 2010
KÖSZÖNETNYILVÁNÍTÁS Köszönetet mondok témavezetőmnek, Dr Baran Ágnesnek a dolgozat készítése során nyújtott segítő munkájáért, azért, hogy mindig irányt mutatott a tudomány útvesztőjében, és nem utolsó sorban végtelen türelméért.
Tartalomjegyzék
1.
Bevezetés ........................................................................................................................2
2.
Parciális differenciálegyenletek......................................................................................4
3.
A Poisson-egyenlet .........................................................................................................5
4.
A variációs feladat ..........................................................................................................7
5.
Szoboljev-terek...............................................................................................................9
6.
A variációs feladat megoldhatósága .............................................................................11
7.
Véges elemek................................................................................................................13
8.
A T6 elem .....................................................................................................................16
9.
A diszkrét egyenlet .......................................................................................................19
10.
A feladat végeselem megoldása ...................................................................................21
11.
IRODALOMJEGYZÉK .......................................................................................................33
Bevezetés Fizikai
jelenségek
matematikai
modellezésekor
gyakran
parciális
differenciálegyenletek felírására kerül sor. Ilyen például a hullámegyenlet, mely homogén közegben terjedő hullámszerű mozgást ír le; a hővezetés egyenlete, mely a hő terjedését írja le homogén közegben; a Schrödinger-egyenlet, mellyel részecske állapotok különböző tulajdonságai határozhatók meg vagy a rezgő húr probléma. Az említett példákból is következtethetünk arra, hogy a műszaki- és természettudományokban fontos szerepet töltenek be az ott gyakran előforduló másodrendű lineáris parciális differenciálegyenletek. Elliptikus egyenletek olyan fizikai jelenségek modellezésénél fordulnak elő, melyek egyensúlyi állapotokat írnak le, illetve az idő nem befolyásoló tényező. A dolgozat elliptikus differenciálegyenletek, elsősorban a Poisson–egyenlet numerikus megoldását ismerteti a többféle megoldási eljárás közül kiválasztott végeselem módszerrel. Az elliptikus egyenletek megoldása – ezen belül a Poisson–egyenleté – igen fontos feladat, hiszen a parabolikus differenciálegyenletek is elliptikus egyenletekre vezethetők vissza azok numerikus megoldásához, mindemellett a Poisson- egyenlettel többféle fizikai jelenség modellezhető (pl: vékony rudak csavarása; mágneses tér vagy elektrosztatikus mező potenciálja; gázok illetve folyadékok áramlása; a hővezetési egyenlet, ha a hőmérséklet nem függ az időtől). A végeselem módszernek a mérnöki munkában nagyon sok alkalmazási területe van (például repülőgép- és autóiparban). Néhány modern végeselem program különböző szimulációs munkakörnyezetet is tartalmaz. Egy rendszer modellezéséhez és elemzéséhez hatékony segítséget nyújthat a szimuláció, a különböző (esetleg extrém) környezeti paraméterek beállítása. Nem beszélve arról, hogy egy ilyen software segítségével tervezett virtuális prototípus tesztelésének, módosításának a költsége, a tér- és az időigénye, valamint a munkaerőigénye jóval kevesebb, mint egy fizikailag megvalósítotté. A végeselem módszer összefoglalva a következő lépésekből áll: − A kitűzött peremérték feladatról egy alkalmasan kiválasztott Hilbert–térben áttérünk a variációs egyenletre.
2
− Ennek létezik egyértelmű, stabil (ún. általánosított) megoldása − A variációs feladat közelítő megoldásához a kiválasztott Hilbert–térről annak egy végesdimenziós alterére térünk át. Ebben az altérben szintén létezik egyértelmű, stabil megoldás, és a Céa-lemma alapján a diszkrét feladat közelítő megoldásának eltérése a pontos megoldástól nem nagyobb mint a legjobb közelítésé. − A diszkrét feladat ekvivalens egy lineáris egyenletrendszerrel, melynek mátrixa szimmetrikus és pozitív definit. A dolgozatban ezen lépések végrehajtásához szükséges elméleti ismeretek feldolgozására kerül sor, ennek gyakorlati megvalósítása a MATLAB rendszer segítségével történik. A dolgozat szerves részét képező ellip_megold.m program elsőfajú peremfeltétellel adott Poisson–egyenlet u megoldását szolgáltatja bármely
f
(adott) függvény esetén az
egységnégyzeten – melynek diszkretizációjához T6 elemeket használunk.
3
Parciális differenciálegyenletek Parciális differenciálegyenlet alatt olyan differenciálegyenletet értünk, melyet egy ismeretlen többváltozós függvény parciális deriváltjai ismeretében –esetleg kezdeti feltételek mellett- írunk fel. Parciális differenciálegyenlet rendjén a benne szereplő parciális deriváltak maximális rendjét értjük. A főrészében lineáris másodrendű parciális differenciálegyenlet általános alakja:
∑ Aij ( x1 ,..., xn ) i, j
∂ 2u ∂u ∂u + F ( x1 ,..., x n , u , ,..., )=0 ∂xi x j ∂x1 ∂x n
ahol u : R n → R , n ∈ N és i, j = 1,..., n Konstans együtthatós az egyenlet, ha a parciális deriváltak együtthatói nem függenek az x1 ,..., x n változóktól. Ebben az esetben a független változók egy homogén lineáris transzformációjával az egyenlet a következő, egyszerűbb alakra hozható:
∑ ci i
∂ 2u ∂u ∂u + F ( x1 ,..., x n ,u , ,..., ) = 0 ( n ∈ N és i = 1,..., n ), 2 ∂x1 ∂x n ∂xi
ahol ci ∈ {− 1,0,1} . A ci együtthatók értékétől függően a következő csoportokba sorolhatjuk a parciális differenciálegyenleteket: 1.
Elliptikus differenciálegyenletről beszélünk, ha bármely i-re ci = 1 .
2.
Hiperbolikus differenciálegyenletről beszélünk, ha bármely i-re ci ≠ 0 és az egyik együttható előjele különbözik az összes többiétől.
3.
Parabolikus differenciálegyenletről beszélünk, ha a ci együtthatók közül pontosan egy nulla, a többi azonos előjelű. [5]
4
A Poisson-egyenlet A Poisson-egyenlet ∆u + f = 0 ,
ahol Ω ⊂ R 2 , és adott f : Ω → R függvény esetén keressük az u : Ω → R függvényt. ∆ a Laplace-operátort jelöli. A kétdimenziós Laplace-operátor:
∂ 2u ∂ 2u ∆u = 2 + 2 ∂x1 ∂x2 Elliptikus egyenleteket legtöbbször valamilyen peremfeltétel mellett oldunk meg. A Poisson-egyenlet peremérték feladata:
∆u + f = 0
x∈Ω
u ( x) = g ( x)
ha x ∈ Γ
ahol Γ az Ω tartomány pereme. Ha g ( x) ≡ 0 , akkor homogén peremfeltételekről beszélünk.
Az adott Ω tartomány Γ peremén meg kell adnunk a keresett megoldás értékét, a tartomány belsejében érvényes a differenciálegyenlet, és ez határozza meg a megoldást.
5
Megfelelő feltételek – az Ω tartomány Γ peremének Lipschitz-folytonossága és a differenciáloperátor egyenletes elliptikussága – mellett egy ilyen feladat megoldása létezik, egyértelmű és folytonosan függ a feladat adataitól. A peremérték feladat korrekt kitűzésű, ha rendelkezik az egzisztenciához, az unicitáshoz és a stabilitáshoz szükséges tulajdonságokkal.
Definíció: (Lipschitz-folytonos perem) Legyen Ω ⊂ R n korlátos tartomány. Az Ω tartomány Γ pereme Lipschitzfolytonos, ha van olyan véges lefedése nyílt R n -beli környezetekből, melyek mindegyikében, alkalmasan választott lokális koordinátarendszerben a perem egyenlete x n = φ ( x1 , K , x n−1 ) alakú, ahol φ : R n-1 → R Lipschitz-folytonos. [3]
Megjegyzés: Az n-dimenziós gömb és az n-dimenziós kocka is ilyen peremmel rendelkezik. Így feladatunkban - ahol Ω két dimenziós téglalap tartomány – a peremre teljesül a Lipschitz-folytonossági feltétel.
Definíció: (egyenletes elliptikusság) Egy differenciáloperátort egyenletesen elliptikusnak nevezünk Ω-ban, ha az együtthatóiból összeállított A = A( x ) mátrix egyenletesen pozitív definit, azaz az euklideszi skalárszorzatban pozitív k 0 konstanssal teljesül
( A(x )ξ ,ξ ) ≥ k 0 ξ 2 , minden ξ ∈ R n -re és minden
x ∈ Ω -ra. [3]
A Poisson-egyenletben szereplő ∆ differenciáloperátor esetén az együtthatómátrix: 1 0 = E , A = 0 1
így ( A( x)ξ , ξ ) = (ξ ,ξ ) = ξ , ami k 0 =1 mellett éppen az operátor elliptikusságát eredményezi. 2
Definíció: Azt mondjuk, hogy a ∆u + f = 0
x∈Ω
u ( x) = g ( x)
ha x ∈ Γ
6
peremérték
feladatnak
van
klasszikus
megoldása,
ha
létezik
olyan
( )
u ∈ C 2 (Ω ) ∩ C Ω függvény, mely kielégíti az egyenletet és a hozzá tartozó
peremfeltételt. [3]
( )
( )
Megjegyzés: A fenti definícióban az u ∈ C 2 (Ω ) ∩ C Ω követelmény u ∈ C 2 Ω helyett azt fejezi ki, hogy a második deriváltakra csak a tartomány belsejében van szükség, míg a
( )
peremfeltételek miatt u ∈ C Ω -ra számítunk.
A variációs feladat A homogén Poisson-egyenlet fenti megfogalmazásáról áttérünk az úgynevezett gyenge vagy variációs egyenletre. Ehhez szükséges a parciális integrálás tétele (kétdimenzióban): ∂u
∫ ∂x
Ω
i
v dx = ∫ u ⋅ v ⋅ ni ds − ∫ u Γ
Ω
∂v dx ∂xi
ahol Ω ⊂ R 2 tartomány, u , v : Ω → R , Γ a tartomány pereme, n a perem kifelé mutató normális vektora: n = (n1 , n2 ) és x = ( x1 , x 2 ) ∈ R 2 Tekintsük a Poisson-egyenlet homogén peremfeltételek mellett: − ∆u = f
x∈Ω
u ( x) = 0
ha x ∈ Γ
A
∂ 2u ∂ 2u − 2 + 2 = f ( x) ∂x1 ∂x 2 egyenletet szorozzuk meg egy olyan v(x) úgynevezett tesztfüggvénnyel, mely az Ω tartományon folytonosan differenciálható és kielégíti a homogén peremfeltételeket (azaz v( x) = 0, ha x ∈ Γ teljesül), majd a kapott
∂ 2u ∂ 2u − 2 + 2 ⋅ v( x) = f ( x)v( x) ∂x1 ∂x 2 egyenletet integráljuk Ω felett:
∂ 2u ∂ 2u − ∫ 2 + 2 ⋅ v( x) dx = ∫ f ( x)v( x) dx ∂x 2 Ω ∂x1 Ω
7
Az egyenlet bal oldala:
∂ 2u ∂ 2u ∂ 2u ∂ 2u − ∫ 2 + 2 ⋅ v( x) dx = − ∫ 2 v( x) dx + ∫ 2 v( x) dx ∂x 2 Ω ∂x1 Ω ∂x 2 Ω ∂x1 A parciális integrálás tételét alkalmazva: ∂ 2u ∂u ∂u ∂v ∫Ω ∂x12 v( x) dx = ∫Γ ∂x1 ⋅ v( x) ⋅ n1 ds − Ω∫ ∂x1 ⋅ ∂x1 dx ∂ 2u ∂u ∂u ∂v ∫Ω ∂x22 v( x) dx = ∫Γ ∂x2 ⋅ v( x) ⋅ n2 ds − Ω∫ ∂x2 ⋅ ∂x 2 dx
Mivel v( x) = 0, ha x ∈ Γ , a peremen vett integrálok eltűnnek: ∂ 2u ∂u ∂v ∫Ω ∂x12 v( x) dx = − Ω∫ ∂x1 ⋅ ∂x1 dx ∂ 2u ∂u ∂v ∫Ω ∂x22 v( x) dx = − Ω∫ ∂x2 ⋅ ∂x 2 dx
Így
∂ 2u ∂ 2u ∂u ∂v ∂u ∂v − ∫ 2 + 2 ⋅ v( x) dx = ∫ ⋅ + ⋅ dx ∂x1 ∂x1 ∂x 2 ∂x 2 ∂x 2 Ω ∂x1 Ω Ezért az egyenlet a következő alakot ölti: ∂u
∫ ∂x
Ω
1
⋅
∂v ∂u ∂v + ⋅ dx = ∫ f ( x)v( x) dx ∂x1 ∂x 2 ∂x 2 Ω
Vezessük be a következő jelöléseket: ∂u ∂v ∂u ∂v ⋅ + ⋅ dx ∂x1 ∂x1 ∂x 2 ∂x 2 Ω
a (u , v) := ∫
ϕ (v) := ∫ f ( x)v( x) dx Ω
Egy alkalmasan kiválasztott V Hilbert-térben a fenti jelölések mellett a peremérték feladatról áttérhetünk a következő egyenletre: Keresünk olyan u ∈ V -t, hogy a (u , v) = ϕ (v) , minden v ∈ V esetén
Ezt nevezzük gyenge vagy variációs megfogalmazásnak.
8
Szoboljev-terek L p -terek: Legyen 1 ≤ p < ∞ valós szám, L p ( X , A , µ ) az X mérhető halmazon a µ mértékre
nézve p-edik hatványon integrálható függvények halmaza. f , g ∈ L p esetén f = g akkor és csak akkor, ha f ( x) = g ( x) majdnem minden x ∈ X -re. Az L p -beli norma: ha f ∈ L p , akkor 1
f
p
p p = ∫ f dµ X
A p = 2 esetben L2 a négyzetesen integrálható függvények tere, mely az
( f , g ) = ∫ f ⋅ g dµ X
skaláris szorzattal Hilbert-tér.
Hilbert-tér: Egy X halmazt Hilbert-térnek nevezünk, ha az X belsőszorzattal ellátott lineáris tér teljes (a belsőszorzatból származó norma által indukált metrikával). Legyen Ω ⊂ R n olyan nyílt összefüggő halmaz, melynek Γ = ∂Ω határa szakaszonként sima. A Soboljev-terek az L2 (Ω, λ ) terekből származtathatók (itt λ a Lebesgue-mérték).
L2 (Ω, λ ) az Ω fölött négyzetesen integrálható valós függvények tere. Két függvényt,
f , g ∈ L2 -t azonosnak tekintünk, ha egy Lebesgue-szerint nullmértékű
halmaztól eltekintve minden x ∈ X esetén f ( x) = g ( x) teljesül.
L2 (Ω, λ ) egy Hilbert-tér a következő skalárszorzattal
( f , g )0 := ( f , g )L
2
= ∫ f ( x) ⋅ g ( x) d x Ω
és a megfelelő, belsőszorzatból származó normával: u
0
=
(u,u )0
.
9
Definíció: Azt mondjuk, hogy az f ∈ L2 (Ω ) függvény α-adrendű általánosított deriváltja, vagy Szoboljev-féle deriváltja a g = ∂ α f , g ∈ L2 (Ω ) függvény, ha teljesül
(ϕ , g )0 = (− 1) α (∂ α ϕ , f )0
∀ϕ ∈ C0∞ (Ω) . [5]
Megjegyzés: 1.
Itt C ∞ (Ω ) jelöli az akárhányszor differenciálható függvények terét, C 0∞ (Ω ) pedig azon C ∞ (Ω ) -beli függvények, melyek csak Ω egy kompakt részhalmazán nem tűnnek el.
2.
Ha egy függvény differenciálható a klasszikus értelemben, akkor létezik a Szoboljev-deriváltja is, és a két derivált megegyezik.
Definíció: Legyen m ≥ 0 egész szám adott, és jelölje H m (Ω) azon f ∈ L2 (Ω ) függvények halmazát, melyeknek létezik a ∂ α f Szoboljev-féle derivátja minden α ≤ m esetén. H m (Ω) -n a következő skalárszorzatot definiálhatjuk:
( f , g )m := ∑ (∂ α f , ∂ α g )0 α ≤m
amelyből a következő norma származik: f
m
:=
( f , f )m
=
∑ α
∂α f
≤m
2 L2 ( Ω )
.
Az u m :=
∑ α
=m
∂α f
2 L2 ( Ω )
függvény félnormát definiál H m (Ω) -n. H m (Ω) -t Szoboljev-térnek nevezzük. [6]
Megjegyzés: H m (Ω) teljes a ⋅
m
normával, így Hilbert-tér.
Definíció: Azon H m (Ω) -beli függvények halmazát, amelyeknek a nyoma 0 Γ -n H 0m (Ω) -val jelöljük.
Megjegyzés: Ha Ω korlátos, akkor ⋅
m
norma H 0m (Ω) -n (amely ekvivalens ⋅
m
-val).
10
A variációs feladat megoldhatósága A továbbiakban a H m (Ω) , illetve a H 0m (Ω) terekre lesz szükségünk. Az előző fejezetben leírtak alapján: u
0
=
∫u
2
dx ,
Ω
2
∂u ∂u + u 1 = ∫ ∂ x 1 ∂x 2 Ω
u1= Ekkor ⋅
1
2
u1 + u
normát definiál H 1 (Ω ) -n, míg ⋅
1
2 0
2
dx ,
.
félnorma. Ha Ω korlátos, akkor ⋅
1
norma
H 01 (Ω) -n. A variációs feladatban válasszuk a V = H 01 (Ω ) Hilbert-teret. Így olyan u ∈ H 01 ( Ω ) függvényt keresünk, melyre
a (u , v) = ϕ (v) teljesül ∀v ∈ H 01 esetén. A variációs egyenletben szereplő a(u ,v ) kifejezés V × V -n definiált funkcionál:
a : H 01 × H 01 → R .
Az a bilineáris:
− mindkét változójában additív és homogén: ∂(u + α ⋅ v ) ∂w ∂ (u + α ⋅ v ) ∂w ⋅ + ⋅ dx = ∂x1 ∂x1 ∂x 2 ∂x 2 Ω
a (u + α ⋅ v, w) = ∫
∂u ∂v ∂w ⋅ = ∫ +α ∂x1 ∂x1 ∂x1 Ω
∂u ∂v ∂w ⋅ + +α dx = ∂x 2 ∂x 2 ∂x 2
∂v ∂w ∂v ∂w ∂u ∂w ∂u ∂w + α dx = ⋅ + ⋅ ⋅ + ⋅ = ∫ ∂x1 ∂x1 ∂x 2 ∂x 2 ∂x 2 ∂x 2 ∂x 2 ∂x 2 Ω = a(u , w) + α ⋅ a(v , w) a(u ,v + α ⋅ w) = a(u ,v ) + α ⋅ a(u , w) minden α ∈ R -re, és minden u , v, w ∈ H 01 -re.
11
Az a korlátos (vagyis folytonos): létezik olyan M ∈ R , hogy
a(u, v ) ≤ M u 1 v 1 Az a szimmetrikus:
∂u ∂v ∂u ∂v ∂v ∂u ∂v ∂u ⋅ + ⋅ dx = ∫ ⋅ + ⋅ dx = a (v, u ) ∂x1 ∂x1 ∂x 2 ∂x 2 ∂x1 ∂x1 ∂x 2 ∂x 2 Ω Ω
a (u , v ) = ∫
és V-elliptikus: létezik olyan m > 0 konstans, hogy 2
2
∂u ∂u + dx ≥ m v 1 minden v ∈ H 01 esetén. a (u ,u ) = ∫ ∂x1 ∂x 2 Ω A ϕ (v ) kifejezés lineáris funkcionál: −
ϕ additív és homogén a v argumentumában:
ϕ (v + α ⋅ w) = ∫ f ( x ) ⋅ [v( x ) + α ⋅ w( x )] dx = ∫ f ( x )v( x ) dx + α ⋅ ∫ f ( x )w(x ) dx = Ω
Ω
Ω
= ϕ (v ) + α ⋅ ϕ (w) minden v , w ∈ H 01 -re és α ∈ R -re. A ϕ korlátos:
ϕ (v ) ≤ f
0
v
0
Tétel: (variációs feladat megoldásának egzisztenciája, unicitása, stabilitása)
Legyen ⋅
V
=
az
( ⋅ ,⋅ )V
a (u , v) = ϕ (v)
( ∀ v ∈V )
variációs
feladatban
V
Hilbert-tér
a
normával, az a bilineáris forma szimmetrikus és V-elliptikus, ϕ
lineáris funkcionál V-n. Ekkor a variációs feladatnak pontosan egy általánosított megoldása van, mely eleget tesz a következő stabilitási becslésnek: u
v
≤
ϕ m
. [2]
12
Véges elemek A variációs feladat numerikus megoldása során az Ω ⊂ R 2 tartományt fel fogjuk osztani véges sok résztartományra, és minden résztartomány fölött egy adott fokszámú polinommal közelítjük a keresett u : Ω → R 2 függvényt. Definíció: Legyen Th = {Ω i , i = 1,..., M } az Ω tartomány egy felosztása háromszögekre vagy
négyszögekre. A Th felosztást Ω egy triangulációjának hívjuk, ha teljesülnek a következő feltételek: M
a) Ω = ∪ Ω i i =1
b) ha Ω i ∩ Ω j pontosan egy pontot tartalmaz, akkor az Ω i és Ω j közös csúcspontja, c) ha i ≠ j és Ω i ∩ Ω j egynél több pontot tartalmaz, akkor Ω i ∩ Ω j az Ω i és Ω j közös oldala. [6] Megjegyzés:
1. Ilyen felosztás akkor létezik, ha az Ω ⊂ R 2 tartomány Γ pereme töröttvonal. 2. A Th = {Ω i , i = 1,..., M } felosztást triangulációnak nevezzük, még akkor is, ha Ω i -k nem háromszögek. A továbbiakban csak háromszöges elemekkel foglakozunk. A fenti definíció kizárja, hogy a felosztás során valamelyik háromszög csúcspontja egy másik háromszög oldalának belső pontja legyen:
Miután az Ω poligonális tartományt felosztottuk háromszögekre, minden háromszög fölött egy adott fokszámú – esetünkben másodfokú – polinommal fogjuk az u függvényt közelíteni.
13
Egy kétváltozós k-adfokú polinomnak
(k + 1) ⋅ (k + 2) 2
hogy egy háromszög fölött alkalmasan megválasztott
szabad paramétere van. Ez azt jelenti,
(k + 1) ⋅ (k + 2) 2
darab pontban előírva
egy k-adfokú polinom értékét az egyértelműen megadható. Egy ilyen alkalmas pontrendszert kapunk, ha a háromszög minden oldalát k + 1 ponttal k egyenlő részre osztjuk úgy, hogy a szélső pontok a háromszög csúcsai; majd a szomszédos oldalak megfelelő osztópontjait összekötjük. Az oldalakon felvett pontok, illetve a szakaszok metszéspontjai
(k + 1) ⋅ (k + 2) 2
darab pontot határoznak meg.
Például k = 4 esetén:
Ezen pontokban előírva egy kétváltozós k-adfokú polinom értékeit a polinom egyértelműen definiált a háromszög fölött. A polinom a háromszög tetszőleges oldalára leszűkítve egy egyváltozós, k-adfokú polinom lesz, amelyet az adott oldalon felvett (k + 1) darab pontban előírt értéke egyértelműen meghatároz. Ez azt jelenti, hogy két szomszédos háromszög közös oldalán ugyanazokat a helyettesítési értékeket előírva az egyes háromszögek fölött adott polinomok számára a két polinom meg fog egyezni a közös oldalon. Így elérhető, hogy a trianguláció háromszögein elemenként definiált polinom folytonos legyen az egész tartomány fölött. Definíció: Azokat a lineáris funkcionálokat, melyek megadják az egyes csomópontokban a
függvényértékeket, szabadsági fokoknak nevezzük. [3]
14
{
}
Az Ω i -hez tartozó szabadsági fokok halmaza: ψ i = ψ k , x k ∈ Ω i , illetve a lokális számozás
{
}
esetén: ψ i = ψ i ,1 ,...,ψ i , mi [3] Megjegyzés:
1.
Egy x k csomóponthoz hozzárendelt ψ k lineáris funkcionál, mint szabadsági fok esetén:
ψ k (u ) = u (x k ) 2.
Szabadsági fok lehet függvényérték és a deriváltak értékei is. Lagrange-típusú elemről beszélünk, ha a szabadsági fokok mind függvényértékek. Mi Lagrange-típusú elemekkel fogunk foglalkozni.
Minden ψ k szabadsági fokhoz rendeljünk hozzá egy wk ∈ V függvényt a következők szerint: −
wk leszűkítése Ω i -re legyen pi ,k polinom
− minden egyes Ω i -n P (Ω i ) a {pi ,k } polinomok lineáris halmaza, és ha az Ω i -hez tartozó csomópontok halmaza
{x
i ,1
} {
,..., x i ,mi = x k1 ,..., x
k mi
}
akkor a mi
pi ( x ) = ∑ p i ,k j ( x )α ij j =1
által definiált P (Ω i ) -beli polinom α ij együtthatói egyértelműen meg legyenek határozva − ha x k1 ,..., x k s az Ω i egy tetszőleges rögzített oldalán lévő csomópontok, akkor a s
qi ( x ) = ∑ pi ,k j ( x )α ij j =1
{
}
polinom egyértelműen meg legyen határozva Ω i oldalán a ψ k1 (q i ),...,ψ k s (qi ) értékek megadása után. Ha ψ m (wk ) = δ mk , akkor a Vh végeselem tér bázisa éppen {w1 ,..., wl }. Definíció: Véges elemnek nevezünk egy {Ωi ,ψ i , P (Ωi )} hármast. [3]
15
A T6 elem Osszuk fel az Ω ⊂ R 2 tartományt véges sok Ω i = ∆ i háromszögre, úgy hogy a háromszögek teljesítsék a trianguláció definíciójában szereplő feltételeket. A csomópontok a háromszög csúcsai: x i ,1 , x i ,2 , x i ,3 , valamint az oldalfelező pontok: x i ,4 =
(
)
(
)
(
1 i ,1 1 1 x + x i ,2 , x i ,5 = x i ,2 + x i ,3 , x i ,6 = x i ,1 + x i ,3 2 2 2
)
P (Ω i ) = P2 (∆ i ) , a ∆ i felett definiált másodfokú polinomok tere ( azaz a 6 szabadsági fokkal
( )
rendelkező, a i ,1 + a i ,2 x1 + a i ,3 x 2 + ai ,4 x1 x 2 + a i ,5 x12 + ai ,6 x 22 = u x i alakú polinomok). A lineáris funkcionálok: ψ i = ψ i,j ψ i,j ( p) = p xi,j ,
{
( )
}
p ∈ P2 (∆ i ), 0 ≤ j ≤ 6 .
Egy ilyen másodrendű háromszögelemet T6-tal jelölünk.
A wl bázisfüggvények leszűkítése ∆ i -re, azaz a pi ,l polinomok a legegyszerűbben a háromszög λ1 , λ 2 , λ3 baricentrikus koordinátáinak segítségével adható meg. Definíció: Legyen ∆ egy háromszög, jelölje S i = ( xi , y i ) , i = 1, 2, 3 a háromszög csúcsait.
Legyen továbbá P = ( x, y ) a háromszög egy pontja. Ha λ1 , λ 2 , λ3 -ra
x1 y1 1
x2 y2 1
x3 λ1 x y 3 λ 2 = y 1 λ3 1
teljesül, akkor λ1 , λ 2 , λ3 -at a P pont baricentrikus koordinátáinak nevezzük a ∆ háromszögben. [1]
16
A defínicióban szereplő mátrix determinánsa egyenlő a ∆ háromszög területének kétszeresével: ∆ =
1 (x1 y 2 − y1 x2 + x2 y3 − y 2 x3 + x3 y1 − y3 x1 ) 2
x1
x2
x3
y1 1
y2 1
y 3 = x1 y 2 + x 2 y 3 + x3 y1 − y 2 x3 − y 3 x1 − y1 x 2 = 2 ∆ 1
Megjegyzés: A baricentrikus koordináták használatával a háromszög fölött definiált
polinomok a háromszögek helyzetétől függetlenül írhatók föl. A fentieknek megfelelően a T6 elem lokális bázisfüggvényei: x i ,1 -hez p1 = λ1 (2λ1 − 1) x i ,2 -höz p 2 = λ 2 (2λ 2 − 1) x i ,3 -hoz p3 = λ3 (2λ3 − 1) x i ,4 -hez p 4 = 4λ1λ 2 x i ,5 -höz p5 = 4λ 2 λ3 x i ,6 -hoz p 6 = 4λ1λ3 tartozik. Megjegyzés:
1.
Ezen függvényekre teljesül, hogy a nekik megfelelő x i , j pontban értékük 1, míg a többi x i , j pontban nulla.
2.
A Vh végeselem tér wl bázisfüggvényei folytonosak a szomszédos háromszögek közös oldalán keresztül is, hiszen Vh ⊂ H 01 (Ω ) .
17
A ∆ háromszög fölött definiált polinomok integráljainak kiszámítását a következő Lemma írja le. Lemma: ( [8] ) Legyen ∆ ∈ Th egy tetszőleges háromszög, és jelölje λ1 , λ 2 , λ3 a ∆-beli
baricentrikus koordinátákat. Ekkor
∫λ λ λ k 1
∆
l 2
m 3
dxdy = 2 ∆
k! l! m! . (k + l + m + 2)!
18
A diszkrét egyenlet A variációs feladat közelítő megoldásához a V = H 01 helyett, annak egy Vh ⊂ V véges n-dimenziós alterét választjuk. Vh egy bázisa legyen {v1 ,..., v n } , ekkor u h ∈ Vh pontosan akkor, ha: n
u h = ∑ y i vi i =1
teljesül, ahol y1 ,..., y n keresett. Ezek után a variációs feladat diszkrét megfogalmazása a következő: Olyan u h ∈ Vh függvényt keresünk, melyre a (u h , v h ) = ϕ (v h ) teljesül minden v h ∈ Vh esetén. Megjegyzés: Miután Vh a V Hilbert-tér altere, maga is Hilbert-tér a ⋅
V
normával, így az a
és ϕ funkcionálok megtartják V -beli tulajdonságaikat Vh -ban. Ezért a variációs feladat megoldására vonatkozó tétel érvényes Vh -ban, azaz a közelítő u h ∈ Vh megoldás létezik, egyértelmű és stabil. A diszkrét feladat megoldása legyen u h , a variációs feladaté pedig u . Ekkor minden v h ∈ Vh ⊂ V esetén:
a(u ,v h ) = ϕ (v h ) a(u h ,v h ) = ϕ (v h ) azaz a (u − u h , v h ) = 0 . Ez azt jelenti, hogy az u − u h hiba ortogonális a Vh altérre az
a( ⋅ , ⋅ ) skalárszorzatban,
pontosabban az u h közelítő megoldás az u általánosított megoldás legjobb a-normájú közelítése Vh -ban.
19
Megjegyzés: Egy Hilbert-térben egyértelműen létezik a legjobb közelítés, és annak
meghatározása a fent említett ortogonalitási tulajdonsággal történhet. Az a-norma és a V-norma ekvivalensek, de általában a ( ⋅ , ⋅ )V skalárszorzatban nincs szó ortogonalitásról. Lemma: (Céa-lemma)
Legyen az a(u, v ) V × V -n definiált bilineáris funkcionál V-ben pozitív definit (azaz szimmetrikus és V-elliptikus). Ekkor a végeselem módszer hibája egy konstans szorzótól eltekintve ugyanakkora, mint a legjobb közelítésé. [2] A numerikus megoldáshoz a {v1 ,..., v n } Vh -beli bázis segítségével a a (u h , v h ) = ϕ (v h ) (minden v h ∈ Vh esetén) diszkrét feladat átfogalmazható a következővé: n
a (∑ y i vi , v j ) = ϕ (v j ) , ahol j = 1,..., n i =1
amiből n
∑y
i
⋅ a (vi , v j ) = ϕ (v j ) , ahol j = 1,..., n
i =1
azaz az Ay = b , y1 ,..., y n ismeretlenekre felírt lineáris egyenletrendszerhez jutunk, ahol
y1 y= M y n
ϕ (v1 ) b= M ϕ (v ) n
A = (a i j )
a ij = a (vi , v j )
Megjegyzés: A lineáris egyenletrendszer A = (a i j ) mátrixa szimmetrikus és pozitív definit,
amely tulajdonságokat az a bilineáris funkcionáltól – és így az eredeti feladattól – örököl. [2]
20
A feladat végeselem megoldása Tekintsük a − ∆u = f
x∈Ω
u(x ) = 0
ha x ∈ Γ
homogén Poisson-egyenletet az egységnégyzeten. A gyakorlati megvalósításhoz: − Kiválasztjuk a Vh tér háromszögenként másodfokú polinomokból álló bázisát − Felépítjük az A mátrixot − A jobboldali vektor komponenseit kiszámítjuk numerikus integrálással − Megoldjuk a keletkezett lineáris egyenletrendszert − A kapott y megoldási vektor szolgáltatja a diszkrét feladat közelítő megoldását. Végezzük el a tartomány diszkretizációját a trianguláció definíciójában szereplő feltételeknek eleget tevő T6 háromszög elemekkel. Adjuk meg a T6 elemhez tartozó bázist. Ehhez vizsgáljuk meg a B lineáris leképezést, mely a következő ∆ 0 , úgynevezett standard háromszöget transzformálja tetszőleges ∆ háromszöggé:
Ekkor ξ x B = η y ahol
x − x1 B = 2 y 2 − y1
x3 − x1 y 3 − y1
21
B -1 =
− ( x3 − x1 ) 1 y 3 − y1 x 2 − x1 det (B) − ( y 2 − y1 )
ahol det (B) = 2 ∆ , a háromszög területének kétszerese. Az A mátrix elemeinek kiszámításához szükséges az integráltranszformációs formula:
∂v(ξ , η ) ∂ξ ∂v(ξ , η ) ∂η ∂u (x, y ) q( x, y ) dxdy = det (B) ∫ + ⋅ p(ξ , η ) dξdη ∂x ∂ξ ∂x ∂η ∂x ∆ ∆0
∫ illetve
∂v(ξ , η ) ∂ξ ∂v(ξ , η ) ∂η ∂u ( x, y ) ⋅ p(ξ , η ) dξdη , q( x, y ) dxdy = det (B) ∫ + ξ η ∂ y ∂ ∂ y ∂ ∂ y ∆ ∆0
∫ ahol
v(ξ , η ) = u ( x(ξ , η ), y (ξ , η )) p(ξ , η ) = q( x(ξ , η ), y (ξ , η )) . Az egységnégyzet standard triangulációja: Felosztjuk az egységnégyzetet N 2 darab egybevágó kisebb négyzetre, majd minden kapott négyzetet az DNY-ÉK irányú átlójával két háromszögre:
N=4 esetén a standard trianguláció Ekkor kétféle háromszöggel van dolgunk: 1 1 Az első típus a (0, 0) , 0, , − , 0 csúcspontú háromszögből eltolással kapható N N háromszögek.
22
Itt a B leképezés mátrixa:
0 B= 1 N
−
1 N = 1 ⋅ 0 − 1 0 N 1 0
Ekkor det (B) =
1 N2
és 0 1 B −1 = N ⋅ −1 0 1 1 A második típusú háromszögek a (0, 0) , , 0 , 0, − csúcspontú háromszögből N N eltolással kapható háromszögek.
Ebben az esetben a B leképezés mátrixa:
0 B= − 1 N
1 N = 1 ⋅ 0 1 0 N − 1 0
23
továbbá det (B) =
1 N2
és 0 − 1 B −1 = N ⋅ 1 0 Ha ∆ egy első típusú háromszög, akkor az integráltranszformációs formula alapján:
∂u ∂u j ∂u i ∂u j dx1 dx 2 = a (u i , u j ) = ∫ i + ∂ x ∂ x ∂ x ∂ x 1 1 2 2 ∆ ∂v ∂ξ ∂vi ∂η ∂v j ∂ξ ∂v j ∂η + ⋅ = det (B) ∫ i + + ξ η ξ η ∂ ∂ x ∂ ∂ x ∂ ∂ x ∂ ∂ x 1 1 1 1 ∆0
∂v ∂ξ ∂vi ∂η ∂v j ∂ξ ∂v j ∂η dξdη = ⋅ + i + + ∂ξ ∂x 2 ∂η ∂x 2 ∂ξ ∂x 2 ∂η ∂x 2 =
1 N2
∂vi
∂vi
∂v j
∫ ∂ξ ⋅ 0 + ∂η ⋅ (− N ) ⋅ ∂ξ
⋅0+
∆0
∂v j
⋅ (− N ) + ∂η
∂v j ∂v ∂v ∂v j + i ⋅ N + i ⋅ 0 ⋅ ⋅N+ ⋅ 0 dξdη = ∂η ∂ξ ∂η ∂ξ =
1 N2
∫N
∆0
2
∂v ∂v j ∂vi ∂v j ∂vi ∂v j ∂v ∂v j dξdη . + + N2 i dξdη = ∫ i ∂η ∂η ∂ξ ∂ξ ∂ξ ∂ξ ∂η ∂η ∆0
Ha ∆ egy második típusú háromszög, akkor:
∂u ∂u j ∂u i ∂u j dx1 dx 2 = a (u i , u j ) = ∫ i + ∂x1 ∂x1 ∂x 2 ∂x 2 ∆ ∂v ∂ξ ∂vi ∂η ∂v j ∂ξ ∂v j ∂η + ⋅ + + = det (B) ∫ i ∂ξ ∂x1 ∂η ∂x1 ∂ξ ∂x1 ∂η ∂x1 ∆0
∂v ∂ξ ∂vi ∂η ∂v j ∂ξ ∂v j ∂η dξdη = ⋅ + i + + ξ η ξ η ∂ ∂ x ∂ ∂ x ∂ ∂ x ∂ ∂ x 2 2 2 2 =
1 N2
∂vi
∂vi
∂v j
∫ ∂ξ ⋅ 0 + ∂η ⋅ N ⋅ ∂ξ
∆0
⋅0 +
∂v j
⋅ N + ∂η
24
∂v j ∂v ∂v ∂v j + i ⋅ (− N ) + i ⋅ 0 ⋅ ⋅ (− N ) + ⋅ 0 dξdη = ∂η ∂ξ ∂η ∂ξ =
1 N2
2 ∫N
∆0
∂v ∂v j ∂vi ∂v j ∂vi ∂v j ∂v ∂v j dξdη . dξdη = ∫ i + N2 i + ∂η ∂η ∂ξ ∂ξ ∂ ∂ ∂ ∂ ξ ξ η η ∆0
Ez azt jelenti, hogy a lineáris egyenletrendszer A alapmátrixának felírásához azt az Aloc lokális mátrixot kell meghatározni, melynek a ij elemei ∂p ∂p j ∂pi ∂p j dξ dη ( i, j = 1, K , 6 ), a ij = ∫ i + ∂ ∂ ∂ ∂ ξ ξ η η ∆0 ahol a p i ( i = 1, K , 6 ) függvények a referencia háromszög fölött definiált másodfokú polinomok. A ∆ 0 referencia háromszög baricentrikus koordinátái:
0 1 0 λ1 ξ 0 0 1 λ 2 = η , azaz λ1 = 1 − ξ − η , λ 2 = ξ , λ3 = η . 1 1 1 λ 1 3 Így a ∆ 0 fölött definiált polinomok:
p1 = λ1 (2λ1 − 1) = (1 − ξ − η )(1 − 2ξ − 2η ) p 2 = λ 2 (2λ 2 − 1) = ξ (2ξ − 1) p3 = λ3 (2λ3 − 1) = η (2η − 1)
p 4 = 4λ1λ 2 = 4(1 − ξ − η )ξ p5 = 4λ 2 λ3 = 4ξη p 6 = 4λ1λ3 = 4(1 − ξ − η )η
25
A referencia háromszög fölött adott bázisfüggvények parciális deriváltjai:
[
]
[
]
∂p1 ∂[(1 − ξ − η )(1 − 2ξ − 2η )] ∂ 1 − 3ξ − 3η + 2ξ 2 + 4ξη + 2η 2 = = = ∂ξ ∂ξ ∂ξ = 4ξ + 4η − 3 = 4λ 2 + 4λ3 − 3
∂p1 ∂[(1 − ξ − η )(1 − 2ξ − 2η )] ∂ 1 − 3ξ − 3η + 2ξ 2 + 4ξη + 2η 2 = = = ∂η ∂η ∂η = 4ξ + 4η − 3 = 4λ 2 + 4λ3 − 3
[
]
[
]
∂p 2 ∂[ξ (2ξ − 1)] ∂ 2ξ 2 − ξ = = = 4ξ − 1 = 4λ 2 − 1 ∂ξ ∂ξ ∂ξ ∂p 2 ∂[ξ (2ξ − 1)] = =0 ∂η ∂η ∂p 3 ∂[η (2η - 1)] = =0 ∂ξ ∂ξ
∂p3 ∂[η (2η − 1)] ∂ 2η 2 − η = = = 4η − 1 = 4λ3 − 1 ∂η ∂η ∂η
[
]
[
]
[
]
[
]
∂p 4 ∂[4(1 − ξ − η ) ξ ] ∂ 4ξ − 4ξ 2 − 4ξη = = = 4 − 8ξ − 4η = 4λ1 − 4λ 2 ∂ξ ∂ξ ∂ξ ∂p 4 ∂[4(1 − ξ − η ) ξ ] ∂ 4ξ − 4ξ 2 − 4ξη = = = −4ξ = −4λ 2 ∂η ∂η ∂η ∂p 5 ∂[4ξη ] = = 4η = 4λ3 ∂ξ ∂ξ ∂p 5 ∂[4ξη ] = = 4ξ = 4λ 2 ∂η ∂η
∂p 6 ∂[4(1 − ξ − η )η ] ∂ 4η − 4ξη − 4η 2 = = = −4η = −4λ3 ∂ξ ∂ξ ∂ξ ∂p 6 ∂[4(1 − ξ − η )η ] ∂ 4η − 4ξη − 4η 2 = = = 4 − 4ξ − 8η = 4λ1 − 4λ3 ∂η ∂η ∂η Az Aloc lokális mátrix szimmetrikus, így elegendő a felső háromszöghöz tartozó a ij elemeket kiszámolni: a ij = a ( p i , p j ) ( i, j = 1, K , 6 )
26
a11 =
∂p1 ∂p1 ∂p1 ∂p1 2 2 + dξdη = ∫ (4λ 2 + 4λ3 − 3) + (4λ 2 + 4λ3 − 3) dξdη = ∂ξ ∂η ∂η ∆0 ∆0
∫ ∂ξ
= ∫ 32λ22 + 32λ32 + 64λ 2 λ3 − 48λ 2 − 48λ3 + 18 dξdη = ∆0
2! 2! 1!1! 1! 1! 0! = 2 ⋅ ∆ 0 32 + 32 + 64 − 48 − 48 + 18 = 2! (2 + 2)! (1 + 1 + 2)! (1 + 2)! (1 + 2)! (2 + 2 )! 1 3 = 2⋅ ⋅ =1 2 3 a12 =
∂p1 ∂p 2 ∂p1 ∂p 2 + dξdη = ∫ (4λ 2 + 4λ3 − 3) ⋅ (4λ 2 − 1) +(4λ 2 + 4λ3 − 3) ⋅ 0 dξdη = ∂ ∂ ∂ ξ η η ∆0 ∆0
∫ ∂ξ
= ∫ 16λ22 − 16λ 2 + 16λ 2 λ3 − 4λ3 + 3 dξdη = ∆0
2! 1! 1!1! 1! 0! 1 = = 2 ⋅ ∆ 0 16 ⋅ − 16 ⋅ + 16 ⋅ − 4⋅ + 3⋅ (2 + 2)! (1 + 2)! (1 + 1 + 2)! (1 + 2)! (0 + 2 )! 6 a13 =
∫ (4λ
2
+ 4λ3 − 3) ⋅ 0 + (4λ 2 + 4λ3 − 3) ⋅ (4λ3 − 1) dξdη =
∆0
= ∫ 16 λ32 − 16 λ3 + 16 λ2 λ3 − 4 λ2 + 3 dξ d = ∆0
2! 1! 1!1! 1! 0! 1 = = 2 ⋅ ∆ 0 16 ⋅ − 16 ⋅ + 16 ⋅ − 4⋅ + 3⋅ (2 + 2)! (1 + 2)! (1 + 1 + 2)! (1 + 2)! (0 + 2 )! 6 a14 =
∫ (4λ
2
+ 4λ3 − 3) ⋅ (4λ1 − 4λ 2 ) + (4λ 2 + 4λ3 − 3) ⋅ (− 4λ 2 ) dξdη =
∆0
= ∫ 16λ1λ 2 + 16λ 1 λ3 − 32λ 2 λ3 − 32λ22 − 12λ1 + 24λ 2 dξdη = ∆0
1!1! 1!1! 1!1! 2! = 2 ⋅ ∆ 0 16 ⋅ + 16 ⋅ − 32 ⋅ − 32 ⋅ − (1 + 1 + 2)! (1 + 1 + 2)! (1 + 1 + 2)! (2 + 2)! − 12 ⋅ a15 =
1! 1! 2 + 24 ⋅ =− (1 + 2)! (1 + 2)! 3
∫ (4λ
2
+ 4λ3 − 3) ⋅ 4λ3 + (4λ 2 + 4λ3 − 3) ⋅ 4λ 2 dξdη =
∆0
= ∫ 16λ22 + 16λ32 + 32λ 2 λ3 − 12λ 2 −12λ3 dξdη = ∆0
27
2! 2! 1!1! 1! 1! =0 = 2 ⋅ ∆ 0 16 ⋅ + 16 ⋅ + 32 ⋅ − 12 ⋅ − 12 ⋅ (2 + 2)! (2 + 2)! (1 + 1 + 2)! (1 + 2)! (1 + 2)! a16 =
∫ (4λ
2
+ 4λ3 − 3) ⋅ (− 4λ3 ) + (4λ 2 + 4λ3 − 3) ⋅ (4λ1 − 4λ3 ) dξdη =
∆0
= ∫ (−32λ32 + 16λ1λ 2 + 16λ1λ3 − 32λ 2 λ3 − 12λ1 + 24λ3 ) dξdη = ∆0
2! 1!1! 1!1! = 2 ⋅ ∆ 0 − 32 ⋅ + 16 ⋅ + 16 ⋅ − (2 + 2)! (1 + 1 + 2)! (1 + 1 + 2)! − 32 ⋅ a 22 =
1!1! 1! 1! 2 = − − 12 ⋅ + 24 ⋅ (1 + 1 + 2)! (1 + 2)! (1 + 2)! 3
∫ (4λ
− 1) + 0 2 dξdη = ∫ 16λ22 − 8λ 2 + 1 dξdη = 2
2
∆0
∆0
2! 1! 0! 1 = = 2 ⋅ ∆ 0 16 ⋅ −8⋅ + (2 + 2)! (1 + 2)! (0 + 2)! 2 a 23 =
∫ (4λ
1
− 1) ⋅ 0 + 0 ⋅ (4λ3 − 1) dξdη = 0
∆0
a 24 =
∫ (4λ
2
∆0
− 1) ⋅ (4λ1 − 4λ 2 ) + 0 ⋅ (− 4λ 2 ) dξdη = ∫ 16λ1λ 2 − 16λ22 − 4λ1 + 4λ 2 dξ dη = ∆0
1!1! 2! 1! 1! 2 = − = 2 ⋅ ∆ 0 16 ⋅ − 16 ⋅ − 4⋅ + 4⋅ (1 + 1 + 2)! (2 + 2)! (1 + 2)! (1 + 2)! 3 a 25 =
∫ (4λ
2
∆0
− 1) ⋅ 4λ3 + 0 ⋅ 4λ 2 dξdη = ∫ 16λ 2 λ3 − 4λ3 dξdη = ∆0
1!1! 1! =0 = 2 ⋅ ∆ 0 16 ⋅ − 4⋅ (1 + 1 + 2)! (1 + 2)! a 26 =
∫ (4λ
2
− 1)(− 4λ3 ) + 0 ⋅ (4λ1 − 4λ3 ) dξdη =
∆0
∫ (− 16λ λ 2
3
+ 4λ3 ) dξdη =
∆0
1!1! 1! =0 = 2 ⋅ ∆ 0 − 16 ⋅ + 4⋅ (1 + 1 + 2)! (1 + 2 )! a 33 = ∫ 0 2 + (4λ3 − 1) dξdη = ∫ 16λ32 − 8λ3 + 1 dξdη = 2
∆0
∆0
2! 1! 0! 1 = = 2 ⋅ ∆ 0 16 ⋅ −8⋅ + (2 + 2)! (1 + 2)! (0 + 2)! 2
28
a 34 = ∫ 0 ⋅ (4λ1 − 4λ 2 ) + (4λ3 − 1) ⋅ (− 4λ 2 ) dξdη = ∆0
∫ (− 16λ λ 2
3
+ 4λ 2 ) dξdη =
∆0
1!1! 1! = 0 = 2 ⋅ ∆ 0 − 16 ⋅ + 4⋅ ( 1 + 1 + 2 ) ! ( 1 + 2 ) ! a 35 = ∫ 0 ⋅ 4λ3 + (4λ3 − 1) ⋅ 4λ 2 dξdη = ∫ 16λ 2 λ3 − 4λ 2 dξdη = ∆0
∆0
1!1! 1! =0 = 2 ⋅ ∆ 0 16 ⋅ − 4⋅ (1 + 1 + 2)! (1 + 2)! a 36 = ∫ 0 ⋅ (− 4λ3 ) + (4λ3 − 1) ⋅ (4λ1 − 4λ3 ) dξdη = ∫ 16λ1λ3 − 16λ32 + 4λ1 − 4λ3 dξdη = ∆0
∆0
1!1! 2! 1! 1! 2 = − = 2 ⋅ ∆ 0 16 ⋅ − 16 ⋅ + 4⋅ − 4⋅ (1 + 1 + 2)! (2 + 2)! (1 + 2)! (1 + 2)! 3 a 44 =
∫ (4λ
1
∆0
− 4λ 2 ) + (− 4λ 2 ) dξdη = ∫ 16λ12 + 32λ22 − 32λ1λ 2 dξdη = 2
2
∆0
2! 2! 1!1! 8 = = 2 ⋅ ∆ 0 16 ⋅ + 32 ⋅ − 32 ⋅ (2 + 2)! (2 + 2 )! (1 + 1 + 2)! 3 a 45 =
∫ (4λ
1
∆0
− 4λ 2 ) ⋅ 4λ3 + (− 4λ 2 ) ⋅ 4λ 2 dξdη = ∫ 16λ1λ3 − 16λ 2 λ3 − 16λ22 dξdη = ∆0
1!1! 1!1! 2! 4 = − = 2 ⋅ ∆ 0 16 ⋅ − 16 ⋅ − 16 ⋅ (1 + 1 + 2)! (1 + 1 + 2)! (2 + 2)! 3 a 46 =
∫ (4λ
1
∆0
− 4λ 2 ) ⋅ (− 4λ3 ) + (− 4λ 2 ) ⋅ (4λ1 − 4λ3 ) dξdη = ∫ 32λ 2 λ3 − 16λ1λ 2 − 16λ1λ3 dξdη = ∆0
1!1! 1!1! 1!1! =0 = 2 ⋅ ∆ 0 32 ⋅ − 16 ⋅ − 16 ⋅ (1 + 1 + 2)! (1 + 1 + 2)! (1 + 1 + 2)!
a55 =
∫ (4λ ) + (4λ ) 2
3
2
2
∆0
2! 2! 8 = dξdη = ∫ 16λ22 + 16λ32 dξ dη = 2 ⋅ ∆ 0 16 ⋅ + 16 ⋅ (2 + 2)! (2 + 2)! 3 ∆0
a 56 = ∫ 4λ3 ⋅ (− 4λ3 ) + 4λ 2 ⋅ (4λ1 − 4λ3 ) dξdη = ∫ 16λ1λ 2 − 16λ 2 λ3 − 16λ23 dξdη = ∆0
∆0
1!1! 1!1! 2! 4 = − = 2 ⋅ ∆ 0 16 ⋅ − 16 ⋅ − 16 ⋅ (1 + 1 + 2)! (1 + 1 + 2)! (2 + 2)! 3 a 66 =
∫ (− 4λ ) + (4λ 2
3
∆0
1
− 4λ3 ) dξdη = ∫ 16λ12 + 32λ32 − 32λ1λ3 dξdη = 2
∆0
29
2! 2! 1!1! 8 = = 2 ⋅ ∆ 0 16 ⋅ + 32 ⋅ − 32 ⋅ (2 + 2)! (2 + 2 )! (1 + 1 + 2)! 3
Aloc
1 1 6 1 6 = − 2 3 0 2 − 3
−
1 6
1 6
−
2 3
0
1 2
0
−
2 3
0
0
1 2
0
0
2 3
0
8 3
0
0
0
−
−
2 3
−
4 3 0
4 3 8 3
−
4 3
2 − 3 0 2 − 3 0 4 − 3 8 3
A lineáris egyenletrendszer felírásához szükséges
ϕ (v j ) = ∫ f ( x )v j ( x ) dx
( j = 1,K, n )
Ω
integrálok meghatározásához a következő kubatúra képletet használjuk:
∫
Ω
m
m
i =1 ∆ i
i =1
f ( x ) dx = ∑ ∫ f ( x ) dx ≈ ∑ ∆ i
f i1 + f i 2 + f i 3 , 3
ahol ∆ 1 , K , ∆ m az Ω egy triangulációja, f i1 , f i 2 , f i 3 az f függvény helyettesítési értékei a ∆ i háromszög oldalfelező pontjaiban. [1] Ha az egységnégyzet standard triangulációjakor a négyzet oldalait N egyenlő részre osztjuk, akkor a Vh tér leírásához n = (2 N + 1)
2
függvényt használunk. Ezek a kis háromszögek
csúcspontjaihoz és oldalfelező pontjaihoz rendelt függvények, minden vi függvény ezen pontok közül pontosan egyben vesz fel nullától különböző értéket: 1-et.
30
N= 4 esetén az ábrán látható pontokhoz tartozó vi függvényekre van szükség. Az
∫ f (x )v (x ) dx j
integrálok közelítésénél a kubatúra képlet alapján az
f ( x )v j ( x )
Ω
függvényeknek a háromszögek oldalfelező pontjaiban felvett helyettesítési értékeivel dolgozunk.
∫ f (x )v (x ) dx ≈ j
∆i
∆i 3
(f
)
v + f i2 vi2 + f i3 vi3 ,
i1 i1
ahol f i1vi1 , f i2 vi2 , f i3 vi3 a ∆ i háromszög oldalfelező pontjaiban az f ( x )v j ( x ) helyettesítési értéke. Ha a v j függvény valamelyik háromszög csúcspontjához tartozik, akkor értéke az összes oldalfelező pontban 0. Ha v j valamelyik háromszög egyik oldalfelező pontjához tartozik, akkor itt az értéke 1, a másik két oldalfelező pontban pedig 0. Így az integrálközelítő összeg értéke egy ilyen háromszög fölött
∆ 3
f j , ahol f j az f függvény azon oldalfelező pontbeli
értéke, melyben v j nem tűnik el., ∆ a kis háromszög területe. A négyzet belsejében elhelyezkedő oldalfelező pontok két háromszöghöz is tartoznak, így az Ω fölötti integrálközelítő összeg ilyen v j függvény esetén felezőpontok esetén
∆ 3
2⋅ ∆ 3
f j , míg a peremen elhelyezkedő
fj.
Az inhomogén peremfeltétellel adott feladatot a diszkrét megoldás során homogén peremfeltétellel adott feladatra vezetjük vissza. Ha a
v1 , v 2 , K , v n
( n = (2 N + 1) ) 2
függvényeket két csoportra osztjuk aszerint, hogy belső- ( vi , i ∈ I b ) vagy perempontokhoz
31
{
}
( vi , i ∈ I p ) tartoznak – ahol I b ∪ I p = 1, 2, K , (2 N + 1) , akkor a keresett 2
n
u h = ∑ y i vi i =1
függvény
u h = ∑ y i vi + ∑ y i v i i∈I b
i∈I p
alakba írható, ahol a perempontokban ismert az u h függvény értéke. Az
Ay = b
lineáris egyenletrendszer felírásánál figyelmen kívül hagyjuk a
peremfeltételt, majd az y vektor perempontbeli ( i ∈ I p ) koordinátái helyett a g függvény azokban felvett helyettesítési értékeit írjuk. Ekkor az
y ismeretlen vektor néhány
koordinátája már adott, az ezeknek megfelelő tagokat az egyenletrendszer jobboldalára rendezzük, majd az i ∈ I p indexű sorokat elhagyjuk a rendszerből. A kapott egyenletrendszer megoldása után a keresett y vektort kapjuk.
32
IRODALOMJEGYZÉK [1] Stoyan Gisbert-Takó Galina: Numerikus módszerek I., ELTE-TypoTEX, Budapest, 2002 [2] Stoyan Gisbert-Takó Galina: Numerikus módszerek II., ELTE-TypoTEX, Budapest, 1995 [3] Stoyan Gisbert-Takó Galina: Numerikus módszerek III., ELTE-TypoTEX, Budapest, 1997 [4] Stoyan Gisbert: MATLAB, Typotex, Budapest, 2008 [5] J. N Bronstein, K. A Szemangyejev, G. Musiol, H. Mühlig: Matematikai kézikönyv, TypoTEX, Budapest, 2000 [6] Braess, Dietrich: Finite elements, Theory, fast solvers, and applications in solid mechanics, Cambridge University Press, 2007 [7] S. Brenner, R. L. Scott: The Mathematical Theory of Finite Element Methods, 2nd edition, Springer, 2005 [8] Grundman, Maller: Invariant integration formulas for the n-simplex by combinatorical methods, SIAM Journal on Numerical Analysis 15, 282-290 [9] http://www.mathworks.com
33
