Eötvös Loránd Tudományegyetem Természettudományi Kar
Takács Sára Matematika BSc, Matematikai elemz® szakirány
Peremérték-feladatok numerikus megoldása Szakdolgozat
Témavezet®: Kurics Tamás, egyetemi tanársegéd Alkalmazott Analízis és Számításmatematikai Tanszék
Budapest, 2009
Tartalomjegyzék 1. Bevezetés
3
1.1. Rövid összefoglaló a dolgozat témájáról . . . . . . . . . . . . .
3
1.2. Motiváció . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2. Peremérték-feladat megoldása egydimenzióban 2.1. Egyértelm¶ megoldás . . . . . . . . . . . . . . . . . . . . . . .
6 7
2.2. Létezik a megoldás . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3. Gyenge megoldások . . . . . . . . . . . . . . . . . . . . . . . . 12
3. Numerikus módszerek
14
3.1. Belövéses módszer . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2. Véges dierencia módszer . . . . . . . . . . . . . . . . . . . . 19 3.3. Végeselem módszer . . . . . . . . . . . . . . . . . . . . . . . . 30 3.4. Összefoglalás . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4. A függelék
37
4.1. M-mátrixok tulajdonságai . . . . . . . . . . . . . . . . . . . . 37 4.2. Funkcionálanalízis alapok
. . . . . . . . . . . . . . . . . . . . 38
4.3. Elméleti alapok a végeselem módszerhez . . . . . . . . . . . . 43
5. B függelék
44
5.1. Belövéses módszer algoritmus . . . . . . . . . . . . . . . . . . 44 5.2. Véges dierencia algoritmus . . . . . . . . . . . . . . . . . . . 46
2
1. fejezet Bevezetés 1.1. Rövid összefoglaló a dolgozat témájáról Jelen dolgozat a peremérték-feladatok numerikus módszereibe ad betekintést. Az els® részben a peremérték-feladatok megoldhatóságáról és a megoldások egyértelm¶ségér®l adok elméleti bevezetést, illetve levezetem a peremértékproblémák gyenge megoldásait. A második részben bemutatom a három numerikus módszert. Els®ként a az ún. belövéses módszert (Shooting method), aminek a lényege, hogy a peremérték-feladatot kezdetiérték-feladatra vezeti vissza, és a kezdeti deriváltat lényegében próbálgatással becsli. Azaz egy tetsz®leges becsült kiindulási értéket addig javít, amíg az eredményben a másik végpont tetsz®leges közelségbe kerül a megadotthoz. Második módszer a véges dierencia módszer. Ez a módszer egy sz¶kebb problémakörre ad közelít® megoldást. Adott egy peremérték feladat, olyan jobboldallal ami csak az id®t®l függ® változótól függ. Ebben az esetben a folytonos megoldást egy numerikus rácsfüggvénnyel közelítjük, amit pedig egy egyszer¶, bár kell®en nom közelítés esetén meglehet®sen nagy méret¶ lineáris egyenletrendszerrel számítunk. A harmadik módszer a végeselem módszer, melyr®l csak röviden lesz szó, elmélete magasszint¶ funkcionálanalízis ismereteket igényel. Mindhárom módszert egy-egy részletesen levezetett példával illusztrálom. 3
Az A függelékben a fenti módszerekhez f¶z®d® elméleti háttérismeretek mutatom be, a véges dierencia módszerhez tartozó M-mátrixok elméletét, illetve a végeselem módszerhez tartozó lineáris terek elméletét vázlatosan. A B függelékben pedig a dolgozatban bemutatott példák kiszámításához használt programok forráskódját ismertetem.
1.2. Motiváció Az alábbiakban bemutatok két, a gyakorlati életben el®forduló zikai példát, melyek megoldását másodrend¶ dierenciálegyenletek peremérték-problémájára tudunk visszavezetni, kiszámításukat pedig a kés®bbiekben tárgyalt numerikus megoldások megkönnyítik.
1. Célbalövési probléma A számegyenes 0 pontjában állunk, adott egy L távolság. Kérdés: milyen szögbe kell az ágyút beállítani, hogy a golyó pontosan L-ben csapódjon be? A modellben feltételezzük, hogy az x irányú sebesség állandó és az y irányú elmozdulás pedig csak a gravitációtól függ. Vagyis a légellenellást és egyéb tényez®ket nem vettünk gyelembe. Jelölje x = x(t),
y = y(t) az x illetve
y irányú elmozdulásokat, u pedig az x irányú sebességet. Ekkor x(t) ˙ =u
x(0) = 0
y¨(t) = −g y(0) = 0 A második sorra alkalmazzuk a láncszabályt1 :
−g = y¨(t) =
dy˙ dx dy˙ d d dy d2 y d y˙ = = u = u y˙ = u (u ) = u2 2 dt dx dt dx dx dx dx dx
Átrrendezve az alábbi közönséges dierenciálegyenletet kapjuk, homogén peremfeltétellel:
y 00 (x) = − ug2 y(0) = 0 y(L) = 0 1 d (y(x(t))) dt
=
dy dx dx dt
4
A pontos megoldás:
y(x) =
gx (L − x) 2u2
2. Gerenda meghajlása Az alábbi feladat az [5] könyvb®l származik. Adott egy, a két végén megtámasztott vízszintes gerenda, amire függ®leges irányú er® hat. Feladatunk meghatározni a meghajlott gerenda maximális behajlását, illetve adott maximális behajlás esetén a gerenda hosszát. Jelölje
P a legnagyobb terhelést, L a gerenda hosszát, J a gerenda metszetének tehetetlenségi nyomatékát, E a rugalmassági modulust, a pedig a gerenda vége és az er® támaszpontja közötti távolságot. Tegyük fel, hogy a P, L, J, a értékek mind ismertek az adott gerendára. Mindkét feladathoz szükséges meghatároznunk a gerenda neutrális szálának függvényét. Azt az u(x) függvényt, amely mentén a gerenda rostjai er®hatásra nem deformálódnak. A neutrális szál görbületi sugarára vonatkozó képlet:
EJ , M ahol M a hajlítónyomaték, M = −P a. A görbületi sugárra ismert képlet: R=
3
(1 + u0 2 ) 2 R= , u00 ahol u0 2 értékét elhanyagolhatjuk, mert a gerenda meghajlásánál az érint® du dx
iránytangense minden pontban nagyon kicsi. A neutrális szál dieren-
ciálegyenlete ekkor az alábbi alakú, a végpontokban pedig nincs elmozdulás, tehát homogén peremfeltétellel: Pa u00 (x) = − EJ = állandó
u(0) = 0,
u(L) = 0
Ha nem hanyagoljuk el az elhajlás érint® iránytengensét, akkor pedig az 3
Pa u00 (x) = f (x, u0 (x), u(x)) = − EJ (1 + u0 2 ) 2
u(0) = 0,
u(L) = 0
peremérték-feladatot kapjuk, amit numerikus módszerekkel meg tudunk oldani. 5
2. fejezet Peremérték-feladat megoldása egydimenzióban Jelen fejezet a [2] könyv második fejezete alapján készült. Megmutatom, hogy egy másodrend¶ közönséges dierenciálegyenleten alapuló peremértékfeladatnak milyen feltételek mellett létezik megoldása, illetve egyértelm¶-e ez a megoldás. Egy inhomogén másodrend¶ lineáris dierenciálegyenletet Dirichletperemfeltétellel:
(
Au := −(au0 )0 + bu0 + cu = f u(0) = u0
Ω = (0, 1)-en
u(1) = u1
(2.1)
ahol az a, b, c együtthatók a(x), b(x), c(x) sima függvények, kell®en sokszor folytonosan dierenciálhatóak, valamint
¯ = [0, 1]-en. a(x) ≥ a0 > 0 és c(x) ≥ 0 ∀x ∈ Ω
(2.2)
Tekintsük el®sz®r is azt a speciális esetet, amikor a ≡ 1 és b, c ≡ 0. Ekkor a kiindulási feladat
(
−u00 (x) = f (x) x ∈ Ω-n u(0) = u0 6
u(1) = u1
(2.3)
feladatra redukálódik. Ennek a feladatnak egyértelm¶ a megoldása, ui. kétszer integrálva megkapjuk u(x)-et: Z u(x) = −
x
Z
y
f (s)dsdy + αx + β 0
0
ahol β = u(0) = u0 , α-t pedig Ω másik végpontjában megadott érték felhasználásával számoljuk a következ® integrállal: Z 1Z y f (s)dsdy + α + u0 = u(1) = u1 0
0
Ha speciálisan f ≡ 0 is teljesül, akkor ( −u00 (x) = 0 x ∈ Ω-n
u(0) = u0
u(1) = u1
(2.4)
a (2.4) alakban felírt egyenlet megoldása lineáris, azaz
u(x) = (u1 − u0 )x + u0 = u0 (1 − x) + u1 x, vagyis u(0) ≤ u(x) ≤ u(1), a függvény minden bels® pontban az intervallum végpontjaiban felvett függvényértékek között van, tehát u maximumát és a minimumát az intervallum végpontjaiban veszi fel. Ezt általánosítja a következ® tétel.
2.1. Egyértelm¶ megoldás ¯ 1. Tétel (Maximum-elv). Au := −(au0 )0 +bu0 +cu = f Legyen u ∈ C 2 (Ω) és Au ≤ 0 Ω-n. Ekkor
© ª 1. ha c = 0 =⇒ max u = max u(0), u(1) ¯ Ω
© ª 2. ha c ≥ 0 =⇒ max u ≤ max u(0), u(1), 0 ¯ Ω
Ha Au ≥ 0, akkor a függvény minimumára állnak fenn hasonló állítások.
7
Bizonyítás: 1. (a) El®ször vegyünk egy szigorúbb feltételt, mint amit a tétel eredetileg kimond tegyük fel, hogy Au < 0
Au ≤ 0 helyett. Ha u-nak
van maximuma egy x0 ∈ Ω bels® pontban, akkor u0 (x0 ) = 0 és
u00 (x0 ) ≤ 0. Au := −(au0 )0 + bu0 + cu = −a0 u0 − au00 + bu0 + 0u Helyettesítsük be az egyenletbe x0 maximumhelyet:
0 > Au(x0 ) = −a0 (x0 )u0 (x0 ) − a(x0 )u00 (x0 ) + bu0 (x0 ) = ahol az els® és harmadik tag kiesik, mert u00 (x0 ) = 0, így:
= −a(x0 )u00 (x0 ) ≥ 0, mert a(x0 ) ≥ a0 szigorúan pozitív számnál és u00 (x0 ) ≤ 0, mert x0 maximum, vagyis szorzatuk nemnegatív, ami ellentmondás. (b) Tegyük fel, hogy Au ≤ 0. Legyen ϕ olyan függvény, hogy ϕ ≥ 0
Ω-on és Aϕ < 0 Ω-n. (Például ϕ = eλx jó, ha λ elég nagy.) Aϕ = (−aλ2 + (b − a0 )λ)ϕ < 0, ha λ elég nagy. Indirekten tegyük fel, hogy u függvény felveszi a maximumát egy x0 bels® pontban, de a végpontokban nem. Ekkor v = u + εϕ is ilyen függvény, ha ε elég kicsi. Av = Au + εAϕ < 0
∀x ∈ Ω-ra, ami ellentmondás.
2. c ≥ 0 Au ≤ 0. Kell: maxΩ u ≤ max{u(0), u(1), 0}. Ha u ≤ 0 Ω-n, akkor készen vagyunk. Ellenkez® esetben tegyük fel, hogy
maxΩ u = u(x0 ) > 0. De ez ellentmondás ugyanis x0 6= 0 és x0 6= 1. Vegyük Ω-ának egy (α, β) maximális részintervallumát, melyre x0 ∈ (α, β) ˜ := Au − cu < 0 (α, β)-n. Az 1. pontés u > 0. Ekkor legyen Au ban beláttottak alapján u(x0 ) = max{u(α), u(β)} valamint α és β nem lehetnek mindketten Ω bels® pontjai, így ellentmondásba kerültünk az
(α, β) intervallum megválasztásával. Tehát u(x0 ) = max{u(0), u(1)}. ¤ 8
¯ , akkor ∃C ≥ 0, melyre 2. Tétel. Ha u ∈ C 2 (Ω) kukC(Ω) ¯ ≤ max{|u(0)|, |u(1)|} + CkAukC(Ω) ¯ . A 2. tétel bizonyítása a [2] könyvben megtalálható.
1. Következmény. A peremértékprobléma megoldása egyértelm¶. Bizonyítás: Indirekten tegyük fel, hogy u és v is megoldásai
Au = f
u(0) = u0
u(1) = u1
feladatnak. Legyen w = u − v . Ekkor
Ω = (0, 1)-en és
Aw = Au − Av = f − f = 0 w(0) = 0 w(1) = 0.
Az el®z® tétel alapján: kwkC(Ω) ¯ ≤ max{|w(0)|, |w(1)|} + CkAwkC(Ω) ¯ . Vagyis ¯ = [0, 1]-en, tehát u = v . w(x) ≡ 0 Ω
2. Következmény. Ha (
Au = f
Ω-n
u(0) = u0 akkor ku − vkC(Ω) ¯
( és
Av = g
Ω-n
u(1) = u1 v(0) = v0 v(1) = v1 © ª ≤ max |u0 − v0 |, |u1 − v1 | + Ckf − gkC(Ω) ¯ , azaz a feladat
stabil.
3. Következmény. Au = f Ω-n u(0) = u0 u(1) = u1 és f ≤ 0 u0 , u1 ≤ 0, akkor u ≤ 0. Bizonyítás: Ha feltesszük, hogy Au = f ≤ 0, és alkalmazzuk a maxi© ª mum elvet, mely szerint maxΩ¯ u = max u(0), u(1), 0 = 0, tehát u legfeljebb 0 lehet.
4. Következmény (Monotonitás). (
Au = f u(0) = u0
Ω-n u(1) = u1
és
(
Av = g v(0) = v0
Ω-n v(1) = v1
és f ≤ g, u0 ≤ v0 , u1 ≤ v1 . Ekkor u ≤ v . Bizonyítás: Legyen w = u−v ismét a különbségfüggvény. A feltételek szerint Aw = f − g ≤ 0. w(0) = u0 − v0 ≤ 0 w(1) = u1 − v1 ≤ 0. Alkalmazva az el®z® következményt megkapjuk az állítást. 9
2.2. Létezik a megoldás Konkrét megoldás el®állítása a b ≡ 0 speciális esetben. Mostantól
Au := −(au0 )0 + cu és keressük az Au = f
u(0) = u(1) = 0
(2.5)
u(0) = 0 = u(1) feladat megoldását.
Legyenek U0 és U1 megoldásai az
AU0 = 0 Ω-n
U0 (0) = 1, U0 (1) = 0
AU1 = 0 Ω-n
U1 (0) = 0, U1 (1) = 1
feladatoknak. Ilyen U0 , U1 függvények léteznek, tekintsük ugyanis az
Au = 0 u(0) = 0, u0 (0) = 1 kezdetiérték feladatot, aminek tudjuk, hogy létezik megoldása. És u(1) 6= 0, mert ha u(1) = 0 lenne Au = 0 u(0) = 0 mellett, akkor u ≡ 0 függvény lenne, de ez ellentmond annak, hogy u0 (0) = 1. Elérhet® (például alkalmas konstanssal megszorozva), hogy Au = 0 u(0) = 0 u(1) = 1 legyen, ez lesz U1 .
U0 -t is hasonlóan állíthatjuk el®.
3. Tétel. Legyen b=0. A (2.5)-beli alakú Au = f, u(0) = 0, u(1) = 0 Ω-n feladat megoldása:
Z
1
u(x) = G(x, y) =
0 1 U (x)U1 (y) κ 0
ha
G(x, y)f (y)dy, ahol 0≤y≤x≤1 Green-függvény, és ahol
1 U (x)U0 (y) κ 1
ha
0≤x≤y≤1
κ = a(x)(U0 (x)U10 (x) − U00 (x)U1 (x)) ≡ konstans. Bizonyítás: El®sz®r belátjuk, hogy κ 6= 0 konstans. Felhasználjuk, hogy
AU0 = −(aU00 )0 + cU0 = 0, 10
amib®l következik, hogy
(aU00 )0 = cU0 és (aU10 )0 = cU1 . Felírjuk κ deriváltját, és ha ez 0 lesz akkor κ konstans:
κ0 (x) = U0 (aU10 )0 + U00 aU10 − U1 (aU00 )0 − U10 aU00 = 0 Behelyettesítve x = 0-t
¢ ¡ κ(0) = a(0) U10 (0) − 0 6= 0, mert a(0) ≥ a0 pozitív számnál, és U10 (0) sem lehet nulla, ugyanis ellenkez® esetben U1 (0) = 0 lenne, így U1 azonosan nulla függyvény lenne, ami ellentmond U1 (1) = 1-nek. AU1 = 0
U1 (0) = 0,
U1 (1) = 1 egyenletek-
ben mindhárom egyenl®ség helyettesíhet® "≥"-vel, így alkalmazhatjuk a 3. következményt, miszerint U1 ≥ 0, így U10 (0) > 0, ami ekvivalens azzal, hogy
κ > 0. Másodszor belátjuk, hogy a fent de niált u függvény kielégíti a homogén peremfeltételeket, (azaz u(0) = 0 és u(1) = 0-t): Z 1 1 u(0) = G(0, y)f (y)dy = U1 (0)U0 (y)f (y)dy = 0 κ 0 Z 1 1 u(1) = G(1, y)f (y)dy = U0 (1)U1 (y)f (y)dy = 0 κ 0 u(x)-be behelyettesítve a homogén peremfeltételeket Z x Z 1 u(x) = G(x, y)f (y)dy + G(x, y)f (y)dy = 0
x
= κ1 U0 (x) deriváljuk
Rx 0
µ 0
u (x) =
U1 (y)f (y)dy + κ1 U1 (x)
Z
x
1 κ
U00 (x)
+ κ1
µ Z 0 U1 (x)
0
R1 x
U0 (y)f (y)dy
¶ U1 (y)f (y)dy + U0 (x)U1 (x)f (x) + 1
x
¶ U0 (y)f (y)dy − U1 (x)U0 (x)f (x) 11
Kell még, hogy (au0 )0 = cu
µ 0 0
−(au ) (x) =
− κ1
¶0 Z a(x)U00 (x)
µ − κ1
x 0
¶0 Z a(x)U10 (x)
1 x
U1 (y)f (y)dy −
U0 (y)f (y)dy−
− κ1 a(x)(U0 (x)U10 (x) − U00 (x)U1 (x))f (x) = Z =
− κ1 c(x)U0 (x)
x 0
Z − κ1 c(x)U1 (x) Z
1
= −c(x) 0
1 x
U1 (y)f (y)dy −
U0 (y)f (y)dy + f (x) =
G(x, y)f (y)dy + f (x)
Azaz (−au0 )0 + cu = f , ezzel a tételt bebizonyítottuk. ¤
2.3. Gyenge megoldások A fentiekben láttuk, hogy egyértelm¶ megoldása létezik egy peremértékfeladatnak. Ebben a részben az egyváltozós peremérték-problémát az L2 =
L2 (Ω) Hilbert-térben tárgyaljuk, és levezetjük az ún. gyenge megoldásokat, ami a következ® fejezetben tárgyalt végeselem módszer hez lesz fontos. A felhasznált funkcionálanalízisbeli eszközöket az A függelékben ismertetem. Tekintsük az alábbi peremérték-problémát homogén peremfeltétellel: ( Au := −(au0 )0 + bu0 + cu = f Ω = (0, 1)-en
u(0) = u(1) = 0
(2.6)
Feltesszük, hogy az a, b, c együtthatók kell®en sima függvények és a (2.2) feltételek helyett az
a(x) ≥ a0 > 0
és
c(x) − 12
b0 (x) ≥0 2
¯ -n. feltételek teljesülnek ∀x ∈ Ω Szorozzuk be a dierenciálegyenletet egy ϕ ∈ C01 = C01 (Ω) un. tesztfügg-
vénnyel. C01 (Ω) a folytonosan dierenciálható és az értelmezési tartomány szélén 0 értéket felvev® függvények tere. Integráljuk a (2.6)-beli dierenciálegyenletet az Ω intervallumon Z 1 Z 0 0 0 (−(au ) + bu + cu)ϕdx = 0
1 0
f ϕ dx
Parciálisan integráljunk és használjuk ki, hogy ϕ(0) = ϕ(1) = 0:
Z
1 0
Z 0
0
1
0
(au ϕ + bu ϕ + cuϕ)dx =
0
f ϕ dx
∀ϕ ∈ C01
(2.7)
Vezessük be az alábbi bilineáris formát: Z 1 a(v, w) = (av 0 w0 + bv 0 w + cvw)dx 0
és L(w) lineáris funkcionált:
Z
1
L(w) = (f, w) = 0
f w dx
Kihasználva, hogy C01 s¶r¶ H01 = H01 (Ω)-n a (2.7)-t felírhatjuk a következ®képpen:
a(u, ϕ) = L(ϕ),
∀ϕ ∈ H01
(2.8)
A (2.6) feladat u megoldását gyenge megoldás nak nevezzük, ha u ∈ H01 és a (2.8) fennáll.
13
3. fejezet Numerikus módszerek A fejezetben tárgyalt három numerikus módszert [1] mutatja be.
3.1. Belövéses módszer Az
(
u00 = f (x, u, u0 ) u(a) = α,
a≤x≤b u(b) = β
(3.1)
peremértékproblémát szeretnénk megoldani, úgy hogy azt a következ® kezdetiérték feladattá alakítjuk
u00 = f (x, u, u0 ) u(a) = α
(3.2)
u0 (a) = s ahol s ∈ R tetsz®leges paraméter. Azaz megpróbáljuk s-et úgy megválasztani, hogy s az ismeretlen függvény a-beli érint®ének meredeksége legyen. Ez általában nem szokott sikerülni, így a következ® algoritmus segítségével közeledünk a megfelel® s-hez. Mivel a kiinduló feladat peremérték probléma, olyan s-et keresünk, amire a (3.2) u(x, s) megoldása teljesíti u(b, s) = β -t, azaz s paraméterrel is a megoldásfüggvény b-ben β kell legyen. Ezért de niálom a következ® F : R → R függvényt:
F (s) := u(b, s) − β 14
Azt az s∗ -ot keressük, amire
F (s∗ ) = 0. A kereséshez egy gyökkeres® algoritmust kell segítségül hívni. A belövéses módszer a Newton-módszert használja, ami a függvény deriváltja segítségével közelíti a gyököt. (sn → s∗ ):
sn+1 = sn −
F (sn ) . F 0 (sn )
Azonban problémát jelent F deriváltjának meghatározása. F 0 kiszámításához de niáljuk az alábbi v segédfüggvényt:
v :=
∂u ∂s
és
u00 (x, s) = f (x, u(x, s), u0 (x, s)) s-szerinti deriváltját írjuk fel a v függvény segítségével: ∂ 00 u (x, s) ∂s
= ∂1 f (x, u(x, s), u0 (x, s)) · 0+ +∂2 f (x, u(x, s), u0 (x, s)) · v+ +∂3 f (x, u(x, s), u0 (x, s)) ·
∂u0 (x,s) . ∂s
Ahhoz, hogy a Newton-módszerbeli F 0 (sn ) meghatározható legyen, ki kell számolnunk az alábbi segédfüggvény által meghatározott kezdetiérték feladatot: ( v 00 (x, s) = fu (x, u(x, s), u0 (x, s))v(x, s) + fu0 (x, u(x, s), u0 (x, s))v 0 (x, s)
v(a, s) = 0,
v 0 (a, s) = 1
Tehát F 0 (s) = v(b, s), ahol v a fenti feladat megoldása.
15
Algoritmus 1. s ∈ R választása (kezdeti meredekség találomra) 2. Az u00 = f (x, u, u0 ), u(a) = α, u0 (a) = s kezdetiérték feladatot oldjuk meg numerikusan 3. Oldjuk meg a
v 00 = ∂2 f (x, u, u0 ) · v + ∂3 f (x, u, u0 ) · v 0 v(a) = 0 v 0 (a) = 1 feladatot is. 4. Ha |u(b) − β| < ² → STOP Ha |u(b) − β| ≥ ², akkor sn+1 := sn −
u(b)−β v(b)
és folytassuk a 2. lépést®l.
1. Példa. Vegyük az alábbi peremértékproblémát: √ 1√ u(1) = 2, u(2) = 2, 2 √ pontos megoldása az u(x) = 2/x. Numerikusan is megoldjuk a fenti algou00 = u3 ,
ritmus használatával. El®ször vegyük az
u00 = u3 ,
u(1) =
√
2,
u0 (1) = s,
kezdeti érték feladatot. Oldjuk meg numerikusan az Euler-módszerrel. Ehhez a másodrend¶ kezdetiérték-feladatot els®rend¶ dierenciaegyenlet rendszerré alakítjuk az alábbiak szerint: !0 Ã ! Ã u w , = w u3
Ã
u w
!Ã
1 1
! =
à √ ! 2 s
A Newton-iteráció elkezdéséhez válasszuk s-et 0-nak. (A pontos kezdeti s √ érték s = − 2 = −1.414214 lenne.) Az F 0 kiszámolásához szükséges
v 00 = ∂u f · v + ∂u0 f · v 0 = 3u2 v 16
másodrend¶ dierenciélegyenletet is els®rend¶ rendszerré kell alakítani:
Ã
v
!0
à =
z
!
z
à ,
3u2 v
v
!Ã
z
1 1
!
à =
0
!
1
A fenti dierenciálegyenletrendszert ugyancsak az Euler-módszerrel oldjuk meg. Így tovább tudunk lépni az s∗ -ot, mint gyököt keres® ciklusban. A program akkor áll meg, ha F (s) kisebb lesz, mint egy el®re beállított ², jelen esetben ² = 10−6 , vagy ha elér bizonyos lépésszámot. Példánkban már hat lépés elég volt, hogy F (s) értéke ²-nál kisebb legyen. A numerikus eredményeket az alábbi táblázatok mutatják, lefuttatva három különböz® lépéshosszra. (A program forráskódja a B függelékben található.)
h=0.1
h=0.01
h=0.001
s
F(s)
s
F(s)
s
F(s)
0
2.5018
0
3.6034
0
3.8177
-1.0873
0.4810
-0.7975
1.1191
-0.7513
1.2493
-1.4054
0.0226
-1.3107
0.1636
-1.2851
0.2065
-1.4218
0.0037
-1.4126
0.0039
-1.4101
0.0065
-1.4218
0.0001
-1.4151
0.0000
-1.4143
0.0000
-1.4218
0.0000
-1.4151
0.0000
-1.4143
0.0000
A megoldásokat ábrázoltam is az alábbi ábrán. A folytonos vonal a pontos megoldást jelöli, a szaggatottak pedig a belövéses módszer algoritmusának megoldásfüggvényei. Mint látható, a megoldások egyre közelítik a pontos függvénygörbét, az utolsó három lépésben már nincs is szemmel látható eltérés. 17
Belövéses módszer 4.5 4 s(0)
3.5 3 2.5 2
s(1) 1.5 s(2)
1
s(3)−s(5) 0.5
1
1.2
1.4
1.6
1.8
2
Többszörös belövéses módszer Rosszul kondícionáltságból fakadó numerikus problémák azt jelentik, hogy
s nagyon kicsi megváltoztatása nagyon nagy megváltozást idéz el® a végpontbeli megoldásban. Rosszul kondícionált feladatot az egyszer¶ belövéses módszerrel nem tudunk megoldani. Az ilyen esetek orvoslására való a több-
szörös belövéses módszer. Az [a,b] intervallumot felosztjuk n db részintervallumra
a = x0 < x 1 < . . . < x n = b u és s Rn -beli vektorok, u := (u0 , . . . , un−1 )T , 18
s = (s0 , . . . , sn )T .
Nézzük az alábbi kezdetiérték problémát 00 0 u = f (x, u, u ) [xj , xj+1 ]-en,
u(xj ) = uj u0 (x ) = s j j minden j = 0, . . . , n − 1-re.
x1 , . . . , xn−1 -ben az u(·, uj , sj ) megoldás egybeesik az els® deriváltjával (a függvény kétszer folytonosan dierenciálható [a,b]-n és u(b) = β is teljesül). Kapunk 2n-1 db nemlineáris egyenletet 2n-1 db ismeretlennel:
u1 , . . . , un−1 , s0 . . . , sn−1 . A Newton-módszer itt is alkalmazható.
3.2. Véges dierencia módszer A peremérték-problémák közelítését megvalósító véges dierencia módszerek lényege, hogy a dierenciálegyenletbeli deriváltakat dierenciákkal helyettesítjük. Az egyszer¶ség kedvéért csak lineáris dierenciálegyenleteket vizsgálunk. Tekintsük az alábbi másodrend¶ lineáris dierenciálegyenletet, inhomogén peremfeltétellel Ω = (0, 1)-en: ( Au := −au00 + bu0 + cu = f
u(0) = u0 ,
u(1) = u1
a = a(x) > 0,
b = b(x),
Ω = (0, 1)
(3.3)
ahol
c = c(x) ≥ 0
(3.4)
¯ -n értelmezett folytonos függvények, és tegyük fel, hogy u ∈ C 4 (Ω) ¯ . Ω
Diszkretizáció Osszuk fel a (0, 1) intervallumot M egyenl® részre. x0 , x1 , . . . , xM darab rácspont,
x0 = 0,
xM = 1,
xj = jh, 19
h=
1 , M
j = 0, . . . , M.
M +1
A bels® osztópontokban a második deriváltat az alábbi dierenciával közelítjük:
¯ j≈ ∂ ∂U
Uj+1 − 2Uj + Uj−1 h2
j = 1, . . . , M − 1.
(3.5)
Ha pedig b nem azonosan nulla függvény, akkor az els® deriváltakat a bels® osztópontokban az alábbi három dierencia közül az egyikkel becsülhetjük: 1. haladó dierencia: ∂Uj ≈
Uj+1 −Uj h
¯ j≈ 2. retrográd dierencia: ∂U ˆ j≈ 3. centrális dierencia: ∂U
Uj −Uj−1 h
Uj+1 −Uj−1 2h
Jelen dolgozatban a harmadikat, a centrális dierenciát használom az els® derivált közelítésére. Új jelöléseket vezetünk be:
aj := a(xj ),
bj := b(xj ),
cj := c(xj ),
fj := f (xj )
csak a rácspontokban értelmezett függvények. Ekkor a (3.3) helyett az alábbi formulával dolgozunk: ( ¯ j + bj ∂U ˆ j + c j Uj = fj ; Ah Uj = −aj ∂ ∂U
U 0 = u0 ,
U 1 = u1
Keressük az
U = (U0 , U1 , . . . , UM −1 , UM ) vektort, ahol U1 , . . . , UM −1 ismeretlenek, U0 és UM pedig a peremértékek. Ha xj egy bels® osztópont, akkor a rá vonatkozó egyenlet:
−aj
Uj+1 − Uj−1 Uj+1 − 2Uj + Uj−1 + b + c j Uj = fj j h2 2h
h2 -tel megszorozva az egyenlet mindkét oldalát és elvégezve az egyszer¶sítéseket (2aj + h2 cj )Uj − (aj −
hbj hbj )Uj+1 − (aj + )Uj−1 = h2 fj 2 2 20
alakú egyenleteket kapunk
j = 1, 2, . . . , M −1-re, ami az alábbi lineáris egyenl®ségrendszerre vezet: (3.6)
AU = g,
ahol keresend® az U vektor; A (M − 1) × (M − 1)-es mátrix, U és g M −1 dimenziós vektorok, ahol g annyiban különbözik f -t®l, hogy dimenziója kett®vel kisebb, azaz f els® és utolsó elemét nem tartalmazza, illetve g els® és utolsó koordinátáit úgy módosítottuk, hogy a peremfeltételek is szerepeljenek az egyenletrendszerben. Tehát M + 1 darab osztópont esetén a lineáris egyenletrendszer M − 1 dimenziós. Az A mátrix általános alakjában a j-edik sor:
..
..
.
..
.
A·,j = . . . 0 −aj −
hbj 2
.
2aj + h2 cj −aj + .. .. . .
0... .. .
hbj 2
Részletesebben nézzük az alábbi példát. Legyen M=5, vagyis a (0,1) intervallumot öt egyenl® részre osszuk fel. Így a végpontokat is beleértve hat osztópontunk van, a lépésköz h = 1/5. Keressük az U = (U0 , U1 , U2 , U3 , U4 , U5 ) vektort ahol U0 = u0 és pedig U5 = u1 . Az egyenletrendszerünk tehát négy soros lesz:
(2a1 + (2a + 2 (2a3 + (2a4 +
c1 )U1 9 c2 )U2 9
− (a1 −
c3 )U3 9 c4 )U4 9
− (a3 −
Mátrixos alakja pedig: 2a1 + c91 −a1 + b61 −a2 − b62 2a2 + c92 0 −a3 − b63
0
0
− (a2 − − (a4 −
b1 )U2 6 b2 )U3 6
− (a1 +
b3 )U4 6 b4 )U5 6
− (a3 +
− (a2 + − (a4 +
0
0
b2 6 c3 2a3 + 9 −a4 − b64
0
−a2 +
21
−a3 + 2a4 +
b1 )U0 6 b2 )U1 6
=
b3 )U2 6 b4 )U3 6
=
b3 6 c4 9
U1
= =
f1 9 f2 9 f3 9 f4 9
g1
U2 g2 = U3 g3 g4 U4
,
ahol a jobboldalon g 6= h2 f , de csak az els® és utolsó komponenseiben különbözik t®le, melyek a ki vannak egészítve peremfeltételekkel. Vagyis
g=
f1 9 f2 9 f3 9 f4 9
+ (a1 +
b1 )u0 6
+ (a4 −
b4 )u1 6
A B függelékben bemutatok egy a fenti azámításokkal dolgozó példát Matlab forráskóddal. Az alábbiakban a következ® kérdések megválaszolására kerül sor: 1. Egyértelm¶en megoldható-e a (3.6)-beli egyenletrendszer? 2. Mekkora a hiba az uj approximált megoldás és az u(xj ) pontos megoldások között? 3. Konvergál-e a közelít® megoldás a pontos megoldáshoz, ha h → 0?
4. Tétel. Ha h elég kicsi, akkor a (3.6) egyenletrendszernek egyértelm¶ megoldása van. Bizonyítás: Az A mátrix tridiagonális, mint láttuk. Diagonálisan domináns, ha h elég kicsi, ugyanis:
¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ − aj − hbj ¯ + ¯ − aj + hbj ¯ ≤ 2aj + h2 cj > 0 ¯ ¯ ¯ 2 2 ¯
¤
A következ® állításban belátjuk, hogyha a fentiek és még egyéb feltételek amik az A mátrixra szintén fennállnak teljesülnek, akkor A invertálható, azaz a (3.6) lineáris egyenletrendszer megoldható.
1. Állítás. Ha A mátrix tridiagonális, diagonálisan domináns, a mellékátlókban nincs 0 elem és |a11 | > |a12 |, akkor A-nak létezik inverze. 22
Bizonyítás:
a1
b1
c1 a2 A = 0 c2 .. .
0... b2 .. . .. .
0... .. . .. . bn−1
. . . 0 cn−1
an
Tegyük fel, hogy létezik A = LU felbontás. Ekkor L és U a következ®képpen néznek ki, ugyanis a Gauss-elimináció meghagyja a háromátlósságot: u1 v 1 1 l1 1 u2 v2 .. . l2 1 L= U = u3 .. .. .. . . . vn−1 ln−1 1 un Az A = LU i-edik sorából az alábbiakat jelenti: ui−1 vi−1 0 ai = li−1 vi−1 + ui (ci−1 , ai , bi ) = (li−1 , 1, 0) ui vi 0 =⇒ bi = vi c =l u 0 0 ui+1 i i−1 i−1 (3.7) Kihasználjuk, hogy |a1 | > |b1 | és |ai | ≥ |bi |+|ci−1 |. Ezért a rendszer megoldása:
a1 = u1 , b1 = v1 ⇒ |u1 | > |v1 | ⇒ u1 6= 0. Teljes indukcióval belátjuk, hogy ui sem nulla: ha ui 6= 0 és |ui | > |vi |, akkor ui+1 6= 0 és kui+1 k > kvi+1 k. A (3.7)-ban kapott egyenletrendszert átrendezve, és a háromszög-egyenl®tlenséget felhasználva:
|ui+1 | = |ai+1 − li vi | ≥ |ai+1 | − |li ||vi | = |ai+1 | −
|ci | |v | |ui | i
>
> |ai+1 | − |ci | ≥ |bi+1 | = |vi+1 | =⇒ ui+1 6= 0 Így U f®átlójában csupa nemnulla elem áll, ezért létezik A−1 . 23
¤
Speciális eset Nézzük az olyan hiányos másodrend¶ dierenciálegyenletet, ahol
a ≡ 1,
b ≡ 0,
c ≡ 0.
Ekkor a (3.3)-beli dierenciálegyenletünk így írható: ( −u00 (x) = f (x) x ∈ (0, 1)
u(0) = u0 ,
u(1) = u1
A (3.6) egyenletrendszer pedig 2 −1 0 . . . U1 −1 2 −1 0 . . . U2 . .. .. .. . . −1 . . . 0 −1 2 UM −1
=
g1 g2 .. .
gM −1
alakúra egyszer¶södik. Egészítsük ki az egyenletrendszert M +1 dimenziósra: h2 0 0 . . . 0 U0 u0 −1 2 −1 . . . U1 f1 1 .. .. .. .. = . . . . . h2 U f −1 2 −1 M −1 M −1 ...0 0 h2 UM u1 Jelölje a mátrixot A˜h , a szorzóvektort Uh , a jobboldalt pedig fh . Egyenletrendszerünk tehát
A˜h Uh = fh , ahol Uh = (U0 , . . . , UM ),
fh = (u0 , f1 , . . . , fM −1 , u1 ) és fi = f (xi ).
24
2. Állítás. A˜h M-mátrix1 . Bizonyítás: A de níció alapján kell egy g ∈ RM +1 , hogy egyrészt g pozitív legyen, másrészt A˜h g > 0 legyen. Legyen
gi = 2(xi − x2i ). Ekkor
g(x) = 2 + (x − x2 ) > 0 x ∈ (0, 1)-en −g 00 (x) = 2 > 0
Tehát g > 0, az els® feltétel teljesül. Másodszor kell: A˜h g > 0.
(A˜h g)0 =
1 2 h g0 h2
(A˜h g)M =
1 2 h gM h2
= 2 + x0 − x20 = 2 + 0 − 0 = 2 > 0 = 2 + xM − x2M = 2 + 0 − 0 = 2 > 0
£
¤ (−1)gi−1 + 2gi − 1gi+1 =
£
¤ − 2 − xi−1 + x2i−1 + 4 + 2xi − 2x2i − 2 − xi+1 + x2i+1 =
£
¤ − (i−1)h + ((i−1)h)2 + 2ih − 2(ih)2 − (i+1)h + ((i+1)h)2 =
£
¤ −ih+h+i2 h2 −2ih2 +h2 +2ih−2i2 h2 −ih−h+i2 h2 +2ih2 +h2 =
(A˜h g)i =
1 h2
=
1 h2
=
1 h2
=
1 h2
=
1 2h2 h2
=2
Vagyis
2 . . A˜h g = . >0 2
¤ 1 Az M-mátrixok tulajdonságairól szóló de níciókat és állításokat a függelék M-mátrixok
cím¶ fejezetében taglalom.
25
Mivel A˜h M-mátrix, így invertálható és kgk∞ = 2 +
kA˜−1 h k∞ ≤
2+ kgk∞ = mini (Ag)i 2
1 4
=1+
1 4
1 8
2. Példa. Bemutatok egy konkrét példát, ami megmutatja hogyan is számolunk a véges dierencia módszer segítségével. A peremértékfeladat a következ®:
(
−u00 (x) = f (x) =
cos x (sin x+1)2
Ω = (0, 1)-en
u(0) = u(1) = 0 Legyen az osztópontok száma 5, ekkor négy részre osztjuk fel az intervallumot, a lépésköz pedig 1/4. Azaz M = 4,
h = 1/4.
x0 = 0, x1 = 0.25, x2 = 0.5, x3 = 0.75, x4 = 1 32
A= −16 0
−16
0
−16 −16 32 32
j = 0, . . . , 4.
0.6227
fh = 0.4010 0.2587
Mert
cos 0.25 = 0.6227, (sin 0.25 + 1)2
cos 0.5 = 0.4010 (sin 0.5 + 1)2
cos 0.75 = 0.2587 (sin 0.75 + 1)2
Megoldjuk az AUh = fh egyenletrendszert, majd kiírom a kapott értékeket. Uh közelít®, mellette uh pontos értékek, majd ábrázolom:
0
0
0.0466 uh = 0.0535 0.0350 0
0.0458 Uh = 0.0526 0.0344 0 26
Függvény közelítése véges differencia módszerrel 5 osztópont esetén 0.06
0.05
0.04
0.03
0.02
0.01
0
0
0.2
0.4 0.6 (0,1) intervallum
0.8
1
A lenti ábrán pedig ugyanennek a feladatnak a megoldása egy nagyobb (h = 50-es) lépésközzel. (Mindkét ábrán a folytonos vonal jelöli a pontos megoldást, a keresztek pedig az approximációt.) Függvény közelítése véges differencia módszerrel 50 osztópont esetén 0.06
0.05
0.04
0.03
0.02
0.01
0
0
0.2
0.4 0.6 (0,1) intervallum
27
0.8
1
Hibabecslés 1. Lemma (A második deriváltat közelít® dierencia hibája). Tegyük fel, hogy a (3.3) megoldása u ∈ C 4 [a, b]. Ekkor
¯ ¯ ¯ ¯ ¯ ¯ ¯∂ ∂U (xj ) − u00 (xj )¯ ≤ Ch2 ku(4) kC(Ω) Bizonyítás:
u(xj+1 ) − 2u(xj ) + u(xj−1 ) 00 u(xj + h) − 2u(xj ) + u(xj − h) 00 −u (xj ) = −u (xj ) = 2 h h2 Fejtsük Taylor-sorba, xj−1 < η < xj < ξ < xj+1 :
µ =
1 h2
u(xj ) + hu0 (xj ) +
h2 00 u (xj ) 2
+
h3 000 u (xj ) 6
+
h4 (4) u (ξ) 24
− 2u(xj )+
¶ 0
+u(xj ) − hu (xj ) +
h2 00 u (xj ) 2
µ =
1 h2
h4 (4) u (ξ) 24
+
h4 (4) u (η) 24
−
h3 000 u (xj ) 6
¶ + O(h ) = 5
+
h4 (4) u (η) 24
h2 u(4) (ξ)+u(4) (η) 12 2
Alkalmazzuk a Darboux-Bolzano tételt, mely szerint
=
− u00 (xj ) =
+ O(h5 ) =
f (a)+f (b) 2
= f (ξ), így
h2 (4) u (ρ), 12
vagyis 2 ¯ (xj ) − u00 (xj )| ≤ h2 ku(4) kC(Ω) |∂ ∂U ¯ h
1 12 ¤
Ahogy a kezdetiérték feladat esetében itt is átörökl®dik a lokális diszkretizációs hiba rendje a globális hibáéba. A fenti hibabecslés általában nem praktikus (a megoldás négyszer folytonosan dierenciálhatóságának megkövetelése miatt), azért a gyakorlatban a hibát h és h/2 lépésközök eredményeib®l becslik. 28
5. Tétel. A (3.3) feladat megoldása centrális dierencia használata esetében konvergens és
kUh − uh k∞ ≤ Ch2 , ahol uh = (u0 , u1 , . . . , uM ) és Uh = (U0 , U1 , . . . , UM ) illetve ui = u(xi ). (Azaz ha h → 0, akkor kuh − Uh k∞ → 0.) Bizonyítás:
A˜h Uh = fh
Uh = A˜−1 h fh ,
=⇒
ahol fh = (u0 , f1 , . . . , fM −1 , u1 )T
˜−1 ˜ kUh − uh k∞ = kA˜−1 h fh − uh k∞ = kAh (fh − Ah uh )k∞ ≤ 9 ˜ ˜ ≤ kA˜−1 h k∞ kf − Ah uh k∞ ≤ kΨh k∞ , 8 ˜ h = A˜h uh − fh képlethiba. ahol Ψ ( ¯ i − u00 (xi ) (∂ ∂u) i = 1, 2, . . . , M − 1 ˜h = Ψ 0 i = 0, M ¤
Véges dierencia módszer elliptikus dierenciálegyenletekre Vegyük a
(
−4u + qu = r
D-n
u=0
∂D-n
D = (0, 1) × (0, 1)
elliptikus peremérték problémát, ahol 4 =
2 ∂2u + ∂∂xu2 ∂x21 2
(3.8)
Laplace-operátort jelöli.
Ekvidisztáns rácshálót választunk: xij = (ih, ij). A Laplace-operátort az alábbi dierenciasémával közelítjük:
4u(xij ) ≈
1 [u(xi+1,j ) + u(xi−1,j ) + u(xi,j+1 ) + u(xi,j−1 ) − 4u(xi,j )] h2
ami kielégíti az
1 [(4 + h2 qij )uij − ui+1,j − ui−1,j − ui,j+1 − ui,j−1 ] = rij 2 h 29
i, j = 1, . . . , n (3.9)
egyenletrendszert, ahol uij közelít® érték, u(xij ) pontos megoldás, és qij =
q(xij ), rij = r(xij )-t jelöli. Kiegészítve a peremfeltételekkel u0,j = un+1,j = 0 j = 0, . . . , n + 1 ui,0 = ui,n+1 = 0
(3.10)
i = 1, . . . , n
6. Tétel. A (3.9) és a (3.10) egyneletrendszernek egyértelm¶ megoldása van minden h > 0 esetén.
7. Tétel. Tegyük fel, hogy a (3.8)-ben adott feladat megoldása négyszer folytonosan dierenciálható. Ekkor a véges dierenciás approximáció hibabecslését a következ® formulával számítjuk: ° 4 ° ¸ ·° 4 ° ° ° ° h2 ° ∂ u ° ° + ° ∂ u ° , i, j = 1, . . . , n |u(xi j) − ui j| ≤ 4° 4° ° ° 96 ∂x1 ∞ ∂x2 ∞
3.3. Végeselem módszer A függelék 4.3. fejezetében kimondott Riesz és Lax-Milgram tételekre épül az alábbi módszer. A tételek bizonyítása az [1] könyvben található.
(
−(pu0 )0 + qu = r
(3.11)
u(a) = u(b) = 0
A végeselem a (3.11) dierenciálegyenletre a homogén peremfeltétellel a Galjorkin módszer alkalmazásán alapszik, a spline tereket közelít® alterekként használva. Megjegyezzük, hogy a polinomok nem alkalmasak közelít® altérnek, mert rosszul kondícionált teli-mátrixos lineáris rendszerekhez vezetnek. Vegyük a lineáris spilne-ok esetét. Az
xj := a + jh,
j = 0, . . . , n + 1,
h = (b − a)/(n + 1),
n∈N
lépésközzel válasszuk a folytonos és darabonként lineáris függvények terét
Xn -nek. Azaz ∀u ∈ Xn :
u ∈ C[a, b], 30
u(a) = u(b) = 0
és minden egyes [xj−1 , xj ] részintervallumon egybeesik egy P1 -beli polinommal j = 0, . . . , n-re. E spline tér beli függvények a H01 [a, b] térhez tartoznak darabonként konstans gyenge deriváltakkal.
Xn báziselemeinek az un. wk (x) :=
kalapfüggvények et választjuk: 1 (x h
− xk−1 ),
1 (xk+1 h
x ∈ [xk−1 , xk ],
− x),
x ∈ [xk , xk+1 ],
0
(3.12)
x∈ / [xk−1 , xk+1 ].
Minden u ∈ Xn felírható
u=
n X
αk wk ,
k=1
alakban, ahol αk = u(xk )
k = 1, . . . , n. Így Z
b
S(wj , wk ) = a
{pwj0 wk0 + qwj wk }dx = 0
ha (xj−1 , xj+1 ) ∩ (xk−1 , xk+1 ) = ∅, azaz ha |j − k| > 2. Ekkor az S(wj , wk ) mátrix a (3.12) kalapfüggvények, mint báziselemek mellet tridiagonális és ritkamátrix. A mátrix elemeit következ®képpen számítjuk: ½ Z xj ¾ Z xj+1 Z xj+1 1 1 2 2 S(wj , wj ) = 2 q(x)(x−xj−1 ) dx+ p(x)dx+ 2 q(x)(xj+1 −x) dx h xj−1 h xj−1 xj és
1 S(wj , wj+1 ) = S(wj+1 , wj ) = − 2 h
Z
xj+1 xj
1 p(x)dx+ 2 h
Z
xj+1 xj
q(x)(xj+1 −x)(x−xj )dx,
és az egyneletrendszer jobb oldala: ½ Z xj ¾ Z xj+1 1 F (wj ) = r(x)(x − xj−1 )dx + r(x)(xj+1 − x)dx . h xj−1 xj A fenti képletek illusztrálják a végeselem módszer két általános tulajdonságát: 31
1. Minden részintervallumon ugyanazzal a képlettel számítjuk ki az együtthatókat 2. A Galjorkin módszert szeretnénk diszkrétté tenni, ezért használjuk a numerikus kvadratúraformulákat. Lineáris spline-okkal approximálva
p, q, r-et: S(wj , wj ) ≈
1 h (pj−1 + 2pj + pj+1 ) + (qj−1 + 6qj + qj+1 ) 2h 12
és
S(wj , wj+1 ) ≈
−1 h (pj + pj+1 ) + (qj + qj+1 ) 2h 12
a mátrixelemekre és
F (wj ) ≈
h (rj−1 + 4rj + rj+1 ) 6
a jobb oldalakra, ahol pj = p(xj ), qj = q(xj ), rj = r(xj ) ahogy a korábbi fejezetekben is. A (3.6) véges dierencia módszerbeli lineáris egyenletrendszerhez hasonlóan a fenti rendszer is irreducibilis és gyengén diagonálisan domináns, így a Jacobi-iteráció konvergenciája gyorsítható relaxációs módszerekkel. A lineáris spilne-okkal végzett végeselem módszer hibájának levezetéséhez szükségünk van a lineáris spline interpoláció hibájának becslésére a H 1 normára tekintettel.
2. Lemma. Legyen f ∈ C[a, b]. Legyen R1 f := f − L1 f a és b végpontokra illeszked® lineáris interpoláció. Ekkor az alábbi becslések igazak R1 f -re:
kR1 f kL2 ≤ (b − a)2 kf 00 kL2 , k(R1 f )0 kL2 ≤ (b − a)2 kf 00 kL2 .
8. Tétel. A lineáris spline-okkal végzett végeselem approximáció hibája a (3.11) dierenciálegyenletre a homogén peremfeltétellel az alábbi képlettel számítható ki:
kun − ukH 1 ≤ Cku00 kL2 h, ahol C tetsz®leges pozitív konstans. 32
Bizonyítás:
Felhasználva a fenti lemmabeli egyenl®tlenségeket és feltéve,
hogy minden részintervallum hossza h, a wn ∈ Xn lináris interpoláló splinera
wn (xj ) = un (xj ),
j = 0, . . . , n feltételekkel azt kapjuk, hogy kwn0 − u0 kL2 ≤ ku00 kL2 h
és
kwn − ukL2 ≤ ku00 kL2 h így
inf kv − ukH 1 ≤ kwn − ukH 1 ≤ (1 + b − a)ku00 kL2 h.
v∈Xn
3. Példa.
(
−u00 (x) = 1 u(−1) = 0
¤
(3.13)
u(1) = 0
A fenti dierenciálegyenletet szeretnénk megoldani a Dirichlet-peremfeltétellel a végeselem módszer alkalmazásával. Szorozzuk be a (3.13) dierenciálegyenletet egy ϕ(x) tesztfüggvénnyel, ami az intervallum két végpontjában 0-t vesz fel és folytonosan dierenciálható:
−u00 ϕ = ϕ Majd integráljuk mindkét oldalt a [-1,1] intervallumon. Z 1 Z 1 00 −u ϕ = ϕ −1
Az egyenlet baloldala: Z Z 1 00 0 1 −u ϕ = [−u ϕ]−1 + −1
1
−1
Z 0
0
1
uϕ = −1
−1
Z 0
0
u (x)ϕ (x)dx =
1 −1
Mivel C01 [−1, 1] s¶r¶ H01 [−1, 1]-ben, ezért Z 1 Z 1 0 0 uϕ = ϕ ∀ϕ ∈ H01 [−1, 1]-re −1
−1
ϕ(x)dx
∀ϕ ∈ C01 [−1, 1]
(3.14)
Ezt úgy oldjuk meg, hogy H01 helyett egy véges dimenziós altérben keresünk megoldást. Legyen például dim=3. 33
1
1
1
0.9
0.9
0.9
0.8
0.8
0.8
0.7
0.7
0.7
0.6
0.6
0.6
0.5
0.5
0.5
0.4
0.4
0.4
0.3
0.3
0.3
0.2
0.2
0.1
0.1
0 −1
−0.5
0
0.5
0.2 0.1
0 −1
1
−0.5
0
w1
0.5
0 −1
1
−0.5
0
w2
0.5
1
w3
Keressük a megoldást abban az altérben, amelyet a w1 , w2 , w3 ∈ H01 [−1, 1] függvények feszítenek ki.
3 X
u=
ci wi
i=1
Összeget tudunk tagonként deriválni: 3 X
u0 =
ci wi0
i=1
Folytatva a (3.14)-beli egyenlet bal oldalát Z Z 1X Z 1 3 3 X 0 0 0 0 ci wi ϕ = ci uϕ = −1 i=1
−1
Összekapcsolva a jobboldallal: Z Z 1 3 X 0 0 wi ϕ = ci i=1
−1
i=1
1
1 −1
wi0 ϕ0
∀ϕ ∈ H01 [−1, 1]-re
ϕ −1
De ezt most csak az altér elemire követeljük meg: ∀ϕ ∈ hw1 , w2 , w3 i. Elég ha a báziselemekre teljesül: Z 3 X ci i=1
1 −1
Z wi0 wj0
1
=
wj · 1
(3.15)
j = 1, 2, 3,
−1
ahol a jobb oldalon az "1"-es a (3.13)-beli jobboldal, vagyis f (x) függvény. Írjuk mátrixos alakba a (3.15)-beli egyenletrendszert: Z 1 Z 1 Z 1 Z 1 0 0 0 0 0 0 −1 w1 w1 −1 w2 w1 −1 w3 w1 c1 −1 w1 Z 1 Z 1 Z 1 Z 1 0 0 0 0 0 0 = w2 w3 w2 w2 w2 w1 w2 c −1 −1 −1 Z−11 2 Z 1 Z 1 Z 1 0 0 0 0 0 0 w3 w3 w3 w2 w3 w1 w3 −1 −1 −1 −1 c3 34
Behelyettesítve a mátrixba a konkrét értékeket c1 1/2 4 −2 0 c2 = 1/2 −2 4 −2 1/2 0 −2 4 c3 Kiszámolva a fenti egyenletrendszert 3/8 c= 1/2 3/8
3 1 3 U (x) = c1 w1 + c2 w2 + c3 w3 = w1 + w2 + w3 8 2 8 Végeselem módszer közelítése 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 −1
−0.5
0
0.5
1
Az ábrán a töröttvonal a végeselem módszer U (x) eredménye, a görbe pedig az u(x) =
u2 2
+
1 2
pontos megoldás három dimenziós altérben. (A pon-
tosság a dimenziószám növelésével javítható.)
35
3.4. Összefoglalás Három numerikus módszert mutattam be különböz® közönséges dierenciálegyenletek peremérték-feladataira. A belövéses módszer esetén nincs semmilyen megkötés sem a dierenciálegyenletre (az expliciten kifejezett u00 =
f (x, u, u0 ) egyenlet jobboldala bármilyen lehet), sem a peremfeltételekre, kifejezetten nem-lineáris feladatok közelítésere alkották. A véges dierencia és végeselem módszerek is alkalmazhatóak nem-lineáris esetben, de ott nem m¶ködnek jól, ellenben lineáris esetben egyszer¶ek és gyorsak. A végeselem módszerre jelen dolgozatban csak els®-derivált tagot nem tartlamazó lineáris dierenciálegyenletre használható módszert mutattam be, de valójában ez a legelterjedtebb módszer. Ritkán fordulnak el® olyan feladatok, amiket végeselem módszerrel nem lehet megoldani. Ellenben a belövéses módszer túl sok számolást igényel, és az F (s) = 0 egyenlet megoldása is bizonytalan, ui. a Newton-módszer nem mindig konvergens. (Rosszul kondícionált feladatok esetén pedig csak a többszörös belövéses módszer használható.) A [3] könyv egy (rosszulkondícionált példán2 ) hasonlítja össze a különböz® módszerek hibáit. A bemutatott példában: módszer
átlagos hiba
Belövéses módszer
∞
Többszörös belövéses módszer
≈ 10−13
Véges dierencia módszer
10−2 10−10
Végeselem módszer
10−3 10−6
Megjegyzés: ahol intervallumot adtam meg hibára, ott értéke arányos valamilyen a módszerben választható tényez® nagyságával. A véges dierencia módszer esetén a hiba nagysága h lépésköz megválasztásától függ, ami itt
h = 2−4 2−10 , a végeselem módszer esetén pedig a dimenziómérett®l, ami itt n = 10100 közötti.
2 −u00
+ 400u = −400 cos2 πx − 2π 2 cos 2πx
36
4. fejezet A függelék A függelékben szerepl® állításokat és tételeket bizonyítás nélkül közlöm. A bizonyítások megtalálhatók [1]-ben.
4.1. M-mátrixok tulajdonságai 1. De níció. Egy A ∈ Rn×n mátrix M-mátrix, ha • az átlóbeli elemek pozitívak, azaz aii > 0 (i = 1, . . . , n); • az átlón kívüli elemek nem-pozitívak, azaz aij ≤ 0 (i 6= j); • létezik egy 0 < g ∈ Rn minden koordinátájában pozitív vektor, amire 0 < Ag (szintén elemenként).
1. Megjegyzés. • A fenti de nícióban szerepl® g vektort A dominálóvektorának nevezzük. • A fenti de níció els® feltétele akövetkezik a másik kett®b®l, máskülönben a
akk ≤ 0
=⇒
0 < (Ag)k =
X i
ellentmondásra jutunk. 37
akj gj ≤ 0 |{z} |{z} ≤0
>0
9. Tétel. Legyen A egy M-mátrix, ekkor • A reguláris, azaz létezik A−1 • A−1 ≥ 0 elemenként.
10. Tétel. Legyen A egy M-mátrix g domináló vektorral, ekkor
kA−1 k∞ ≤
kgk∞ . mini (Ag)i
4.2. Funkcionálanalízis alapok Normált terek 2. De níció. Legyen X egy valós vagy komplex vektortér a k · k : X → R függvény a következ® tulajdonságokkal:
(N 1) kxk ≥ 0
(pozitivitás)
(N 2) kxk = 0 ⇔ x = 0
(meghatározottság)
(N 3) kαxk = |α|kxk
(homogenitás)
(N 4) kx + yk ≤ kxk + kyk (háromszög-egyenl®tlenség) minden x, y ∈ X és ∀α ∈ C( vagy R)-ra norma X -en. X Normált tér, ha értelmezve van rajta egy norma.
2. Megjegyzés. A másik háromszög-egyenl®tlenség is érvényes a normákra: ¯ ¯ ¯kxk − kyk¯ ≤ kx − yk ∀x, y ∈ X
3. De níció. Egy normált téren vett (xn ) vektorsorozat konvergens, ha ∃x ∈ X amire lim kxn − xk = 0. Divergens, ha nem konvergens. n→∞
11. Tétel. Egy konvergens sorozat határértéke pontosan meghatározott. 4. De níció. Egy lineáris téren vett bármely két norma ekvivalens, ha ugyanazok a sorozatok konvergálnak hozzájuk. 38
12. Tétel. k · ka és k · kb normák pontosan akkor ekvivalensek, ha ∃c, C > 0:
ck · ka ≤ k · kb ≤ Ck · ka minden x ∈ X -re.
13. Tétel. Véges dimenziós vektortereken minden norma ekvivalens. 5. De níció. Egy X normált tér U alterét korlátosnak nevezzük, ha létezik olyan C pozitív szám amire kxk ≤ C minden x ∈ U -ra. (A konvergens sorozatok korlátosak.)
14. Tétel. X véges-dimenziós normált téren bármely korlátos sorozatnak van konvergens részsorozata.
Skalárszorzat 6. De níció. Legyen X egy komplex (vagy valós) lineáris tér. Ekkor a (·, ·) : X × X → C( vagy R) függvényt az alábbi tulajdonságokkal (H1) (x, x) ≥ 0
(pozitivitás)
(H2) (x, x) = 0 ⇔ x = 0
(meghatározottság)
(H3) (x, y) = (y, x)
(szimmetria)
(H4) (αx + βy, z) = α(x, z) + β(y, z) (linearitás) minden x, y, z ∈ X
α, β ∈ C( vagy R)-ra X -en vett skalárszorzatnak vagy
bels®-szorzatnak nevezzük.
7. De níció. Egy skalárszorzattal ellátott teret pre-Hilbert térnek nevezünk. (H3) és (H4) triviális következménye az antilinearitás:
¯ z) (H40 ) (x, αy + βz) = α ¯ (x, y) + β(x,
15. Tétel (Cauchy-Bunyakovszky-Schwarz egyenl®tlenség). |(x, y)|2 ≤ (x, x)(y, y) ∀x, y = X -re és egyenl®ség pontosan akkor áll fenn, ha x és y linaárisan függetlenek. 39
16. Tétel. Az X lineáris téren vett (·, ·) skalárszorzat meghatároz egy normát, amit az alábbi képlet ad meg:
kxk := (x, x)1/2 minden x ∈ X -re, azaz a pre-Hilbert tér mindig normált tér is egyben.
Korlátos lineáris operátorok Az alábbiakban A : X → X kifejezésre a leképezés, függvény és operátor azavakat szinonimaként használjuk.
8. De níció. U ⊆ X, Y normált terek, A : X → X leképezés folytonos ∀x ∈ U -ban, ha minden (xn ) : lim xn = x
U -beli sorozat esetén lim Axn = Ax.
n→∞
n→∞
9. De níció. Egy A : X → X leképezés lineáris, ha A(αx + βy) = αAx + βAy
∀x, y ∈ X,
α, β ∈ C vagy R.
17. Tétel. Egy lineáris operátor folytonos, ha egy elemében folytonos. 10. De níció. A : X → Y lineáris operátor korlátos, ha létezik C pozitív szám, hogy kAxk ≤ Ckxk minden x ∈ X.
18. Tétel. A : X → Y lineáris operátor akkor és csak akkor korlátos, ha kAk := sup kAxk < ∞. kxk=1
Az kAk azám A legalsó korlátja, és A normájának nevezzük.
3. Megjegyzés. X, Y, Z normált terek, A : X → Y, B : Y → Z korlátos leképezések. Szorzatuk BA : X → Z,
(BA)x := B(Ax) ∀x ∈ X szintén
korlátos leképezések a kBAk ≤ kAkkBk korláttal. 40
Mátrixnormák 19. Tétel. Legyen (ajk ) ∈ Rn×n vagy Cn×n -beli mátrix. Ekkor az alábbi képlettel megadott A : Rn → Rn és A : Cn → Cn lineáris operátor
(Ax)j :=
n X
ajk xk ,
j = 1, . . . , n,
k=1
megfelel egyes normáknak Rn -en illetve Cn -en. Konkrét esetben: n X
kAk1 = max
k=1,...,n
kAk∞ = max
n X
j=1,...,n
à kAk2 ≤
n X
|ajk |,
j=1
|ajk |,
k=1
!1/2 |ajk |2
.
j,k=1
Ekkor a fenti normákat mátrixnormáknak nevezzük.
20. Tétel. Minden A mátrixhoz létezik egy Q unitér mátrix, hogy Q∗ AQ egy fels®-háromszög mátrix.
21. Tétel. Az n × n-es Hermite-mátrix1 sajátértékei valósak, és a sajátvektorok ortogonális bázist alkotnak C n -ben.
11. De níció. %(A) := max{|λ| : λ sajátértéke A-nak} A spektrálsugara. 22. Tétel. A n × n-es mátrix. kAk2 =
p %(A∗ A)
Ha A Hermite-mátrix, akkor kAk2 = %(A)
23. Tétel. Minden Cn -en értelmezett normára és minden n×n-es A mátrixra %(A) ≤ kAk. 1A
Hermite-mátrix, ha A∗ = A, vagyis, ha önadjungált
41
Teljesség 12. De níció. Az X normált téren egy (xn ) sorozat Cauchy-sorozat, ha ∀ε > 0 ∃N (ε) ∈ Z : Máshogy:
kxn − xm k < ε
∀n, m ≥ N (ε) esetén.
lim kxn − xm k = 0
n,m→∞
24. Tétel. Minden konvergens sorozat Cauchy-sorozat. 13. De níció. Az U ⊆ X altér teljes, ha U minden elemének Cauchysorozata U egy eleméhez konvergál.
14. De níció. Ha egy normált tér teljes Banach-térnek nevezzük. Ha egy pre-Hilbert2 tér teljes, Hilbert-térneknevezzük.
4. Példa. Az alábbi módon de niált maximum- vagy sornormával ellátott [a, b] intervallumon folytonos függvények tere Banach-tér. kf k∞ := max |f (x)| [a,b]
5. Példa. Az L1 normával ellátott [a, b] intervallumon folytonos függvények tere, ahol
Z
b
kf k1 := a
f (x)dx
nem teljes.
6. Példa. Az L2 normával ellátott [a, b] intervallumon folytonos függvények tere, ahol
µZ
b
kf k2 := a
¶1/2 2
|f (x)| dx
sem teljes.
25. Tétel. Minden véges dimenziós tér Banach-tér. 4. Megjegyzés. A teljes terek zártak és minden teljes tér zárt altere is teljes. 2 Ismétlés:
a pre-Hilbert tér a skalárszorzattal ellátott lineáris tér
42
4.3. Elméleti alapok a végeselem módszerhez A végeselem módszer bevezetéséhez megemlítek néhány funkcionálanalízisbeli tételt és állítást:
26. Tétel (Riesz-tétel). Legyen H Hilbert-tér. Ekkor minden korlátos lineáris F : H → C függvényhez ∃!f ∈ H amire F (u) = (u, f ) ∀u ∈ H esetén. És f illetve F lineáris függvények normái egyenl®ek: kf k = kF k.
15. De níció. Egy A : H → H lineáris operátor az H pre-Hilbert téren szigorúan koercitív, ha ∃c > 0 konstans: Re(au, u) ≥ ckuk2
∀u ∈ H -re.
27. Tétel (Lax-Milgram-tétel). Legyen H Hilbert-tér, a(·, ·) egy korlátos(1), koercitív(2) bilineáris forma, l egy folytonos(3) lineáris funkcionál, azaz 1. ka(u, v)k ≤ M kukkvk 2. ka(u, v)k ≥ αkvk2 , 3. |l(v)| ≤ ckvk Ekkor az a(u, v) = l(v)
∀u, v ∈ H ahol α > 0
∀v ∈ H
∀v ∈ H ∀v ∈ H feladatnak egyértelm¶en létezik megoldása.
43
5. fejezet B függelék 5.1. Belövéses módszer algoritmus %A belövéses módszerrel oldjuk meg az %u''(x)=(u(x))^3
[a,b]-n
%u(a)=alpha=gyök(2), u(b)=beta=gyök(2)/2 %peremérték-feladatot alpha = sqrt(2); beta = 0.5*sqrt(2); xspan = [1, 2]; aeps = 1e-6; h = 0.01; error = Inf; counter = 0; s = 0; while ((error > aeps) & (counter < 10)) U = []; f1 = @(x,u)(u(2)); f2 = @(x,u)(u(1)^3); [x, U] = ode_Euler(f1, f2, xspan, [alpha, s], U, h, 1); 44
g1 = @(x,z)(z(2)); g2 = @(x,z)(3*z(1)); [x, Z] = ode_Euler(g1, g2, xspan, [0, 1], U(:, 1), h, 0); plot(x,U(:,1),'k:'); hold on; F = U(end, 1) - beta; error = abs(F); s = s - F/Z(end, 1); counter = counter + 1; end counter; U1=U(:,1); Y=sqrt(2)*ones(length(x),1)./x; plot(x,Y,'k'); title('Belövéses módszer'); text(1.85,3.5,'s(0)'); text(1.9,1.85,'s(1)'); text(1.85,1,'s(2)'); text(1.8,0.65,'s(3)-s(5)'); A fenti program futás közben az alábbi ode_Euler.m fájlt hívja meg:
function [x Y] = ode_Euler(f1, f2, xspan, Y0, U, h, fuggveny) n = (xspan(2)-xspan(1))/h; x = linspace(xspan(1), xspan(2), n+1)'; k = zeros(1, 2); Y = Y0; for i=1:n if fuggveny == 1 ertek1 = feval(f1, x(i), Y(i, :)); ertek2 = feval(f2, x(i), Y(i, :)); 45
k(1, :) = [ertek1, ertek2]; Y(i+1, :) = Y(i, :) + h*k(1, :); else ertek1 = feval(f1, x(i), Y(i, :)); ertek2 = feval(f2, x(i), Y(i, :)); k(1, :) = [ertek1, ertek2*(U(i, 1))^2]; Y(i+1, :) = Y(i, :) + h*k(1, :); end end
5.2. Véges dierencia algoritmus Speciális eset A 3.2. fejezet 2. példájában levezetett példában használt algoritmus forráskódja:
%Véges differencia módszer %-u''(x)=f(x) u(1)=y(1) u(N)=y(N) %a példa: f(x)=cos(x)/(sin(x)+1)^2 M=4; N=M+1;
%N az osztópontok száma beleértve a végpontokat is
x=linspace(0,1,N)'; h=(x(N)-x(1))/(N-1); %elkészül az A=(1/h^2)tridiag(-1,2,-1) mátrix B1=-ones(N-3,1); B2=2*ones(N-2,1); A1=diag(B1,-1); A2=diag(B2,0); A3=diag(B1,1); A=(1/h^2)*(A1+A2+A3); %a lineáris egyenletrendszer jobb oldala 46
g=zeros(N,1); for i=2:N-1, g(i)=cos(x(i))/(sin(x(i))+1)^2; end g=g(2:N-1) % mivel a feladatot homogén peremfeltétellel adtuk meg nincs szükség g % további módosítására y1=A\g y=zeros(N,1); y(2:N-1)=y1; %a pontos megoldás (u(x)=-2/(1+tan(x)/2))-.70659*x+2) xx=linspace(0,1,100); u=-2./(1+tan(xx/2))-.70659.*xx+2; plot(xx,u,'k'); hold on; %közelít® megoldás plot(x,y,'k+','LineWidth',2);
Általános eset Bemutatok egy másik programot, amely az általánosabb ( Au := −au00 + bu0 + cu = f Ω = (0, 1)
u(0) = u0 ,
u(1) = u1
alakú peremérték-feladatokon is végrehajtja a véges dierencia módszert.
7. Példa. A feladat a következ®: −(1 + x2 )u00 (x) − xu0 (x) + u(x) = 0 Ω = (0, 1)-en √ x(0) = 1, x(1) = 2 47
%Véges differencia módszer %függvényegyütthatós lineáris másodrend¶ differenciálegynletre %a(x)u''(x)+b(x)u'(x)+c(x)u(x)=f(x) %u(1)=y(1) u(N)=y(N) %(a példa: -(1+x^2)*u''(x)-x*u'(x)+u(x)=0 Omega=(0, 1)-en %
x(0)=1, x(1)=gyok(2) peremfeltétellel
N=10; x=linspace(0,1,N)'; y=zeros(N,1); Y(1)=1; Y(N)=sqrt(2); h=(x(N)-x(1))/(N-1); a=1+x.^2; b=-x; c=ones(N,1); f=zeros(N,1); %elkészül az A tridiagonális mátrix aj=a(2:N-1); aj1=a(3:N-1); aj2=a(2:N-2); bj1=b(3:N-1); bj2=b(2:N-2); cj=c(2:N-1); B0=2*aj+h^2*cj; B1=-aj1-0.5*h*bj1; B2=-aj2+0.5*h*bj2; A0=diag(B0,0); A1=diag(B1,-1);
A2=diag(B2,1); A=(A0+A1+A2);
for i=2:N-1, f(i)=0; end g=f(2:N-1); g(1)=(aj(1)+0.5*h*bj2(1))*Y(1); g(N-2)=(aj(N-2)-0.5*h*bj1(N-3))*Y(N); 48
y1=A\g; y=zeros(N,1); y(2:N-1)=y1; y(1)=Y(1); y(N)=Y(N); %a pontos megoldás: u(x)=(1+x^2)^(1/2) xx=linspace(0,1,100)'; uu=sqrt(1+xx.^2); plot(xx,uu); hold on; %közelít® megoldás plot(x,y,'+'); %a pontos megoldás értékei az osztópontokban: u=sqrt(1+x.^2); [y u] Megoldása pedig 10 osztópont esetén a fels® sorban a közelít®, az alsóban a pontos megoldásokkal: uu 1.0000 1.0061 1.0243 1.0540 1.0942 1.1439 1.2018 1.2668 1.3379 1.4142 y
1.0000 1.0062 1.0244 1.0541 1.0943 1.1440 1.2019 1.2669 1.3380 1.4142 1.45 1.4 1.35 1.3 1.25 1.2 1.15 1.1 1.05 1
0
0.2
0.4
0.6
0.8
49
1
Irodalomjegyzék [1] R. Kress: Numerical analysis, Springer (1998) [2] S. Larsson, V. Thomee: Partial dierential equations with numerical
methods, Springer (2003) [3] J. Stoer, R. Bullirsch: Introduction to numerical analysis, Springer (1992) [4] Stoyan Gisbert: MATLAB 4. és 5. verzió, TYPOTEX (1999) [5] K. K. Ponomaljov: Dierenciálegyenletek felállítása és megoldása, Tankönyvkiadó, Budapest (1969)
Köszönetnyílvánítás Ezúton is szeretnék köszönetet mondani témavezet®mnek, Kurics Tamásnak szakdolgozatom elkészítésében nyújtott segítségéért. Köszönöm segít®kész támogatását, dolgozatom alapos és kritikus átnézését, javaslatait.
50
