A Poisson-egyenletre alkalmazott multigrid mo dszer

A Poisson-egyenletre alkalmazott multigrid ´ dszer mo

Szakdolgozat ´Irta: Klimaj Bettina Matematika BSc Alkalmazott matematikus szakir´any

T´emavezet˝o: Dr. G´asp´ar Csaba egyetemi tan´ar Numerikus Anal´ızis Tansz´ek

E¨ otv¨ os Lor´and Tudom´anyegyetem Termeszettudom´anyi Kar 2012

K¨ osz¨ onetnyilv´ an´ıt´ as Szeretn´ek k¨osz¨onetet mondani t´emavezet˝omnek, G´asp´ar Csab´anak, ami´ert felkeltette ´erdekl˝od´esemet a t´ema ir´ant, hasznos tan´acsokkal, ´eszrev´etelekkel ´es seg´edanyagokkal l´atott el, ´es hogy k´erd´eseimmel b´armikor bizalommal fordulhattam hozz´a. K¨osz¨onettel tartozom a csal´adomnak, akikt˝ol rengeteg t´amogat´ast ´es biztat´ast kaptam. V´eg¨ ul szeretn´em megk¨osz¨onni Herczeg Bonif´acnak, hogy matematikai hozz´a´ert´es´evel seg´ıtette a munk´amat, a rengeteg t¨ urelmet, ´es hogy mindig mindenben mellettem ´allt.

2

Tartalomjegyz´ ek Bevezet´ es

5

1. Parci´ alis differenci´ alegyenletek ´ es megold´ asi m´ odszereik

7

1.1. A parci´alis differenci´alegyenletek fogalma, oszt´alyoz´asa. . . . . . . .

7

1.2. Fizikai p´eld´ak elliptikus egyenletekre ´es peremfelt´etelekre (3D-ben) .

9

1.3. Numerikus megold´asi m´odszerek . . . . . . . . . . . . . . . . . . . . . 13 1.3.1. V´eges differencia m´odszer . . . . . . . . . . . . . . . . . . . . 13 1.3.2. V´egeselem-m´odszer . . . . . . . . . . . . . . . . . . . . . . . . 16 2. A multigrid m´ odszer

18

2.1. Diszkretiz´aci´o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2. A m´odszer alap¨otlete . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3. Jav´ıt´as durv´abb r´acson . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4. A sim´ıt´o iter´aci´ok ´es sim´ıt´o tulajdons´agaik . . . . . . . . . . . . . . . 24 2.5. K´etr´acsos m´odszer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.6. A multigrid ciklus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.7. A teljes multigrid m´odszer . . . . . . . . . . . . . . . . . . . . . . . . 32 3. A multigrid m´ odszer implement´ al´ asa MATLAB-ban

36

3.1. Multigrid program . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2. A teljes multigrid m´odszer hat´ekonys´ag´anak vizsg´alata . . . . . . . . 38 3.3. A kaszk´ad m´odszer megval´os´ıt´asa . . . . . . . . . . . . . . . . . . . . 42 3.4. Az FMG ´es a Seidel-iter´aci´o konvergenciasebess´eg´enek o¨sszehasonl´ıt´asa 43 ¨ Osszefoglal´ as

48

F¨ uggel´ ek

49

Irodalomjegyz´ ek

61

3

´ ak jegyz´ Abr´ eke 2.1. Jav´ıt´as durva r´acson . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2. Ml iter´aci´os m´atrix saj´at´ert´ekei µhl f¨ uggv´eny´eben . . . . . . . . . . . 27 2.3. K´etr´acsos m´odszer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.4. Egy multigrid iter´aci´os l´ep´es az Ll ul = fl (l ≥ 1) megold´as´ara

. . . . 33

3.1. A k¨ozel´ıt˝o megold´as hib´aja ´es a kisz´am´ıt´as´ahoz sz¨ uks´eges id˝o a szintek f¨ uggv´eny´eben . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2. A k¨ozel´ıt´esek hib´aj´anak cs¨okken´ese - 2. tesztf¨ uggv´eny (l=3, 4) . . . . 40 3.3. A k¨ozel´ıt´esek hib´aj´anak cs¨okken´ese - 2. tesztf¨ uggv´eny (l=5, 6) . . . . 40 3.4. Els˝o tesztfeladat - sim´ıt´o iter´aci´o n´elk¨ ul: mgteszt(6,0,0,1,1) . . . . . . 41 3.5. M´asodik tesztfeladat - simit´o iter´aci´o n´elk¨ ul: mgteszt(6,0,0,1,2) . . . 41 3.6. A multigrid ´es a kaszk´ad m´odszer a´ltal nyert k¨ozel´ıt˝o megold´asok hib´aja ´es a kisz´am´ıt´asukhoz sz¨ uks´eges id˝o a szintek f¨ uggv´eny´eben . . 43 3.7. A Seidel- ´es a multigrid m´odszer a´ltal nyert k¨ozel´ıt˝o megold´asok l´enyeg´eben azonos hib´aja ´es a kisz´am´ıt´asukhoz sz¨ uks´eges id˝o a szintek f¨ uggv´eny´eben . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.8. A Seidel- ´es a multigrid m´odszer a´ltal nyert k¨ozel´ıt˝o megold´asok hib´aja nagyj´ab´ol azonos m˝ uvelet eset´en, a szintek f¨ uggv´eny´eben . . . 44 3.9. Az FMG ´es a Seidel-iter´aci´o o¨sszehasonl´ıt´asa - 1. tesztf¨ uggv´eny (l=2) 45 3.10. Az FMG ´es a Seidel-iter´aci´o o¨sszehasonl´ıt´asa - 1. tesztf¨ uggv´eny (l=3) 46 3.11. Az FMG ´es a Seidel-iter´aci´o o¨sszehasonl´ıt´asa - 1. tesztf¨ uggv´eny (l=4) 46 3.12. Az FMG ´es a Seidel-iter´aci´o o¨sszehasonl´ıt´asa - 1. tesztf¨ uggv´eny (l=5) 47

4

Bevezet´ es A diszkr´et elliptikus perem´ert´ek-probl´em´ak megold´as´ara szolg´al´o iterat´ıv multigrid algoritmus a k¨ovetkez˝ok¨on alapszik: A klasszikus iter´aci´ok seg´ıts´eg´evel el´erhet˝o, hogy egy tetsz˝oleges k¨ozel´ıt´es hib´aj´anak magas frekvenci´as komponensei kell˝oen kicsik legyenek, ebb˝ol ad´od´oan lehet˝os´eg ny´ılik az adott k¨ozel´ıt´es jav´ıt´as´ara a hiba egy durv´abb r´acson sz´am´ıtott k¨ozel´ıt´es´evel. Ezen ¨otlet rekurz´ıv m´odon t¨ort´en˝o alkalmaz´asa egyre durv´abb ´es durv´abb r´acsokon egy aszimptotikusan optim´alis, azaz a diszkr´et ismeretlensz´ammal ar´anyos m˝ uveletig´eny˝ u iter´aci´os elj´ar´ast eredm´enyez. A multigrid m´odszer vizsg´alat´aban el´ert els˝o eredm´enyek Fedorenko nev´ehez f˝ uz˝odnek, aki az egys´egn´egyzeten ´ertelmezett Poisson-egyenleten tanulm´anyozta a m´odszer konvergenci´aj´at. Ut´ana Bakhvalov vizsg´alta a m´odszer viselked´es´et m´asodrend˝ u, nem a´lland´o egy¨ utthat´os elliptikus perem´ert´ek-probl´em´ak eset´eben. Az algoritmus gyakorlatban mutatott hat´ekonys´ag´anak felismer´ese Achi Brandt [1] ´es a t˝ole f¨ uggetlen¨ ul alkot´o Hackbush [5, 6] ´erdeme. Megeml´ıtend˝o, hogy a multigrid elj´ar´as sz´amos egy´eb ter¨ uleten haszn´alatos, p´eld´aul nem-line´aris perem´ert´ek-probl´em´ak, saj´at´ert´ek-probl´em´ak, bifurk´aci´os probl´em´ak, parabolikus ´es hiperbolikus probl´em´ak, s˝ot integr´alegyenletek megold´as´ara is alkalmazz´ak. A dolgozat els˝o fejezet´eben a´ttekintj¨ uk, hogy mik is azok a parci´alis differenci´alegyenletek, ´es hogy milyen gyakorlati alkalmaz´asok ¨oszt¨on¨ozt´ek a bevezet´es¨ uket. Konkr´et fizikai p´eld´akat l´athatunk, melyeket a Poisson-egyenlet (egy a´ltal´anosabb v´altozata) ´ır le, mely egy speci´alis elliptikus t´ıpus´ u parci´alis differenci´alegyenlet. A szakdolgozat f˝o c´elja egy ezen egyenlet k¨ozel´ıt˝o megold´as´anak megkeres´es´ere szolg´al´o algoritmus bemutat´asa. Ennek a numerikus probl´em´anak a megold´as´ahoz a´ttekintj¨ uk a v´eges differenci´ak m´odszer´et, melynek seg´ıts´eg´evel diszkretiz´alhatjuk a Poisson-egyenletet, azaz helyette egy v´eges sok ismeretlent tartalmaz´o k¨ozel´ıt˝oegyenlet megold´as´ara szor´ıtkozunk. A m´asodik fejezetben betekint´est nyer¨ unk a multigrid m´odszer rejtelmeibe. R´eszletesen ´attekintj¨ uk az algoritmus ´ep´ıt˝ok¨oveit, alapegys´egeit, melyek teh´at a hiba-

5

sim´ıt´o iter´aci´ok, a k¨ozel´ıt´es korrekci´oj´anak kisz´am´ıt´asa egy durv´abb h´al´on, valamint a m´odszer ¨onmegh´ıv´o mivolta. Ennek egy m´eg hat´ekonyabb tov´abbfejleszt´ese a m´ar nem iterat´ıv teljes multigrid m´odszer, melynek alap¨otlete az, hogy az egyre jobb k¨ozel´ıt´esek el´er´es´ehez multigrid algoritmust alkalmazunk az el˝oz˝o k¨ozel´ıt´est haszn´alva az iter´aci´o kezdeti k¨ozel´ıt´es´enek. A multigrid m´odszer elm´elet´enek megismer´ese ut´an a harmadik fejezetben a gyakorlatban val´o alkalmaz´as´ara l´athatunk p´eld´at. H´arom p´eldaprogram hivatott szeml´eltetni az algoritmus rendk´ıv¨ ul gyors konvergenci´aj´at a Poisson-egyenlet megold´asa sor´an. Az els˝on´el megfigyelhetj¨ uk, hogy egyre finomabb h´al´ok alkalmaz´asa eset´en hogy cs¨okken az egyre javul´o k¨ozel´ıt´es hib´aja amellett, hogy m´ar a m´asodperc t¨ored´eke alatt is szert tehet¨ unk egy el´eg j´o k¨ozel´ıt˝o megold´asra. A klasszikus Seideliter´aci´oval szembe´all´ıtva a m´odszer hat´ekonys´aga m´eg szeml´eletesebb, melyet a m´asodik ´es a harmadik program seg´ıts´eg´evel k¨ovethet¨ unk nyomon. El˝obbiben azt m´erj¨ uk, mennyi id˝ore van sz¨ uks´ege a Seidel-iter´aci´onak olyan pontoss´ag´ u k¨ozel´ıt´es el˝o´all´ıt´as´ahoz, mint amilyet a multigrid m´odszer szolg´altat, ut´obbiban pedig megfigyelj¨ uk, hogy nagyj´ab´ol azonos m˝ uvelet v´egrehajt´as´aval melyik m´odszer mennyire pontos k¨ozel´ıt´est eredm´enyez. V´eg¨ ul a negyedik program az egyszer˝ u kaszk´ad algoritmus ´es a j´oval bonyolultabb teljes multigrid m´odszer gyorsas´ag´at ´es pontoss´ag´at hasonl´ıtja o¨ssze. A megfelel˝o programk´odok a f¨ uggel´ekben olvashat´ok.

6

1. fejezet Parci´ alis differenci´ alegyenletek ´ es megold´ asi m´ odszereik 1.1.

A parci´ alis differenci´ alegyenletek fogalma, oszt´ alyoz´ asa.

Vegy¨ unk egy Ω ⊂ Rn tartom´anyt (azaz ¨osszef¨ ugg˝o, ny´ılt halmazt), melyre n ≥ 2. Az α = (α1 , α2 , . . . , αn ) u ´gynevezett multiindexekre (αj ≥ 0 eg´esz) vezess¨ uk be a k¨ovetkez˝o jel¨ol´eseket: |α| :=

Pn

j=1

αj

´es

∂ α u := ∂1α1 ∂2α2 . . . ∂nαn u, tetsz˝oleges u

n-v´altoz´os f¨ uggv´enyre. Legyen N azon multiindexek sz´ama, melyek kiel´eg´ıtik az |α| ≤ m egyenl˝otlens´eget, ahol m ∈ N adott sz´am. Tekints¨ unk egy G ⊂ RN tartom´anyt ´es egy F : Ω × G → R adott f¨ uggv´enyt. Keress¨ unk olyan u : Ω → R m-szer folytonosan differenci´alhat´o f¨ uggv´enyt, melyre F (x, u(x), ∂1 u(x), . . . , ∂ α u(x), . . .) = 0 teljes¨ ul minden x ∈ Ω ´es |α| ≤ m eset´en, azaz fenn´all, hogy F ◦ (id, u, ∂1 u, . . . , ∂ α u, . . .) = 0.

(1.1)

Az (1.1)-et parci´alis differenci´alegyenletnek, az m sz´amot pedig a rendj´enek nevezz¨ uk. A egyik legfontosabb speci´alis t´ıpus´ u egyenlet a f˝or´esz´eben line´aris parci´alis differenci´alegyenlet, mely a´ltal´anos alakja a k¨ovetkez˝o: X

aα ∂ α u = f ◦ (id, u, ∂1 u, . . . , ∂ β u, . . .),

|α|=m

ahol aα , f : Ω → R adott f¨ uggv´enyek ´es |β| ≤ m − 1. f = 0 eset´en homog´en, egy´ebk´ent inhomog´en egyenletr˝ol besz´el¨ unk. Amennyiben aα a´lland´o, a differenci´al7

egyenletet ´alland´o egy¨ utthat´osnak nevezz¨ uk. (A parci´alis differenci´alegyenletekr˝ol, illetve ezek tov´abbi t´ıpusair´ol r´eszletesebb anyag [2]-ben tal´alhat´o.) A m´asodrend˝ u f˝or´esz¨ ukben line´aris parci´alis differenci´alegyenletek a k¨ovetkez˝o alakban ´ırhat´ok:

n X

ajk ∂j ∂k u = g ◦ (id, u, ∂1 u, . . . , ∂n u),

j,k=1

ahol ajk : Ω → R, g : Ω × G → R adott f¨ uggv´enyek (itt G ⊂ Rn+1 tartom´any). ajk egy¨ utthat´okr´ol feltehet˝o, hogy val´osak ´es ajk = akj . Egy x0 ∈ Ω r¨ogz´ıtett pontban tekints¨ uk az A val´os szimmetrikus egy¨ utthat´om´atrixot: A := A(x0 ) = [ajk (x0 )]nj,k=1 . Ismeretes, hogy ezen m´atrix minden saj´at´ert´eke val´os. Jel¨olje n+ , n− , illetve n0 az A m´atrix pozit´ıv, negat´ıv, illetve null´aval egyenl˝o saj´at´ert´ekeinek (multiplicit´assal vett) sz´am´at. Ekkor n+ + n− + n0 = n. 1.1.1. Defin´ıci´ o. Azt mondjuk, hogy a f˝or´esz´eben line´aris m´asodrend˝ u parci´alis differenci´alegyenlet az Ω tartom´anyon 1. elliptikus, ha n+ = n vagy n− = n, 2. hiperbolikus, ha n+ = 1 ´es n− = n − 1 vagy n+ = n − 1 ´es n− = 1, 3. parabolikus, ha n0 = 1 ´es n+ = n − 1 vagy n− = n − 1, ´es a megfelel˝o felt´etel a tartom´any minden pontj´aban teljes¨ ul. P´ eld´ ak: 1.

a ∆u :=

n X

∂j2 u = 0 Laplace-egyenlet elliptikus

j=1

2.

a ∂02 u −

n X

∂j2 u = 0 egyenlet hiperbolikus

j=1

3.

a ∂0 u −

n X

∂j2 u = 0 egyenlet parabolikus

j=1

Sz´amunkra a k´es˝obbiek sor´an a Laplace-egyenlet inhomog´en a´ltal´anos´ıt´asa, azaz a ∆u =

n X

∂j2 u = f,

j=1

u ´gynevezett Poisson-egyenlet lesz fontos – annak is az egy- illetve k´etdimenzi´os v´altozata –, mely teh´at egy elliptikus t´ıpus´ u m´asodrend˝ u line´aris parci´alis differenci´alegyenlet.

8

1.2.

Fizikai p´ eld´ ak elliptikus egyenletekre ´ es peremfelt´ etelekre (3D-ben)

N´ezz¨ unk meg most n´eh´any a val´os ´eletb˝ol vett fizikai folyamatot, melyek le´ır´asa az elliptikus parci´alis differenci´alegyenletek, azon bel¨ ul is a Poisson-egyenlet egy a´ltal´anosabb v´altozat´anak seg´ıts´eg´evel t¨ort´enik. Mivel egy parci´alis differenci´alegyenletnek a´ltal´aban v´egtelen sok megold´asa van, a sz´obanforg´o jelens´egek egy´ertelm˝ u le´ır´as´ahoz sz¨ uks´eges megadnunk bizonyos mell´ekfelt´eteleket is. A tov´abbiakban n´eh´any p´eld´an kereszt¨ ul a´ttekintj¨ uk, hogyan vezethet valamely fizikai folyamat parcia´lis differenci´alegyenletre, ´es hogy milyen mell´ekfelt´eteleket sz¨ uks´eges el˝o´ırnunk a megoldhat´os´ag ´erdek´eben.

Stacion´ arius h˝ ovezet´ es Legyen Ω ⊂ R3 adott tartom´any Γ hat´arral, melyben stacion´arius h˝ovezet´es j´atsz´odik le, azaz a h˝om´ers´ekleteloszl´as az id˝ot˝ol f¨ uggetlen. Jel¨olje u(x) az id˝oben a´lland´osult h˝om´ers´ekletet (m´ert´ekegys´ege: K (kelvin)) az x ∈ Ω pontban. Ekkor az u f¨ uggv´eny az Ω tartom´anyban kiel´eg´ıti a 3 X

− α · ∆u = −α ·

∂j2 u = f

(1.2)

j=1

egyenletet, ahol f : Ω → R f¨ uggv´eny az Ω-ban l´ev˝o h˝oforr´asok intenzit´as´at mutatja (m´ert´ekegys´ege: W/m3 (watt/m3 )), α : Ω → R pedig az Ω egyes pontjaihoz tartoz´o h˝odiff´ uzi´os t´enyez˝oket adja meg (m´ert´ekegys´ege: W/mK (watt/m·kelvin)). Integr´alva az egyenlet¨ unket az Ω tartom´anyon, alkalmazva az Z

div~v dx =

Z

~v · ~n ds

Γ

Ω

Gauss-f´ele divergencia-t´etelt (ahol ~n a Γ perem Ω-b´ol kifel´e mutat´o norm´alisa), valamint felhaszn´alva a ∆u = div(∇u) ´es ∇u · ~n = o¨sszef¨ ugg´eseket az Z ∂u ds = f dx n Γ ∂~ Ω egyenletet nyerj¨ uk, melyet ´atrendezve

−

Z Γ

Z

α

α

Z ∂u ds + f dx = 0 ∂~n Ω

9

∂u ∂~n

ad´odik. Ennek fizikai jelent´ese az, hogy az Ω tartom´anyban a h˝oforr´asok keltette energia ´es a peremen (id˝oegys´eg alatt) be´araml´o energia ¨osszege megegyezik a peremen kereszt¨ ul ki´araml´o h˝oenergi´aval (energiamegmarad´as). Speci´alisan, ha a tartom´anyban nincsenek h˝oforr´asok (f ≡ 0), akkor azt kapjuk, hogy a peremen kereszt¨ ul id˝oegys´eg alatt a tartom´anyba be´araml´o illetve az onnan ki´araml´o energia megegyezik. Ahhoz, hogy az Ω tartom´any h˝om´ers´ekleteloszl´as´at egy´ertelm˝ uen meg tudjuk hat´arozni, ismern¨ unk kell Ω hat´ar´an az al´abbi peremfelt´etelek egyik´et: 1. peremfelt´etel: u|Γ = ϕ 2. peremfelt´etel: α

∂u = ϕ (~ n a Γ perem Ω-b´ol kifel´e mutat´o norm´alisa) ∂~n Γ



3. peremfelt´etel: β α

 ∂u + δu|Γ = ϕ (β, δ adott f¨ uggv´enyek Γ peremen) ∂~n Γ

´Igy kapjuk az a´ltal´anosabb Poisson-egyenletre vonatkoz´o Dirichlet-feladatot ((1.2)+ +1. peremfelt´etel), Neumann-feladatot ((1.2)+2. peremfelt´etel), valamint a harmadik t´ıpus´ u perem´ert´ek-feladatot ((1.2)+3. peremfelt´etel). Ezen peremfelt´etelek fizikai jelent´ese: 1. Dirichlet-peremfelt´etel: a tartom´any perem´en a h˝om´ers´eklet el˝o´ırt 2. Neumann-peremfelt´etel: a tartom´any perem´en az α

∂u mennyis´eg, vagyis az u ´gynevezett norm´alis ir´a∂~n

ny´ u h˝oa´ram-s˝ ur˝ us´eg (m´ert´ekegys´ege: W/m2 ) el˝o´ırt. (Speci´alisan, ha ez z´erus, az azt jelenti, hogy Ω h˝oszigetelt.) 3. Harmadik t´ıpus´ u, vagy Robin-peremfelt´etel: a tartom´any perem´en a h˝om´ers´eklet ´es a h˝oa´ram-s˝ ur˝ us´eg egy line´aris kombin´aci´oja el˝o´ırt. (K¨ozel´ıt˝oleg ilyen peremfelt´etel ´erv´enyes, ha a testet egy rossz h˝ovezet´es˝ u, de a test m´ereteihez k´epest v´ekony anyagba burkoljuk, ´es a h˝om´ers´eklet ennek a k¨ uls˝o fel¨ ulet´en lesz adott.)

Stacion´ arius elektromos ´ aram kiterjedt vezet˝ oben Legyen Ω ⊂ R3 adott tartom´any Γ hat´arral, melyben stacion´arius (id˝ot˝ol f¨ uggetlen) a´ramterjed´es folyik. Jel¨olje u(x) az (id˝oben a´lland´osult) elektromos potenci´alt (m´ert´ekegys´ege: V (volt)) az x ∈ Ω pontban. Ekkor az u f¨ uggv´eny az Ω tartom´anyban kiel´eg´ıti a −α · ∆u = −α ·

3 X j=1

10

∂j2 u = f

egyenletet, ahol f : Ω → R f¨ uggv´eny az Ω-ban l´ev˝o ´aramforr´asok keltette a´rams˝ ur˝ us´eget mutatja (m´ert´ekegys´ege: A/m3 (amper/m3 )), α : Ω → R pedig Ω pontjainak fajlagos vezet˝ok´epess´eg´et adja meg (m´ert´ekegys´ege: 1/mΩ ( 1/m·ohm )). Vegy¨ uk a m´ ultkori meggondol´asnak megfelel˝oen nyert: Z ∂u ds + f dx = 0 α n Ω Γ ∂~

Z

egyenletet. Ennek az a fizikai jelent´ese, hogy az Ω tartom´anyban az a´ramforr´asok keltette ´aram valamint a peremen bel´ep˝o a´ram o¨sszege megegyezik a peremen kil´ep˝o a´rammal, vagyis az id˝oegys´eg alatt a´t´araml´o t¨olt´esmennyis´eggel (t¨olt´esmegmarad´as). Speci´alisan, ha a tartom´anyban nincsenek a´ramforr´asok (f ≡ 0), akkor a peremen a tartom´anyba bel´ep˝o illetve az onnan kil´ep˝o a´ramok megegyeznek. A perem´ert´ek-feladatok megad´asa ugyan´ ugy t¨ort´enik, mint az el˝oz˝o p´elda eset´eben. Ahhoz, hogy az Ω tartom´any potenci´aleloszl´as´at egy´ertelm˝ uen meg tudjuk hat´arozni, ismern¨ unk kell Ω hat´ar´an valamelyik peremfelt´etelt. N´ezz¨ uk meg enn´el a p´eld´an´al milyen fizikai jelent´esek b´ ujnak meg a peremfelt´etelek m¨og¨ott: 1. Dirichlet-peremfelt´etel: a tartom´any perem´en az elektromos potenci´al el˝o´ırt 2. Neumann-peremfelt´etel: a tartom´any perem´en az α

∂u mennyis´eg, vagyis a norm´alis ir´any´ u elektromos ∂~n

a´rams˝ ur˝ us´eg (m´ert´ekegys´ege: A/m2 ) el˝o´ırt. (Speci´alisan, ha ez z´erus, az azt jelenti, hogy Ω elektromosan szigetelt.) 3. Harmadik t´ıpus´ u, vagy Robin-peremfelt´etel: a tartom´any perem´en az elektromos potenci´al ´es az a´rams˝ ur˝ us´eg egy line´aris kombin´aci´oja el˝o´ırt. (K¨ozel´ıt˝oleg ilyen peremfelt´etel ´erv´enyes, ha a testet egy elektromosan rossz vezet˝ok´epess´eg˝ u, de a test m´ereteihez k´epest v´ekony anyagba burkoljuk, ´es az elektromos potenci´al ennek a k¨ uls˝o fel¨ ulet´en lesz adott.)

Stacion´ arius sziv´ arg´ as por´ ozus k¨ ozegen kereszt¨ ul Legyen Ω ⊂ R3 adott tartom´any (por´ozus k¨ozeg) Γ hat´arral, melyben stacion´arius (id˝ot˝ol f¨ uggetlen) sziv´arg´as t¨ort´enik.

Jel¨olje u(x) a folyad´ek (id˝oben

a´lland´osult) sziv´arg´asi potenci´alj´at (m´ert´ekegys´ege: m) az x ∈ Ω pontban. Ezen u f¨ uggv´enyre ´erv´enyes az u=h+

11

p γ

o¨sszef¨ ugg´es, ahol h egy r¨ogz´ıtett alapszintt˝ol m´ert geodetikus magass´ag, p jel¨oli a nyom´ast, γ pedig a fajs´ ulyt. Ekkor az u f¨ uggv´eny az Ω tartom´anyban kiel´eg´ıti a −α · ∆u = −α ·

3 X

∂j2 u = f

j=1

egyenletet, ahol f : Ω → R f¨ uggv´eny az Ω-ban l´ev˝o bels˝o v´ızforr´asok intenzit´as´at mutatja (azaz, hogy a v´ızforr´as egys´egnyi t´erfogatban h´any m3 /sec o¨sszhozam´ u, m´ert´ekegys´ege enn´elfogva: (m3 /sec)/m3 = 1/sec), α : Ω → R pedig az Ω egyes pontjaihoz tartoz´o sziv´arg´asi t´enyez˝oket adja meg (m´ert´ekegys´ege: m/sec), tov´abb´a α · ∇u a sziv´arg´asi sebess´eget jelenti. Az Z ∂u ds + f dx = 0 n Γ ∂~ Ω egyenlet jelen esetben is egy fontos fizikai jelens´eget ´ır le. Jelent´ese az, hogy az Ω Z

α

tartom´anyban a bels˝o v´ızforr´asok ¨osszhozama ´es a peremen bel´ep˝o v´ızhozam ¨osszege megegyezik a peremen kil´ep˝o v´ızhozammal, vagyis az id˝oegys´eg alatt a´t´araml´o v´ızt´erfogat-mennyis´eggel (t¨omegmegmarad´as). Speci´alisan, ha a tartom´anyban nincsenek v´ızforr´asok (f ≡ 0), akkor a peremen a tartom´anyba bel´ep˝o illetve az onnan kil´ep˝o v´ızhozamok megegyeznek. N´ezz¨ uk meg ebben az esetben milyen fizikai tartalommal b´ırnak a peremfelt´etelek: 1. Dirichlet-peremfelt´etel: a tartom´any perem´en a sziv´arg´asi potenci´al el˝o´ırt. (Ilyen peremfelt´etel ´erv´enyes, ha a perem szabad v´ızt´errel (foly´o, t´o) ´erintkezik, ekkor a sziv´arg´asi potenci´al a szabad v´ızt´er felsz´ın´enek geodetikus magass´ag´aval egyezik meg.) 2. Neumann-peremfelt´etel: ∂u mennyis´eg, vagyis a norm´alis ir´any´ u sziv´arg´asi ∂~n sebess´egkomponens (m´ert´ekegys´ege: m/sec) el˝o´ırt. (Speci´alisan, ha ez z´erus, a tartom´any perem´en az α

az azt jelenti, hogy Ω sziv´arg´asi tartom´anyt v´ızz´ar´o anyag hat´arolja.) 3. Harmadik t´ıpus´ u, vagy Robin-peremfelt´etel: a tartom´any perem´en a sziv´arg´asi potenci´al ´es az norm´alis ir´any´ u sebess´egkomponens egy line´aris kombin´aci´oja el˝o´ırt. (K¨ozel´ıt˝oleg ilyen peremfelt´etel ´erv´enyes, ha a sziv´arg´asi tartom´anyt egy majdnem v´ızz´ar´o, de a tartom´any m´ereteihez k´epest v´ekony r´eteg bor´ıtja, ´es a sziv´arg´asi potenci´al ennek a k¨ uls˝o fel¨ ulet´en lesz adott. Ez a helyzet p´eld´aul foly´ok vagy tavak fenek´en, ha az id˝ok folyam´an le¨ ulepedett finomszemcs´es anyag r´eszben elt¨om´ıti a talaj p´orusait, melyeken kereszt¨ ul a sziv´arg´as v´egbemegy. Ez a jelens´eg a partk¨ozeli kutak m˝ uk¨od´es´ere kedvez˝otlen hat´assal van, cs¨okkenti a kutak v´ızhozam´at.) 12

Tov´ abbi ´ altal´ anos megjegyz´ esek: El˝ofordulhat – s˝ot, ´altal´aban ez a jellemz˝o – hogy Ω hat´ar´an egyszerre t¨obb, k¨ ul¨onb¨oz˝o t´ıpus´ u peremfelt´etel is adott, azaz Γ diszjunkt m´odon felbonthat´o v´eges sok peremszakaszra, melyeken m´as-m´as t´ıpus´ u peremfelt´etelt (de csak egyf´el´et) ´ırunk el˝o. Ilyenkor kevert peremfelt´etelr˝ol besz´el¨ unk. Amennyiben a vizsg´alt fizikai folyamatban az anyagra jellemz˝o α t´enyez˝o nem konstans, hanem a helyt˝ol f¨ ugg, a folyamatot az −α · ∆u = f egyenlet helyett a bonyolultabb −div(α · ∇u) = f egyenlet ´ırja le. Ha α konstans, akkor term´eszetesen beolvaszthat´o a jobb oldali f¨ uggv´enybe, ´ıgy a formailag egyszer˝ ubb ∆u = f Poisson-egyenletet nyerj¨ uk. Hangs´ ulyozzuk azonban, hogy ekkor a megold´asnak, a jobb oldalnak valamint a peremfelt´eteleknek m´ar nem a fentebb v´azolt term´eszetes” ” fizikai jelent´ese van.

1.3.

Numerikus megold´ asi m´ odszerek

Mint l´attuk sz´amos a gyakorlatban felmer¨ ul˝o probl´ema megold´asa sor´an parci´alis differenci´alegyenletekbe u ¨tk¨oz¨ unk. Ezen egyenletek numerikus megold´asa a diszkretiz´aci´o seg´ıts´eg´evel t¨ort´enik, amikoris a parci´alis differenci´alegyenletet egy m´asik, v´eges sok ismeretlent tartalmaz´o egyenlettel k¨ozel´ıtj¨ uk. Ilyenkor az ´ertelmez´esi tartom´anyt egy v´eges sok pontb´ol ´all´o h´al´oval helyettes´ıtj¨ uk. A k´et legismertebb diszkretiz´aci´os technika a v´eges differencia m´odszer ´es a v´egeselem-m´odszer. A k¨ovetkez˝o r´eszben ezen m´odszerek r¨ovid ismertet´es´ere ker¨ ul sor (a r´eszletek´ert lsd. [7]).

1.3.1.

V´ eges differencia m´ odszer

A v´eges differencia m´odszer a parci´alis deriv´altak lok´alis approxim´aci´oj´an alapul. A k¨ozel´ıt´eseket a parci´alis deriv´altak Taylor-sorba fejtett alakj´ab´ol nyerj¨ uk. Egy adott u f¨ uggv´eny x pontbeli els˝orend˝ u deriv´altj´anak legegyszer˝ ubb k¨ozel´ıt´ese a k¨ovetkez˝o:

∂u u(x + h) − u(x) (x) ≈ . ∂x h 13

Ha u n´egyszer folytonosan differenci´alhat´o x egy k¨ornyezet´eben, akkor a Taylorformula felhaszn´al´as´aval ∂u h2 ∂ 2 u h3 ∂ 3 u h4 ∂ 4 u u(x + h) = u(x) + h · (x) + · · · (x) + (x) + (ξ1 ), ∂x 2 ∂x2 6 ∂x3 24 ∂x4

(1.3)

illetve u(x − h) = u(x) − h ·

h2 ∂ 2 u h3 ∂ 3 u h4 ∂ 4 u ∂u (x) + · 2 (x) − · 3 (x) + · (ξ2 ) ∂x 2 ∂x 6 ∂x 24 ∂x4

(1.4)

ad´odik, ahol ξ1 ∈ (x, x + h) ´es ξ2 ∈ (x − h, x). (1.3) ´atrendez´es´eb˝ol ∂u u(x + h) − u(x) h ∂ 2 u (x) = − · 2 (x) + O(h2 ) ∂x h 2 ∂x ad´odik, melyb˝ol az els˝orend˝ u deriv´alt el˝obb l´atott k¨ozel´ıt´ese ered. (1.3)-b´ol ´es (1.4)b˝ol – alkalmazva a negyedrend˝ u parci´alis deriv´altakra a k¨oz´ep´et´ekt´etelt – a k¨ovetkez˝o formul´at nyerj¨ uk: u(x + h) − 2u(x) + u(x − h) h2 ∂ 4 u ∂ 2u (x) = − · (ξ), ∂x2 h2 12 ∂x4 ahol ξ1 ≤ ξ ≤ ξ2 . Differencias´ ema a Laplace-oper´ ator k¨ ozel´ıt´ es´ ere Vegy¨ unk u k´etv´altoz´os f¨ uggv´eny ´ertelmez´esi tartom´any´aban egy t´eglalap alak´ u r´acsot, melyben h1 az x ir´any´ u, h2 az y ir´any´ u r´acsm´eret. Ha a m´asodrend˝ u deriv´alt el˝obb ismertetett k¨ozel´ıt´es´et a Laplace oper´atorban szerepl˝o ∂x2 u ´es ∂y2 u kifejez´esekre alkalmazzuk, a k¨ovetkez˝ot kapjuk eredm´eny¨ ul: ∆u(x, y) ≈

u(x+h1 , y)−2u(x, y)+u(x−h1 , y) u(x, y+h2 )−2u(x, y)+u(x, y−h2 ) + , h21 h22

Abban a speci´alis esetben, amikor h := h1 = h2 , az el˝oz˝o k¨ozel´ıt´es a al´abbi alakra egyszer˝ us¨odik: ∆u(x, y) ≈

1 [u(x + h, y) + u(x − h, y) + u(x, y + h) + u(x, y − h) − 4u(x, y)], h2

melyet a Laplace oper´ator k¨oz¨ons´eges ¨otpontos approxim´aci´oj´anak nevez¨ unk. Ezen k¨ozel´ıt´es egyszer˝ us´ıtett jel¨ol´es´ere a k¨ovetkez˝o ”stencil” haszn´alatos:     

−1 −1

4 −1

14

 

 −1  

V´ eges differencia m´ odszer egydimenzi´ os probl´ em´ ara N´ezz¨ uk a k¨ovetkez˝o egydimenzi´os perem´ert´ek-probl´em´at: −u00 (x) = f (x) x ∈ (0, 1), u(x) = 0 x ∈ {0, 1}. A [0, 1] intervallumot u ´gy diszkretiz´aljuk, hogy vessz¨ uk azon n + 2 pontj´at, melyekre xi = i · h (i = 0, . . . , n + 1) ´es h = 1/(n + 1). Alkalmazzuk most a v´eges differenci´ak m´odszer´et! Eszerint az u(xi ) pontok ui k¨ozel´ıt´eseit (i = 1, ..., n) a −ui−1 + 2ui − ui+1 = h2 fi

(i = 1, . . . , n)

rendszer megold´as´aval nyerj¨ uk (ahol fi = f (xi ) ´es u0 = u(x0 ), un+1 = u(xn+1 )). Mivel u0 = un+1 = 0, az eredeti probl´ema diszkretiz´as´aval ´es a peremfelt´etelek elimin´al´as´aval nyert line´aris rendszer a k¨ovetkez˝o: Lh uh = fh , ahol 



         

       ,      

2 −1    −1 2 −1  1 Lh = 2 h

−1

2 ..

−1 . . . .. .. −1

2

−1

−1

2

uh = (u1 , . . . , un ),

fh = (f1 , . . . , fn ).

V´ eges differencia m´ odszer k´ etdimenzi´ os probl´ em´ ara Vegy¨ uk az al´abbi k´etdimenzi´os perem´ert´ek-probl´em´at: ∂ 2u ∂ 2u − + = f (Ω-ban), ∂x2 ∂y 2 u = 0 (Γ ment´en), !

ahol Ω egy (0, A1 ) × (0, A2 )-es t´eglalap ´es Γ a hat´ara. Ω diszkretiz´al´as´aval egy Ωh r´acsot kapunk, melyre az x ir´any´ u r´acsm´eret h1 , az y ir´any´ u h2 (h := (h1 , h2 )) ´es a r´acsot (n1 + 2) · (n2 + 2) pont alkotja, azaz Ωh = {(x, y) = (κ1 h1 , κ2 h2 ) : κj = 0, . . . , nj + 1, hj =

Aj (j = 1, 2)}, nj + 1

N´ezz¨ uk most a h := h1 = h2 esetet, valamint legyen uh az ismeretlen vektor, melyet a bels˝o r´acspontokban felvett u(x, y) ´ert´ekek k¨ozel´ıt´esei alkotnak (a r´acspontokat 15

balr´ol jobbra, lentr˝ol felfele, sorfolytonosan szok´as sorba rendezni), fh pedig az f f¨ uggv´eny bels˝o r´acspontokban felvett ´ert´ekeib˝ol a´ll´o oszlopvektor. Alkalmazzuk most a v´eges differencia m´odszert! Eszerint az u(x+h, y), u(x−h, y), u(x, y+h), u(x, y−h) t´ıpus´ u pontn´egyesek uE , uW , uN , uS k¨ozel´ıt´eseit a uE + uW + uN + uS − 4uC = h2 fC t´ıpus´ u egyenletekb˝ol a´ll´o rendszer megold´as´aval nyerj¨ uk. Mivel u a hat´aron 0, az eredeti probl´ema diszkretiz´as´aval ´es a peremfelt´etelek elimin´al´as´aval nyert line´aris rendszer a k¨ovetkez˝o: Lh uh = fh , ahol 

Lh =

    1   h2     

1.3.2.

B −I

−I B

−I .. .

.

..

−I

B

−I

−I

B

..

.

     ,     





         

       .      

4 −1    −1 4 −1 



B=

−1

4 −1 .. .. .. . . . −1

4

−1

−1

4

V´ egeselem-m´ odszer

Vegy¨ uk a k´etdimenzi´os Poisson-egyenletet Dirichlet-peremfelt´etellel: −∆u = f

(Ω-ban),

u = 0 (Γ-n), A v´egeselem-m´odszer meg´ert´es´ehez sz¨ uks´eg¨ unk lesz a Green-formul´ara, mely szerint: Z

∇v · ∇u dx = −

Ω

Z

v∆u dx +

Ω

Z Γ

v

∂u ds. ∂~n

A fenti perem´ert´ek-probl´ema gyenge alakj´anak defini´al´as´ahoz vezess¨ uk be a k¨ovetkez˝o jel¨ol´eseket: a(u, v) :=

Z

∇u · ∇v dx

Ω

(f, v) :=

Z

f v dx.

Ω

Szorozzuk meg az egyenlet¨ unk mindk´et oldal´at egy tetsz˝oleges v ∈ H01 (Ω) tesztf¨ uggv´ennyel, majd integr´aljunk Ω-n. Ekkor azon f¨ uggv´enyekre, melyek teljes´ıtik a peremfelt´etelt, a Green-formula felhaszn´al´as´aval a(u, v) = −(∆u, v). 16

ad´odik. Mostm´ar megfogalmazhatjuk az eredeti probl´ema gyenge alakj´at: Keress¨ uk azt uggv´enyt, melyre az az u ∈ H01 (Ω) f¨ a(u, v) = (f, v) minden v ∈ H01 (Ω) vari´aci´os egyenlet teljes¨ ul, ahol a(u, v) = −(∆u, v). Ezen u f¨ uggv´enyt az eredeti perem´ert´ek-feladat gyenge megold´as´anak nevezz¨ uk (u egy´ertelm˝ u). A v´egeselem m´odszer a gyenge megold´as egy uh k¨ozel´ıt´es´et a´ll´ıtja el˝o. Sz¨ uks´eg¨ unk lesz egy Vh ⊂ H01 (Ω) v´eges dimenzi´os alt´erre, melyen az uh k¨ozel´ıt˝o megold´ast keress¨ uk. K´esz´ıts¨ uk el Ω egy felbont´as´at h´aromsz¨ogel´essel, ´es jel¨olj¨ uk Ωh -val a Kj h´aromsz¨ogek uni´ojak´ent el˝o´all´o tartom´anyt. Ekkor Vh alt´er a k¨ovetkez˝o f¨ uggv´enyekb˝ol ´ep¨ ul fel: Vh = {v : v|Ωh folytonos, v|∂Ωh = 0, v|Kj line´aris ∀j}. Tov´abb´a az al´abbi vi f¨ uggv´enyek Vh egy b´azis´at alkotj´ak: vi (xj ) =

 

1, ha xj = xi



0, ha xj 6= xi ,

ahol xj pontok a h´aromsz¨ogel´es csom´opontjai. Ekkor uh ezen Vh -beli b´azisf¨ uggv´enyek line´aris kombin´aci´ojak´ent a´ll el˝o. Vagyis keress¨ uk azokat a ci egy¨ utthat´okat, melyekre uh =

X

ci vi ´es a(uh , v) = (f, v) minden v ∈ Vh

i

teljes¨ ul, melyb˝ol egy line´aris egyenletrendszer ´ırhat´o fel, melynek megold´as´ara valamilyen numerikus megold´o m´odszert alkalmazunk. A hat´ekony numerikus probl´emamegold´ashoz egy megfelel˝oen pontos diszkretiz´aci´o ¨onmag´aban kev´es, ezt m´eg egy alkalmas megold´o algoritmussal is ki kell eg´esz´ıteni. Elliptikus parci´alis differenci´alegyenletek megold´as´ara az egyik leggyorsabb ilyen algoritmus a multigrid m´odszer. Ez alapvet˝oen a v´eges differencia m´odszerhez kapcsol´odik, de alkalmazhat´o a v´egeselem m´odszerrel kombin´alva is, b´ar ennek megval´os´ıt´asa komplik´altabb feladat.

17

2. fejezet A multigrid m´ odszer 2.1.

Diszkretiz´ aci´ o

Egy line´aris perem´ert´ek-probl´ema a k¨ovetkez˝o alakban ´ırhat´o fel: LΩ u = f Ω (x) (x ∈ Ω) LΓ u = f Γ (x) (x ∈ Γ). Itt Ω ⊂ Rn adott ´ertelmez´esi tartom´any Γ := ∂Ω hat´arral, LΩ line´aris differenci´aloper´ator Ω-n, LΓ line´aris oper´ator Γ-n, f Ω ´es f Γ adott f¨ uggv´enyek Ω-n ´es Γ-n, tov´abb´a u = u(x) jel¨oli a perem´ert´ek-probl´ema megold´as´at. Az el˝oz˝o differenci´alegyenlet ´es peremfelt´etel diszkretiz´aci´oj´aval nyert rendszer: Ω LΩ h uh = fh (x) (x ∈ Ωh )

LΓh uh = fhΓ (x) (x ∈ Γh ). Itt h a r´acsm´eret, uh a diszkr´et megold´as, mely egy Ωh ∪ Γh r´acson ´ertelmezett f¨ uggv´eny, fhΩ ´es fhΓ r´acsf¨ uggv´enyek (f Ω ´es f Γ diszkr´et anal´ogjai), valamint LΩ es LΓh h ´ r´acsoper´atorok, melyek r´acsf¨ uggv´enyeket r´acsf¨ uggv´enyekbe visznek. LΩ h -t differencia oper´atornak nevezz¨ uk. A peremfelt´etel elimin´al´as´aval az el˝oz˝o rendszer egyszer˝ ubb alakba ´ırhat´o: Lh uh = fh

(Ωh -ban).

(2.1)

Ez egy r´acsegyenlet, melyben uh ´es fh az Ωh r´acson ´ertelmezett f¨ uggv´enyek, tov´abb´a Lh : G(Ωh ) → G(Ωh ) line´aris oper´ator, ahol G(Ωh ) jel¨oli az Ωh -n ´ertelmezett r´acsf¨ uggv´enyek line´aris ter´et. Amennyiben Ω t´eglalap alak´ u, vagyis Ω = (0, A1 ) × (0, A2 ) ⊂ R2 , jel¨olje Gh := {(x, y) = (κ1 h1 , κ2 h2 ) : (κ1 , κ2 ) ∈ Z2 , hj :=

Aj , Nj ∈ N (j = 1, 2)}, Nj

a v´egtelen r´acsot, ahol h := (h1 , h2 ). Ekkor Ωh = Ω ∩ Gh teljes¨ ul. 18

A Poisson-egyenlet diszkretiz´ aci´ oja A legalapvet˝obb p´elda perem´ert´ek-feladatra a k´etdimenzi´os Poisson-egyenlet Dirichlet-peremfelt´etellel, az egys´egn´egyzeten. Ez a k¨ovetkez˝ok´eppen n´ez ki: LΩ u := −∆u = f Ω (x) (x ∈ Ω := (0, 1)2 ) LΓ u := u = f Γ (x) (x ∈ Γ), ahol ∆ := ∂x2 + ∂y2 . Diszkretiz´aljuk az egyenletet ´es a peremfelt´etelt egy h-n´egyzetr´acson a k¨oz¨ons´eges o¨tpontos approxim´aci´o seg´ıts´eg´evel. Ekkor a k¨ovetkez˝o ad´odik a differenciaoper´atorra: 

LΩ h = −∆h =

1 h2

   

−1



4

   

−1

−1

−1

h 2

Itt Ωh = Ω ∩ Gh , Gh := {(x, y) = (κ1 h, κ2 h) : κ ∈ Z , h =

1 ,N N

∈ N}. A diszkr´et

peremfelt´etel elimin´al´as´aval (2.1)-nek megfelel˝o r´acsegyenletet kapunk.

2.2.

A m´ odszer alap¨ otlete

A k¨ovetkez˝okben r´at´er¨ unk a multigrid m´odszer t´argyal´as´ara, melynek legnagyobb el˝onye, hogy seg´ıts´eg´evel egy N ismeretlen˝ u, diszkr´et elliptikus rendszer k¨ozel´ıt˝o megold´as´at O(N ) m˝ uvelettel megkaphatjuk, azaz a m´odszer nagys´agrend szerint optim´alis. Tekints¨ uk ´at most ezen m´odszer alap¨otlet´et. A differenci´alegyenlet ´es a peremfelt´etel diszkretiz´aci´oj´aval nyert egyenlet: Lh uh = fh . Legyen az uh diszkr´et megold´as tetsz˝oleges k¨ozel´ıt´ese ujh . Jel¨olje ezen j-edik k¨ozel´ıt´es hib´aj´at: vhj := uh − ujh , ´es a marad´ekot: djh := fh − Lh ujh . Ekkor az eredetivel ekvivalens egyenletet kapunk (marad´ekegyenlet): Lh vhj = djh . Amennyiben vhj -t csak k¨ozel´ıt˝oleg a´ll´ıtjuk el˝o (jel¨olje ezt vˆhj ) egy olyan egyenletrendszerb˝ol, melyet az eredeti feladat egy durv´abb r´acson val´o diszkretiz´al´as´aval nyer¨ unk, kisz´am´ıt´asa a kevesebb ismeretlen miatt kevesebb m˝ uveletet ig´enyel. Ekkor uj+1 = ujh + vˆhj , h 19

jav´ıtott megold´as m´ar egy jobb k¨ozel´ıt´ese uh -nak. Az elj´ar´ast iterat´ıv m´odon ism´etelj¨ uk, a finom Ωh ´es a durva ΩH r´acson dolgozva. S˝ot, egy´altal´an nem sz¨ uks´eges az j LH vˆH = djH egyenletet sem pontosan megoldani, ΩH -r´ol egy m´eg durv´abb r´acsra

t´erhet¨ unk a´t ´es ´ıgy tov´abb, eljutva ezzel a t¨obbr´acsos m´odszerhez. Ez a m´odszer m´eg o¨nmag´aban a´ltal´aban nem konvergens. Vezess¨ uk be a k¨ovetkez˝o jel¨ol´est: nH :=dimG(ΩH ) ´es nh :=dimG(Ωh ). Rendezz¨ uk az Lh ´es LH m´atrixok saj´atvektorait a hozz´ajuk tartoz´o saj´at´ert´ekek nagys´aga szerint n¨ovekv˝o sorrendbe. Ekkor elmondhat´o, hogy LH saj´atvektorai k¨ozel´ıtik Lh els˝o nH saj´atvektor´at, viszont az utols´o nh − nH saj´atvektor´anak nincsen megfeleltet´ese a durva r´acson. ´Igy, b´ar az alacsonyabb sorsz´am´ u saj´atvektoroknak megfelel˝o megold´asi komponensek ΩH -n gyorsabban el˝oa´ll´ıthat´ok, a magas sorsz´am´ u saj´atvektorok sz´am´ara nem szerezhet˝o hasznos inform´aci´o a durva r´acson. Teh´at az elj´ar´as konvergenci´aj´aban csak akkor rem´enykedhet¨ unk, ha a magas frekvenci´aj´ u saj´atvektorokkal kapcsolatos megold´asi komponensek l´enyeg´eben m´ar ujh -ben jelen vannak, azaz uh − ujh -b˝ol hi´anyoznak. M´ask´epp fogalmazva, ha a j-edik k¨ozel´ıt´es hib´aja sima abban az ´ertelemben, hogy l´enyeg´eben nincs benne magas frekvenci´aj´ u saj´atvektor. A k´es˝obbiekben megbizonyosodhatunk arr´ol, hogy erre a probl´em´ara kiel´eg´ıt˝o megold´ast jelentenek a sim´ıt´o iter´aci´ok. A m´odszer alapja teh´at a k´et egym´ast kieg´esz´ıt˝o m˝ uvelet: 1. hiba k¨ozel´ıt´ese durva r´acson 2. sim´ıt´o iter´aci´ok finom r´acson.

2.3.

Jav´ıt´ as durv´ abb r´ acson

Az el˝oz˝o r´eszben bemutatott iter´aci´o l´enyege teh´at az, hogy az eredetivel ekvivalens marad´ekegyenlet k¨ozel´ıt˝o megold´as´at u ´gy kapjuk, hogy a´tt´er¨ unk egy durv´abb r´acsra ´es az u ´jonnan kapott r´acsegyenletet oldjuk meg, mely j´oval kevesebb m˝ uveletbe ker¨ ul. A marad´ekegyenlet k¨ozel´ıt˝o megold´as´aval jav´ıtva az el˝oz˝o k¨ozel´ıt´est, az eredeti egyenletrendszer uh megold´as´anak egy jobb k¨ozel´ıt´es´et kapjuk. Teh´at az Lh vhj = djh egyenletr˝ol ´att´er¨ unk a durva r´acson diszkretiz´alt j LH vˆH = djH

20

egyenletre, ahol az LH : G(ΩH ) → G(ΩH ) oper´ator Lh alkalmas k¨ozel´ıt´ese a durv´abb j r´acson, valamint nH  nh teljes¨ ul, vˆH ´es djH pedig r´acsf¨ uggv´enyek ΩH -n.

A legfontosabb ´es leggyakrabban alkalmazott v´alaszt´as ΩH -ra az adott r´acsm´eret megdupl´az´as´ab´ol ad´odik, azaz H := 2h v´alaszt´assal. Sz¨ uks´eg¨ unk lesz k´et line´aris transzfer oper´atorra, melyek lehet˝ov´e teszik az a´tt´er´est a finom r´acsr´ol a durv´abbra ´es ford´ıtva: IhH : G(Ωh ) → G(ΩH ),

h IH : G(ΩH ) → G(Ωh ).

h IhH az u ´gynevezett restrikci´os vagy lesz˝ uk´ıt˝o oper´ator, melyre djH := IhH djh , ´es IH az h j vˆH . interpol´aci´os vagy m´asn´even prolong´aci´os oper´ator, melyre vˆhj := IH

N´ezz¨ unk most meg a leggyakrabban haszn´alatos restrikci´os ´es prolong´aci´os oper´atorokat egy-, illetve k´etdimenzi´os perem´ert´ek-probl´ema eset´en, H := 2h v´alaszt´assal. Egy prolong´aci´os oper´atort legegyszer˝ ubben a line´aris interpol´aci´o felhaszn´al´as´aval defini´alhatunk. Vegy¨ uk el˝osz¨or az egydimenzi´os esetet, amikoris a durva h´al´ot nH bels˝o pont ´es k´et hat´arpont alkotja. Ekkor a finomabb h´al´o a szomsz´edos ponj tok k¨ozti felez˝opontok hozz´av´etel´evel j¨on l´etre. Amennyiben a hiba vˆH k¨ozel´ıt´es´et

szeretn´ek ´attranszform´alni a durv´ar´ol a finomabb r´acsra, akkor a vˆhj k´epvektort ´ıgy defini´aljuk a finom h´al´o pontjaiban: vˆhj (x) Legyen nh :=

:= 1 h

j )(x) (pˆ vH

 

j vˆH (x), x ∈ ΩH =  ˆ[v j (x − h) + v j (x + h)]/2, x ∈ Ω \Ω h H H H

− 1, nH :=

1 H

− 1, ekkor az el˝oz˝o prolong´aci´ot le´ır´o nh × nH -as

m´atrix a k¨ovetkez˝o:

1 p= 2



1

                      

2



0

1 1 2 ..

. 2 1 1

0

2 1

           ,           

melyet egyszer˝ ubb alakban ´ıgy szok´as jel¨olni: p=

h 1 i 1 2 1 . 2

K´etdimenzi´oban ennek anal´ogj´ara a vektor egy adott finom r´acspontban felvett ´ert´eke a szomsz´edos durva r´acspontokban felvett ´ert´ekekb˝ol ker¨ ul kisz´am´ıt´asra 21

line´aris interpol´aci´oval: j vˆhj (0, 0) := vˆH (0, 0);

j vˆhj (0, H) := vˆH (0, H);

j vˆhj (H, 0) := vˆH (H, 0);

j vˆhj (H, H) := vˆH (H, H);

1 j 1 j (0, 0) + vˆH (0, H); vˆhj (0, h) := vˆH 2 2 1 j 1 j vˆhj (h, 0) := vˆH (0, 0) + vˆH (H, 0); 2 2 1 j 1 j (H, 0) + vˆH (H, H); vˆhj (H, h) := vˆH 2 2 1 j 1 j vˆhj (h, H) := vˆH (0, H) + vˆH (H, H); 2 2 1 j 1 j 1 j 1 j vˆhj (h, h) := vˆH (0, 0) + vˆH (0, H) + vˆH (H, 0) + vˆH (H, H); 4 4 4 4 A t¨obbi ilyen, 9 pontb´ol a´ll´o cell´an (azaz [k · H, k · H + H] × [l · H, l · H + H]-n, ahol k, l = 0, 1, 2, . . . nH ) ugyan´ıgy sz´am´ıtjuk ki a vektor ´ert´ekeit. Ezt kilenc pontos prolong´aci´onak nevezz¨ uk ´es a k¨ovetkez˝o egyszer˝ us´ıtett jel¨ol´es (stencil) seg´ıts´eg´evel reprezent´aljuk: 

p=

1 2 1

1   2 4 2 4 1 2 1

    

A prolong´aci´o megford´ıt´asa teh´at a restrikci´o, mely az adott vektort a a finomabbr´ol a durv´abb r´acsra k´epezi. Term´eszetesen a legegyszer˝ ubb p´elda restrikci´ora az identit´as (rinj ), hiszen minden durva r´acspont egyben finom r´acspont is: djH (x) := (rinj djh )(x) = djh (x) minden x ∈ ΩH ⊂ Ωh -ra. Egy gyakrabban alkalmazott p´elda a k¨ovetkez˝o s´ ulyozott restrikci´o, az u ´gynevezett ”full weighting”(FW), melynek egydimenzi´os v´altozata a marad´ekvektorra alkalmazva a k¨ovetkez˝o: 1 djH (x) := (rdjh )(x) = [djh (x − h) + 2djh (x) + djh (x + h)] minden x ∈ ΩH -ra. 4 Az el˝obb defini´alt restrikci´ot az al´abbi nH × nh -as m´atrix ´ırja le: 

 1 2 1   1 2 1 

r=

1  4    

0 ..

. 1 2 1

0

1 2 1

22

      ,     

mely az el˝oz˝oekben bemutatott prolong´aci´o m´atrix´anak transzpon´altj´anak 1/2-szerese. A fenti oper´atort le´ır´o stencil: r=

i 1h 1 2 1 . 4

Teh´at a p prolong´aci´o pT transzpon´altja, mely az ellent´etes ir´anyba k´epez, azaz a finomabb r´acsr´ol visz a durv´abbra, alkalmazhat´o restrikci´os oper´atork´ent. A FW restrikci´o k´etdimenzi´os anal´ogj´at ´eppen a m´ar bevezet´esre ker¨ ult kilenc pontos prolong´aci´o transzpon´altjak´ent kapjuk (konstans szorz´ot´ol eltekintve). Ezt a k¨ovetkez˝o egyszer˝ us´ıtett jel¨ol´essel adhatjuk meg: 

r=

1 2 1

1    2 4 2 16  1 2 1

   . 

Mostm´ar minden r´eszlet´eben ´attekintett¨ uk a durva r´acs alkalmaz´as´aval jav´ıt´o ¨ iter´aci´o m´odszer´et. Osszefoglalva, egy iter´aci´os l´ep´es (azaz hogyan kapjuk uj -b˝ol h

uj+1 h -t)

a k¨ovetkez˝ok´eppen fest:

1. A marad´ek kisz´am´ıt´asa: djh := fh − Lh ujh . 2. A marad´ek transzform´al´asa a durv´abb r´acsra: djH := IhH djh . j = djH . 3. A pontos megold´as kisz´am´ıt´asa a durv´abb r´acson: LH vˆH h j vˆH . 4. A hiba k¨ozel´ıt´es´enek transzform´al´asa a finomabb r´acsra: vˆhj := IH

5. Az u ´j k¨ozel´ıt´es kisz´am´ıt´asa: uj+1 := ujh + vˆhj . h Az elj´ar´ast a k¨ovetkez˝o ´abra szeml´elteti:

ujh → djh = fh − Lh ujh

vˆhj

↓ IhH djH

→

h ↑ IH

→

j LH vˆH = djH

2.1. ´abra. Jav´ıt´as durva r´acson 23

uj+1 = ujh + vˆhj h

Mint azt kor´abban meggondoltuk, a durva h´al´on jav´ıt´o iter´aci´o o¨nmag´aban a´ltal´aban m´eg csak nem is konvergens. Tudjuk, a probl´ema abb´ol ad´odik, hogy amikor a´tt´er¨ unk a durv´abb r´acsra, a vhj hiba egyes komponenesei ezen nem reprezent´alhat´ok, ´ıgy nem is cs¨okkenthet˝ok. Az iter´aci´o ezen h´atr´any´anak kik¨ usz¨ob¨ol´ese v´egett alkalmazzuk a sim´ıt´o iter´aci´okat, melyek nagyon hat´akonyan lesimitj´ak a hib´at, azaz a hiba magas frekvenci´aj´ u komponenseit gyorsan cs¨okkentik. ´Igy a k´et iter´aci´o kombin´aci´oja rendk´ıv¨ ul gyorsan konverg´al annak ellen´ere, hogy k¨ ul¨on-k¨ ul¨on igen lassan, illetve egy´altal´an nem. Ezen kombin´aci´ot k´etr´acsos iter´aci´onak nevezz¨ uk.

2.4.

A sim´ıt´ o iter´ aci´ ok ´ es sim´ıt´ o tulajdons´ agaik

Az iter´ aci´ ok alapjai Vegy¨ uk a k¨ovetkez˝o r´acsegyenletet: Ll ul = fl

(Ωl -ben).

(2.2)

Ennek megold´asa egy ϕl line´aris lek´epez´es a´ltal gener´alt iterat´ıv elj´ar´assal nyerhet˝o: u0l 7→ u1l 7→ . . . 7→ ulj−1 7→ ujl 7→ . . . , melyre uj+1 = ϕl (ujl , fl ) l teljes¨ ul. A klasszikus iter´aci´okat ´es a multigrid iter´aci´ot is ´ıgy reprezent´aljuk. Mivel ϕl line´aris lek´epez´es, ez´ert ϕl (ul , fl ) = Ml ul + Nl fl form´aban is magadhat´o, ´ıgy az iter´aci´o = Ml ujl + Nl fl uj+1 l

(2.3)

alakba ´ırhat´o. Ml -t iter´aci´os m´atrixnak nevezz¨ uk. Tov´abb´a k´ezenfekv˝o felt´etel, hogy a (2.2)-es egyenlet ul megold´asa az iter´aci´onak fixpontja legyen, vagyis ul = Ml ul + Nl fl

(2.4)

teljes¨ ulj¨on. Ebb˝ol azt kapjuk, hogy I = Ml + Nl Ll , ahol I az identit´as m´atrix. Amennyiben Ll nem szingul´aris, Nl egyszer˝ uen kifejezhet˝o: Nl = (I − Ml )L−1 . ´Igy l

az iter´aci´o megadhat´o csak az iter´aci´os m´atrix´aval is: uj+1 = Ml ujl + (I − Ml )L−1 l l fl . Kivonva (2.4)-b˝ol (2.3)-at ul − uj+1 = Ml (ul − ujl ) j = 0, 1, 2, . . . l 24

ad´odik, ´ıgy a j-edik k¨ozel´ıt´es hib´aj´ara ul − ujl = Mlj (ul − u0l ) teljes¨ ul. Teh´at a hiba cs¨okken´ese csak az iter´aci´os m´atrixt´ol f¨ ugg. Ismeretes, hogy egy A m´atrix spektr´alsugar´at ´ıgy defini´aljuk: %(A) := max{|λ| : λ saj´at´ert´eke A-nak}. Ennek tudat´aban megeml´ıtend˝o az al´abbi, ismert lemma, mely a k´es˝obbiek folyam´an m´eg fontos lesz a sz´amunkra. 2.4.1. Lemma. A (2.3)-as iter´aci´o pontosan akkor konverg´al b´armely kezdeti u0l k¨ozel´ıt´esre, ha %(Ml ) < 1. Az iter´aci´os m´atrix spektr´alsugar´at ez´ert a konvergencia rendj´enek nevezz¨ uk.

N´ eh´ any klasszikus iter´ aci´ o A (2.3) t´ıpus´ u iter´aci´ok szerkezet´et a´ltal´anosan a k¨ovetkez˝ok´eppen is megadhatjuk: Vegy¨ uk az Al ul = fl egyenletrendszert. Tegy¨ uk fel, hogy Al = Bl − Cl , ahol Bl nem szingul´aris. Legyen Bl uj+1 − Cl ujl = fl , ebb˝ol a k¨ovetkez˝o iter´aci´o ad´odik: l uj+1 = Bl−1 Cl ujl + Bl−1 fl . l Teh´at (2.3)-ban Ml = Bl−1 Cl ´es Nl = Bl−1 . Vegy¨ uk Al egy m´asik felbont´as´at: Al = Dl −Ll −Ul , ahol Dl =diag(a11 , a22 , ..., ann ) az Al diagon´alelemeib˝ol alkotott m´atrix, −Ll ´es −Ul pedig Al szigor´ u als´o ´es fels˝o h´aromsz¨ogm´atrixai. Ekkor Bl = ω1 Dl ´es Cl = ω1 [(1 − ω)Dl + ω(Ll + Ul )], 0 < ω ≤ 1 v´alaszt´assal ´eppen a csillap´ıtott Jacobi iter´aci´ot kapjuk: uj+1 = (I − ωDl−1 Al )ujl + ωDl−1 fl , l valamint a Bl = Dl − Ll ´es Cl = Ul v´alaszt´as a Gauss-Seidel iter´aci´ot eredm´enyezi: uj+1 = (Dl − Ll )−1 Ul ujl + (Dl − Ul )−1 fl . l

A klasszikus iter´ aci´ ok sim´ıt´ o tulajdons´ aga Sim´ıt´o iter´aci´oknak megfelel˝ok lesznek a klasszikus iter´aci´ok. B´ar lassan konverg´alnak, de az adott k¨ozel´ıt´es hib´aj´at gyorsan lesim´ıtj´ak, azaz a magasabb sorsz´am´ u saj´atvektorokkal kapcsolatos megold´asi komponenseket j´oval gyorsabban pontos´ıtj´ak, mint az alacsony sorsz´am´ uak´et. 25

Az iter´aci´ok sim´ıt´o tulajdons´ag´anak analiz´al´as´ahoz vegy¨ uk p´eldak´ent (lsd. [6]) az egydimenzi´os Poisson-egyenletet a legegyszer˝ ubb Dirichlet-peremfelt´etellel: −u00 (x) = f (x) Ω = (0, 1) u(x) = 0 x ∈ Γ = {0, 1}. Defini´aljuk r´acsm´eretek egy sorozat´at: 1 1 1 h0 > h1 > h2 > ... > hl > ..., ahol h0 = , h1 = , . . . , hl = l+1 , . . . 2 4 2 l-et szintsz´amnak nevezz¨ uk. Legyen az l-edik szinthez tartoz´o r´acs: Ωl := {khl : k = 1, 2, ..., 2l+1 − 1} (l = 0, 1, 2, . . .). Diszkretiz´aljuk a perem´ert´ek-feladatot a v´eges differenci´ak m´odszer´evel. Ekkor l ul := (ul (khl ))nk=1 ´es

l fl := (fl (khl ))nk=1

(nl =

1 − 1, az ismeretlenek sz´ama) hl

vektorokkal a k¨ovetkez˝o r´acsegyenleteket kapjuk (l = 0, 1, 2, . . .): 



         

        ul      

2 −1    −1 2 −1  1 h2l

−1

2 ..

−1 . . . .. .. −1

2

−1

−1

2

= Ll ul = fl

Vizsg´aljuk most a csillap´ıtott Jacobi iter´aci´o sim´ıt´o tulajdons´ag´at ezen rendszer seg´ıts´eg´evel. A Jacobi iter´aci´o, ahogy azt m´ar l´attuk, a k¨ovetkez˝o alakban ´ırhat´o: uj+1 = (I − ΘDl−1 Ll )ujl + ΘDl−1 fl l

Θ ∈ (0, 1).

Ezt ´atalak´ıtva uj+1 = ujl − ΘDl−1 (Ll ujl − fl ) Θ ∈ (0, 1) l ad´odik. Jelen esetben Dl−1 =

h2l I, 2

´ıgy ω = Θ/2-vel azt kapjuk, hogy

1 uj+1 = ujl − ωh2l (Ll ujl − fl ) ω ∈ (0, ). l 2 Jel¨olj¨ uk ezen Jacobi iter´aci´ohoz tartoz´o iter´aci´os m´atrixot Ml -el. A konvergencia vizsg´alat´ahoz sz¨ uks´eg¨ unk van ennek saj´atvektoraira ´es saj´at´ert´ekeire, melyekre: Ml vlµ = λµl vlµ

(µ = 1, . . . nl ). 26

El˝osz¨or n´ezz¨ uk Ll saj´atvektorait, ezek a k¨ovetkez˝ok: vlµ =

√

l 2hl (sin(kµπhl ))nk=1

(µ = 1, . . . nl ).

Ekkor Ll vlµ =

4 sin2 (µπhl /2)vlµ h2l

(µ = 1, . . . nl )

miatt ´es abb´ol, hogy Ml = I − ωh2l Ll , azt kapjuk, hogy Ml saj´atvektorai megegyeznek Ll saj´atvektoraival ´es a megfelel˝o saj´at´ert´ekek: λµ (ω) = 1 − 4ω sin2 (µπhl /2) (µ = 1, . . . nl ).

2.2. ´abra. Ml iter´aci´os m´atrix saj´at´ert´ekei µhl f¨ uggv´eny´eben A (2.2)-es a´br´ar´ol leolvashat´o, hogy ω =

1 2

eset´en a legnagyobb abszol´ ut´ert´ek˝ u

saj´at´ert´ek λ1 , ´ıgy a konvergencia rendje: %(Ml ) = λ1 ( 12 ) = 1 − 2 sin2 (πhl /2) = 1 − 12 π 2 h2l + O(h4l ), azaz ebben az esetben a Jacobi iter´aci´o nagyon lassan konverg´al. K¨onnyen ad´odik, hogy ω ∈ (0, 21 ) v´alaszt´as m´eg lassabb konvergenci´at eredm´enyez. Vegy¨ uk most az ω = gas frekvenci´aj´ u

vlµ

1 4

esetet. A (2.2)-es a´bra alapj´an elmondhatjuk, hogy a ma-

saj´atvektorok (µ ≥ 1/(2hl )) minden iter´aci´os l´ep´esben legal´abb

1/2-szeres´ere cs¨okkennek. Teh´at ha csak a magas frekvenci´aj´ u saj´at´ert´ekeket vessz¨ uk figyelembe, a csillap´ıtott Jacobi iter´aci´o nagyon gyorsan konverg´al. Ebb˝ol azt a 27

k¨ovetkeztet´est vonhatjuk le, hogy a lass´ u konvergenci´a´ert az alacsony frekvenci´aj´ u komponensek a felel˝osek. Vizsg´aljuk meg ezut´an a hiba cs¨okken´es´et. ´Irjuk fel a kezdeti k¨ozel´ıt´es hib´aj´at a saj´atvektorok seg´ıts´eg´evel: u0l − ul =

X

αµ vlµ .

Ebb˝ol j iter´aci´os l´ep´es elv´egz´ese ut´an a j-edik k¨ozel´ıt´es hib´aj´ara a k¨ovetkez˝ot kapjuk: ujl − ul = Mlj (u0l − ul ) =

P

αµ Mlj vlµ =

P

αµ [λµ ( 41 )]j vlµ =

P

βµ vlµ .

Az el˝oz˝o megfontol´as szerint αµ ≈ βµ teljes¨ ul az alacsony, |βµ |  |αµ | a magas frekvenci´as komponensekre. Ekkor azt mondjuk, hogy az ujl − ul hiba sim´abb, mint az u0l − ul . Azaz a csillap´ıtott Jacobi iter´aci´o sim´ıt´o iter´aci´ok´ent funkcion´al. (Megjegyezend˝o, hogy a t´ ulrelax´alt Jacobi-iter´aci´o (ω > 1) m´ar nem rendelkezik a k´ıv´ant sim´ıt´asi tulajdons´aggal.) A klasszikus iter´aci´ok sim´ıt´o tulajdons´ag´anak a´ltal´anos esetben t¨ort´en˝o igazol´as´a´ert lsd. [3].

2.5.

K´ etr´ acsos m´ odszer

¨ Osszegezve teh´at az eddig tapasztaltakat, azt kaptuk, hogy ´erdemes kombin´alnunk a k´et egym´ast j´ol kieg´esz´ıt˝o elj´ar´ast, az iter´aci´os sim´ıt´ast ´es a durva r´acson t¨ort´en˝o hibajav´ıt´ast. Ez adja az alapj´at a k´etr´acsos m´odszernek. K´et r´acson dolgozunk, egy finomabb Ωh ´es durv´abb ΩH r´acson Az elj´ar´ast iterat´ıv m´odon ism´etelj¨ uk. Egy iter´aci´os l´ep´es (azaz hogy kapjuk ujh -b˝ol uj+1 ovetkez˝o h -t) a k¨ algoritmussal adhat´o meg: 1. Iter´aci´os sim´ıt´as I. – ujh kisz´am´ıt´asa ujh -b˝ol az adott sim´ıt´o iter´aci´o ν1 -szer t¨ort´en˝o alkalmaz´as´aval: ujh := Shν1 (ujh , Lh , fh ). 2. Hibajav´ıt´as a durva r´acson j

– A marad´ek kisz´am´ıt´asa: dh := fh − Lh ujh . j

j

– A marad´ek transzform´al´asa a durv´abb r´acsra: dH := IhH dh . j

j – A pontos megold´as kisz´am´ıt´asa a durv´abb r´acson: LH vˆH = dH . h j – A hiba k¨ozel´ıt´es´enek transzform´al´asa a finomabb r´acsra: vˆhj := IH vˆH .

– A jav´ıtott k¨ozel´ıt´es kisz´am´ıt´asa: ujh + vˆhj . 28

3. Iter´aci´os sim´ıt´as II. – uj+1 kisz´am´ıt´asa ujh + vˆhj -b˝ol az adott sim´ıt´o iter´aci´o ν2 -szer t¨ort´en˝o alh kalmaz´as´aval: uj+1 := Shν2 (ujh + vˆhj , Lh , fh ). h

ν1 ujh

→

ν2 ujh

→

j dh

= fh −

Lh ujh

vˆhj

j

+

vˆhj

→

uj+1 h

h ↑ IH

↓ IhH

j

j LH vˆH = dH

→

dH

→

ujh

2.3. ´abra. K´etr´acsos m´odszer Az elj´ar´ast a (2.3)-as a´bra szeml´elteti. Err˝ol k¨onnyen leolvashatjuk az iter´aci´os m´atrixot: MhH = Shν2 KhH Shν1 ,

ahol

h −1 H KhH = Ih − IH LH Ih Lh ,

´es Sh a sim´ıt´o iter´aci´onak megfelel˝o opre´ator. ¨ Osszegezz¨ uk most k´etr´acsos m´odszer azon komponenseit, melyek nincsenek el˝ore meghat´arozva (ezen komponensek megv´alaszt´asa mindig a megoldand´o feladatt´ol f¨ ugg): – a sim´ıt´o iter´aci´o (Sh ) – az el˝o- ´es ut´osim´ıt´asok v´egrehajt´asainak sz´ama (ν1 , ν2 ) – a durva r´acs (ΩH ) – a restrikci´os oper´ator (IhH ) – az Lh oper´ator megfeleltet´ese a durva r´acson (LH ) h – az interpol´aci´os oper´ator (IH )

A c´elunk a m´odszer ezen komponenseit u ´gy megv´alasztani, hogy az algoritmus optim´alis legyen. A v´alaszt´asok jelent˝osen befoly´asolj´ak az algoritmus hat´ekonys´ag´at. Sajnos nincsen a´ltal´anos ´erv´eny˝ u, minden feladatra optim´alisan m˝ uk¨od˝o elj´ar´as, ezen komponensek megv´alaszt´asa az optimalit´as ´erdek´eben mindig az adott feladat f¨ uggv´eny´eben kell hogy t¨ort´enjen. 29

A m´odszer hat´ekonys´ag´anak vizsg´alat´ahoz sz¨ uks´eg van a %(MhH ) ´ert´ek, azaz az aszimptotikus konvegencia faktor vagy a megfelel˝o kMhH k norma meghat´aroz´as´ara. A m´odszer egyes komponenseinek megv´alaszt´asakor azt kell figyelembe venni, hogy azok milyen hat´assal vannak %(MhH ) ´es kMhH k ´ert´ekekre, melyekb˝ol teh´at a m´odszer hat´ekonys´ag´ara vonatkoz´o k¨ovetkeztet´eseket vonhatunk le. A k´etr´acsos m´odszer alapj´aul szolg´al a multigrid m´odszernek, mely az el˝oz˝o algoritmus azon alapul´o tov´abbfejleszt´ese, hogy nem csak egy, hanem t¨obb durv´abb r´acsot alkalmazunk. Azaz nem sz¨ uks´eges a durva r´acson diszkretiz´alt marad´ekegyenletet pontosan megoldanunk, hanem elegend˝o el´ern¨ unk egy k¨ozel´ıt˝o megold´ast, egy u ´jabb k´etr´acsos m´odszer, egy m´eg durv´abb r´acson t¨ort´en˝o alkalmaz´as´aval.

2.6.

A multigrid ciklus

Eddig csak a multigrid m´odszer k´etr´acsos verzi´oj´at l´attuk, melynek kev´esb´e gyakorlati, de ann´al ink´abb elm´eleti jelent˝os´ege van, hiszen a val´odi multigrid elj´ar´as alapj´aul szolg´al.

A m´odszerhez az alap¨otletet az a t´eny adja, hogy egy adott

k´etr´acsos elj´ar´as sor´an keletkez˝o j

j = dH LH vˆH

marad´ekegyenletet nem sz¨ uks´eges pontosan megoldanunk, ehelyett egy megfelel˝o k¨ozel´ıt˝o megold´as megad´asa is elegend˝o. Szerencs´ere ez kivitelezhet˝o a kovergenciasebess´eg l´enyeges cs¨okken´ese n´elk¨ ul. Term´eszetes m´odon ad´odik az a lehet˝os´eg, hogy ezt a k¨ozel´ıt˝o megold´ast egy u ´jabb k´etr´acsos algoritmussal a´ll´ıtsuk el˝o, ahol egy m´eg durv´abb r´acsot haszn´alunk. Vil´agos, hogy amennyiben a k´etr´acsos m´odszer¨ unk konvergencia faktora el´eg kicsi, elegend˝o mind¨ossze n´eh´any (γ) iter´aci´os l´ep´est alkalmaznunk egy megfelel˝o k¨ozel´ıt´eshez el´er´es´ehez. K´ezenfekv˝o, hogy ez az o¨tlet rekurz´ıvan alkalmazhat´o egyre durv´abb ´es durv´abb r´acsokat haszn´alva, eljutva ezzel egy olyan r´acsegyenlethez, melynek pontos megold´asa m´ar elegend˝oen kev´es m˝ uveletbe ker¨ ul. N´ezz¨ uk meg most a m´odszer form´alis le´ır´as´at. Ehhez sz¨ uks´eg¨ unk van egyre finomabb r´acsok egy Ωhl sorozat´ara, ahol hl (l = 0, 1, 2, . . .) a r´acsm´eretek sorozata. Az egyszer˝ us´eg kedv´e´ert az elk¨ovetkez˝okben hl helyett l-et ´ırunk. Minden Ωl -re legyenek adottak a k¨ovetkez˝o line´aris r´acsoper´atorok: Ll : G(Ωl ) → G(Ωl ),

Sl : G(Ωl ) → G(Ωl ),

30

Ill−1 : G(Ωl ) → G(Ωl−1 ),

l Il−1 : G(Ωl−1 ) → G(Ωl ),

´es a diszkr´et r´acsegyenletek: Ll ul = fl

(Ωl -ben),

ahol Ll invert´alhat´o, G(Ωl ) az Ωl -en ´ertelmezett r´acsf¨ uggv´enyek tere, Sl pedig az adott simit´o iter´aci´onak megfelel˝o oper´atort jel¨oli. γ jel¨oli teh´at, hogy h´any multigrid iter´aci´os l´ep´est kell v´egrehajtanunk az egyes durv´abb r´acsokon. A tapasztalatok azt mutatj´ak, hogy nem ´erdemes γ-t 2-n´el nagyobbra v´alasztani, mert b´ar pontosabb k¨ozel´ıt´est kapunk az adott szinten, de egy u ´jabb iter´aci´os l´ep´es m´ar annyival n¨oveli a l´ep´essz´amot, hogy nem nyer¨ unk, hanem vesz´ıt¨ unk a pontosabb k¨ozel´ıt´es el˝oa´ll´ıt´as´aval, ´ıgy a konvergencia sebess´ege romlik. A γ = 1 esetet V-ciklusnak, a γ = 2 esetet W-ciklusnak nevezz¨ unk. Most pedig n´ezz¨ uk meg a multigrid iter´aci´o -vagyis az (l + 1)-r´acsos iter´aci´o egy l´ep´es´et (multigrid ciklus) az Ll ul = fl

(Ωl -ben),

r´acsegyenlet megold´as´ara, ahol l ≥ 1 fix. Ehhez sz¨ uks´eg¨ unk lesz Ωk r´acsokra, k oper´atorokra (k = Lk oper´atorokra (k = l, l − 1, . . . , 0), valamint Sk , Ikk−1 , Ik−1

l, l − 1, . . . , 1). Legyenek tov´abb´a ν1 , ν2 , γ fix param´eterek. A multigrid ciklust egy o¨nmag´at megh´ıv´o, rekurz´ıv algoritus defini´al´as´aval kapjuk. k¨ozel´ıt´es kisz´am´ıt´asa a ´j uj+1 Az ul megold´as egy adott ujl k¨ozel´ıt´es´eb˝ol az u l k¨ovetkez˝o algoritmussal t¨ort´enik: Ha l = 1: a k´etr´acsos m´odszer alkalmaz´asa Ω1 , Ω0 r´acsokra Ha l > 1: 1. Iter´aci´os sim´ıt´as I. – ujl kisz´am´ıt´asa ujl -b˝ol az adott sim´ıt´o iter´aci´o ν1 -szer t¨ort´en˝o alkalmaz´as´aval: ujl := Slν1 (ujl , Ll , fl ). 2. Hibajav´ıt´as a durva r´acson j

– A marad´ek kisz´am´ıt´asa: dl := fl − Ll ujl . j

j

– A marad´ek transzform´al´asa a durv´abb r´acsra: dl−1 := Ill−1 dl .

31

– Az j

j Ll−1 vˆl−1 = dl−1

(Ωl−1 )

j egyenlet v˜l−1 k¨ozel´ıt˝o megold´as´as´anak kisz´am´ıt´asa az l-r´acsos m´odszerrel,

γ iter´aci´os l´ep´essel az azonosan 0 r´acsf¨ uggv´enyt v´eve kezdeti k¨ozel´ıt´esnek (a Ωl−1 , Ωl−2 , . . . , Ω0 r´acsok ´es a megfelel˝o r´acsf¨ uggv´enyek haszn´alat´aval). j l v˜l−1 . – A hiba k¨ozel´ıt´es´enek transzform´al´asa a finomabb r´acsra: v˜lj := Il−1

– A jav´ıtott k¨ozel´ıt´es kisz´am´ıt´asa: ujl + v˜lj . 3. Iter´aci´os sim´ıt´as II. – uj+1 kisz´am´ıt´asa ujl + v˜lj -b˝ol az adott sim´ıt´o iter´aci´o ν2 -szer t¨ort´en˝o all kalmaz´as´aval: uj+1 := Slν2 (ujl + v˜lj , Ll , fl ). l Az algoritmus lefut´as´at nyomon k¨ovethetj¨ uk a (2.4)-es a´bra seg´ıts´eg´evel. Itt 0 ≤ C(k) ≤ γ sz´aml´alja, hogy az egyes szinteken eddig h´any multigrid iter´aci´os l´ep´es t¨ort´ent. ν1 -et, ν2 -t ´es γ-t fix param´etereknek v´alasztottuk, de ez nem volt felt´etlen¨ ul sz¨ uks´eges. Sok esetben ugyanis γ-t k f¨ uggv´eny´eben v´alasztjuk, p´eld´aul γ = 1, ha k p´aratlan, γ = 2, ha k p´aros. A most defini´alt multigrid elj´ar´as iter´aci´os m´atrixa: γ l−1 l )L−1 Ll )Slν1 . Ml = Slν2 (Il − Il−1 (Il−1 − Ml−1 l−1 Il

Erre a k´es˝obbiek sor´an, a teljes multigrid m´odszer konvergenci´aj´anak t´argyal´asakor lesz sz¨ uks´eg¨ unk. Ezzel el is jutottunk az utols´o l´ep´eshez, a teljes multigrid m´odszer k¨ovetkezik, mely a multigridt˝ol elt´er˝oen m´ar nem iterat´ıv elj´ar´as.

2.7.

A teljes multigrid m´ odszer

A diszkr´et, elliptikus perem´ert´ek-probl´em´ak megold´as´ara alkalmazott iterat´ıv multigrid elj´ar´ast m´ar defini´altuk. Hab´ar ezen m´odszer megfelel˝o verzi´oi is nagyon hat´ekonyak lehetnek, hat´ekonys´aguk m´eg tov´abb n¨ovelhet˝o, amennyiben teljes multigrid m´odszerben form´aj´aban ker¨ ulnek felhaszn´al´asra. A teljes multigrid algoritmus azon az o¨tleten alapszik, hogy a megold´as u ´j k¨ozel´ıt´es´et n´eh´any multigrid iter´aci´os l´ep´es elv´egz´es´evel kapjuk, az (interpol´alt) el˝oz˝o k¨ozel´ıt´est v´alasztva a multigrid iter´aci´o kezdeti k¨ozel´ıt´es´enek, melyet szint´en n´eh´any hasonl´o, de az eggyel durv´abb r´acson alkalmazott multigrid iter´aci´os l´ep´essel kaptunk. Teh´at az alacsonyabb szinteken nagyon j´o kezdeti k¨ozel´ıt´esre tehet¨ unk szert 32

2.4. ´abra. Egy multigrid iter´aci´os l´ep´es az Ll ul = fl (l ≥ 1) megold´as´ara

33

a k¨ovetkez˝o magasabb szinten elv´egzend˝o iter´aci´ohoz. Adott az al´abbi LΩ u = f Ω Γ

L u=f

Γ

(Ω-ban), (Γ := ∂Ω ment´en),

elliptikus perem´ert´ek-probl´ema korl´atos Ω tartom´anyon. Ezen perem´ert´ek-feladat diszkr´et k¨ozel´ıt´esei a k¨ovetkez˝o r´acsegyenletek: Ll ul = fl

(Ωl -ben) (l = 0, 1, 2, ...)

Minden l ≥ 1-re jel¨olj¨ uk ´ıgy: MGIµ (·, l, Ll , fl ) : G(Ωl ) → G(Ωl ) a µ iter´aci´os l´ep´esb˝ol a´ll´o, Ω0 , . . . Ωl r´acsokon alkalmazott, adott iterat´ıv multigrid elj´ar´ast. A elk¨ovetkez˝okben legyen µ fix (mint eml´ıtett¨ uk, µ a´ltal´aban 1 vagy 2, ´es esetleg f¨ ugghet l-t˝ol is). Legyen l ≥ 0 fix. Ekkor az Ll ul = fl (Ωl -ben) r´acsegyenlet k¨ozel´ıt˝o megold´as´anak el˝o´all´ıt´as´ara szolg´al´o teljes multigrid m´odszer (FMG) algoritmusa a k¨ovetkez˝o: 1. k := 0,

u˜0 : az L0 u0 = f0 pontos megold´asa

2. k := k + 1 k 3. uˆk := Ik−1 u˜k−1

4. u˜k := MGIµ (ˆ uk , k, Lk , fk ) 5. ha k < l : → 2. ha k = l :

u˜l : az Ll ul = fl megold´as´anak FMG elj´ar´assal nyert k¨ozel´ıt´ese

k Itt Ik−1 : G(Ωk−1 ) → G(Ωk ) egy megfelel˝o interpol´aci´os elj´ar´as.

Az FMG algoritmus rendelkezik a k¨ovetkez˝o tulajdons´agokkal: • k˜ uh − uh k kisebb, mint a ku − uh k diszkretiz´aci´os hiba, ahol u˜h a diszkr´et uh megold´as FMG m´odszerrel nyert k¨ozel´ıt´ese. • A m˝ uveletig´eny ar´anyos a Ωh r´acs pontjainak sz´am´aval: O(Nk ) (optim´alis l´ep´essz´am).

34

A teljes multigrid m´ odszer konvergenci´ aja Az el˝oz˝oek alapj´an legyen Lk uk = fk a k-adik szinthez tartoz´o egyenletrendszer k az az interpol´aci´o, mely az adott k¨ozel´ıt´est a k − 1-edik szintr˝ol (0 ≤ k ≤ l), Ik−1

a k-adikra viszi, valamint µ az MGI elj´ar´as h´ıv´asainak sz´ama a mindenkori legmagasabb szinten. Jel¨olje tov´abb´a κ a diszkretiz´aci´o rendj´et, mely az alkalmazott differencias´ema pontoss´ag´at mutatja, illetve legyen ζk := kMk k, ahol Mk a k-adik szinthez tartoz´o multigrid elj´ar´as iter´aci´os m´atrixa ´es ζ := max0≤k≤l ζk . 2.7.1. T´ etel. (Hackbusch): Minden 0 ≤ k ≤ l-re legyen hk = 1/2k+1 . Legyen µ olyan, hogy ζ µ ≤ 1/(1 + 2κ ), ´es tegy¨ uk fel, hogy teljes¨ ulnek az al´abbi felt´etelek: k 1. kIk−1 k≤1 k 2. kIk−1 uk−1 − uk k ≤ c0 hκk , ahol c0 = c0 (u) konstans az eredeti feladat u pontos

megold´as´at´ol f¨ ugg 3. a k-adik szinten az MGIµ (·, k, Lk , fk ) elj´ar´as az azonosan 0 kezdeti k¨ozel´ıt´essel ζk -ra cs¨okkenti a hib´at, ´es ζ < 1 teljes¨ ul. Ekkor a teljes multigrid m´odszer konverg´al, ´es minden u˜k k¨ozel´ıt´es hib´aj´ara k˜ uk − uk k ≤ c0 hκk , ad´odik, ahol uk az Lk uk = fk egyenletrendszer megold´asa. Ez azt jelenti, hogy a k-adik szinten O(Nk ) m˝ uvelettel el˝oa´ll´ıtott u˜k k¨ozel´ıt´es hib´aja O(hκk ) (ahol Nk = nk · nk a diszkr´et ismeretlensz´am a k-adik szinten). A t´etel bizony´ıt´as´a´ert lsd. [3]. A multigrid elj´ar´as elm´eleti h´atter´enek a´ttekint´es´evel v´egezt¨ unk, a k¨ovetkez˝o fejezetben r´at´er¨ unk a m´odszer gyakorlati hat´ekonys´ag´anak bemutat´as´ara. A multigrid m´odszerr˝ol a tov´abbi r´eszletek´ert lsd. [1, 6, 5, 3, 7].

35

3. fejezet A multigrid m´ odszer implement´ al´ asa MATLAB-ban 3.1.

Multigrid program

Az fmgcycle.m f´ajlban egy teljes multigrid m´odszert defini´altunk, mely az al´abbi k´et perem´ert´ek-feladat megold´as´ara alkalmas: − ∆u = 2π 2 sin(πx) sin(πy) x ∈ Ω := (0, 1)2 ,

(3.1)

u(x) = 0 x ∈ Γ := ∂Ω,

− ∆u = 13π 2 sin(2πx) sin(3πy) x ∈ Ω := (0, 1)2 ,

(3.2)

u(x) = 0 x ∈ Γ := ∂Ω. Teh´at k´et speci´alis, az egys´egn´egyzeten ´ertelmezett, k´etdimenzi´os Poisson-egyenlet k¨ozel´ıt˝o megold´as´at keress¨ uk a megadott Dirichlet-peremfelt´etel teljes¨ ul´ese eset´en. A programban a r´acsm´eretek a szintek n¨oveked´es´evel felez˝odnek, azaz a k-adik szinthez tartoz´o r´acsm´eret: hk = 1/2k+1 . Megadhatjuk a programnak, hogy h´any szintet haszn´aljon, h´anyszor sim´ıtson (el˝o- ´es ut´osim´ıt´as), valamint h´any multigrid ciklust alkalmazzon az egyes szinteken, ´es hogy melyik tesztfeladatot szeretn´enk megoldani ((3.1) vagy (3.2)). Azaz input param´eterei a k¨ovetkez˝ok: – l: szintsz´am – nu1: el˝osim´ıt´asok sz´ama – nu2: ut´osim´ıt´asok sz´ama – gamma: multigrid iter´aci´ok sz´ama az egyes szinteken 36

– t: tesztf¨ uggv´eny sorsz´ama Az u output a megfelel˝o tesztfeladathoz tartoz´o k¨ozel´ıt˝o megold´as, w pedig egy seg´edoutput, melyre az csak egyik tesztfeladatban lesz sz¨ uks´eg (arra, hogy mi is ez pontosan, majd k´es˝obb t´er¨ unk ki). A program a m´odszer r´eszegys´egeinek megfelel˝oen a k¨ovetkez˝o alprogramokb´ol ´ep¨ ul fel: – laplace.m: az adott szinthez tartoz´o Laplace-oper´atort a´ll´ıtja el˝o a k¨oz¨ons´eges o¨tpontos approxim´aci´o seg´ıts´eg´evel – input: l: szintsz´am – output: L: Laplace-oper´ator – rhs.m: a r´acsegyenlet jobboldali vektor´at a´ll´ıtja el˝o a szintsz´am ´es a tesztfeladat f¨ uggv´eny´eben – input: l: szintsz´am, t: tesztf¨ uggv´eny sorsz´ama – output: f : jobboldali vektor – pre smooth.m: n´eh´any Jacobi-iter´aci´os l´ep´est v´egez (el˝osim´ıt´as) – input: l: szintsz´am, u: bemen˝o vektor, f : jobboldali vektor, nu: sim´ıt´asok sz´ama – output: u: a lesim´ıtott vektor – post smooth.m: n´eh´any Gauss-Seidel iter´aci´os l´ep´est v´egez (ut´osim´ıt´as) – input: l: szintsz´am, u: bemen˝o vektor, f : jobboldali vektor, nu: sim´ıt´asok sz´ama – output: u: a lesim´ıtott vektor – restri.m: az adott vektort a el˝oz˝o durv´abb szintre viszi (l-edikr˝ol az l-1-edikre) – input: l: szintsz´am, r in: bemen˝o vektor – output: r out: a bemen˝o vektorb´ol kapjuk ”full weighting” restrikci´oval – prolong.m: az adott vektort a k¨ovetkez˝o finomabb szintre viszi (l-1-edikr˝ol az l-edikre) – input: l: szintsz´am, p in: bemen˝o vektor – output: p out: a bemen˝o vektor megszorozva az adott szinthez tartoz´o kilencpontos prolong´aci´o m´atrix´aval 37

– mgcycle.m: n´eh´any multigrid iter´aci´os l´ep´est hajt v´egre az adott szinten – input: l: szintsz´am, u in: r´egi k¨ozel´ıt´es, f : jobboldali vektor, nu1: el˝osim´ıt´asok sz´ama, nu2: ut´osim´ıt´asok sz´ama, gamma: multigrid iter´aci´ok sz´ama – output: u out: u ´j k¨ozel´ıt´es, w: seg´edoutput – alprogramjai: laplace.m, rhs.m, pre smooth.m, post smooth.m, restri.m, prolong.m

3.2.

A teljes multigrid m´ odszer hat´ ekonys´ ag´ anak vizsg´ alata

Az mgteszt.m f´ajlban meg´ırt program a (3.1)-es ´es a (3.2) perem´ert´ek-feladatok pontos megold´asait hasonl´ıtja ¨ossze az fmgcycle.m f´ajl haszn´alat´aval kapott k¨ozel´ıt˝o megold´asokkal. A (3.1)-es feladat pontos megold´asa az u(x) = sin(πx) sin(πy), a (3.2)-es´e pedig az u(x) = sin(2πx) sin(3πy). A program a pontos ´es a multigrid m´odszer a´ltal szolg´altatott k¨ozel´ıt˝o megold´as k¨ozti elt´er´est a k´et f¨ uggv´eny k¨ ul¨onbs´eg´enek a´br´azol´as´aval szeml´elteti, egyre finomod´o r´acsok eset´eben. Erre az´ert van sz¨ uks´eg, mert a m´odszer olyan pontos, hogy a megfelel˝o f¨ uggv´enyek a´br´azol´asakor szabad szemmel nem is ´eszlelhet˝o az elt´er´es. Minden r´acsm´eret eset´en ki is sz´am´ıtjuk, hogy pontosan mennyi az adott k¨ozel´ıt´es relat´ıv hib´aja, vagyis az ku − u˜k/kuk ar´any (%-os alakban), ahol u a pontos megold´as, u˜ pedig a k¨ozel´ıt˝o megold´as vektora. Emellett a program ki´ırja, hogy az egyes k¨ozel´ıt´esek kisz´am´ıt´as´ahoz mennyi id˝ore volt sz¨ uks´ege. Input param´eterei ugyanazok, mint az fmgcycle.m alprogramj´anak. N´ezz¨ uk meg mit eredm´enyez az mgteszt() f¨ uggv´eny, ha input param´etereit p´eld´aul ´ıgy v´alasztjuk: 1. l = 9, azaz 10 r´acsot alkalmazunk 2. el˝osim´ıt´asok sz´ama: ν1 = 4 3. ut´osim´ıt´asok sz´ama: ν2 = 4 4. multigrid iter´aci´ok sz´ama az egyes szinteken: γ = 1 38

5. tesztf¨ uggv´eny sorsz´ama: 2 Az mgteszt(9,4,4,1,2) parancs be´ır´as´aval az al´abbi eredm´eny sz¨ uletik: Id˝ o

Relat´ıv hiba

1. szint

0.002239 sec

48.1111%

2. szint

0.006548 sec

10.1226%

3. szint

0.011695 sec

2.194%

4. szint

0.040206 sec

0.50671%

5. szint

0.147134 sec

0.1222%

6. szint

0.641772 sec

0.030129%

7. szint

2.879955 sec

0.0074966%

8. szint

12.128315 sec

0.0018713%

9. szint

51.297938 sec

0.00046762%

3.1. a´bra. A k¨ozel´ıt˝o megold´as hib´aja ´es a kisz´am´ıt´as´ahoz sz¨ uks´eges id˝o a szintek f¨ uggv´eny´eben A kapott eredm´enyekb˝ol azt a k¨ovetkeztet´est vonhatjuk le, hogy a multigrid m´odszer rendk´ıv¨ ul r¨ovid id˝o alatt nagyon pontos megold´ast szolg´altatott, mivel a h´al´ok finom´ıt´as´aval a k¨ozel´ıt´esek hib´aja igen gyorsan lecs¨okkent. A k¨ovetkez˝o k´epsoron ((3.2)-es ´es (3.3)-as ´abr´ak) megfigyelhetj¨ uk a k¨ozel´ıt˝o megold´as hib´aj´anak cs¨okken´es´et a szintek el˝orehaladt´aval (l=3, 4, 5, 6). Az a´br´akon a k¨ozel´ıt˝o ´es a pontos megold´as k¨ ul¨onbs´ege l´athat´o, valamint az adott k¨ozel´ıt´es hib´aja is fel van t¨ untetve. Az x ´es y tengelyek a r´acspontok szerint vannak felosztva, a z tengelyen a f¨ uggv´eny´ert´ekek minden a´br´an azonosan -0.02 ´es 0.02 k¨oz¨ott vannak a´br´azolva. Amint l´athattuk sim´ıt´o iter´aci´ok haszn´alat´ara az´ert van sz¨ uks´eg, mert a hiba magas frekvenci´aj´ u komponenseit – melyek jav´ıt´asa a durv´abb r´acsokon m´ar nem ´ lehets´eges – gyorsan lecs¨okkentik. Erdekes megfigyelni, hogyan viselkedik a m´odszer a sim´ıt´asok elv´egz´ese n´elk¨ ul, azaz hogyan torzul el a k¨ozel´ıt˝o megold´as, ha a magas sorsz´am´ u saj´atvektorokhoz tartoz´o megold´asi komponensek jav´ıt´asa nem t¨ort´enik meg a szintv´alt´asokkor. A sim´ıt´as n´elk¨ uli multigriddel nyert k¨ozel´ıt˝o megold´asokat – szembe´all´ıtva a pontos megold´asokkal – a (3.4)-es ´es (3.5)-¨os a´br´akon tekinthetj¨ uk meg. Amint l´athat´o, m´ar az els˝o tesztf¨ uggv´eny eset´en is jelent˝os az elt´er´es, de a m´asodik esetben a k¨ozel´ıt´esek m´eg csak nem is hasonl´ıtanak a pontos megold´asokra. 39

3.2. ´abra. A k¨ozel´ıt´esek hib´aj´anak cs¨okken´ese - 2. tesztf¨ uggv´eny (l=3, 4)

3.3. ´abra. A k¨ozel´ıt´esek hib´aj´anak cs¨okken´ese - 2. tesztf¨ uggv´eny (l=5, 6)

40

3.4. ´abra. Els˝o tesztfeladat - sim´ıt´o iter´aci´o n´elk¨ ul: mgteszt(6,0,0,1,1)

3.5. ´abra. M´asodik tesztfeladat - simit´o iter´aci´o n´elk¨ ul: mgteszt(6,0,0,1,2)

41

3.3.

A kaszk´ ad m´ odszer megval´ os´ıt´ asa

Ebben a r´eszben megvizsg´aljuk az u ´gynevezett kaszk´ad m´odszer hat´ekonys´ag´at, m´eghozz´a a j´oval bonyolultabb teljes multigrid m´odszerhez viszony´ıtva. A kaszk´ad algoritmus l´enyege a k¨ovetkez˝o: Ugyan´ ugy, mint a multigrid m´odszer eset´eben egy adott perem´ert´ek-probl´ema megold´as´ahoz diszkretiz´alnunk kell a feladatot l + 1 r´acson. Els˝o l´ep´esben megoldjuk a feladatot a legdurv´abb szinten. Az ´ıgy nyert megold´ast transzform´aljuk az eggyel magasabb szintre egy megfelel˝o prolong´aci´o seg´ıts´eg´evel. Ezut´an v´egrehajtunk ν db sim´ıt´o iter´aci´os l´ep´est (ν adott, ´altal´aban 5 ´es 10 k¨oz¨otti sz´am), mely lehet p´eld´al Jacobi- vagy Seidel-iter´aci´o. A kapott jav´ıtott k¨ozel´ıt´est szint´en emelj¨ uk az eggyel finomabb szintre, ´es sim´ıtsuk le a hib´aj´at. V´egezz¨ uk el ezeket a l´ep´eseket minden szinten eg´eszen addig, m´ıg eljutunk a legfinomabb r´acsig. ´Igy az eredeti feladat egy k¨ozel´ıt˝o megold´as´at nyerj¨ uk. (A m´odszerr˝ol r´eszletesebb anyag [8]-ban tal´alhat´o.) Az al´abbi programok seg´ıts´eg´evel meg fogjuk mutatni, hogy ez az egyszer˝ u m´odszer is k´epes igen j´o eredm´enyek el´er´es´ere. A kaszkad.m f´ajlban egy olyan kaszk´ad m´odszert defini´altuk, mely a multigrid program anal´ogj´ara k´esz¨ ult, szint´en a (3.1)-es ´es a (3.2)-es perem´ert´ek-feladatok megold´as´ara alkalmas, ugyanazokkal a r´acsokkal dolgozik ´es Gauss-Seidel iter´aci´ot v´egez ut´osim´ıt´ask´ent. Input param´eterei a k¨ovetkez˝ok: – l: szintsz´am – nu2: ut´osim´ıt´asok sz´ama – t: tesztf¨ uggv´eny sorsz´ama Az u output a megfelel˝o tesztfeladathoz tartoz´o k¨ozel´ıt˝o megold´as. A k´et m´odszer ¨osszehasonl´ıt´as´ara szolg´al´o, a kaszkadteszt.m f´ajlban tal´alhat´o program inputparam´eterei megegyeznek a t¨obbi tesztprogram´eval, azaz a multigrid ´es jelen esetben a kaszk´ad algoritmus speci´alis t´enyez˝oit a´ll´ıtj´ak be, vagyis azonos param´eterekkel b´ır´o elj´ar´asokat fogunk o¨sszevetni (ugyanannyi szint- ´es ut´osim´ıt´assz´am, ugyanaz a tesztfeladat). A program futtat´as´anak v´egezt´evel a megfelel˝o m´odszerek a´ltal szolg´altatott k¨ozel´ıt˝o megold´asok hib´ait ´es a kisz´am´ıt´asukhoz sz¨ uks´eges id˝oket tekinthetj¨ uk ´at a szintek sz´am´anak f¨ uggv´eny´eben. A kaszkadteszt(9,4,4,1,1) utas´ıt´as kiad´as´aval kapott eredm´enyeket a (3.6)-os ´abr´an ¨osszegezz¨ uk. L´athat´o, hogy a kaszk´ad m´odszer, egyszer˝ us´ege dac´ara, eg´esz j´o k¨ozel´ıt´eseket eredm´enyez, m´ıg a full multigriddel, b´ar pontosabb megold´asokra tett¨ unk szert, kisz´am´ıt´asuk m´ar j´oval t¨obb m˝ uveletbe ker¨ ult. A kaszk´ad k¨ozel´ıt´es hib´aja, a multigrid´evel ellent´etben, a szintek el˝orehaladt´aval n¨ovekszik, de ez a n¨oveked´es el´eg kism´ert´ek˝ u. 42

Multigrid m´ odszer

Kaszk´ ad m´ odszer

id˝ o

id˝ o

relat´ıv hiba

relat´ıv hiba

1. szint

0.002661 sec

5.1281%

0.001247 sec

3.8369 %

2. szint

0.004144 sec

1.0874%

0.001885 sec

1.5777%

3. szint

0.012083 sec

0.24749%

0.007339 sec

2.9006%

4. szint

0.040443 sec

0.059298%

0.027537 sec

3.2342%

5. szint

0.153938 sec

0.014586%

0.105074 sec

3.3178%

6. szint

0.666729 sec

0.0036261%

0.458630 sec

3.3387%

7. szint

2.941658 sec

0.0009049%

1.901906 sec

3.3439%

8. szint

12.645947 sec

0.0002261%

7.857439 sec

3.3452%

9. szint

52.440102 sec

0.000056515%

32.077767 sec

3.3455%

3.6. a´bra. A multigrid ´es a kaszk´ad m´odszer ´altal nyert k¨ozel´ıt˝o megold´asok hib´aja ´es a kisz´am´ıt´asukhoz sz¨ uks´eges id˝o a szintek f¨ uggv´eny´eben

3.4.

Az FMG ´ es a Seidel-iter´ aci´ o konvergenciasebess´ eg´ enek ¨ osszehasonl´ıt´ asa

A teszt.m ´es a teszt2.m f´ajlokban megtal´alhat´o programok a multigrid m´odszer hat´ekonys´ag´at h´ıvatottak igazolni, melyet a Gauss-Seidel iter´aci´o ´es a multigrid elj´ar´as konvergenciasebess´egeinek ¨osszehasonl´ıt´as´aval ´erz´ekeltet¨ unk. A programok input param´eterei megegyeznek az mgteszt.m f´ajlban meg´ırt program´eval, azaz most is a multigrid m´odszer speci´alis t´enyez˝oit hivatottak be´all´ıtani. N´ezz¨ uk el˝osz¨or a teszt.m programj´at, melynek seg´ıts´eg´evel azt fogjuk o¨sszevetni, hogy mennyi id˝ore van sz¨ uks´ege a multigrid m´odszernek ´es a Seidel-iter´aci´onak majdnem ugyanolyan pontoss´ag´ u k¨ozel´ıt´es el˝oa´ll´ıt´as´ahoz. Pontosabban mindig annyiszor futtatjuk a Seidel-iter´aci´ot, hogy legal´abb olyan j´o k¨ozel´ıt´est adjon, mint az adott szinthez tartoz´o multigrid k¨ozel´ıt´es. Tesztelj¨ uk ezt a f¨ uggv´enyt is a megszokott inputon (teszt(6,4,4,1,1) parancs). Mint azt a kapott eredm´enyek is mutatj´ak (lsd. (3.7)-es ´abra), a Gauss-Seidel iter´aci´onak nagys´ag-rendekkel t¨obb id˝ore van sz¨ uks´ege egy adott pontoss´ag´ u k¨ozel´ıt˝o megold´as el˝oa´ll´ıt´as´ahoz, mint a multigrid m´odszernek. A teszt2.m f´ajl programja a m´asik szemsz¨ogb˝ol hasonl´ıtja ¨ossze a Seidel-iter´aci´ot ´es a multigrid m´odszert, vagyis azt vizsg´alja, hogy k¨or¨ ulbel¨ ul azonos sz´am´ u m˝ uvelet 43

Multigrid m´ odszer id˝ o

relat´ıv hiba

Seidel-iter´ aci´ o id˝ o

relat´ıv hiba

1. szint

0.002375 sec

5.1281%

0.001651 sec

4.2227%

2. szint

0.006400 sec

1.0874%

0.016854 sec

1.0405%

3. szint

0.011440 sec

0.24749%

0.129532 sec

0.24609%

4. szint

0.039755 sec

0.059298%

2.404275 sec

0.058976%

5. szint

0.148012 sec

0.014586%

49.214101 sec

0.014585%

6. szint

0.648248 sec

0.0036261%

1191.548492 sec

0.0036243%

(≈ 20 perc) 3.7. ´abra.

A Seidel- ´es a multigrid m´odszer a´ltal nyert k¨ozel´ıt˝o megold´asok

l´enyeg´eben azonos hib´aja ´es a kisz´am´ıt´asukhoz sz¨ uks´eges id˝o a szintek f¨ uggv´eny´eben elv´egz´es´evel melyik elj´ar´as milyen relat´ıv hib´at produk´al. Amennyiben a Seidelben annyi iter´aci´ot hajtunk v´egre, amennyi a teljes multigridben a legfinomabb szinten van (el˝o- ´es ut´osim´ıt´asok o¨sszesen), a k´et m´odszer m˝ uveletsz´ama l´enyeg´eben megegyezik (ez annak k¨osz¨onhet˝o, hogy a multigridben az egyre durv´abb r´acsokon az iter´aci´ok sz¨ uks´eges m˝ uveletsz´ama sokkal kisebb, mint a Seidel-iter´aci´o´e). Ez´ert van sz¨ uks´eg az mgcycle.m ´es az fmgcycle.m f´ajlokban a w v´altoz´ora, mellyel o¨sszesz´aml´aljuk a multigridben az egyes szinteken elv´egzett iter´aci´okat, hogy megkapjuk az FMG-ben a legfinomabb szinten v´egrehajtott iter´aci´ok sz´am´at. Vizsg´aljuk meg mit kapunk eredm´eny¨ ul, ha a multigrid m´odszer param´etereit az el˝oz˝oekhez hasonl´oan v´alasztjuk (teszt2(6,4,4,1,1) parancs kiad´asa): Relat´ıv hiba Multigrid m´ odszer Seidel iter´ aci´ o 1. szint

5.1281%

7.1377%

2. szint

1.0874%

35.7734%

3. szint

0.24749%

68.8501%

4. szint

0.059298%

88.372%

5. szint

0.014586%

96.2163%

6. szint

0.0036261%

98.8499%

3.8. ´abra. A Seidel- ´es a multigrid m´odszer a´ltal nyert k¨ozel´ıt˝o megold´asok hib´aja nagyj´ab´ol azonos m˝ uvelet eset´en, a szintek f¨ uggv´eny´eben 44

Amint tapasztalhattuk, nagyj´ab´ol az a l´ep´essz´am, mellyel a multigrid m´odszer egy nagyon pontos k¨ozel´ıt´est ´er el, a Gauss-Seidel iter´aci´o sz´am´ara igen kev´esnek bizonyul, az eredm´eny¨ ul kapott Seidel-k¨ozel´ıt´es minden szinten sokkal rosszabb lesz a multigridesn´el. R´aad´asul a szintek el˝orehaladt´aval a Seidel-iter´aci´o egyre nagyobb ´es nagyobb relat´ıv hib´at eredm´enyez, ellenben a multigriddel, mely hib´aja cs¨okken˝o tendenci´at mutat. A k¨ovetkez˝o a´br´ak j´ol szeml´eltetik a a szintek el˝orehaladt´aval n¨ovekv˝o k¨ ul¨onbs´eget a k´et m´odszer a´ltal nyert k¨ozel´ıt˝o megold´asok k¨oz¨ott:

3.9. ´abra. Az FMG ´es a Seidel-iter´aci´o o¨sszehasonl´ıt´asa - 1. tesztf¨ uggv´eny (l=2)

45

3.10. ´abra. Az FMG ´es a Seidel-iter´aci´o o¨sszehasonl´ıt´asa - 1. tesztf¨ uggv´eny (l=3)

3.11. ´abra. Az FMG ´es a Seidel-iter´aci´o o¨sszehasonl´ıt´asa - 1. tesztf¨ uggv´eny (l=4)

46

3.12. ´abra. Az FMG ´es a Seidel-iter´aci´o o¨sszehasonl´ıt´asa - 1. tesztf¨ uggv´eny (l=5)

47

¨ Osszefoglal´ as Mint l´athattunk, a fizika sz´amos ter¨ ulet´en elker¨ ulhetetlen, hogy perem´ert´ekprobl´em´ak numerikus megold´as´ahoz kelljen folyamodnunk. Az el˝oz˝oekben egy olyan megold´o m´odszert ismerhett¨ unk meg, mely rendk´ıv¨ ul gyors konvergenciasebess´eg´enek k¨osz¨onhet˝oen emelkedik ki a hasonl´o numerikus megold´o algoritmusok k¨oz¨ ul. M˝ uveletig´enye a diszkr´et ismeretlensz´ammal ar´anyos (r´aad´asul a konstans szorz´o igen kicsi), vagyis a m´odszer nagys´agrend szerint optim´alis. Fontos el˝onye, hogy konvergenciasebess´ege a r´acsm´eret cs¨okkent´es´evel, vagyis a h´al´o finom´ıt´as´aval nem lassul le, ahogyan az a klasszikus iter´aci´okn´al megfigyelhet˝o (ezt mi is tapasztalhattuk a harmadik tesztprogram seg´ıts´eg´evel). Mindez persze nem azt jelenti, hogy lenne egy konkr´et multigrid algoritmusunk, melyet minden perem´ert´ek-probl´ema megold´asa sor´an alkalmazhatunk. Adva van egy multigrid technika, az algoritmus alapszerkezete, de minden egyes megoldand´o probl´ema eset´en a m´odszer speci´alis t´enyez˝oit a feladat f¨ uggv´eny´eben kell meghat´aroznunk. ´Igy egy adott probl´ema megold´as´ara szolg´al´o multigrid algoritmus hat´ekonys´aga a speci´alis komponenseinek megv´alaszt´as´an m´ ulik. Mindazon´altal a m´odszer igazi jelent˝os´ege nem a Poisson-egyenlet megold´as´aban felmutatott gyorsas´aga, hanem hogy hat´ekonys´aga az elliptikus perem´ert´ek-probl´em´akn´al j´oval bonyolultabb feladatok megold´asa sor´an is megmarad.

48

F¨ uggel´ ek Ebben a f¨ uggel´ekben olvashat´ok azok a programok, melyek bemutat´as´aval, elemz´es´evel a harmadik fejezet foglalkozott. Az egyes MATLAB parancsok jelent´es´e´ert lsd. [4].

laplace.m function [L]=laplace(l) % az l-edik szinthez tartoz´o diszkr´et Laplace-oper´ator el˝o´all´ıt´asa h=1/2ˆ (l+1); % r´acsm´eret n=(1/h)-1; N = n*n; % bels˝o r´acspontok sz´ama elemek=zeros(5*N, 1); sor=zeros(5*N, 1); oszlop=zeros(5*N, 1); index=0; sorszam=0; for j=1:n for i=1:n sorszam=sorszam+1; if j>1 index=index+1; elemek(index)=-1; sor(index)=sorszam; oszlop(index)=sorszam-n; end if i>1 index=index+1; elemek(index)=-1; 49

sor(index)=sorszam; oszlop(index)=sorszam-1; end index=index+1; elemek(index)=4; sor(index)=sorszam; oszlop(index)=sorszam; if i
rhs.m function [f] = rhs(l, t) % a t-edik tesztf¨ uggv´enyhez tartoz´o, l-edik szinten diszkretiz´alt egyenlet jobboldali vektora h=1/2ˆ (l+1); % r´acsm´eret n=(1/h)-1; 50

N=n*n; % bels˝o r´acspontok sz´ama [X,Y]=meshgrid([1:n]/(n+1), [1:n]/(n+1)); if t==1 B=2*(piˆ 2)*sin(pi*X).*sin(pi*Y); elseif t==2 B=13*(piˆ 2)*sin(pi*2*X).*sin(pi*3*Y); end f=reshape(B, N, 1); % jobboldali vektor end

pre smooth.m function [u]=pre smooth(l, u, f, nu) % nu db Jacobi-iter´aci´os l´ep´es az Ll ul = fl egyenlet u k¨ozel´ıt´es´enek jav´ıt´as´ahoz D=4/5*(1./spdiags(laplace(l), 0)); for i=1 : nu u=u+D.*(f-laplace(l)*u); end end

post smooth.m function [u] = post smooth(l, u, f, nu) % nu db Seidel iter´aci´os l´ep´es az Ll ul = fl egyenlet u k¨ozel´ıt´es´enek jav´ıt´as´ahoz L=tril(laplace(l)); for i=1:nu u=u+L\(f-laplace(l)*u); end end

restri.m function [r out] = restri(l, r in) % a finomabb, l-edik szintr˝ol a durv´abb (l-1)-edik szintre transzform´al´as h f=1/2ˆ (l+1); % r´acsm´eret a finomabb szinten n f=(1/h f)-1; h c=1/2ˆ (l); % r´acsm´eret a durv´abb szinten n c=(1/h c)-1; N c=n c*n c; % bels˝o r´acspontok sz´ama a durv´abb szinten 51

n0 f=n f+2; uls˝o r´acspontok sz´ama a finomabb szinten N0 f=n0 f*n0 f; % bels˝o + k¨ dx=1; dy=n0 f; % r in kiterjeszt´ese a peremre r0=zeros(N0 f, 1); for j=1:n f for i=1:n f r0(n0 f+1 + i + n0 f*(j-1))=r in(i + n f*(j-1)); end end % a durva r´acspontok kijel¨ol´es´ehez elk´esz´ıtj¨ uk az I vektort I=zeros(N c, 1); for j=1:n c for i=1:n c I(i + n c*(j-1))=2*i + 2*j*n0 f + 1; end end %full weighting restrikci´o r out=0.25*r0(I) + 0.125*(r0(I+dx) + r0(I-dx) + r0(I+dy) + r0(I-dy)) + . . . . . . + 0.0625*(r0(I+dx+dy) + r0(I-dx+dy) + r0(I+dx-dy) + r0(I-dx-dy)); end

prolong.m function [ p out ] = prolong( l, p in ) % a durv´abb, (l-1)-edik szintr˝ol a finomabb, l-edik szintre transzform´al´as h c=1/2ˆ (l); % r´acsm´eret a durv´abb szinten n c=(1/h c)-1; N c=n c*n c; % bels˝o r´acspontok sz´ama a durv´abb szinten h f=1/2ˆ (l+1); % r´acsm´eret a finomabb szinten n f=(1/h f)-1; N f=n f*n f; % bels˝o r´acspontok sz´ama a finomabb szinten sor=zeros(9*N c, 1); oszlop=zeros(9*N c, 1); 52

elemek=zeros(9*N c, 1); % a prolong´aci´os oper´ator elk´esz´ıt´ese v=[1:9]’; for j=1:n c for i=1:n c k=n f + (j-1)*2*n f + 2*i; sor(v)=[k-n f-1, k-n f, k-n f+1, k-1, k, k+1, k+n f-1, k+n f, k+n f+1]’; k=(j-1)*n c + i; oszlop(v)=k*ones(9, 1); elemek(v)=[1; 2; 1; 2; 4; 2; 1; 2; 1]; v=v+9; end end P=0.25*sparse(sor, oszlop, elemek, N f, N c); % a prolong´aci´os oper´ator p out=P*p in; % prolong´aci´o end

mgcycle.m function [u out, w] = mgcycle(l, u in, f, nu1, nu2, gamma) % a multigrid ciklus % input param´eterek: l: szintsz´am; u in: r´egi k¨ozel´ıt´es; f: jobboldali vektor; nu1: el˝osim´ıt´asok sz´ama; nu2: ut´osim´ıt´asok sz´ama; gamma: multigrid iter´aci´ok sz´ama a durv´abb szinteken % output param´eterek: u out: u ´j k¨ozel´ıt´es, w: elv´egzett iter´aci´ok sz´ama w=0; % w-ben o¨sszegzem az iter´aci´ok sz´am´at if l==1 L=laplace(l); u=pre smooth(l, u in, f, nu1); w=w+nu1; d=f-L*u; d=restri(l, d); v=laplace(l-1)\d; v=prolong(l, v); u out=post smooth(l, u+v, f, nu2); w=w+nu2; else L=laplace(l); u=pre smooth(l, u in, f, nu1); w=w+nu1; 53

d=f-L*u; d=restri(l, d); v=zeros(length(d), 1); for k=1 : gamma [v,ww]=mgcycle(l-1, v, d, nu1, nu2, gamma); w=w+ww; end v=prolong(l, v); u out=post smooth(l, u+v, f, nu2); w=w+nu2; end end

fmgcycle.m function [u,w] = fmgcycle(l, nu1, nu2, gamma, t) % a teljes multigrid m´odszer % input param´eterek: l: szintsz´am; nu1: el˝osim´ıt´asok sz´ama; nu2: ut´osim´ıt´asok sz´ama; gamma: multigrid iter´aci´ok sz´ama az egyes szinteken; t: tesztf¨ uggv´eny sorsz´ama % output param´eterek: u: fmg k¨ozel´ıt´es; w: a legmagasabb szinten megtett iter´aci´ok sz´ama u=laplace(0)\rhs(0, t); for k=1: l u=prolong(k, u); [u,w]=mgcycle(k, u, rhs(k, t), nu1, nu2, gamma); end end

mgteszt.m function [ ] = mgteszt( l , nu1 , nu2, gamma, t) % a teljes multigrid m´odszer hat´ekonys´ag´anak vizsg´alata disp(’A k¨ozel´ıt˝o megold´as hib´aja ´es a kisz´am´ıt´as´ahoz sz¨ uks´eges id˝o ’); for k=1: l h=1/2ˆ (k+1); % r´acsm´eret n=(1/h)-1; 54

N=n*n; % bels˝o r´acspontok sz´ama n0=n+2; N0=n0*n0; % bels˝o + k¨ uls˝o r´acspontok sz´ama disp([’a(z) ’ int2str(k) ’. szinten:’]); tic u=fmgcycle(k, nu1, nu2, gamma, t); % multigrides megold´as toc % a megold´as kiterjeszt´ese a peremre u0=zeros(N0, 1); for j=1:n for i=1:n u0(n0+1 + i + n0*(j-1))=u(i + n*(j-1)); end end % pontos megold´as [X,Y]=meshgrid([0:n+1]/(n+1), [0:n+1]/(n+1)); if t==1 U=sin(pi*X).*sin(pi*Y); elseif t==2 U=sin(2*pi*X).*sin(3*pi*Y); end % relat´ıv hiba hiba=[num2str(norm(reshape(U, N0, 1)-u0)/norm(reshape(U, N0, 1))*100) ’%’], % a k¨ozel´ıt˝o ´es a pontos megold´as k¨ ul¨onbs´eg´enek a´br´azol´asa mesh(reshape(u0, sqrt(length(u0)), sqrt(length(u0)))-U); axis([1, n0, 1, n0, -0.02, 0.02]); end end

kaszkadteszt.m function [ ] = kaszkadteszt( l, nu1, nu2, gamma, t ) % a kaszk´ad m´odszer bemutat´asa, ¨osszehasonl´ıt´asa a FMG-del

55

disp(’A k¨ozel´ıt˝o megold´asok hib´aja ´es a kisz´am´ıt´asukhoz sz¨ uks´eges id˝o: ’); for k=1: l h=1/2ˆ (k+1); % r´acsm´eret n=(1/h)-1; N=n*n; % bels˝o r´acspontok sz´ama n0=n+2; N0=n0*n0; % bels˝o + k¨ uls˝o r´acspontok sz´ama % pontos megold´as [X,Y]=meshgrid([0:n+1]/(n+1), [0:n+1]/(n+1)); if t==1 U=sin(pi*X).*sin(pi*Y); elseif t==2 U=sin(2*pi*X).*sin(3*pi*Y); end tic % multigrides megold´as u=fmgcycle(k, nu1, nu2, gamma, t); disp([’Multigrid - a(z) ’ int2str(k) ’. szinten:’]); toc % a megold´as kiterjeszt´ese a peremre u0=zeros(N0, 1); for j=1:n for i=1:n u0(n0+1 + i + n0*(j-1))=u(i + n*(j-1)); end end % multigrid m´odszer hib´aja hiba mg=[num2str(norm(reshape(U,N0,1)-u0)/norm(reshape(U,N0,1))*100) ’%’], tic % kaszkados megold´as v=kaszkad(k, nu2, t); disp([’Kaszk´ad - a(z) ’ int2str(k) ’. szinten:’]); toc % a megold´as kiterjeszt´ese a peremre

56

v0=zeros(N0, 1); for j=1:n for i=1:n v0(n0+1 + i + n0*(j-1))=v(i + n*(j-1)); end end % kaszk´ad m´odszer hib´aja hiba k=[num2str(norm(reshape(U,N0,1)-v0)/norm(reshape(U,N0,1))*100) ’%’], end

teszt.m function [ ] = teszt( l, nu1, nu2, gamma, t ) % a Seidel ´es a multigrid sz´am´ara sz¨ uks´eges id˝ok adott pontoss´ag´ u k¨ozel´ıt´es el´er´es´ehez disp(’A k¨ozel´ıt˝o mo.-k kb. azonos hib´aja ´es a kisz´am´ıt´asukhoz sz¨ uks´eges id˝o:’); for k=1: l h=1/2ˆ (k+1); % r´acsm´eret n=(1/h)-1; N=n*n; % bels˝o r´acspontok sz´ama n0=n+2; N0=n0*n0; % bels˝o + k¨ uls˝o r´acspontok sz´ama %pontos megold´as [X,Y]=meshgrid([1:n]/(n+1), [1:n]/(n+1)); if t==1 U=sin(pi*X).*sin(pi*Y); elseif t==2 U=sin(2*pi*X).*sin(3*pi*Y); end % pontos megold´as a hat´aron is megadva [X0,Y0]=meshgrid([0:n+1]/(n+1), [0:n+1]/(n+1)); if t==1 U0=sin(pi*X0).*sin(pi*Y0); elseif t==2 U0=sin(2*pi*X0).*sin(3*pi*Y0); 57

end tic u=fmgcycle(k, nu1, nu2, gamma, t); % multigrides megold´as disp([’Multigrid - a(z) ’ int2str(k) ’. szinten:’]); toc % a megold´as kiterjeszt´ese a peremre u0=zeros(N0, 1); for j=1:n for i=1:n u0(n0+1 + i + n0*(j-1))=u(i + n*(j-1)); end end % a multigrid m´odszer hib´aja hiba mg=[num2str(norm(reshape(U0,N0,1)-u0)/norm(reshape(U0,N0,1))*100) ’%’], v=zeros(length(rhs(k, t)), 1); D=4/5*(1./spdiags(laplace(k), 0)); %annyi Seidel-iter´aci´os l´ep´es, am´ıg el´er olyan k¨ozel´ıt´esi hib´at, mint a multigrid %az eltelt id˝o ker¨ ul m´er´esre ´es o¨szehasonl´ıt´asra tic while norm(reshape(U, N, 1)-u)<norm(reshape(U, N, 1)-v) v=v+D.*(rhs(k, t)-laplace(k)*v); end disp(’Seidel:’); toc % a megold´as kiterjeszt´ese a peremre v0=zeros(N0, 1); for j=1:n for i=1:n v0(n0+1 + i + n0*(j-1))=v(i + n*(j-1)); end end % Seidel-iter´aci´o hib´aja hiba S=[num2str(norm(reshape(U0,N0,1)-v0)/norm(reshape(U0,N0,1))*100) ’%’], end

58

teszt2.m function [ ] = teszt2( l, nu1, nu2, gamma, t ) % Seidel ´es multigrid hib´aj´anak ¨osszevet´ese, kb. azonos m˝ uvelet eset´en disp(’A k¨ozel´ıt˝o megold´asok hib´aja nagyj´ab´ol azonos m˝ uvelet eset´en: ’); for k=1: l h=1/2ˆ (k+1); % r´acsm´eret n=(1/h)-1; N=n*n; % bels˝o r´acspontok sz´ama n0=n+2; N0=n0*n0; % bels˝o + k¨ uls˝o r´acspontok sz´ama % pontos megold´as [X,Y]=meshgrid([0:n+1]/(n+1), [0:n+1]/(n+1)); if t==1 U=sin(pi*X).*sin(pi*Y); elseif t==2 U=sin(2*pi*X).*sin(3*pi*Y); end % multigrides megold´as % w: a legmagasabb szinten megtett iter´aci´ok sz´ama [u, w]=fmgcycle(k, nu1, nu2, gamma, t); % a megold´as kiterjeszt´ese a peremre u0=zeros(N0, 1); for j=1:n for i=1:n u0(n0+1 + i + n0*(j-1))=u(i + n*(j-1)); end end % annyi Seidel-iter´aci´os l´ep´es, amennyi a multigridben a legfinomabb szinten van v=zeros(length(rhs(k, t)), 1); D=4/5*(1./spdiags(laplace(k), 0)); for m=1:w v=v+D.*(rhs(k, t)-laplace(k)*v); end 59

% a megold´as kiterjeszt´ese a peremre v0=zeros(N0, 1); for j=1:n for i=1:n v0(n0+1 + i + n0*(j-1))=v(i + n*(j-1)); end end % a multigrid m´odszer hib´aja disp([’Multigrid - a(z) ’ int2str(k) ’. szinten:’]); hiba mg=[num2str(norm(reshape(U,N0,1)-u0)/norm(reshape(U,N0,1))*100) ’%’], % a Seidel-iter´aci´o hib´aja disp(’Seidel:’); hiba S=[num2str(norm(reshape(U,N0,1)-v0)/norm(reshape(U,N0,1))*100) ’%’], % a k¨ozel´ıt˝o megold´asok a´br´azol´asa subplot(1, 2, 1) mesh(reshape(u0, sqrt(length(u0)), sqrt(length(u0)))); if t==1 axis([1, n0, 1, n0, 0, 1]); elseif t==2 axis([1, n0, 1, n0, -1, 1]); end subplot(1, 2, 2) mesh(reshape(v0, sqrt(length(v0)), sqrt(length(v0)))); if t==1 axis([1, n0, 1, n0, 0, 1]); elseif t==2 axis([1, n0, 1, n0, -1, 1]); end end

60

Irodalomjegyz´ ek [1] A. Brandt, Multigrid Techniques: 1984 guide with applications to fluid dynamics. GMD-Studien, Nr. 85, Bonn, 1984. [2] Simon L., M´asodrend˝ u line´aris parci´alis differenci´alegyenletek, Tank¨onyvkiad´o, Budapest, 1983. [3] Stoyan G., Tak´o G., Numerikus M´odszerek III, Elte Informatikai Kar, 2008. [4] Stoyan G., Matlab, Typotex, 2005. [5] W. Hackbusch, U. Trottenberg (eds.), Multigrid Methods. Lecture Notes in Mathematics 960. Springer-Verlag, Berlin, 1982. [6] W. Hackbusch, Multi-Grid Methods and Applications. Springer-Verlag, Berlin, 1985. [7] Y. Saad, Iterative Methods for Sparse Linear Systems, Society for Industrial and Applied Mathematics, United States of America, 2003. [8] F. A. Bornemann, P. Deuflhard, The cascadic multigrid method for elliptic problems, Numerische Mathematik 75 (1996) 135-152. http://pcbornemann4.ma.tum.de/foswiki/pub/M3/Allgemeines/ FolkmarBornemannPublications/num75.pdf

61