Numerikus módszerek 2 (MSc)

Numerikus m´odszerek 2 (MSc) Horv´ath R´ obert Budapesti M˝ uszaki ´ es Gazdas´ agtudom´ anyi Egyetem Term´ eszettudom´ anyi Kar Matematika Int´ ezet Anal´ızis Tansz´ ek

2015. tavasz

1. Tant´argyismertet˝o

Hasznos inform´aci´ok

I

El´erhet˝os´egeim: e-mail: [email protected], Iroda: H.24/b

I

Fogad´o´ora: H´etf˝o 9-10 ´ ora (az irod´amban).

I

Gyakorlat: Beadand´ o feladatok, programok. A feladatok teljes´ıt´ese ´es jelenl´et sz¨ uks´eges az al´a´ır´ashoz.

I

El˝oad´as: k¨ovetelm´enyek ismertet´ese

I

Vizsg´ak: vizsga m´ odja, vizsg´ara bocs´at´as felt´etele

I

Jegyzet: ´orai jegyzet, Horv´ath, Izs´ak, Kar´atson, Parci´alis differenci´alegyenletek numerikus m´ odszerei, elektronikus jegyzet, 2013, az el˝ oad´as di´ai, interneten szerepl˝o anyagok: www.math.bme.hu/~rhorvath

Felhaszn´alt irodalom - Arent, Urban, Partielle Differenzialgleichungen, Spektrum, 2010. - Larry J. Segerlind, Applied finite element analysis, John Wiley and sons, 1984. - Morton, Mayers, Numerical solutions of partial differential equations, Cambridge University Press, 2005. - Alfio Borzi, Introduction to multigrid methods, Institut f¨ ur Mathematik und Wissenschaftliches Rechnen, Graz (http://www.kfunigraz.ac.at/imawww/borzi/mgintro.pdf) - C. Johnson, Numerical solutions of PDEs by the finite element methods, Dover Publications 2009. - Stoyan Gisbert, Matlab, Typotex 2005. - Stoyan Gisbert, Tak´o Galina, Numerikus m´ odszerek II-III, Typotex. - Farag´o, Horv´ath, Numerikus m´ odszerek BSc, egyetemi jegyzet (honlapomr´ol let¨olthet˝ o).

2. Bevezet´es

Mivel foglalkozik a numerikus matematika?

A numerikus matematika a folytonos matematika probl´em´ainak megold´as´ahoz konstru´al algoritmusokat ´es elemzi azokat (pontoss´ag, hat´ekonys´ag, viselked´es a sz´am´ıt´ og´epes megval´os´ıt´asok sor´an). A folytonos matematikai probl´em´ak ´altal´aban valamilyen m´as tudom´any´ag probl´em´ainak matematikai modelljei (fizika, biol´ogia, k¨ozgazdas´agtan, stb.).

Modellalkot´as

A val´ odi probl´ema ↓ Tudom´anyos modell ↓ Matematikai modell → PDE ↓ Numerikus modell ↓ Sz´am´ıt´ og´epes modell

N´eh´any p´elda (tud. ´es mat. modellek) ”A boldog csal´adok hasonl´ ok egym´ashoz; minden boldogtalan csal´ad, a maga m´odj´an az.” (Tolsztoj)

Transzportegyenletek: Pl. A keresztmetszet˝ u cs˝ oben folyad´ek ´aramlik, oldott anyag s˝ ur˝ us´ege u(t, x). Egy [x, x + ∆x] cs˝ odarabot tekint¨ unk a [t, t + ∆t] id˝ointervallumban. M´erlegegyenlet (megmarad´asi t¨ orv´eny): T¨omegv´altoz´as = Be´araml´ o - Ki´araml´ o + Forr´as ψ(t, x): fluxus: egys´egnyi keresztmetszeten, egys´egnyi id˝o alatt ´at´araml´o anyag. f (t, x): forr´as: id˝oegys´eg alatt egys´egnyi t´erfogatban keletkez˝o anyag.

N´eh´any p´elda (tud. ´es mat. modellek) Z

x+∆x

A(u(t + ∆t, y) − u(t, y)) dy = x

Z

t+∆t

Z

t+∆t

Z

x+∆x

A(ψ(τ, x) − ψ(τ, x + ∆x)) dτ +

= t

Af (τ, y) dydτ t

x

⇓ u0t (t, x) + ψx0 (t, x) = f (t, x) Speci´alis esetek: ψ(t, x) = cu(t, x) → line´aris transzportegyenlet u0t (t, x) + cu0x (t, x) = f (t, x), u0t (t, x) + cu0x (t, x) − λu(t, x) = 0 (advekci´o-reakci´o egyenlet)

N´eh´any p´elda (tud. ´es mat. modellek) ψ(t, x) = αu(β − u) → Burgers-egyenlet (k¨ ozleked´esi dug´o modell). A v(t, x) = β − 2u(t/α, x) v´altoz´ ocser´evel vt0 (t, x) + v(t, x)vx0 (t, x) = 0.

Diff´ uzi´ o, h˝ ovezet´ es Fourier-t¨orv´eny, Fick-t¨ orv´eny: ψ(t, x) = −ku0x (t, x) → u0t (t, x) = ku00xx (t, x) + f (t, x).

Hull´ amegyenlet Newton 2. t¨orv´enye: u00tt = c2 u00xx + f (t, x).

N´eh´any p´elda (tud. ´es mat. modellek) Ugyanezek magasabb t´erdimenzi´ oban: ~ x) = 0. - Transzportegyenlet: u0t (t, x) + ∇ψ(t, ~ x) = c~bu(t, x) → u0 (t, x) + c~b∇u(t, x) = 0. - ψ(t, t - Diff´ uzi´o, h˝ovezet´es u0t (t, x) = k∆u + f (t, x).

Spec.: (∆u = f Poisson-egyenlet, ∆u = 0 Laplace-egyenlet) - Hull´amegyenlet u00tt (t, x) = c2 ∆u + f (t, x).

N´eh´any p´elda (tud. ´es mat. modellek)

Black–Scholes-egyenlet Opci´o´araz´as egyenlete: 1 00 + rSVS0 − rV = 0 Vt0 + σ 2 S 2 VSS 2

Navier–Stokes-egyenlet Nem¨osszenyomhat´o folyad´ek ´araml´asa:
% ~vt0 + ~v · ∇~v = −∇p + µ∆~v + f~.

N´eh´any p´elda (tud. ´es mat. modellek) Maxwell-egyenletek Elektrom´agneses t´er le´ır´asa: ∂(εE) = ∇ × H − σE − Je , ∂t ∂(µH) = −∇ × E, ∂t ∇(εE) = 0, ∇(µH) = 0.

Schr¨ odinger-egyenlet ~2 ∆Ψ + V (t, x)Ψ. 2m A PDE-ket kezdeti- ´es peremfelt´etelek teszik korrekt kit˝ uz´es˝ uv´e. i~Φ0t = −

3. M´odszert´ıpusok

M´odszert´ıpusok ´attekint´ese I. V´ eges differenci´ ak m´ odszere II. Integr´ alm´ odszerek I

Vari´ aci´ os m´ odszer (Ritz-m´ odszer)

I

S´ ulyozott reziduum m´ odszerek - Kollok´aci´os m´odszer - V´eges t´erfogat m´ odszer (vagy r´esztartom´any m´odszer) - Galjorkin-m´odszer - Momentumok m´ odszere - Legkisebb n´egyzetek m´ odszere

Megj.: Ha az integr´alm´ odszerekben a tesztf¨ uggv´enyeket kis tart´oj´ u, szakaszonk´ent polinomi´alis f¨ uggv´enyeknek v´alasztjuk, akkor v´egeselem m´odszerr˝ ol besz´el¨ unk.

Egy egydimenzi´os perem´ert´ekfeladat

Lehajl´o r´ ud egyenlete: y 00 =

M (x) ≡ M0 , y(0) = y(H) = 0. EI

A pontos megold´as: y(x) =

M0 2 (x − Hx), 2

a r´ ud k¨ozep´en a lehajl´as y(H/2) = −

M0 H 2 = −0.125M0 H 2 . 8

V´eges differenci´ak m´odszere Az ismeretlen f¨ uggv´enyt r´acsf¨ uggv´ennyel k¨ ozel´ıtj¨ uk, a deriv´altakat pedig v´eges differenci´akkal. Az ´ıgy nyert egyenletrendszert megoldjuk az ismeretlen r´acsf¨ uggv´eny´ert´ekekre. Legyenek Q = tridiag [−1, 2, −1] ∈ Rn×n , h = H/(n + 1), n ∈ N, xk = kh (k = 0, . . . , n + 1). yk k¨ ozel´ıtse az y(xk ) f¨ uggv´eny´ert´eket. Ekkor a −1 Qy = M0 e h2 egyenletrendszert kell megoldani. n=101 eset´en a megold´asvektor k¨ oz´eps˝ o eleme: y51 = −0.124999999999999M0 H 2 ≈ −M0 H 2 /8. (l´asd: Numerikus m´odszerek 1.)

Integr´alm´odszerek

Tesztf¨ uggv´enyek Vh N dim. vektorter´eb˝ ol v´alasztjuk ki valamilyen szempontb´ol a legmegfelel˝ obbet a tesztf¨ uggv´enyeket megfelel˝o integr´alokba helyettes´ıtve. Vari´aci´os m´odszer: A differenci´alegyenletet egy vari´aci´os feladat Euler–Lagrange-egyenletek´ent fogjuk fel, majd a kapott funkcion´alt minimaliz´aljuk a tesztf¨ uggv´enyek halmaz´an.

Integr´alm´odszerek S´ ulyozott reziduum m´ odszer: Az Lu? = f (u? a pontos megold´as) ´altal´anos alak´ u feladat eset´en a reziduum: r = f − Lu. ´Igy az Z w(x)r(x) dx = 0 Ω

egyenl˝os´eg minden w(x) u ´n. s´ ulyf¨ uggv´enyre fenn´all, ha a reziduumot a pontos megold´assal sz´am´ıtjuk. A s´ ulyozott reziduum m´odszerben a fenti egyenl˝ os´eget k¨ ovetelj¨ uk meg N adott s´ ulyra (wi (x)) a tesztf¨ uggv´enyekre. Olyan φ(x) tesztf¨ uggv´enyt keres¨ unk teh´at, amire Z Ii (φ(x)) = wi (x)r(φ(x)) dx = 0, i = 1, . . . , N. Ω

Vari´aci´os m´odszer

Walter Ritz (1878-1909) sv´ajci elm´eleti fizikus.  Z H 1 0 2 (φ (x)) + M0 φ(x)) dx I(φ(x)) = 2 0 φ(x) = y(x)-re minimum a V = {φ : C 2 [0, H], φ(0) = φ(H) = 0} halmazon m y(x) megold´asa a φ00 = M0 , y(0) = y(H) = 0 perem´ert´ekfeladatnak (Euler–Lagrange-egyenlet).

Vari´aci´os m´odszer Legyen a pr´obaf¨ uggv´enyek halmaza most Vh = {φ : [0, H] → R | φ(x) = A sin(πx/H), A ∈ R}. Minimaliz´aljuk az I funkcion´alt Vh -n. Legyen φA ∈ Vh . Ekkor Z H   πx 2  πx  1 π cos + M0 A sin I(φA (x)) = dx 2 H H H 0 =

π 2 2 2M0 H A + A. 4H π

Minim´alis, ha

4M0 H 2 , π3 azaz a  πx  4M0 H 2 φ(x) = − sin π3 H f¨ uggv´enyre. Maxim´alis lehajl´as A=−

−

4M0 H 2 = −0.129006M0 H 2 . π3

Kollok´aci´os m´odszer

N tetsz˝oleges pontban a tartom´anyon bel¨ ul legyen a reziduum nulla. Integr´alalakban is fel´ırhat´ o a wi (x) = δxi (x) Dirac-f´ele f¨ uggv´eny seg´ıts´eg´evel. (xi -ben 1, m´ashol nulla az ´ert´eke, ´es egy adott f¨ uggv´ennyel vett szorzat´anak integr´alja a f¨ uggv´eny xi -beli ´ert´eke.) Z Ii (φ(x)) = δxi (x)r(x) dx = r(xi ) = 0, i = 1, . . . , N. Ω

Kollok´aci´os m´odszer A p´eldabeli reziduum: r(x) = M0 − y 00 (x). A tesztf¨ uggv´enyeken az alakja  πx   π 2 sin . r(x) = M0 + A H H Az x = H/2 felez˝opontban akkor ad ez null´at, ha A=−

H 2 M0 . π2

Ekkor a k¨ozel´ıt˝o megold´as teh´at φ(x) = −

 πx  H 2 M0 sin , π2 H

Maxim´alis lehajl´as −0.1013M0 H 2 .

V´eges t´erfogat m´odszer (vagy r´esztartom´any m´odszer) A tartom´anyt diszjunkt r´esztartom´anyokra osztjuk: Ω = Ω1 ∪ . . . ∪ ΩN ´es a wi (x) = χΩi (x) (i = 1, . . . , N ) karakterisztikus f¨ uggv´enyeket v´alasztjuk s´ ulyoknak: Z Z Ii (φ(x)) = χΩi (x)r(x) dx = r(x) dx = 0, i = 1, . . . , N. Ω

Ωi

A p´eld´aban csak egy r´esztartom´any van: Z H Z H  πx   π 2 2Aπ sin dx = M0 H+ r(x) dx = M0 +A = 0. H H H 0 0 A ´ert´eke teh´at A = −M0 H 2 /(2π) = −0.1592M0 H 2 , ez a max. lehajl´as. K¨ozel´ıt˝o megold´as φ(x) = −

 πx  M0 H 2 sin . 2π H

Galjorkin-m´odszer

Borisz Grigorjevics Galjorkin (1871-1945) orosz m´ern¨ok / matematikus. Vh b´aP zisa legyen φ1 , . . . , φN , ´es a k¨ ozel´ıt˝ o megold´ast keress¨ uk φ = k ck φk (x) alakban. A s´ ulyf¨ uggv´enyeket v´alasszuk a b´azisf¨ uggv´enyeknek! Z X Ii ( ck φk (x)) = φi (x)r(x) dx = 0, i = 1, . . . , N. k

Ω

Ez a m´odszer a v´egeselem m´ odszer alapja.

Galjorkin-m´odszer A p´eld´aban: φ1 (x) = sin Z

 πx  H

.

H

φ1 (x)r(x) dx

I1 (Aφ1 (x)) = 0

Z = 0

H

 πx    π 2  πx  sin sin M0 + A dx H H H =

Aπ 2 2M0 H + = 0. 2H π

A ´ert´eke teh´at A = −4M0 H 2 /π 3 (mint a Ritz-m´ odszern´el), ez a max. lehajl´as. K¨ozel´ıt˝o megold´as φ(x) =

 πx  −4M0 H 2 sin . π3 H

Momentumok m´odszere

wi (x) = xi−1 , i = 1, 2, 3, . . . A p´eld´aban most wi (x) = 1 ´es egy intervallum van csak, ´ıgy ugyanazt kapjuk, mint a v´eges t´erfogat m´ odszer eset´en. A ´ert´eke teh´at A = −M0 H 2 /(2π) = −0.1592M0 H 2 , ez a max. lehajl´as. K¨ozel´ıt˝o megold´as φ(x) = −

 πx  M0 H 2 sin . 2π H

Legkisebb n´egyzetek m´odszere

Az

Z I(φ(x)) =

(r(x))2 dx → min.

Ω

feladat megold´as´at keress¨ uk. Ha Vh aP φi b´azisf¨ uggv´enyek line´aris kombin´aci´oja ´es a megold´ast φ(x) = k ck φk (x) alakban keress¨ uk, akkor a minimumhelyen a ck ismeretlenek szerinti deriv´altaknak null´anak kell lenni¨ uk. Inn´et kapjuk a Z ∂r(x) dx = 0, (i = 1, . . . , N ) 2 r(x) ∂ci Ω egyenl˝os´egeket. Azaz a reziduum deriv´altjai a s´ ulyok ebben a m´odszerben (wi (x) = ∂r/∂ci ).

Legkisebb n´egyzetek m´odszere A feladatban: Z I(φA (x)) =

H

 M0 + A

0

=

 π 2 H

sin

 πx 2 H

dx

A2 π 4 4M0 πA + + M02 H → min. 2H 3 H

Minimum

4M0 H 2 π3 eset´en. Ez a max. lehajl´as (mint a Ritz-m´ odszern´el), a k¨ozel´ıt˝o megold´as pedig A=−

φ(x) = −

 πx  4M0 H 2 sin . π3 H

Megjegyz´esek

I

¨ Onadjun´ alt oper´atorok eset´en a Ritz- ´es Galjorkin m´odszerek ugyanazt az eredm´enyt adj´ak.

I

V´egeselem m´odszernek h´ıvunk egy integr´alm´ odszert, ha a tesztf¨ uggv´enyek szakaszonk´ent folytonos f¨ uggv´enyek (a b´azis tesztf¨ uggv´enyek tart´ oja kicsi az eg´esz tartom´any m´eret´ehez ´ aban a Galjorkin-m´ k´epest). Altal´ odszert szokt´ak haszn´alni ebben az esetben.

I

Term´eszetesen eddig csak k¨ ozel´ıt˝ o megold´asokat ´all´ıtottunk el˝o, de nem foglalkoztunk azzal, hogy mekkora a k¨ozel´ıt˝o megold´as hib´aja. Ezt a k´erd´est majd k´es˝ obb a konkr´et feladatokn´al fogjuk vizsg´alni.

4. Az advekci´os egyenlet megold´asa (v´eges differencia)

Az advekci´os egyenlet 1D-ban Advekci´os egyenlet u0t (t, x) + cu0x (t, x) = f (t, x), kieg´esz´ıtve a megfelel˝o kezdeti- ´es peremfelt´etelekkel; c az ´araml´as sebess´ege. Egzakt megold´asa a karakterisztikus g¨ orb´ek seg´ıts´eg´evel t¨ort´enik. A megold´as grafikonja ¨ osszerakhat´ oa     t(s) 1 d  =  x(s) c ds u(t(s), x(s)) f (t(s), x(s)) KDER megold´asg¨orb´eib˝ ol ([t(0), x(0), u(t(0), x(0))]T adott vektor).

Karakterisztik´ak Legyen u(0, x) = u0 (x) adott kezdeti f¨ uggv´eny. Ha c = c(t, x) x-ben Lipschitz-folytonos ´es t-ben folytonos, akkor a karakterisztik´ak nem metszik egym´ast. Ha f ≡ 0, akkor u konstans a karakterisztik´akon. Ekkor, ha c = c(u), akkor egyenesek a karakterisztik´ak (u(t, x) = u0 (x − c(u(t, x))t)), c konstans eset´en p´arhuzamosak (u(t, x) = u0 (x − ct)).

Els˝o pr´ob´alkoz´as

Numerikus megold´asi m´ odszerek lehetnek: I

V´eges t´erfogat m´ odszer (nem vizsg´aljuk),

I

a KDER numerikus megold´asa (t¨ obb v´altoz´ oban komplik´alt),

I

v´eges differencia m´ odszer.

Vizsg´aljuk az al´abbi feladatot: Konstans sebess´eges advekci´ o, periodikus peremmel u0t (t, x) + cu0x (t, x) = 0, u0 (x) adott, u(t, 0) = u(t, 1).

Els˝o pr´ob´alkoz´as R´acsh´al´ot defini´alunk a megold´asi tartom´anyon (Ujk ≈ u(j∆t, k∆t)).

FTCS-s´ema: (k = 0, 1, . . . ; j = 1, . . . , n) Ujk+1 − Ujk ∆t

+c

k k Uj+1 − Uj−1 = 0. 2∆x

Els˝o pr´ob´alkoz´as k

A q = c∆t/∆x ´es Uh = [U1k , . . . , Unk ]T v´alaszt´assal a al´abbi vektoriter´aci´os alakot ¨ olti:  1 −q/2 q/2  q/2 1 −q/2  k+1  .. .. .. Uh =  . . .   q/2 1 −q/2 −q/2 q/2 1

s´ema az     k  Uh  

Pr´ob´aljuk ki u0 (x) = sin(2πx) kezdeti felt´etellel (c = 1)! → transeq(n,dt,kmax)=transeq(30,0.05,400) L´athat´o, hogy ”valami nem stimmel” a m´ odszerrel. Mi lehet a rossz viselked´es oka?

Lax ekvivalencia t´etele

Konvergencia: A r´acs finom´ıt´as´aval a numerikus megold´as tart a pontos megold´ashoz. Stabilit´as: Kicsit v´altoztatva a kezdeti felt´etelen a megold´as nem v´altozhat b´armekkor´at a r´acs finom´ıt´as´aval. Konzisztencia: A s´ema k´eplethib´aja null´ahoz tart, ha a r´acst´avols´agok null´ahoz tartanak. Lax ekvivalencia t´etele: Korrekt kit˝ uz´es˝ u feladatokra konstru´alt konzisztens numerikus s´em´ak pontosan akkor konvergensek, ha stabilak.

CFL-felt´etel

I I I

Richard Courant (1888-1972) n´emet sz´armaz´as´ u (amerikai) matematikus Kurt Otto Friedrichs (1901-1982) n´emet sz´armaz´as´ u (amerikai) matematikus Hans Lewy (1904-1988) n´emet sz´armaz´as´ u (amerikai) Wolf-d´ıjas matematikus

¨ Courant, R.; Friedrichs, K.; Lewy, H. (1928), ”Uber die partiellen Differenzengleichungen der mathematischen Physik” (in German), Mathematische Annalen 100 (1): 32-74.

CFL-felt´etel

Piros pontok: A P pontbeli numerikus megold´as f¨ ugg˝os´egi halmaza. K´ek vonal: A P pontbeli pontos megold´as f¨ ugg˝ os´egi halmaza. CFL-felt´etel: A numerikus s´ema, csak akkor lehet konvergens, ha numerikus f¨ ugg˝os´egi halmaza tartalmazza a pontos megold´as f¨ ugg˝os´egi halmaz´at (sz¨ uks´eges a konvergenci´ahoz).

CFL-felt´etel A mintap´eld´aban: A karakterisztika meredeks´ege: 1/c A numerikus f¨ ugg˝os´egi tartom´any oldalainak meredeks´ege ∆t/∆x. A CFL-felt´etel teljes¨ ul, ha c∆t ≤ 1. |q| = ∆x (q elnevez´ese: CFL-sz´am.) A p´eld´aban q = 1.5, azaz nem teljes¨ ul a felt´etel. → transeq(n,dt,kmax)=transeq(30,0.01,1200) Most m´ar teljes¨ ul, de m´eg mindig valami nem j´ o. Tudjuk, hogy a CFL-felt´etel csak sz¨ uks´eges felt´etel a konvergenci´ahoz.

Diszkr´et Fourier-anal´ızis A numerikus s´em´at kiel´eg´ıtik az al´abbi komplex sz´amok: Ujk = g k eijl∆x , ha g ≡ g(l∆x) = 1 − qi sin(l∆x) (g: n¨oveked´esi faktor, amplification factor). Vegy¨ uk ´eszre, hogy |g|2 = 1 + q 2 sin2 (l∆x), amely ´altal´aban nagyobb lesz, mint 1, azaz ´altal´aban a Fourier-¨osszetev˝ok amplit´ uld´ oja n¨ ovekszik. Azaz r¨ ogz´ıtett q mellett finom´ıtott r´acsh´al´on a sz´am´ıt´asokra rak´ od´ o hiba n¨ovekszik, ´ıgy nem lehet a m´odszer stabil. A FTCS-s´ema nem haszn´alhat´ o, mert nem stabil. A CFL-felt´etelt viszont kiel´eg´ıti.

Upwind s´ema Egy egyszer˝ u ´es konvergens s´ema (el˝ oz˝ o id˝ or´etegen val´o line´aris interpol´aci´ob´ol):

Ujk+1

( k k − Ujk ) = (1 + q)Ujk − qUj+1 Ujk − q(Uj+1 = k ) = (1 − q)U k + qU k Ujk − q(Ujk − Uj−1 j−1 j

, c < 0, , c ≥ 0.

CFL-felt´etel: teljes¨ ul, ha |q| ≤ 1. Diszkr´et Fourier: g ≡ g(l∆x) = 1 − |q|(1 − e−il∆x ), azaz |g|2 = 1 − 4|q|(1 − |q|) sin2 (l∆x/2). Azaz a CFL-felt´etel mellett stabil is a s´ema (most ez el´egs´eges is), azaz konvergens is (konzisztencia miatt).

Az upwind s´ema konvergenci´aja (c > 0 esete) Legyen a pontos megold´as jel¨ ol´ese: ukj = u(k∆t, j∆x) (hasonl´oan a deriv´altakra is). Tegy¨ uk fel, hogy a megold´as m´asodrend˝ u deriv´altjai folytonosak, azaz ebben az esetben korl´atosak is a megold´asi tartom´anyon: [0, tmax ] × [0, 1]. K´eplethiba: uk+1 − ukj j

ukj − ukj−1 +c ∆t ∆x    0 k 1    00 k 1 00 k? 0 k = ut j + ∆t[ utt ]j + c ux j − ∆x uxx j ? , |{z} 2 2 Tjk :=

c2 u00 xx

teh´at

1 |Tjk | ≤ c∆x(1 − q)Mxx , 2 ahol Mxx fels˝o becsl´es a megold´as x szerinti m´asodik deriv´altj´ara a tartom´anyon.

Az upwind s´ema konvergenci´aja (c > 0 esete) Legyen a hiba jel¨ol´ese: ekj = ukj − Ujk . Ekkor a numerikus s´ema ´es a k´eplethiba k´eplet´enek kivon´as´aval kapjuk, hogy ek+1 − ekj j ∆t

+c

ekj − ekj−1 = −Tjk , ∆x k

azaz az ekh = [ek1 , . . . , ekn ]T ´es Th = [T1k , . . . , Tnk ]T jel¨ol´esekkel   1−q q  q  1−q   k   k .. .. ek+1 =   eh − ∆tTh . . h     q 1−q

Az upwind s´ema konvergenci´aja (c > 0 esete)

Az iter´aci´os m´atrix maximumnorm´aja 1, ´ıgy ad´ odik a 0

k−1

kekh k∞ ≤ ∆t(kTh k∞ + . . . + kTh

k∞ )

1 1 ≤ ∆tk c∆x(1 − q)Mxx ≤ tmax c∆x(1 − q)Mxx 2 2 becsl´es, ami mutatja a s´ema els˝ orend˝ u konvergenci´aj´at. → transeq upwind(50,0.01,800) Megj.: A bizony´ıt´ashoz felhaszn´altuk a megold´as kell˝o simas´ag´at. Ha a kezdeti felt´etel nem el´eg sima, akkor ez nem teljes¨ ul. Ilyenkor m´as m´odszerek kellenek a konvergencia vizsg´alat´ahoz.

Az upwind s´ema Fourier-anal´ızise

λ: hull´amhossz, λ = 2π/l (l: hull´amsz´am) f : frekvencia, ω = 2πf : k¨ orfrekvencia Hull´am sebess´ege: v = λf = ω/l ´Igy egy szinuszhull´am fel´ırhat´ o az u(t, x) = A0 ei(lx−ωt) alakban. Ez a hull´am kiel´eg´ıti az advekci´ os egyenletet, ha c = ω/l (halad´asi sebess´ege ´eppen c). Amplit´ ud´ o ´alland´ o. ∆t alatt a f´azisv´altoz´as −ω∆t = −cl∆t = −ql∆x.

Az upwind s´ema Fourier-anal´ızise

A numerikus s´ema f´aziseltol´ od´asa: Ujk = g k eijl∆x = eln(|g|)k ei(jl∆x+arg(g)k) L´athat´o, hogy |g| = 6 1 eset´en v´altozik az amplit´ ud´ o ´es a f´azisv´altoz´as ∆t id˝o alatt arg(g). Lemma. arctg(c1 p + c2 p2 + c3 p3 + c4 p4 + . . .) = c1 p + c2 p2 + (c3 − c31 /3)p3 + (c4 − c21 c2 )p4 + . . . .

Az upwind s´ema Fourier-anal´ızise Az upwind s´ema eset´en: g = 1 − |q|(1 − e−il∆x ), melyre q |g| = 1 − 4|q|(1 − |q|) sin2 (l∆x/2). √ A 1 − x = 1 − x/2 − x2 /8 − x3 /16 − . . . sorfejt´es miatt ez azt mutatja, hogy az amplit´ ud´ o hib´aja (l∆x)2 -tel ar´anyos. Mivel g = (1 − q) + q cos(l∆x) − iq sin(l∆x), emiatt az argumentum arg(g) = −arctg

q sin(l∆x) (1 − q) + q cos(l∆x)

= −ql∆x(1 − (1 − q)(1 − 2q)(l∆x)2 /6 + . . .).

Lax–Wendroff-s´ema

Lax P´eter (1926, Budapest), magyar sz´armaz´as´ u amerikai matematikus. Burton Wendroff (1930) amerikai matematikus.

Lax–Wendroff-s´ema Kvadratikus interpol´aci´ o seg´ıts´eg´evel sz´armaztathat´o. 1 1 k k Ujk+1 = q(1 + q)Uj−1 + (1 − q 2 )Ujk − q(1 − q)Uj+1 . 2 2 → transeq laxwendroff(50,0.01,1000) N¨oveked´esi faktor: 2

2

g(l∆x) = 1 − iq sin(l∆x) − 2q sin ahonn´et 2

2

2

4

|g| = 1 − 4q (1 − q ) sin





 1 l∆x , 2

1 l∆x 2



(CFL el´egs´eges is), arg(g) = −ql∆x(1 − (1 − q 2 )l2 ∆x2 /6 + . . .).

Leap-frog s´ema

Ujk+1 − Ujk−1

k k Uj+1 − Uj−1 =0 2∆t 2∆x (ind´ıt´as pl. a Lax-Wendroff s´em´aval).

+c

N¨oveked´esi faktor: q g(l∆x) = −iq sin(l∆x) ± 1 − q 2 sin2 (l∆x), ha a CFL-felt´etel teljes¨ ul, akkor nincs csillap´ıt´as. arg(g) = −ql∆x(1 − (1 − q 2 )l2 ∆x2 /6 + . . .).

5. Elliptikus egyenletek numerikus megold´asa (v´eges differencia)

A Poisson-egyenlet t´eglalapon Poisson-egyenlet 2D-ban, t´eglalapon, Dirichlet-peremmel ∂2u ∂2u − 2 = f (x, y), (x, y) ∈ Ω := (0, a) × (0, b) ∂x2 ∂y u|Γ = g(x, y), Γ := ∂Ω.

−∆u = −

V´eges differencia megold´ashoz r´acsot defini´alunk: Ωh := {(x, y) ∈ Ω | x = ih1 , y = jh2 , | {z } xij

i = 0, . . . , n1 , j = 0, . . . , n2 }, Ωh ∩ Ω =: Ωh , Ωh ∩ Γ =: Γh .

A v´eges differencia s´ema

Alkalmazzuk az Uij ≈ u(ih1 , jh2 ), fij = f (ih1 , jh2 ), gij = g(ih1 , jh2 ) k¨ozel´ıt´eseket, melyekkel a s´ema alakja a k¨ ovetkez˝ o lesz.

A v´eges differencia s´ema V´eges differencia s´ema, Poisson-egyenlet, t´eglalap, Dirichlet-perem

−

j j Ui−1 − 2Uij + Ui+1 Uij−1 − 2Uij + Uij+1 − = fij , xij ∈ Ωh h21 h22

Uij = gij , xij ∈ Γh .

A fenti egyenletek egy (n1 + 1)(n2 + 1) m´eret˝ u line´aris egyenletrendszert eredm´enyeznek Ah Uh = φh , ahol ( fij , (φh )ij = gij ,

xij ∈ Ωh , , kφh k∞ ≤ max{kgkC(Γ) , kf kC(Ω) } xij ∈ Γh ,

´es a m´atrix ill. vektorok alakja att´ ol f¨ ugg, hogy hogy rendezt¨ uk a r´acspontokat. A tov´abbiakban a sorfolytonos sorsz´amoz´ast alkalmazzuk (lentr˝ol felfel´e, balr´ ol jobbra).

Az Ah m´atrix szerkezete k´ek: 1 z¨old: 2 2 + 2 2 h1 h2 piros: −

1 h21

narancs: −

1 h22

Az Ah m´atrix szerkezete T´etel. Az Ah m´atrix M-m´atrix. Biz. El˝ojelek stimmelnek. Tekints¨ uk a w(x, y) = 4 + x(a − x) + y(b − y) f¨ uggv´enyt ´es defini´aljuk a w, wij = w(ih1 , jh2 ) pozit´ıv vektort, amit ugyan´ ugy rendez¨ unk, mint az Uh vektor elemeit. Erre a vektorra ( 4, xij ∈ Ωh , (Ah w)ij = ≥ 4, xij ∈ Γh , ami pozit´ıv vektor.

Az Ah m´atrix szerkezete

Megj. Ha f1 ≥ f2 ´es g1 ≥ g2 , akkor az ezekkel el˝ o´all´ıtott k´et numerikus megold´asra ´erv´enyes (Uh )1 ≥ (Uh )2 . Speci´alisan ebb˝ol k¨ovetkezik a nemnegativit´as-meg˝orz´esi tulajdons´ag.

A numerikus s´ema konvergenci´aja Vezess¨ uk be az uji jel¨ol´est az u(ih1 , jh2 ) pontos ´ert´ekekre, tov´abb´a legyen E = maxi=0,...,n1 ; j=0,...,n2 |Uij − uji |. Def. A numerikus s´em´at konvergensnek h´ıvjuk, ha h := max{h1 , h2 } → 0 eset´en E → 0. Ha E = O(hr ), akkor azt mondjuk, hogy a konvergencia rendje r. Def. A numerikus s´ema stabil, ha kA−1 ulr˝ ol korl´atos a r´acs h k fel¨ finoms´ag´at´ol f¨ uggetlen¨ ul. Def. A numerikus s´ema konzisztens, ha k´eplethib´aja minden pontban null´ahoz tart, ha h → 0. A konzisztenciarend r, ha a u. k´eplethiba minden pontban O(hr ) nagys´agrend˝ T´etel. Konzisztencia + stabilit´as ⇒ konvergencia.

A numerikus s´ema konvergenci´aja

T´etel. A numerikus s´ema stabil, ugyanis kA−1 h k∞ ≤ 1 +

a2 + b2 . 16

Biz. Az M-m´atrixok inverz´enek maximumnormabeli becsl´es´eb˝ol k¨ovetkezik, hogy kA−1 h k∞ ≤

4 + a2 /4 + b2 /4 a2 + b2 =1+ . 4 16

A numerikus s´ema konvergenci´aja T´etel. Amennyiben a vizsg´alt perem´ert´ekfeladat megold´asa C 4 (Ω)-beli, akkor a vizsg´alt numerikus s´ema m´asodrendben konzisztens. Biz. Tij := −

uji−1 − 2uji + uji+1 uj−1 − 2uji + uj+1 i i − − fij h21 h22 ?

j ∂ 2 u ∂y4 ui h22 ∂ 2 u ∂ 4 uj? h2 − fij , =− 2 − x i 1 − 2 − ∂x 12 ∂y 12

azaz |Tij | ≤

M4,0 h21 M0,4 h22 M h2 + ≤ 12 12 6

a r´acspontt´ol f¨ uggetlen¨ ul. (Az M ´ert´ekek a szokott m´odon a deriv´altak fels˝o becsl´esei.) Ez mutatja a m´asodrend˝ u konzisztenci´at.

A numerikus s´ema konvergenci´aja T´etel. Amennyiben a vizsg´alt perem´ert´ekfeladat megold´asa C 4 (Ω)-beli, akkor a vizsg´alt numerikus s´ema m´asodrendben konvergens. Biz. Legyen a k´eplethib´ak vektora megfelel˝ oen rendezve Th . Ekkor az Ah Uh − φh = 0 ´es Ah uh − φh = Th kivon´as´ab´ol ad´odik, hogy Ah (Uh − uh ) = −Th , melyb˝ol ...

A numerikus s´ema konvergenci´aja

... kUh − uh k∞ ≤ kA−1 h k∞ kTh k∞ ´es a kor´abbi becsl´eseket haszn´alva   a2 + b2 M h2 kUh − uh k∞ ≤ 1 + , 16 6 amely mutatja a m´asodrend˝ u konvergenci´at.

Maximum-elv A folytonos feladatra: Ha f ≤ 0, akkor a megold´as a peremen veszi fel a maximum´at. Igaz-e ez a numerikus megold´asra? Rendezz¨ uk ´at az Ah Uh = φh egyenletrendszer sorait ´es oszlopait u ´gy, hogy el˝ ore ker¨ uljenek Ω bels˝o pontjai, ´es ut´ana k¨ ovetkezzenek a hat´arpontok. Ekkor a line´aris egyenletrendszer az al´abbi alak´ u lesz:      A0 A∂ U0 f = 0 E g U∂ T´etel. A numerikus s´ema kiel´eg´ıti a maximum-elv diszkr´et v´altozat´at, azaz, ha f ≤ 0, akkor max(U0 ) ≤ max(g).

Maximum-elv

Vegy¨ uk ´eszre, hogy A0 M-m´atrix (A−1 es A∂ elemei 0 ≥ 0) ´ nempozit´ıvak, valamint [A0 A∂ ][1, . . . , 1]T = 0, [1, . . . , 1]T =: [eT0 , eT∂ ]T , azaz e0 = −A−1 0 A∂ e∂ . Emiatt −1 −1 −1 U0 = A−1 0 f − A0 A∂ U∂ = A0 f − A0 A∂ g −1 ≤ A−1 0 f − A0 A∂ e∂ max(g) ≤ e0 max(g).

Poisson-egyenlet Robin-peremmel Poisson-egyenlet 2D-ban, t´eglalapon, Robin-peremmel −∆u = −

∂2u ∂2u − 2 = f (x, y), (x, y) ∈ Ω := (0, a) × (0, b) ∂x2 ∂y

∂u = σ(x, y)(u0 (x, y) − u), (x, y) ∈ Γ. ∂n Tegy¨ uk fel, hogy f, σ, u0 olyan f¨ uggv´enyek, melyekkel az 4 egyenletnek van C (Ω)-beli megold´asa. A v´eges differenci´as megold´as diszkretiz´aci´ oja nyilv´anval´oan a bels˝o pontokban ugyanaz, mint a Dirichlet-esetben, azaz csak a tartom´any hat´arpontjaiban kellene megadni (lehet˝oleg m´asodfok´ u) approxim´aci´okat.

Approxim´aci´o az oldal´elekn´el

ujn1 − ujn1 −1 h1 h2 = (∂1 u)jn1 − (∂12 u)jn1 + 1 (∂13 u)jn? , 1 h1 2 6 ahol (∂1 u)jn1 = σn1 ,j ((u0 )jn1 − ujn1 ) a peremfelt´etel miatt, ´es (∂12 u)jn1 = −(∂22 u)jn1 − fn1 ,j ´es −(∂22 u)jn1 = −

j j−1 h22 4 j ? uj+1 n1 − 2un1 + un1 − (∂ u) 12 2 n1 h22

a differenci´alegyenletb˝ol.

Approxim´aci´o az oldal´elekn´el

´Igy a s´ema Unj 1 − Unj 1 −1 = σn1 ,j ((u0 )jn1 − Unj 1 ) h1 h1 − 2

Unj+1 − 2Unj 1 + Unj−1 1 − 1 − fn1 ,j 2 h2

Hasonl´oan j´arunk el a t¨ obbi oldal´an a tartom´anynak. A k´eplethiba O(h21 + h22 ).

! .

Approxim´aci´o a sarkokban

h2 h2 un0 2 − un0 2 −1 −2 n? = (∂2 u)0n2 − (∂22 u)n0 2 + 2 (∂23 u)0 2 /· , h2 2 6 h2 2 un1 2 − un0 2 h2 h1 = (∂1 u)n0 2 + (∂12 u)n0 2 + 1 (∂13 u)n0?2 /· . h1 2 6 h1 A k´et egyenletet ¨osszeadva ´es kihaszn´alva a peremfelt´etelt ´es a differenci´alegyenletet, az al´abbi s´em´at kapjuk (H = h1 h2 /2/(h1 + h2 )): H(−2

U0n2 − U0n2 −1 U1n2 − U0n2 +2 ) = −σ0,n2 ((u0 )n0 2 −U0n2 )−Hf0,n2 h22 h21

Approxim´aci´o a sarkokban A t¨obbi sarokban is ´ıgy j´arunk el. A k´eplethiba most is O(h21 + h22 ) nagys´agrend˝ u. T´etel. Ha σ(x, y) ≥ σ0 > 0 ´es az f, σ, u0 f¨ uggv´enyek el´eg sim´ak 4 ahhoz, hogy a feladat u megold´asa C (Ω)-beli legyen, akkor a differencias´em´ank m´asodrendben konverg´al a pontos megold´ashoz. Biz. A konzisztencia az approxim´aci´ os k´epletekb˝ ol l´athat´o. A stabilit´as a kor´abbi M-m´atrixos technik´aval igazolhat´o az  2  a + b2 4 + a + b + kUh k∞ ≤ kφh k 16 4σ0 alakban.

Nem t´eglalap alak´u tartom´anyok

Mit lehet tenni, ha nem t´eglalap alak´ u az Ω tartom´any? I

A perem´ert´ekeket a r´acspontokban k¨ ozel´ıtj¨ uk (elveszik a m´asodrend˝ u konvergencia).

Nem t´eglalap alak´u tartom´anyok

I

S˝ ur´ıtj¨ uk a r´acsot.

I

Hozz´avessz¨ uk r´acspontnak a r´acs ´es a perem metsz´espontjait (Shortly–Weller-approxim´aci´ o).

I

T´eglalapba transzform´aljuk a tartom´anyt.

Shortly–Weller-approxim´aci´o r´eszletesebben A bels˝o pontokban ugyanaz a k¨ ozel´ıt´es, mint kor´abban, a peremen pedig ! j j Ui+1 − Uij Uij − Ui−1 2 − + h− h+ h− 1 + h1 1 1 ! Uij − Uij−1 Uij+1 − Uij 2 + − − = −fij . h2 + h+ h+ h− 2 2 2 Ezen egyenletek egy egyenletrendszerre vezetnek: AS,h Uh = φh , ´erv´enyes a 1 3 kA−1 S,h k∞ ≤ 1 + 4 (diam(Ω)) stabilit´asi becsl´es. Ha u ∈ C 4 (Ω), akkor a m´ odszer m´asodrendben konvergens lesz.

6. T¨obbr´acsos m´odszerek

A megoldand´o line´aris egyenletrendszer 

A0 A∂ 0 E



U0 U∂



 =

f g



Olyan egyenletrendszert kell tudnunk megoldani, melynek egy¨ utthat´om´atrixa A0 . Ez egy szimmetrikus M -m´atrix. T´etel. A szimmetrikus M-m´atrixok pozit´ıv definitek. Biz. Jel¨olje a m´atrixot M. Szimmetrikus m´atrix minden saj´at´ert´eke val´os. Azt kell megmutatni, hogy mindegyik pozit´ıv. A f˝ o´atl´oban ´all´o elemek pozit´ıvak. Legyen ezek maximuma µ, ´es ´ırjuk fel M-et M = µE − H alakban, ahol H nyilv´anval´oan nemnegat´ıv m´atrix. ´Igy a fenti fel´ır´as egy nemnegat´ıv inverz˝ u m´atrix (az M m´atrixok ilyenek) egy regul´aris felbont´asa. Azaz %((1/µ)H) < 1, ami azt jelenti, hogy %(H) < µ. Ha H saj´at´ert´ekeit λk jel¨oli, akkor M µ − λk alak´ u saj´at´ert´ekei nyilv´an pozit´ıvak.

Szimmetrikus, pozit´ıv definit m´atrix´u line´aris egyenletrendszerek megold´asa Az egyenletrendszer azon t´ ul, hogy szimmetrikus, pozit´ıv definit m´atrix´ u, m´eg ritka m´atrix´ u is, ´ıgy j´ ol m˝ uk¨ odnek az iter´aci´os m´odszerek is r´a. Lehets´eges megold´asi m´ odszerek: I

Gauss-m´odszer (2N03 /3, nem praktikus)

I

Cholesky-felbont´as (N03 /3)

I

Jacobi- ´es Gauss-Seidel ill. relax´alt v´altozatai (N02 ln(1/ε), ε: el´erni k´ıv´ant pontoss´ag a le´all´asi felt´etelben)

I

Konjug´alt gradiens m´ odszer (N0

I

T¨obbr´acsos m´odszer. 1960-as ´evek, Fedorenko (1962), Bakvalov (1966). Haszn´aljuk ki, hogy az egyenletrendszer egy elliptikus feladat diszkretiz´aci´ oj´ab´ ol keletkezett!

3/2

ln(1/ε))

Iter´aci´ok sim´ıt´o tulajdons´aga

Iter´aci´ok sim´ıt´o tulajdons´aga Egydimenzi´os perem´ert´ekfeladat −u00 (x) = f (x), x ∈ (0, 1), u(0) = g1 , u(1) = g2 A szok´asos v´eges differenci´as k¨ ozel´ıt´es sor´an megoldand´o line´aris egyenletrendszer: 

−1   −1 . . . . . .   ..  .   2

        −1  −1 2

U1 U2 .. . Un





      = h2    

vagy t¨om¨orebben Ah Uh = φh .

f1 + g1 /h2 f2 .. . fn−1 fn + g2 /h2

      

Iter´aci´ok sim´ıt´o tulajdons´aga Oldjuk meg az Ah Uh = φh egyenletrendszert a relax´alt Jacobi-iter´aci´oval. Ekkor az (k+1)

Uh

(k)

= ((1 − ω)E + ωD−1 (L + R))Uh + ωD−1 φh (k)

= (E − ωD−1 (D L − R})) Uh + ωD−1 φh | − {z Ah

|

{z

BJ (ω)

}

iter´aci´ot kapjuk. A konvergenci´ahoz a %(BJ (ω)) < 1 felt´etelnek kell teljes¨ ulni. Ah saj´atrendszere: vl = sin(lπhj), j = [1, . . . , n]T , λl (Ah ) = 2−2 cos(lπh) l = 1, . . . , n. ´Igy BJ (ω) saj´atrendszere: vl , λl (BJ (ω)) = 1 − ω(1 − cos(lπh)), azaz, ha ω ∈ (0, 1], akkor a JOR-m´ odszer konvergens lesz.

Iter´aci´ok sim´ıt´o tulajdons´aga Megj. Vegy¨ uk ´eszre, hogy ha pl. ω = 1, akkor %(BJ (ω)) = cos(πh) → 1, ha h → 0, azaz egyre lassabban konverg´al a Jacobi-m´odszer. A t¨obbr´acsos m´odszerek eset´en nem a hat´ar´ert´ek el˝o´all´ıt´asa a c´el, hanem a k¨ozel´ıt˝o megold´as sim´ıt´asa az´ert, hogy az j´ol k¨ozel´ıthet˝o legyen egy durv´abb r´acson. tobbracsos(f,n), f=16,1, n=5,10,30,60,160 tobbracsosJOR(omega,n,kmax), (omega=0.1,2/3,1, n=10,100, kmax=1500) A magas frekvenci´aj´ u saj´atvektorokat kellene gyorsan cs¨okkenteni. Alacsony frekvenci´as ¨osszetev˝ ok (LF): 1 ≤ l < n/2. Magas frekvenci´as ¨osszetev˝ ok (HF): n/2 ≤ l ≤ n. Def. µ legyen az a legnagyobb t´enyez˝ o, amennnyivel a HF o¨sszetev˝ok cs¨okkennek egy iter´aci´ oban. µ-t sim´ıt´asi t´enyez˝onek nevezz¨ uk.

Iter´aci´ok sim´ıt´o tulajdons´aga

µ = max{|λk (BJ (ω))| | n/2 ≤ l ≤ n} = max{|1 − ω|, |1 − 2ω|} Ez akkor lesz minimum, ha ω = 2/3, ´es a minimum ´ert´eke µ = 1/3. Azaz a HF ¨osszetev˝ ok 2 iter´aci´ o alatt a tized¨ ukre cs¨okkennek, annak ellen´ere, hogy a konvergencia lass´ u lehet a finom r´acson. Ha a HF ¨ osszetev˝ ok kicsik m´ar, akkor a hibavektor j´ ol k¨ozel´ıthet˝o egy durv´abb r´acson is.

A k´etr´acsos m´odszer

A marad´ekegyenlet T´erj¨ unk vissza a k´etv´altoz´ os Poisson-egyenletb˝ ol sz´armaztatott Ah Uh = φh egyenletrendszer megold´as´ahoz! V´egezz¨ unk el n´eh´any iter´aci´os l´ep´est valamilyen sim´ıt´ o iter´aci´ oval valamilyen ˜ h vektort. kezd˝ovektorr´ol indulva. Kapjuk ekkor az U ˜ h hibavektor egy sima vektor lesz, ´es Ekkor az eh = Uh − U ´erv´enyes r´a az ˜ h ) = φ − Ah U ˜ h =: rh (marad´ekvektor), Ah eh = Ah (Uh − U h egyenl˝os´eg (´ un. marad´ekegyenlet), ami megoldhat´ o eh -ra. De mivel a jobb oldala ´es a megold´asa is sima az egyenletrendszernek, j´ o k¨ozel´ıt´est kapn´ank akkor is, ha egy durv´abb r´acson (H = 2h) oldan´ank meg egy ugyanilyen egyenletrendszert.

Restrikci´o T´erj¨ unk ´at a durv´abb r´acsra: rH = RhH rh . Ezt egy u ´n. restrikci´os oper´atorral val´os´ıthatjuk meg: RhH . Ez ´altal´aban valamilyen s´ ulyozott vet´ıt´es.

Teljes ´es r´eszleges s´ ulyoz´as:    0 1 0 1 2 1 1 1  1 4 1  2 4 2  8 16 0 1 0 1 2 1 

Megoldjuk az AH eH = rH egyenletet, ahol AH a differenci´aloper´ator k¨ ozel´ıt´ese a durv´abb r´acson.

Prolong´aci´o eH hibavektort visszavet´ıtj¨ uk a finomabb r´acsra. Ezt egy h prolong´aci´os PH oper´atorral tessz¨ uk meg, amely ´altal´aban valamilyen interpol´aci´ot jelent: eh = PhH eH .

Course Grid Correction (CGC) Adunk egy jobb k¨ozel´ıt´est Uh -ra: ˜ h + Ph eH . unew =U H h Megj.: Az interpol´aci´o miatt megjelenhetnek HF o¨sszetev˝o, amiket egy u ´jabb sim´ıt´asi ciklussal lehet elt¨ untetni. A bemutatott elj´ar´as egy V-alak´ u s´em´aban foglalhat´o ¨ossze (V-ciklus).

V-ciklus Megj.: Term´eszetesen a V-ciklusb´ ol t¨ obbet kell v´egrehajtani ahhoz, hogy a pontos megold´ast kell˝ oen megk¨ ozel´ıts¨ uk. Lemma. A V-ciklus iter´aci´ os m´atrixa (Sh a sim´ıt´as): k1 H Skh2 (E − PhH A−1 H Rh Ah )Sh ,

ezen m´atrix spektr´alsugara kisebb 1-n´el, azaz konvergens. Lemma. Csak ¨onmag´aban a CGC nem konvergens. Biz.: M´atrixa H E − PhH A−1 H Rh Ah , H ahol PhH A−1 alhat´ o, ´ıgy egyik saj´at´ert´eke nulla. H Rh Ah nem invert´ ´Igy a m´atrixnak 1 saj´at´ert´eke lesz, azaz spektr´alsugara legal´abb 1.

A t¨obbr´acsos m´odszer (multigrid)

A t¨obbr´acsos m´odszer alapgondolata

¨ Otlet: Az egyenlet durva r´acson val´ o megold´as´ahoz haszn´aljuk egy u ´jabb V-ciklust egy m´eg durv´abb r´accsal. T¨ obb r´acsot defini´alunk: h1 > h2 = h1 /2 > h3 = h2 /2 > . . . hL > 0 (L szintek sz´ama). Legyne a k-adik szinten nk a feladat m´erete. A k-adik szinten az Ak Uk = φk egyenletrendszert kell megoldanunk. A restrikci´os oper´ator legyen Rk−1 ´es a prolong´aci´ o Pkk−1 . k

A t¨obbr´acsos m´odszer algoritmusa

Oldjuk meg az Ak Uk = φk egyenletrendszert. 1. Ha k = 1, akkor oldjuk meg k¨ ozvetlen¨ ul az egyenletrendszert. ˜ 2. El˝osim´ıt´as, k1 darab l´ep´es → Uk . ˜ k = rk . 3. Marad´ek kisz´am´ıt´asa: φk − Ak U 4. A marad´ek restrikci´ oja: rk−1 = Rk−1 k rk . 5. Legyen Uk−1 = 0. 6. H´ıvjuk meg γ-szor a t¨ obbr´acsos m´ odszert az Ak−1 Uk−1 = rk−1 egyenletre. ˜ k + Pk Uk−1 . 7. CGC: Uk = U k−1

8. Ut´osim´ıt´as k2 -szer.

A t¨obbr´acsos m´odszer algoritmusa Def.: γ a ciklusindex, azt adja meg, hogy a durva r´acson h´anyszor ´ aban γ = 1 (V-ciklus) haszn´aljuk a t¨obbr´acsos algoritmust. Altal´ vagy γ = 2 (W-ciklus). Lemma. A t¨obbr´acsos m´ odszer iter´aci´ os m´atrixa k = 1 eset´en 0, k¨ ul¨onben k = 2, . . . , L eset´en k1 k−1 Mk = Skk2 (E − Pkk−1 (E − Mγk−1 )A−1 k−1 Rk Ak )Sk .

K´etr´acsos m´odszer: u ¨res k¨ or: pontos megold´as, teli k¨or: sim´ıt´as, marad´eksz´am´ıt´as, CGC

A t¨obbr´acsos m´odszer algoritmusa H´arom r´acs esete:

N´egy r´acs esete:

M˝uveletsz´amok vizsg´alata 2D Poisson-egyenlet, max. 5 nemnulla elem soronk´ent, sim´ıt´ashoz line´aris vektoriter´aci´o A k. szint m˝ uveletsz´ama: 1. 10k1 nk : el˝osim´ıt´as 2. 10nk : marad´eksz´am´ıt´as 3. 10nk /4: teljes restrikci´ o 4. 6nk /4: biline´aris interpol´aci´ o 5. 10k2 nk : ut´osim´ıt´as ¯k. Ez ¨osszesen k1 + k2 = 2 eset´en 34nk = Q ¯ Q1 a legdurv´abb r´acson az egyenlet pontos megold´as´anak m˝ uveletsz´ama.

M˝uveletsz´amok vizsg´alata ¯2 + Q ¯ 1, Qk legyen a k. r´acs m˝ uveletsz´ama. Ekkor Q2 = Q ¯2 + Q ¯ 1) + Q ¯ 3 , stb. Q3 = γ(Q Qk =

=

k−2 X l=0 k−2 X

¯ k−l + γ k−2 Q ¯1 γlQ ¯1 γ l 34nk−l + γ k−2 Q

l=0

=

k−2 X

γ l 34

l=0

nk 4l |{z}

¯1 +γ k−2 Q

2D probl´ ema

= 34nk

k−2   X γ l l=0

≈ 34nk

4

¯1 + γ k−2 Q

1 = O(nk ) 1 − γ/4

Megjegyz´esek

Megj.: γ = 1 eset´en Qk ≈ 40nk , γ = 2 eset´en Qk ≈ 68nk ´es γ = 3 eset´en Qk ≈ 136nk Megj.: A 2D Poisson-feladatra csak a γ ≤ 3 esetek praktikusak Megj.: K´etr´acsos m´odszer eset´en, vagy t¨ obb r´acsra γ = 1, 2 eset´en a konvergencia t´enyez˝o µk1 +k2 a h r´acst´avols´agt´ ol f¨ uggetlen¨ ul.

7. M´asodrend˝u, id˝of¨ugg˝o, line´aris PDE-ek megold´asa a v´eges differencia m´odszerrel

A hull´amegyenlet Az egydimenzi´os hull´amegyenlet 2 ∂2u 2∂ u = c , ∂t2 ∂x2 u(0, x) = u0 (x), ∂u (0, x) = u1 (x). ∂t

A pontos megold´as 1 1 u(t, x) = (u0 (x − ct) + u0 (x + ct)) + 2 2c

Z

x+ct

u1 (y) dy. x−ct

Gyakorlaton szerepelt: I

V´eges differenci´as megold´as, a kezdeti felt´etelek m´asodrend˝ u approxim´aci´oja. (CFL felt´etel q ≤ 1.)

I

A leap-frog s´ema.

A h˝ovezet´esi vagy diff´uzi´os egyenlet

A h˝ovezet´esi vagy diff´ uzi´ os egyenlet ∂u = ∆u + f ∂t + kezdeti- ´es peremfelt´etel 1D-ban ∂u ∂2u = , x ∈ (0, 1) ∂t ∂x2 u(0, x) = u0 (x), u(t, 0) = u(t, 1) = 0. M´as egyenletek, peremfelt´etelek → gyakorlaton

A h˝ovezet´esi vagy diff´uzi´os egyenlet

A pontos megold´as el˝o´all´ıthat´ o Fourier-sor alakban: u(t, x) =

∞ X

al e−l

2 π2 t

sin(lπx),

l=1

ahol al = 2

R1 0

u0 (x) sin(lπx) dx.

Probl´em´ak: I

al meghat´aroz´as´ahoz ´altal´aban numerikus integr´al´as kell.

I

V´egtelen sort kell ¨ osszegezni.

Az egyenesek m´odszere R´acsh´al´ot defini´alunk a [0,1] intervallumon: 0 = x0 < x1 < · · · < xn+1 = 1,

xi+1 − xi = ∆x.

Vezess¨ uk be az al´abbi f¨ uggv´enyeket uj (t) ≈ u(xj , t). Ezen f¨ uggv´enyekre k¨ozel´ıts¨ uk a helyszerinti deriv´altat v´eges differenci´aval: uj−1 (t) − 2uj (t) + uj+1 (t) duj (t) = . dt ∆x2 Ez egy k¨oz¨ons´eges differenci´alegyenlet rendszer, a kezdeti felt´etel uj (0) = u0 (xj ) ´es u0 (t) = 0, un+1 (t) = 0. B´armely kor´abban tanult k¨oz¨ons´eges differenci´alegyenlet megold´ o m´ odszerrel megoldhatjuk.

Explicit Euler-m´odszer (EE) Oldjuk meg pl. az egyenletrendszert az explicit Euler m´odszerrel. Jel¨olje Ujk az u(xj , k∆t) ´ert´ek k¨ ozel´ıt´es´et, ahol ∆t > 0 egy tetsz˝oleges id˝ol´ep´es. Ujk+1 − Ujk ∆t U0k

=

k Un+1

= 0,

Uj0

=

k k − 2Ujk + Uj+1 Uj−1 , ∆x2

adott u0 seg´ıts´eg´evel.

Explicit Euler-m´odszer (EE)

A szok´asos ´es a q = ∆t/∆x2 jel¨ ol´essel az al´abbi vektoriter´aci´oval ´all´ıthat´o el˝o a numerikus megold´as: k+1

U

k

k

= tridiag [q, 1 − 2q, q] U =: Ah U .

Az alap elv´ar´asunk az el˝ o´all´ıtott numerikus elj´ar´ast´ol, hogy ”tartson a pontos megold´ashoz”. Tegy¨ uk fel, hogy a feladatot a [0, tmax ] id˝ointervallumon kell megoldanunk. Legyen Ejk = Ujk − u(k∆t, j∆x) a numerikus m´odszer hib´aja egy adott r´acspontban ´es k E = [E1k , . . . , Enk ]T .

Konvergencia, konzisztencia, stabilit´as Konvergencia: Azt mondjuk, hogy a numerikus m´ odszer konvergens, ha tetsz˝oleges 0 < t < tmax ´ert´ekre, ha ∆x, ∆t → 0, k k∆t → t, akkor kE k → 0 (konvergenciarend: k kE k = O(∆tr , ∆xs )). Konzisztencia: Azt mondjuk, hogy a numerikus m´odszer konzisztens, ha a Tjk+1 k´eplethib´aja (lok´alis hiba / ∆t) tart null´ahoz, ha ∆t, ∆x → 0 (konzisztenciarend: |Tjk+1 | = O(∆tr , ∆xs )). Stabilit´as: Azt mondjuk, hogy a numerikus m´ odszer stabil, ha v´eges id˝o alatt nem v´altozhat tetsz˝ olegesen nagyot k´et megold´as elt´er´ese a kezdeti elt´er´eshez k´epest a r´acs finom´ıt´as´at´ol f¨ uggetlen¨ ul. Most el´eg lesz az kAh k ≤ 1 felt´etel. Lax-t´etel: Stabilit´as + konzisztencia ⇒ konvergencia (konvergenciarend = konzisztenciarend)

EE-m´odszer konvergenci´aja Konzisztencia: Sz´am´ıtsuk ki a k´eplethib´at! Tjk+1 =

u(xj , tk+1 ) − u(xj , tk ) u(xkj−1 ) − 2u(xkj ) + u(xkj+1 ) − ∆t ∆x2

∆t ∂ 2 u ∆x2 ∂ 4 u (x , ξ) − (η, tk ). j 2 ∂t2 12 ∂x4 Ha a megold´as kell˝oen sima f¨ uggv´eny (peremfelt´etellel konzisztens kezdeti felt´etel is kell hozz´a), akkor vannak olyan pozit´ıv c1 , c2 konstansok, hogy =

|Tjk+1 | ≤ c1 ∆t + c2 ∆x2 . Stabilit´as: K¨onnyen l´athat´o, hogy q ≤ 1/2 eset´en kAh k∞ = 1, ´ıgy ilyenkor a m´odszer stabil lesz.

EE-m´odszer konvergenci´aja Megj.: ∆x2 ∂ 4 u ∆t ∂ 2 u (x , ξ) − (η, tk ) + O(∆t2 , ∆x4 ) j 2 ∂t2 12 ∂x4   ∆t ∂ 4 u 1 = + O(∆t2 , ∆x4 ). (x , t ) 1 − j k 2 ∂x4 6q

Tjk+1 =

q = 1/6 v´alaszt´assal magasabbrend˝ u m´ odszert kapunk (szuperkonvergencia). Megj.: Mi t¨ort´enik akkor, ha a q > 1/2 v´alaszt´ast haszn´aljuk, hiszen csak egy el´egs´eges felt´etelt haszn´altunk a stabilit´asra? A k´ıs´erletek azt mutatj´ak, hogy a q ≤ 1/2 felt´etel sz¨ uks´eges is. Megj.: A konvergenci´ahoz felhaszn´altuk, hogy a megold´as kell˝oen sima. Mi van akkor, ha ez nem ´ıgy van?

EE-m´odszer konvergenci´aja Diszkr´et Fourier-anal´ızis: Legyen a szok´asos m´ odon k k imj∆x Uj = g e . Mikor teljes¨ ul erre a differencias´ema? Akkor, ha a n¨oveked´esi faktor g = 1 − 4q sin2 (m∆x/2). Legyen m = lπ, amivel a numerikus megold´as Ujk =

∞ X

al (g(l))k eilπj∆x

l=1

alakot ¨olti. q > 1/2 eset´en a legnagyobb abszol´ ut´ert´ek˝ u n¨ oveked´esi faktor g ≈ 1 − 4q < −1. Ez az eiπj (1, −1, 1, −1, . . .) ¨ osszetev˝o n¨oveked´esi faktora, azaz korl´atlanul n˝ o, ha k n¨ ovekszik. Ilyen o¨sszetev˝o a kerek´ıt´esi hib´ak miatt biztosan beleker¨ ul az iter´aci´oba.

EE-m´odszer konvergenci´aja A konvergencia igazol´asa abban az esetben, ha u0 Fourier-sora abszol´ ut konvergens ´es q ≤ 1/2 (azaz nem kell, hogy u negyedrend˝ u parci´alis deriv´altjai is folytonosak legyenek a megold´asi tartom´anyon). Legyen ε > 0 tetsz˝ ut konvergencia miatt van P oleges. Az abszol´ olyan l0 , hogy l≥l0 |al | ≤ ε/4. Mivel q ≤ 1/2 eset´en a n¨oveked´esi faktorok 1-n´el kisebb abszol´ ut ´ert´ek˝ uek, ez´ert ∞ X 2 2 |Ejk | = al eilπj∆x ((g(l))k − e−l π k∆t ) l=1

≤

X

|al ||(g(l))k − e−l

2 π 2 k∆t

|+

l
X

|al ||(g(l))k − e−l

2 π 2 k∆t

l≥l0

≤

X l
|al ||(g(l))k − e−l

2 π 2 k∆t

|+2

X l≥l0

|al | ≤ . . .

|

EE-m´odszer konvergenci´aja ... ≤

X

|al ||(g(l))k − e−l

2 π 2 k∆t

|+

l
≤

X

|al |C(q)l4 π 4 ∆t2 k +

l
 = C(q)π 4 tmax 

ε 2

ε 2

 X l
ε |al |l4  ∆t + . 2

Mivel ε tetsz˝oleges volt, ´ıgy ∆t → 0 eset´en a hiba val´oban null´ahoz tart. Kihaszn´altuk, hogy I

|a|, |b| ≤ 1 eset´en |ak − bk | ≤ k|a − b|.

I

g(l) − e−l konstans.

2 π 2 ∆t

≤ C(q)l4 π 4 ∆t2 , ahol C(q) q-t´ol f¨ ugg˝o

Implicit Euler-m´odszer (IE) Ujk+1 − Ujk ∆t

=

k+1 k+1 Uj−1 − 2Ujk+1 + Uj+1

∆x2

,

M´atrixos alak: tridiag [−q, 1 + 2q, −q] U

k+1

k

=U .

Egyenletrendszert kell megoldani, de tridiagon´alis a m´atrix, azaz az inga m´odszerrel megoldhat´ o 8n flop m˝ uvelettel (explicit m´ uveletig´enye 5n flop). N¨oveked´esi faktor: g=

1 , 1 + 4q sin2 (m∆x/2)

ami minden q-ra egyn´el kisebb abszol´ ut ´ert´ek˝ u lesz. Ilyenkor nincs teh´at felt´etel q megv´alaszt´as´ara. Ez egy u ´n. felt´etel n´elk¨ ul stabil s´ema.

θ-m´odszer S´ ulyozzuk egy θ ∈ [0, 1] s´ ullyal az EE ´es IE s´em´akat az al´abbi m´odon: Ujk+1 − Ujk ∆t

k k Uj−1 − 2Ujk + Uj+1 ∆x2 k+1 k+1 k+1 Uj−1 − 2Uj + Uj+1 +θ . ∆x2

= (1 − θ)

M´atrixos alak: tridiag [−θq, 1 + 2θq, −θq] U

k+1 k

= tridiag [(1 − θ)q, 1 − 2(1 − θ)q, (1 − θ)q] U . I I I

θ = 0 : EE-m´odszer, θ = 1 : IE-m´odszer, θ = 1/2 : Crank-Nicolson-m´ odszer.

M˝ uveletsz´ama ´altal´aban 13n flop.

θ-m´odszer N¨oveked´esi faktor: g=

1 − 4(1 − θ)q sin2 (m∆x/2) 1 + 4θq sin2 (m∆x/2)

g ≤ 1 mindig igaz. g ≥ −1 teljes¨ ul´es´ehez 1 − 2q sin2 (m∆x/2)(1 − 2θ) ≥ 0 kell. Ez θ ≥ 1/2 eset´en mindig teljes¨ ul (felt´etel n´elk¨ ul stabil), θ < 1/2 eset´en meg el´eg hozz´a (felt´etelesen stabil) q≤

1 . 2(1 − 2θ)

Konzisztencia (k´eplethiba vizsg´alata). 2 0000 2 000 2 0000 (1/2−θ)∆tu000 xxt −(1/12)∆x uxxxx +(1/24)∆t uttt −(1/8)∆t uxxtt +. . .

θ = 1/2 magasabbrend˝ u lesz, mint a t¨ obbi.

Nemnegativit´as meg˝orz´es Mikor biztos´ıtott az, hogy nemnegat´ıv kezd˝ of¨ uggv´eny eset´en nemnegat´ıv lesz a numerikus megold´as? tridiag [−θq, 1 + 2θq, −θq] U

k+1 k

= tridiag [(1 − θ)q, 1 − 2(1 − θ)q, (1 − θ)q] U . A bal oldalon egy M-m´atrix ´all, aminek inverze nemnegat´ıv, a jobb oldali m´atrix pedig nemnegat´ıv, ha q≤

1 . 2(1 − θ)

Igazolhat´o, hogy ez a felt´etel kell is, ´es hogy a maximum elv is pontosan ekkor ´erv´enyes csak (a t´erbeli feloszt´assz´amt´ol f¨ uggetlen¨ ul). θ = 1/2 eset´en q ≤ 1 a felt´etel!

A h˝ovezet´esi egyenlet k´et dimenzi´oban

∂u = ∆u + kezdeti- ´es peremfelt´etel ∂t A Laplace-oper´ator approxim´aci´ oj´ahoz a Poisson-egyenletn´el tanultakat haszn´aljuk (Ah m´atrix). Alkalmazzuk az egyenesek m´odszer´et: u(t) = −Ah u(t), ∆t ahol u(t) a r´acspontokbeli pontos megold´ast k¨ ozel´ıt˝o f¨ uggv´enyek vektora. EE-m´odszer: U

k+1

k

= (E − ∆tAh )U .

ADI m´odszer a k´etdimenzi´os h˝ovezet´esi egyenletre CN-m´odszer: k+1

(E + (1/2)∆tAh )U

k

= (E − (1/2)∆tAh )U .

A megold´as m˝ uveletig´enyes, mert nem lehet kis s´avsz´eless´eg˝ u m´atrixk´ent fel´ırni Ah − t. Egy lehets´eges megold´as lehet az ADI (alternating directions implicit) m´ odszer. ´Irjuk fel −∆tAh -t −∆tAh = qx Xh + qy Yh alakban. Ezzel a CN-m´ odszer: (E − (qx /2)Xh − (qy /2)Yh )U

k+1

= (E + (qx /2)Xh + (qy /2)Yh )U

k

alak´ u lesz, ami k¨ozel´ıthet˝ o az (E−(qx /2)Xh )(E−(qy /2)Yh )U

k+1

= (E+(qx /2)Xh )(E+(qy /2)Yh )U

formul´aval (a k´eplethiba els˝ o¨ ot tagja azonos).

k

ADI m´odszer a k´etdimenzi´os h˝ovezet´esi egyenletre

Ez az iter´aci´o pedig elv´egezhet˝ o az al´abbi m´ odon: k+1/2

(E − (qx /2)Xh )U

k+1

(E − (qy /2)Yh )U

k

= (E + (qy /2)Yh )U

k+1/2

= (E + (qx /2)Xh )U

.

A m´odszer felt´etel n´elk¨ ul stabil ´es a maximumelv teljes¨ ul´es´enek felt´etele qx , qy ≤ 1.

8. A v´egeselem m´odszer elm´elete (elliptikus feladatok megold´asa)

Bevezet˝o feladatok

Bevezet˝o feladatokk´ent megoldottuk az y 00 = M0 , ´es −∆u = f perem´ert´ekfeladatokat homog´en peremfelt´etellel. Mindegyik esetben egy ´altal´anos´ıtott megold´as k¨ ozel´ıt´es´et ´all´ıtottuk el˝o. A k¨ozel´ıt´eseket a folytonos, szakaszonk´ent folytonosan deriv´alhat´o f¨ uggv´enyek ter´eben kerest¨ uk. Megadtuk a k¨ ozel´ıt´esek hib´ait is. A v´egeselem m´odszer menete m´ar ezekb˝ ol a p´eld´akb´ol j´ol l´atszik.

Bevezet˝o feladatok A v´egeselem m´odszer megval´ os´ıt´as´anak menete. I

Az egyenletb˝ol megkonstru´aljuk a vari´aci´ os egyenleteket (gyenge alak) vagy a minimaliz´aland´ o funkcion´alt.

I

A megold´asi tartom´anyt szakaszokra, h´aromsz¨ogekre, stb. bontjuk.

I

V´alasztunk egy v´egeselem teret a k¨ ozel´ıt´esre (kis tart´oj´ u, szakaszonk´ent sima f¨ uggv´enyek). Hibabecsl´est adunk.

I I

Meghat´arozzuk az elemi merevs´egi m´atrixot. ¨ Ossze´ all´ıtjuk a merevs´egi m´atrixot ´es a terhel´esi vektort.

I

Megoldjuk az egyenletrendszert.

I

Szeml´eltetj¨ uk a megold´ast.

A megold´asi halmaz kiterjeszt´ese Tekints¨ uk az al´abbi perem´ert´ek-feladatot: Reakci´o-diff´ uzi´os egyenlet homog´en Dirichlet peremmel −∆u + u = f, (x, y) ∈ Ω, u|∂Ω = 0 ¯ megold´asf¨ Keress¨ uk az u ∈ C 2 (Ω) ∩ C(Ω) uggv´enyt. Tegy¨ uk fel, hogy u ∈ C 1 (Ω) is igaz. Vegy¨ unk egy tetsz˝oleges 1 v ∈ C0 (Ω) f¨ uggv´enyt (peremen nulla), szorozzunk vele, ´es integr´aljuk a k´et oldalt Ω-n. Z Z Z − ∆uv + uv = f v, Ω

Ω

Ω

alkalmazva a Green-formul´at az al´abbi vari´aci´ os egyenletet nyerj¨ uk: Z Z Z uv = f v, ∀v ∈ C01 (Ω). ∇u∇v + Ω

Ω

Ω

Ennek teljes¨ ul´es´ehez el´eg a f¨ uggv´enyekt˝ ol gyeng´ebb simas´agi felt´eteleket megk¨ovetelni.

Szoboljev-terek

Def.: Tekints¨ uk a C k (Ω) halmazt, ´es l´assuk el az X Z hg, hiH k (Ω) = (∂ α g)(∂ α h) |α|≤k

Ω

skal´aris szorzattal (ez val´ oban skal´aris szorzat). Ez norm´at is defini´al: v uX Z u kgkH k (Ω) = t (∂ α g)2 . |α|≤k

Ω

A fenti teret a fenti norm´aban teljess´e t´eve Hilbert-teret kapunk. Jel¨ol´ese H k (Ω) (Szoboljev-t´er).

Szoboljev-terek T´etel. H k (Ω) izomorf azon o, n´egyzetesen R L2 (Ω)-beli (m´erhet˝ integr´alhat´o f¨ uggv´enyek, f 2 < ∞) g f¨ uggv´enyek ter´evel, melyek minden legfeljebb k-adrend˝ u parci´alis deriv´altja is n´egyzetesen integr´alhat´o, a skal´aris szorzat a fenti m´ odon van ´ertelmezve, ´es k uggv´enysorozat, ami a fenti norm´aban a g van olyan C (Ω)-beli f¨ f¨ uggv´enyhez tart. Hasonl´oan ´ertelmezhet˝ o a H0k (Ω) t´er is a C0k (Ω) teret teljess´e t´eve a kor´abbi norm´aban. T´etel. H0k (Ω) izomorf azon L2 (Ω)-beli g f¨ uggv´enyek ter´evel, melyek minden legfeljebb k-adrend˝ u parci´alis deriv´altja is n´egyzetesen integr´alhat´ o, a skal´aris szorzat a fenti m´odon van ´ertelmezve, ´es van olyan C0k (Ω)-beli f¨ uggv´enysorozat, ami a kor´abbi norm´aban a g f¨ uggv´enyhez tart.

Szoboljev-terek T´erj¨ unk vissza a kor´abbi vari´aci´ os feladathoz! Z Z Z ∇u∇v + uv = f v, ∀v ∈ C01 (Ω). Ω

Ω

Ω

A kor´abbiak alapj´an ez a feladat b˝ ovebb f¨ uggv´enyt´eren is ´ertelmezhet˝o: Keress¨ uk azt az u ∈ H01 (Ω)-beli f¨ uggv´enyt, amelyre 1 a fenti egyenl˝os´eg igaz tetsz˝ oleges v ∈ H0 (Ω) f¨ uggv´eny eset´en (f -r˝ol is el´eg, hogy f ∈ L2 (Ω)). A minimaliz´aci´os feladat megfogalmaz´asa: Keress¨ uk azt az u ∈ H01 (Ω)-beli f¨ uggv´enyt, ami minimaliz´alja az  Z  1 2 2 ((∇u) + u ) − uf I(u) = Ω 2 funkcion´alt.

Neumann-feladat esete Reakci´o-diff´ uzi´os egyenlet Neumann peremmel ∂u = g, (x, y) ∈ ∂Ω ∂n

−∆u + u = f, (x, y) ∈ Ω,

Tegy¨ uk fel, hogy u ∈ C 1 (Ω). Vegy¨ unk egy tetsz˝ oleges v ∈ C 1 (Ω) f¨ uggv´enyt , szorozzunk vele, ´es integr´aljuk a k´et oldalt Ω-n. Z Z Z − ∆uv + uv = f v, Ω

Ω

Ω

alkalmazva a Green-formul´at az al´abbi vari´aci´ os egyenletet nyerj¨ uk: Z Z Z Z − ∇uvn ds + ∇u∇v + uv = f v, ∀v ∈ C 1 (Ω), ∂Ω

Ω

Ω

Ω

azaz Z

Z ∇u∇v +

Ω

Z uv =

Ω

Z fv +

Ω

∂Ω

gv ds, ∀v ∈ C 1 (Ω).

Neumann-feladat esete

Gyenge alak teh´at: Keress¨ uk azt az u ∈ H 1 (Ω) f¨ uggv´enyt, mellyel Z Z Z Z ∇u∇v + uv = fv + gv ds, ∀v ∈ H 1 (Ω). Ω

Ω

Ω

∂Ω

A vari´aci´os feladat: Keress¨ uk azt az u ∈ H 1 (Ω) f¨ uggv´enyt, ami minimaliz´alja az  Z Z  1 2 2 ((∇u) + u ) − uf − gu I(u) = ∂Ω Ω 2 funkcion´alt.

L´enyeges ´es term´eszetes peremfelt´etelek

Def.: Egy peremfelt´etel l´enyeges, ha azt a tesztf¨ uggv´enyeknek is ki kell el´eg´ıteni¨ uk, k¨ ul¨onben term´eszetes peremfelt´etelr˝ol besz´el¨ unk. Vegy¨ uk ´eszre, hogy a Neumann-felt´etel nem ad megk¨ot´est u megv´alaszt´as´ara H 1 (Ω)-b´ ol, a peremfelt´etelt csak implicit m´odon tartalmazza a vari´aci´os egyenlet. A kor´abbi Dirichlet-peremfelt´etel l´enyeges peremfelt´etel volt.

A Neumann-feladat v´egeselemes megold´asa H´aromsz¨ogr´acsot gener´alunk (τh ), most a peremen is defini´alunk cs´ ucsokat. Legyen Vh a szok´asos s´atorf¨ uggv´enyek tere. Keress¨ uk azt az uh ∈ Vh f¨ uggv´enyt, mellyel Z Z Z Z ∇u∇v + uv = fv + gv ds, ∀v ∈ Vh . Ω

Ω

Ω

∂Ω

Kisz´amoljuk a merevs´egi m´atrixot ´es a terhel´esi vektort. Megoldjuk az egyenletrendszert. ´Igy megkapjuk uh -t. T´etel. ku − uh kH 1 (Ω) ≤ ku − vkH 1 (Ω) , ∀v ∈ Vh .

T´etel. ku − uh kH 1 (Ω) ≤ Ch.

A minimaliz´aci´os ´es vari´aci´os feladatok Legyen V val´os vektort´er ´es a : V × V → R szimmetrikus (a(u, v) = a(v, u), ∀u, v ∈ V ), pozit´ıv (a(u, u) > 0, ha u 6= 0) ´es biline´aris (mindk´et v´altoz´ oban line´aris) forma, tov´abb´a legyen l : V → R line´aris funkcion´al. Tekints¨ uk az al´abbi u ´n. minimaliz´aci´ os feladatot: Keress¨ uk meg azt az u ∈ V elemet, ami minimaliz´alja a 1 (M ) J(v) = a(v, v) − l(v) 2 f¨ uggv´enyt V -n! Tekints¨ uk az al´abbi u ´n. vari´aci´ os feladatot: Keress¨ uk meg azt az u ∈ V elemet, mellyel (V ) a(u, v) = l(v), ∀v ∈ V !

(M) ´es (V) ekvivalenci´aja T´etel. u pontosan akkor megold´asa (V)-nek, ha (M)-nek. Az (M) feladatnak maximum egy megold´asa lehet csak. Biz.: (V)⇒(M) Legyen u a vari´aci´ os feladat megold´asa, tov´abb´a 0 6= v ∈ V tetsz˝oleges elem ´es t tetsz˝ oleges val´ os sz´am. Ekkor 1 I(u + tv) = a(u + tv, u + tv) − l(u + tv) 2 1 = (a(u, u) + 2ta(u, v) + t2 a(v, v)) − l(u) − tl(v) 2 t2 = I(u) + t (a(u, v) − l(v)) + a(v, v) . | {z } 2 | {z } (V )⇒0

>0

Azaz u az egyetlen elem, ami minimaliz´alja J-t.

(M) ´es (V) ekvivalenci´aja Biz. (folyt.): (M)⇒(V) Legyen u a minimaliz´aci´ os feladat megold´asa. Az el˝oz˝oek alapj´an a I(u + tv) = I(u) + t(a(u, v) − l(v)) +

t2 a(v, v) 2

t-t˝ol f¨ ugg˝o egyv´altoz´os f¨ uggv´enynek minimuma van t = 0-ban, azaz a deriv´alt itt nulla: (a(u, v) − l(v)) + 0 · a(v, v) = a(u, v) − l(v) = 0. Ezt akartuk igazolni. Megj.: A bizony´ıt´asb´ol kider¨ ult, hogy 1 I(u + v) = I(u) + a(v, v) > I(u), ∀0 6= v ∈ V. 2 Megj.: Hasonl´o t´etelt igazoltunk LER-re NM1-b˝ ol.

A megold´as l´etez´ese Mi biztos´ıtja, hogy van megold´asa (M)-nek vagy (V)-nek? T¨obb felt´etel kell V -r˝ol, a-r´ol ´es l-r˝ ol. I

Legyen V Hilbert-t´er a h., .i skal´aris szorzattal ´es a vele induk´alt k.k norm´aval.

I

a legyen koerc´ıv (V -elliptikus) (szimmetrikus, folytonos (|a(u, v)| ≤ γkukkvk, ∀u, v ∈ V ) ´es ∃α > 0, a(v, v) ≥ αkvk2 , ∀v ∈ V ).

I

l legyen korl´atos line´aris funkcion´al (∃Λ > 0, |l(v)| ≤ Λkvk).

Megj.: A fenti felt´ etelek mellett hu, via = a(u, v) skal´aris szorzat p V -n, ´es a kvka = a(v, v) kifejez´es a k.k norm´aval ekvivalens norm´at defini´al V -n (energianorma).

Lax-Milgram-t´etel T´etel. Legyen V Hilbert-t´er, legyen tov´abb´a a : V × V → R koerc´ıv biline´aris forma, ´es l : V → R korl´atos line´aris funkcion´al. Ekkor az (M) minimaliz´aci´ os feladatnak (´es a kor´abbiak szerint (V)-nek is) egyetlen megold´asa van H-ban. Biz.: Az I(v) ´ert´ekek alulr´ ol korl´atosak V -ben, ugyanis α 1 I(v) = a(v, v) − l(v) ≥ kvk2 − Λkvk 2 2 1 Λ2 Λ2 (αkvk − Λ)2 − ≥− . 2α 2α 2α Legyen c1 = inf{I(v) | v ∈ V }, ´es legyen {vk } ⊂ V olyan sorozat, melyre I(vk ) → c1 . Megmutatjuk, hogy {vk } Cauchy-sorozat. =

Lax-Milgram-t´etel αkvn − vm k2 ≤ a(vn − vm , vn − vm ) = 2a(vn , vn ) + 2a(vm , vm ) − a(vn + vm , vn + vm ) −4l(vn ) − 4l(vm ) + 8l((vn + vm )/2) = 4I(vn ) + 4I(vm ) − 8I((vn + vm )/2) ≤ 4I(vn ) + 4I(vm ) − 8c1 . Mivel vk → c1 , ez´ert minden ε > 0 sz´amhoz van olyan n0 , hogy minden n, m ≥ n0 eset´en αkvn − vm k2 ≤ 4I(vn ) + 4I(vm ) − 8c1 < ε. ´Igy a sorozat val´oban Cauchy. Mivel V Hilbert-t´er, ´ıgy van u ∈ V , hogy vk → u. A folytonoss´ag miatt pedig I(u) = c1 . u teh´at megold´asa az (M) feladatnak (a kor´abbiak szerint (V)-nek is).

Lax-Milgram-t´etel Az egy´ertelm˝ us´eg igazol´asa: Tegy¨ uk fel, hogy az u1 6= u2 elemek is megold´asai (M)-nek. Ekkor I(u1 ) = I(u2 ) = c1 , ´es az u1 , u2 , u1 , u2 , stb. sorozatra nyilv´an I(uk ) → c1 , azaz a sorozat a kor´abbiak szerint Cauchy, de ez csak u ´gy lehet, ha ku1 − u2 k = 0, ami ellentmond´ast ad. Megj.: A t´etel ´altal´anosabban kimondhat´ ou ´gy, hogy V helyett V egy z´art ´es konvex r´eszhalmaz´at tekintj¨ uk. K¨ov.: Legyen a V Hilbert-t´eren a(u, v) = hu, vi. Ez val´oban szimmetrikus, biline´aris, folytonos (C-B-Sch) ´es koerc´ıv. Teh´at minden l korl´atos line´aris funkcion´alhoz egy´ertelm˝ uen l´etezik olyan u ∈ V elem, mellyel hu, vi = l(v), ∀v ∈ V. Riesz reprezent´aci´os t´etele.

Lax-Milgram-t´etel

K¨ov.: Igaz az Λ α stabilit´asi becsl´es, hiszen v = u v´alaszt´assal kuk ≤

αkuk2 ≤ a(u, u) = l(u) ≤ Λkuk, melyb˝ol kuk-val ´es α-val val´ o oszt´as ut´an ad´ odik az ´all´ıt´as.

A Galjorkin-m´odszer Legyen VN a V Hilbert-t´er egy N dimenzi´ os altere, legyen ennek egy b´azisa φ1 , . . . , φN . A Lax–Milgram-t´etel szerint egy´ertelm˝ uen l´etezik olyan uN ∈ VN elem, melyre (V ) a(uN , v) = l(v), ∀v ∈ VN , ill. minimaliz´alja az 1 (M ) I(v) = a(v, v) − l(v) 2 funkcion´alt. Azt a m´odszert, amikor a (V) feladat megold´as´aval adunk k¨ozel´ıt´est u-ra Galjorkin-m´ odszernek nevezz¨ uk. (Nemszimmetrikus biline´aris form´ak eset´en is haszn´alhat´ o.)

A Galjorkin-m´odszer P Keress¨ uk uN -t uN = N oleges VN -beli k=1 ci φi alakban, a tetsz˝ P N elemet pedig ´ırjuk v = k=1 di φi alakban. Ekkor a vari´aci´os feladat alakja ! ! N N N X X X a ci φi , di φi = l di φi , ∀dN ∈ RN , k=1

k=1

k=1

amely az T

T

dN AN cN = dN bN egyenletre vezet, ahol (AN )ij = a(φj , φi ) ´es (bN )i = l(φi ). Mivel ennek minden dN ∈ RN vektorra teljes¨ ulni kell, ez´ert az ismeretlen egy¨ utthat´ok az AN cN = bN egyenletrendszer megold´as´aval nyerhet˝ ok. (AN elnevez´ese merevs´egi (stiffness) m´atrix, bN elnevez´ese terhel´esi (load) vektor.)

A Galjorkin-m´odszer Az T dN AN dN

=a

N X k=1

di φi ,

N X

! di φi

k=1

≥ αk

N X

di φi k2

k=1

egyenl˝os´eg miatt AN pozit´ıv definit, a szimmetrikuss´ag a szimmetri´aj´ab´ol k¨ovetkezik. ´Igy a kor´abbi egyenletrendszer egy´ertelm˝ uen megoldhat´ o. Megj.: A kor´abbiakhoz hasonl´ oan igazolhat´ o, hogy az uN megold´asra ugyanolyan stabilit´asi becsl´es igaz, mint az u-ra, azaz kuN k ≤

Λ . α

Ritz-m´odszer Azt a m´odszert, ahol az (M) feladat megold´as´aval adunk becsl´est u-ra, Ritz-m´odszernek nevezz¨ uk. ! ! ! N N N N X X X X 1 ci φi , ci φi − l ci φi I ci φi = a 2 k=1

k=1

k=1

k=1

1 = cTN AN cN − cTN bN , 2 melynek minimuma szint´en az AN cN = bN egyenletrendszer megold´as´an´al van. Megj.: Szimmetrikus form´ak eset´en a Galjorkin ´es Ritz-m´odszerek ugyanarra az egyenletrendszer megold´as´ara vezetnek. Emiatt a m´odszert ebben az esetben Ritz–Galjorkin-m´ odszernek nevezz¨ uk.

C´ea-lemma Legyen a koerc´ıv biline´aris forma ´es l korl´atos line´aris funkcion´al a V Hilbert-t´eren. Legyen a vari´aci´ os egyenletek megold´asa V -n u, ´es VN -en uN . Ekkor ku − uN k ≤

γ inf ku − vk. α v∈VN

Biz.: Vonjuk ki egym´asb´ ol a a(u, v) = l(v), ∀v ∈ V, a(uN , v) = l(v), ∀v ∈ VN , vari´aci´os egyenleteket VN -beli f¨ uggv´enyekre alkalmazva: a(u − uN , v) = 0, ∀v ∈ VN . ´Irjuk fel v-t vN − uN alakban, ahol vN tetsz˝ oleges VN -beli elem.

C´ea-lemma Ekkor kapjuk az u ´n. Galjorkin-ortogonalit´as k´eplet´et: a(u − uN , vN − uN ) = 0, ∀vN ∈ VN . Tov´abb´a αku − uN k2 ≤ a(u − uN , u − uN ) = a(u − uN , u − vN ) + a(u − uN , vN − uN ) ≤ γku − uN kku − vN k, {z } | =0

amib˝ol k¨ovetkezik a lemma ´all´ıt´asa. A v´egeselem m´odszern´el v-t az u f¨ uggv´eny interpol´aci´os k¨ozel´ıt´es´enek v´alasztva a C´ea-lemma seg´ıts´eg´evel becs¨ ulhet˝o a hiba. Feladat. Gondoljuk ´at a kor´abban megoldott feladatokat, hogy hogy illeszkednek a fent ismertetett elm´eleti keretbe!

A v´egeselem t´er defin´ıci´oja

Def.: A (K, PK , Σ) rendezett h´armast v´egeselem t´ernek h´ıvjuk, ha K valamilyen geometriai alakzat, PK egy K-n ´ertelmezett f¨ uggv´enyekb˝ol ´all´o v´eges dimenzi´ os vektort´er, ´es Σ PK f¨ uggv´enyeinek szabads´agfokainak halmaza (Pk minden eleme egy´ertelm˝ uen meghat´arozott a Σ-beli szabads´agfokokkal). Megj.: 1. K ´altal´aban szakasz, h´aromsz¨ og, t´eglalap, tetra´eder, t´eglatest. 2. PK ´altal´aban polinomokat tartalmaz (Pr (K): legfeljebb r-edfok´ u polinomok halmaza K-n). 3. A szabads´agfokok a polinomok ´es deriv´altjainak ´ert´ekei bizonyos pontokban.

Regularit´asi k¨ovetelm´enyek

A Vh t´er olyan v f¨ uggv´enyek vektortere, melyekre v|K ∈ Pr (K). A K halmazok belsejeinek nincs k¨ oz¨ os pontja, lez´artjaik uni´oja Ω (korl´atos). Egyik K halmaz cs´ ucsa sem eshet egy m´asik ´el´ere vagy lapj´ara. K¨onnyen l´athat´o, hogy ¯ Vh ⊂ H 1 (Ω) ⇔ Vh ⊂ C 0 (Ω) ´es ¯ Vh ⊂ H 2 (Ω) ⇔ Vh ⊂ C 1 (Ω).

N´eh´any konkr´et v´egeselem t´er

I

(1D szakaszok,P1 ,fv.-´ert´ekek a v´egpontokban)

I

(1D szakaszok,P2 ,fv.-´ert´ekek a v´egpontokban ´es k¨oz´epen)

I

(1D szakaszok,P3 ,fv. ´ert´ekek ´es deriv´altak a k´et sz´els˝o pontban)

I

(2D h´aromsz¨ogek,P1 ,fv. ´ert´ekek a cs´ ucspontokban)

I

(2D h´aromsz¨ogek,P2 ,fv. ´ert´ekek a cs´ ucspontokban, oldalfelez˝ok¨on)

I

(2D t´eglalapok,Q1 ,fv. ´ert´ekek a cs´ ucspontokban)

I

(3D tetra´ederek,P1 ,fv. ´ert´ekek a cs´ ucspontokban)

I

stb.

V´egeselem hibabecsl´esek

El˝ozm´enyek Kor´abban l´attuk az egydimenzi´ os u00 = M0 , u(0) = u(H) = 0 perem´ert´ekfeladat eset´en, hogy ku − uh kL2 (0,H) ≤ hM2 H. Teh´at a −u00 = f, u(0) = u(1) = 0 feladat v´egeselemes uh megold´as´ara ´erv´enyes ku − uh kL2 (0,1) ≤ hM2 (els˝orend˝ u konvergencia). Lehet-e jobb (magasabb rend˝ u) becsl´est adni?

L2 (0, 1) hibabecsl´es Gyakorlaton l´attuk, hogy a −u00 = f , u(0) = u(1) = 0 perem´ert´ekfeladat eset´en a r´acspontokban a v´egeselemes megold´as a pontos megold´ast adja. ´Igy teh´at ´erv´enyes az ek (x) = u(x) − uh (x) hib´ara a k. oszt´ointervallumon, hogy ek (x) =

u00 (ξk ) (x − xk−1 )(x − xk ). 2!

Ebb˝ol kapjuk, hogy s Z ke(x)kL2 (0,1) =

1

|e(x)|2 dx ≤

0

(m´asodrend˝ u konvergenci´at mutat).

h2 M2 8

H 1 (0, 1) hibabecsl´es Def.: Vezess¨ uk be az al´abbi u ´n. szeminorm´at H k (Ω)-n: vZ u X u (∂ α v)2 . |v|H k (Ω) = t Ω |α|=k

P´elda. Ezzel a jel¨ol´essel pl. kvk2H 1 (Ω) = kvk2L2 (Ω) + |v|2H 1 (Ω) . ♦ Kor´abban l´attuk, hogy |e0 (x)| ≤ hM2 , azaz most |e0 (x)|H 1 (0,1) ≤ hM2 , ´es ´ıgy ´erv´enyes, hogy ke(x)k2H 1 (0,1) = ku − uh k2H 1 (0,1) = ku − uh k2L2 (0,1) + |u − uh |2H 1 (0,1) 2  2 h M2 + (hM2 )2 . ⇒ ke(x)kH 1 (0,1) = O(h)M2 ≤ 8 (els˝orend˝ u konvergencia).

Hibabecsl´es ´altal´aban

Eml.: C´ea-Lemma: Legyen a vari´aci´ os egyenletek megold´asa V -n u, ´es Vh -n uh . Ekkor ku − uh k ≤

γ inf ku − vk ≤ Cku − Ph uk. α v∈Vh

Amennyiben v-t a pontos u megold´ast interpol´al´ o szakaszonk´ent polinomi´alis f¨ uggv´enynek v´alasztjuk, ´es az interpol´aci´os hib´at becs¨ ulni tudjuk, akkor ezzel a m´ odszerrel hibabecsl´es nyerhet˝o. ´Igy teh´at az ku − Ph uk norm´at kell becs¨ uln¨ unk.

Hibabecsl´es ´altal´aban T´etel. Ha az interpol´aci´ ora szakaszonk´ent r-edfok´ u polinomokat haszn´alunk egy poligon´alis tartom´any egy h´aromsz¨ogr´acs´an (szakasz, h´aromsz¨og, tetra´eder), ´es a r´acsfinom´ıt´as sor´an a (maxim´alis ´elhossz) / (be´ırhat´ o g¨ omb sugara) korl´atos marad, akkor ku − Ph ukL2 (Ω) ≤ Chr+1 |u|H r+1 (Ω) , |u − Ph u|H 1 (Ω) ≤ Chr |u|H r+1 (Ω) , |u − Ph u|H 2 (Ω) ≤ Chr−1 |u|H r (Ω) , stb. Megj.: Ha u csak H s (Ω)-beli (s ≤ r + 1), akkor ku − Ph ukL2 (Ω) ≤ Chs |u|H s (Ω) , |u − Ph u|H 1 (Ω) ≤ Chs−1 |u|H s (Ω) .

Hibabecsl´es ´altal´aban K¨ov.: ku − Ph ukH 1 (Ω) ≤ Chr |u|H r+1 (Ω) , azaz a C´ea-lemma miatt ku − uh kH 1 (Ω) ≤ Chr |u|H r+1 (Ω) , tov´abb´a nyilv´anval´oan ku − uh kL2 (Ω) ≤ Chr |u|H r+1 (Ω) . (1D esetben eggyel nagyobb rend j¨ ott ki. Hogy lehetne ezt kihozni?) Pl.: 1D, folytonos, szakaszonk´ent els˝ ofok´ u v´egeselemek eset´en, r = 1-gyel ´erv´enyesek a kor´abbi becsl´esek. Pl.: 1D, folytonos, szakaszonk´ent m´asodfok´ u v´egeselemek eset´en, r = 2-vel ´erv´enyesek a kor´abbi becsl´esek.

A pontos megold´as regularit´asa Tekints¨ uk a Poisson-egyenletet k´etdimenzi´ oban: −∆u = f, u|∂Ω = 0. Ha Ω pereme sima (sarkok ´es cs´ ucsok n´elk¨ uli) ´es f ∈ H s (Ω), akkor kukH s+2 (Ω) ≤ Ckf kH s (Ω) . Ez a becsl´es nem teljes¨ ul, ha Ω-nak cs´ ucsai is vannak. Pl.: Igazoljuk, hogy egy ω-sz¨ og˝ u k¨ orcikkben az u(r, θ) = rδ sin(δθ) pol´arkoordin´at´akban fel´ırt f¨ uggv´eny megold´asa a ∆u = 0, u = 0 a k¨orcikk k´et sz´el´en feladatnak! A fenti felt´etelek mellett l´athat´ o, hogy u 6∈ H 2 (Ω). ♦ T´etel. Ha a tartom´any egy konvex poligon, akkor igaz, hogy kukH 2 (Ω) ≤ Ckf kL2 (Ω) , azaz a megold´as H 2 (Ω)-beli lesz.

A pontos megold´as regularit´asa Ahhoz, hogy u H s (Ω)-beli legyen, az Z

R

(rδ−s )2 r dr

0

integr´alnak kell v´egesnek lenni. Ez akkor ´es csak akkor teljes¨ ul, ha s < δ + 1. Ha pl. a maxim´alis sz¨og a tartom´any cs´ ucsain´al 3π/2, akkor δ = 2/3, ´es ´ıgy az al´abbi becsl´est nyerj¨ uk: |u − uh |H 1 (Ω) ≤ Chδ = Ch2/3 , azaz nem kapjuk meg a teljes konvergenciarendet.

Hibabecsl´es L2 -ben T´etel. Tegy¨ uk fel, hogy Ω konvex poligon´alis tartom´any, ´es legyen e = u − uh (uh lin. b´azisf¨ uggv´enyes megold´as). Ekkor kekL2 (Ω) ≤ Ch2 |u|H 2 (Ω) . Biz.: Tekints¨ uk az al´abbi u ´n. kieg´esz´ıt˝ o egyenletet: −∆ψ = e, ψ|∂Ω = 0. A kor´abbiak szerint ennek megold´as´ara igaz, hogy kψkH 2 (Ω) ≤ CkekL2 (Ω) . R Azt is l´attuk, hogy h∇e, ∇vh i = Ω ∇e∇vh = 0 minden vh ∈ Vh eset´en. ´Igy kek2L2 (Ω) = he, ei = he, −∆ψi = h∇e, ∇ψi = h∇e, ∇(ψ − Ph ψ)i ≤ kekH 1 (Ω) kψ − Ph ψkH 1 (Ω) ≤ ChkekH 1 (Ω) |ψ|H 2 (Ω) ≤ ChkekH 1 (Ω) kekL2 (Ω) , ahonn´et kekL2 (Ω) -val oszt´assal ad´ odik az ´all´ıt´as.

9. Id˝of¨ugg˝o m´asodrend˝u differenci´alegyenletek (v´egeselem m´odszer)

H˝ovezet´esi egyenlet H˝ovezet´esi egyenlet, homog´en Dirichlet-peremfelt´etel ∂t u − ∆u = f, (t, x) ∈ (0, T ) × Ω, u|(0,T )×∂Ω = 0, u(0, x) = u0 (x) adott. Az elliptikus esethez hasonl´ oan szorozzuk az egyenlet mindk´et uggv´ennyel oldal´at egy tetsz˝oleges, megfelel˝ oen sima v(x) tesztf¨ (peremen nulla), majd integr´aljuk mindk´et oldalt Ω-n! Z Z Z ∂t uv − ∆uv = f v, Ω

Ω

Z

Z

∇u∇v =

∂t uv + Ω

Z

Ω

f v. Ω

Gyenge alak ´Igy teh´at az al´abbi gyenge alakhoz jutunk: Keress¨ uk minden t ∈ (0, T ) eset´en azt az u(t) ∈ H01 (Ω) f¨ uggv´enyt, melyre Z Z Z ∂t uv + ∇u∇v = fv Ω

Ω

Ω

minden v ∈ H01 (Ω) tesztf¨ uggv´eny eset´en, tov´abb´a u(0) = u0 (x). Ebb˝ol a gyenge alakb´ol sz´armaztathatjuk az u ´n. szemidiszkr´et (a t´erbeli diszkretiz´aci´ot a Galjorkin v´egeselem m´ odszer adja) alakj´at az egyenletnek: Legyen Vh a szokott m´ odon egy v´eges dimenzi´ os altere H01 (Ω)-nak a φ1 , . . . , φN b´azisf¨ uggv´enyekkel. Keress¨ uk minden t ∈ (0, T ) eset´en azt az uh (t) ∈ Vh f¨ uggv´enyt, melyre Z Z Z ∂t uh v + ∇uh ∇v = f v, Ω

tov´abb´a

R

Ω uh (0)v

=

Ω

R

Ω u0 v

Ω

minden v ∈ Vh tesztf¨ uggv´eny eset´en.

A szemidiszkr´et megold´as el˝o´all´ıt´asa P ıts¨ uk Keress¨ uk uh (t)-t uh (t) = N i=1 ci (t)φi (x) alakban! Helyettes´ ezt a szemidiszkretiz´alt egyenletbe (legyen v = φj (x) kor´abban l´attuk, hogy el´eg csak a b´azisf¨ uggv´enyekre teljes¨ ulni a gyenge alaknak)! !

Z ∂t Ω

X i

ci (t)φi

!

Z ∇

φj + Ω

X

ci (t)φi

i

Z ∇φj =

f φj , Ω

ahonn´et az al´abbi m´atrixos alakban fel´ırt KDE-rendszert nyerj¨ uk: Bh c0h (t) + Ah ch (t) = bh , ahol (AhR)ij = a(φj , φi ) a kor´abbi merevs´egi (stiffness) R m´atrix, (bh )i = Ω f φi a terhel´esi (load) vektor ´es (Bh )ij = Ω φj φi az u ´n. t¨omeg (mass) m´atrix.

A szemidiszkr´et megold´as el˝o´all´ıt´asa A Bh t¨omeg m´atrix szimmetrikus ´es pozit´ıv definit. ´Igy az al´abbi KDE-rendszert kell megoldanunk: −1 c0h (t) = −B−1 h Ah ch (t) + Bh bh .

Ezt a differenci´alegyenlet-rendszert megoldhatjuk b´armely ´ aban ez egy merev kor´abban tanult numerikus m´ odszerrel is. Altal´ rendszer, ´ıgy ennek megfelel˝ oen implicit m´ odszereket kell a megold´ashoz haszn´alni. T´etel. Ha u megold´asa a k´etdimenzi´ os h˝ ovezet´esi egyenletnek egy Ω konvex, poligon´alis tartom´anyon ´es uh (t) megold´asa a szemidiszkr´et feladatnak line´aris v´egeselemes t´erbeli diszkretiz´aci´o eset´en, akkor van olyan C pozit´ıv konstans, mellyel max ku(t) − uh (t)kL2 (Ω)

t∈(0,T )


≤ C 1 + ln(T /h2 ) max h2 ku(t)kH 2 (Ω) . t∈(0,T )

Teljes (t´er-id˝o) diszkretiz´aci´o Alkalmazzuk a θ-m´odszert a gyenge egyenletben a deriv´altak k¨ozel´ıt´es´ere! Legyen ukh ∈ Vh (k = 0, 1, . . . , kmax , kmax ∆t < T , ∆t > 0) az uh (∆tk) f¨ uggv´eny egy k¨ ozel´ıt´ese! Keress¨ uk azokat az ukh ∈ Vh f¨ uggv´enyeket, melyekre ! Z ukh − uhk−1 v ∆t Ω Z + Ω

∇(θukh + (1 − θ)uk−1 h )∇v

Z (θf (k∆t) + (1 − θ)f ((k − 1)∆t))v,

= Ω

R R uggv´eny eset´en ( Ω u0h v = Ω u0 v) minden v ∈ Vh tesztf¨ (θ ∈ [0, 1]).

Teljes (t´er-id˝o) diszkretiz´aci´o P k Keress¨ uk az ukh f¨ uggv´enyeket ukh = N es ´ırjuk fel i ci φi alakban, ´ az egyenleteket csak a v = φj b´azisf¨ uggv´enyekre! Ekkor a ckh = [ck1 , . . . , ckN ]T vektorokra az al´abbi egyenl˝ os´egeket kapjuk (k = 1, . . . , kmax ): 1 k−1 k Bh (ckh − ck−1 h ) + θAh ch + (1 − θ)Ah ch ∆t θ,k

= θbh (k∆t) + (1 − θ)bh ((k − 1)∆t) =: bh . Ez ugyanazt az iter´aci´ ot adja, amit a szemidiszkretiz´aci´o ut´an nyert KDE-rendszer szok´asos θ-m´ odszeres megold´asa sor´an nyer¨ unk. Eml´ekezz¨ unk, hogy a θ-m´ odszer θ ≥ 1/2 eset´en lesz felt´etel n´elk¨ ul stabil (pl. Crank–Nicolson-m´ odszer, implicit Euler-m´odszer).

Nemnegativit´as meg˝orz´ese

Mi biztos´ıtja, hogy ha f ≥ 0, u0 (x) ≥ 0, akkor a megold´asra is a fizik´ab´ol elv´art ukh ≥ 0 teljes¨ ul? El´eg ehhez, hogy (Bh + θ∆tAh )−1 ≥ 0, ´es (Bh + θ∆tAh )−1 (Bh − (1 − θ)∆tAh ) ≥ 0 teljes¨ ulj¨on.

Nemnegativit´as meg˝orz´ese Pl. az egydimenzi´os feladat line´aris v´egeselemes megold´asa eset´en, amikor ekvidiszt´ans a feloszt´as: Bh =

1 h tridiag[1, 4, 1], Ah = tridiag[−1, 2, −1]. 6 h

´Igy az el˝obbi egyszer˝ u el´egs´eges felt´etellel kapjuk, hogy az 1D line´aris v´egeselemes megold´as nemnegativit´asmeg˝ orz˝o lesz, ha θ ∈ [1/3, 1] ´es 1 1 ≤q≤ . 6θ 3(1 − θ) (Igazolhat´o (Farag´o, 1996), hogy ez a felt´etel sz¨ uks´eges is tetsz˝oleges feloszt´as eset´en.)

Tanuls´ag

I

A m´asodrend˝ u Crank–Nicolson-m´ odszer ”hi´aba” felt´etel n´elk¨ ul stabil, a nemnegativit´asmeg˝ orz´eshez q nem lehet tetsz. nagy (2D ´es 3D-ben sem!)

I

Als´o korl´at is van q-ra! (mass lumping eset´en nincs)

I

Implicit Euler-m´odszer u ´gy nemnegativit´asmeg˝orz˝o, hogy csak als´o korl´at van q-ra.

Hull´amegyenlet megold´asa Tekints¨ uk az al´abbi hull´amegyenletet homog´en peremfelt´etellel: Hull´amegyenlet

u|∂Ω

∂2u − ∆u = f, ∂t2 = 0, u0 (x) = u(0, x) ´es u1 (x) = ∂t u(0, x) adottak.

Keress¨ uk az u megold´ast a (0, T ) intervallumon! Gyenge alak: Keress¨ uk minden t ∈ (0, T ) eset´en azt az u(t) ∈ H01 (Ω) f¨ uggv´enyt, melyre Z Z Z 2 ∂t uv + ∇u∇v = fv Ω

Ω

Ω

minden v ∈ H01 (Ω) tesztf¨ uggv´eny eset´en, tov´abb´a u(0) = u0 (x), ∂t u(0) = u1 (x).

Szemidiszkr´et megold´as

Keress¨ uk minden t ∈ (0, T ) eset´en azt az uh (t) ∈ Vh f¨ uggv´enyt, melyre Z Z Z 2 ∇uh ∇v = f v, ∂t uh v + Ω

Ω

R

tov´abb´a Ω uh (0)v = tesztf¨ uggv´eny eset´en.

R

Ω u0 v,

R

Ω ∂t uh (t)v

Ω

=

R

Ω u1 v

minden v ∈ Vh

A szemidiszkr´et megold´as el˝o´all´ıt´asa P Keress¨ uk uh (t)-t uh (t) = N ıts¨ uk i=1 ci (t)φi (x) alakban! Helyettes´ ezt a szemidiszkretiz´alt egyenletbe (legyen v = φj (x) kor´abban l´attuk, hogy el´eg csak a b´azisf¨ uggv´enyekre teljes¨ ulni a gyenge alaknak)! Z Ω

! ∂t2

X i

ci (t)φi

!

Z ∇

φj + Ω

X

ci (t)φi

i

Z ∇φj =

f φj , Ω

ahonn´et az al´abbi m´atrixos alakban fel´ırt KDE-rendszert nyerj¨ uk: Bh c00h (t) + Ah ch (t) = bh . (A kezdeti felt´etelt az eredeti kezdeti felt´etelekb˝ ol nyerhet¨ unk.)

A szemidiszkr´et megold´as el˝o´all´ıt´asa Meg kell oldani a m´asodrend˝ u KDE-rendszert az adott kezdeti felt´etelekkel. Ez hat´ekonyan megtehet˝ oa 1 k k−1 Bh (dh − dh ) + θAh ckh + (1 − θ)Ah ck−1 h ∆t θ,k

= θbh (k∆t) + (1 − θ)bh ((k − 1)∆t) =: bh , ck − ck−1 k k−1 − (αd + (1 − α)d ) = 0, ∆t 0

c0 ´es d a kezdeti felt´etelekb˝ ol ismert. A m´odszer θ, α ≥ 1/2 eset´en felt´etel n´elk¨ ul stabil, ´es pl. ha mindkett˝o 1/2, akkor m´asodrend˝ u. Megj.: Kell˝oen sima kezdeti f¨ uggv´enyek eset´en j´ ol viselkedik a m´odszer, de nemfolytonos kezdeti f¨ uggv´enyek eset´en hamis oszcill´aci´ok l´epnek fel a megold´asban.

10. Az FDTD m´odszer a Maxwell-egyenletek numerikus megold´as´ara

Maxwell-egyenletek

ε∂t E = ∇ × H − σE − Je , µ∂t H = −∇ × E + divergence equations. E = E(t, x, y, z) elektromos t´erer˝ oss´eg H = H(t, x, y, z) m´agneses t´erer˝ oss´eg ε = ε(x, y, z), µ = µ(x, y, z), σ = σ(x, y, z) anyagi param´eterek Je = Je (t, x, y, z) forr´asf¨ uggv´eny Numerikusan kell megoldani!

FDTD (Finite difference time domain) vagy Yee-m´odszer (1966)

I

K´et eltolt r´acs.

I

T´erbeli l´ep´est´avols´agok: ∆x, ∆y, ∆z

I

Id˝ol´ep´es: ∆t

I

Leapfrog s´ema

Yee-m´odszer (1966) El˝ony¨ok

H´atr´anyok

A divergenciafelt´etel automatikusan teljes¨ ul

I

Numerikus diszperzi´o.

I

L´epcs˝ozetes perem.

I

Explicit m´odszer. (Nem kell megoldani egyenletrendszert.)

I

I

Csak a t´erer˝oss´egkomponenseket kell t´arolni.

Szigor´ u stabilit´asi felt´etel. ∆ ∆t < √ c 3

I

K¨onnyen implement´alhat´ o.

I

Az anyagi param´eterek f¨ ugghetnek a helyt˝ ol ´es az ir´anyt´ol.

I

Vill´amcsap´as hat´asa rep¨ul˝on

Frequency: 10kHz −→ hull´amhossz: 30km ↓ 3km-es r´acst´avols´ag elegend˝ o lenne, de a r´acst´avols´agot kb. ∆x = ∆y = ∆z = 30cm-nek kell v´alasztani.

Forr´ as: www.efieldsolutions.com/Pics/SAAB2000.gif

Eltolt r´acsos szemidiszkretiz´aci´o ε∂t E = ∇ × H, µ∂t H = −∇ × E. √ E(t, x, y, z) −→ εE(t, x, y, z) √ H(t, x, y, z) −→ µH(t, x, y, z) [Kole et al, 2001] Ψ0 (t) = AΨ(t), Ψ(0) = Ψ0 is given. I

Az u ´n. Yee-cell´ak sz´ama: Y .

I

A ∈ R6Y ×6Y , A ferd´en szimmetrikus.

I

A ritka m´atrix, max. 4 nemnulla elem soronk´ent az √ 1/( εµ∆) alakban.

2D p´elda T¨ ok´eletes vezet˝o perem. ε=µ≡1

    −1   −1  1 0 −1 Ψ (t) = AΨ(t) =  ∆  1   1   1 −1 1

−1 −1

1 1 1 −1

 1 −1        Ψ(t)      

Yee-m´odszer A = Ayr + Ayg , Szekvenci´alis szeletel´es T´etel. (Ayr )2 = (Ayg )2 = 0, exp(τ Ayr ) = E + τ Ayr , exp(τ Ayg ) = E + τ Ayg .

exp(τ Ayr ) exp(τ Ayg )Ψn = (I + τ Ayr )(E + τ Ayg )Ψn =: Ψn+1 Els˝orend˝ u m´odszer.

Yee-m´odszer

ΨM = exp

1 2 τ Ayr




exp (τ Ayg ) 

exp(τ Ayr ) exp(τ Ayg ) . . . exp(τ Ayr ) exp(τ Ayg )exp

 1 τ Ayr Ψ0 2

↓  exp

   1 1 τ Ayr exp τ Ayr 2 2

Eltolt kezdeti f¨ uggv´enyekkel m´asodrend˝ u s´em´at kapunk (klasszikus Yee-m´odszer, Strang-Marcsuk szeletel´es). Szigor´ u stabilit´asi felt´etel.

ADI-FDTD-m´odszer [Namiki, Zheng et al 1999] A = Anr + Anb SM-szeletel´est haszn´aljuk. A r´eszfeladatokat az EE-, CN-, IE-m´ odszerekkel oldjuk meg.  τ  Anr exp(τ Anb ) exp Anr Ψn = 2 2 ∆t ∆t ∆t ∆t = (E− Anr )−1 (E− Anb )−1 (E+ Anb )(E+ Anr )Ψn =: Ψn+1 2 2 2 2 exp

τ

∆t: a numerikus m´odszer l´ep´est´avols´aga (∆t = τ ). M´asodrend˝ u, felt´etel n´elk¨ ul stabil.

T´erbeli ir´anyok szerinti szeletel´es [Kole et al, 2001]

A = Akb + Akv + Akr + Akg , Szekvenci´alis vagy Strang-Marcsuk szeletel´esek haszn´alhat´ok.

T´erbeli ir´anyok szerinti szeletel´es

T´etel. A r´eszfeladatok egzaktul megoldhat´ ok az     0 α cos α sin α exp = −α 0 − sin α cos α azonoss´ag felhaszn´al´as´aval. exp(τ Akg ) exp(τ Akr ) exp(τ Akv ) exp(τ Akb )Ψn =: Ψn+1 A m´odszer rendje az alkalmazott szeletel´es rendj´evel egyezik meg. Felt´etel n´elk¨ ul stabil.

11. A Black–Scholes-egyenlet

N´eh´any opci´o´araz´asi alapfogalom Def. Vagyoni ´ert´ek: olyan dolog, melynek ismerj¨ uk az ´ar´at a jelenben, ami v´altozhat az id˝ o folyam´an. (Pl. r´eszv´eny, arany, olaj, elektromos ´aram, deviza, stb. Def. Eur´opai v´eteli opci´ o: (European call option, el˝ov´as´arl´asi jog) birtokos´anak megadja a jogot (de nem a k¨ oteless´eget), hogy megvegye az el˝ore r¨ogz´ıtett vagyoni ´ert´eket az opci´o jegyz˝oj´et˝ol az el˝ore r¨ogz´ıtett leh´ıv´asi ´aron ´es el˝ ore r¨ ogz´ıtett id˝ opontban (lej´arati id˝o). Pl. Alad´ar (a jegyz˝o) k¨ ot egy eur´ opai v´eteli opci´ ot B´el´aval (a birtokos). Ez megadja a jogot B´el´anak arra, hogy vegyen Alad´art´ol 100 IBM r´eszv´enyt 1000$-´ert 3 h´ onap m´ ulva. Mi t¨ort´enik h´arom h´onap m´ ulva? I

I

a) Ha 100 IBM r´eszv´eny ´ert´eke t¨ obb lesz, mint 1000$, akkor B´ela leh´ıvja az opci´ ot, ´es az aktu´alis ´ert´ekn´el olcs´obban kapja meg a r´eszv´enyt, ´ıgy profithoz jut. b) Ha 100 IBM r´eszv´eny ´ara kevesebb, mint 1000$, akkor B´ela nem h´ıvja le az opci´ ot.

N´eh´any opci´o´araz´asi alapfogalom Def. Az eur´opai elad´asi opci´ o (put option) hasonl´ o a v´eteli opci´ohoz, csak elad´asi jogot ad. Az amerikai opci´ ok pedig a lej´arati id˝o el˝ott b´armikor leh´ıvhat´ ok. Probl´ema: Ez ´ıgy nem igazs´agos. Mennyibe ker¨ ulj¨on az opci´o jegyz´ese? El˝ony¨ok: I

a) Az opci´o ´ara meghat´arozhat´ o, ´es mindenki megel´eged´es´ere adhat´o ´es vehet˝o.

I

b) Az opci´o ´ara sokkal kevesebb, mint a vagyoni ´ert´ek ´ara. A befektet´eshez k´epest nagyobb profit ´erhet˝ o el, mint a vagyoni ´ert´ek megv´as´arl´as´aval.

I

c) Opci´okat kombin´alva a vagyoni ´ert´ek k¨ ul¨ onb¨oz˝o v´altoz´asaib´ol nyeres´egre tehet¨ unk szert

Egy p´elda

1 r´eszv´eny ´ara 10$, ´es az ´ara 1 ´ev m´ ulva vagy 8$ lesz, vagy 12$. A banki kamatl´ab 10%. Tegy¨ uk fel, hogy 10$-os lej´arati ´aron lehet eur´opai v´eteli opci´ot jegyezni. Az opci´ o ´ara 1$ r´eszv´enyenk´ent. 1. strat´egia: vegy¨ unk 200 opci´ ot ´es adjunk el 100 r´eszv´enyt 200 opci´o v´etele 100 r´eszv´eny elad´asa Befektet´es Nyeres´eg

-200 1000 -800 0

S(T ) = 12 400 -1200 880 +80

S(T ) = 8 0 -800 880 +80

Arbitr´azs: kock´azat n´elk¨ uli p´enzszerz´es. Nem megengedett.

Egy p´elda

A v´eteli opci´o ´ara legyen 2$. Akkor adjunk el 200 opci´ot ´es vegy¨ unk 100 r´eszv´enyt. Elad´as 200 opci´o V´etel 100 r´eszv´eny K¨olcs¨on Nyeres´eg

400 -1000 600 0

S(T ) = 12 -400 1200 -660 +140

S(T ) = 8 0 800 -660 +140

Megj.: 1.36$-os opci´o´ar mellett nincs arbitr´azs.

A kifizet´esi f¨uggv´eny (payoff function) Jel¨ol´esek: S(t): a vagyoni ´ert´ek ´ara a t id˝ opillanatban. T : lej´arati id˝o. E: leh´ıv´asi ´ar. C(t): a v´eteli opci´o ´ara a t id˝ opontban. P (t): az elad´asi opci´o ´ara a t id˝ opontban. V (t): egy ´altal´anos opci´ o ´ert´eke a t id˝ opontban. σ: vagyoni ´ert´ek ´ar´anak volatilit´asa (v´altoz´ekonys´aga, kb. 0.05 ´es 0.5 k¨oz¨ott). r: kamatt´enyez˝o (pl. 0.1 10% kamatl´ab eset´en). Def. A kifizet´esi f¨ uggv´eny azt adja meg, hogy mekkora lesz a birtokos nyeres´ege az opci´ on a lej´aratkor. Pl. tekints¨ uk az eur´opai elad´asi opci´ ot. Ennek kifizet´esi f¨ uggv´enye: ( E − S(T ), ha S(T ) ≤ E, P (T ) = = max{E − S(T ), 0}. 0, ha S(T ) > E,

A Black–Scholes-egyenlet

Fischer Black, Myron Scholes (1973) Robert Merton az opci´ o´araz´asi technika pontos matematikai megfogalmaz´asa. 1997: Merton ´es Scholes k¨ ozgazdas´agi Nobel-d´ıj Eur´opai opci´okra: 1 ∂2V ∂V ∂V + σ 2 S 2 2 + rS − rV = 0 ∂t 2 ∂S ∂S + v´eg- ´es peremfelt´etelek.

A Black–Scholes-egyenlet Eur´opai elad´asi opci´okra: 1 ∂2P ∂P ∂P + σ 2 S 2 2 + rS − rP = 0 ∂t 2 ∂S ∂S + P (S, T ) = max{E − S(T ), 0}, P (0, t) = Ee−r(T −t) , P (S, t) → 0, ha S → ∞. Kezdeti´ert´ek-feladatot csin´alunk: (t = (T − t)) ∂P 1 ∂2P ∂P − σ 2 S 2 2 − rS + rP = 0 ∂t 2 ∂S ∂S + P (S, 0) = max{E − S(0), 0}, P (0, t) = Ee−rt , P (S, t) → 0, ha S → ∞.

A BSch-egyenlet v´eges differenci´as megold´asa S = ∞ helyett csak Smax -ig vessz¨ uk fel az S ´ert´ekeket. ∆S = Smax /(n + 1), ∆t az id˝ ol´ep´es. Explicit Euler-m´odszer k P k − 2Pik + Pi+1 Pik+1 − Pik 1 − σ 2 (i∆S)2 i−1 ∆t 2 (∆S)2 k − Pk Pi+1 i−1 + rPik = 0, 2∆S amely bevezetve a Q = tridiag[1, −2, 1], D = tridiag[−1, 0, 1], N1 = diag(1, 2, . . . , n), N2 = diag(12 , 22 , . . . , n2 ) m´atrixokat, ´es a pk vektort a r´acspontokbeli k¨ ozel´ıt´esekre a k-adik id˝or´etegen, az al´abbi numerikus s´em´at nyerj¨ uk:

−ri∆S

1 1 k pk+1 = pk + ((−r∆t)E + ∆tσ 2 N2 Q + ∆trN1 D) pk + f . 2 2 | {z } AEE

A BSch-egyenlet v´eges differenci´as megold´asa k

Az f vektor tartalmazza a peremfelt´eteleket az al´abbi m´odon. k

k f = [(1/2)∆t(σ 2 − r)P0k , 0, 0, . . . , 0, (1/2)∆tn(σ 2 n + r)Pn+1 ]T , k ahol P0k ´es Pn+1 a peremfelt´etelb˝ ol ismertek, a p0 vektor pedig a kezdeti felt´etelb˝ol.

Implicit Euler-m´odszer 1 1 k+1 pk+1 = pk + ((−r∆t)E + ∆tσ 2 N2 Q + ∆trN1 D) pk+1 + f . 2 2 | {z } AIE

Crank–Nicolson-m´odszer k

(E − AIE /2)pk+1 = (E + AEE /2)pk + f /2 + f

k+1

/2.

Eur´opai elad´asi opci´o megold´asf¨uggv´eny E = 40, Smax = 100, r = 0.1, T = 1, σ = 0.4

Amerikai opci´ok vizsg´alata V´eteli opci´o eset´en nem ´erdemes kihaszn´alni a kor´abbi leh´ıv´as lehet˝os´eg´et. Ugyanis I

Ha t-ben S(t) > E, akkor leh´ıv´as eset´en S(t) − E lenne a nyeres´eg. (t = T -ben er(T −t) (S(t) − E))

I

Helyette adjunk el egy r´eszv´enyt a t id˝ opontban. (Bev´etel: r(T −t) S(t), t = T -ben e S(t))

I

t = T -ben: ha dr´ag´abb a r´eszv´eny, mint E, akkor h´ıvjuk le az opci´ot ´es vegy¨ uk meg E-´ert, ha nem dr´ag´abb, akkor S(T )-´ert.

I

Bev´etel: er(T −t) S(t), kiad´as maximum E. Nyeres´eg nem kevesebb, mint er(T −t) S(t) − E ≥ S(t) − E > 0. ´Igy a korai leh´ıv´asn´al van jobb strat´egia.

I

Teh´at az amerikai v´eteli opci´ o ´ara ugyanaz mint az eur´opai´e.

Amerikai elad´asi opci´ok Jel¨olje a kifizet´esi f¨ uggv´enyt Λ(S(t)) := max{E − S(t), 0}. ´ ıt´as. P Am (S, t) ≥ Λ(S(t)). All´ Biz.: K¨ ul¨onben vegy¨ unk egy opci´ ot ´es egyszerre h´ıvjuk is le. Bev´etel: −P Am (S, t) + Λ(S(t)) > 0. Ez nem lehet (arbitr´azs). Igazolhat´o, hogy a fenti felt´etelen k´ıv¨ ul teljes¨ ul az al´abbi parci´alis differenci´alegyenl˝otlens´eg: ∂P Am ∂P Am 1 2 2 ∂ 2 P Am + σ S + rS − rP Am ≤ 0. ∂t 2 ∂S 2 ∂S A peremfelt´etelek ´es v´egfelt´etel: P Am (S, T ) = max{E − S(T ), 0}, P Am (0, t) = E, P Am (S, t) → 0, ha S → ∞.

Amerikai elad´asi opci´ok

Minden t ∈ [0, T ]-re van egy S ? (t) ´ert´ek, amin´el kisebb S ´ert´ekekre ´erdemes leh´ıvni az opci´ ot, nagyobbakra ´erdemes megtartani. S ? (T ) = E, tov´abb´a S ? (t) t f¨ uggv´eny´eben n¨ovekszik. Hol van ez a hat´ar egy adott t eset´en? S ? (t)-n´el kisebb S ´ert´ekekre nyilv´an az opci´ o ´ert´eke ´eppen a kifizet´esi f¨ uggv´eny, nagyobbakra pedig kiel´eg´ıti a Black–Scholes-egyenletet az ´ert´eke (mozg´ o perem˝ u perem´ert´ekfeladat).

Amerikai elad´asi opci´ok

Az al´abbi u ´n. line´aris komplementarit´asi feladatot kell megoldani (´att´er¨ unk kezdeti´ert´ek-feladatra): 

 ∂P Am ∂P Am 1 2 2 ∂ 2 P Am Am − σ S − rS + rP (P Am −Λ(S(t))) = 0, ∂t 2 ∂S 2 ∂S ∂P Am 1 2 2 ∂ 2 P Am ∂P Am − σ S − rS + rP Am ≥ 0, ∂t 2 ∂S 2 ∂S P Am − Λ(S(t)) ≥ 0.

A line´aris komplementarit´asi feladat numerikus megold´asa 

k

(E − AIE /2)pk+1 − (E + AEE /2)pk − f /2 − f

k+1

 /2 .·(pk+1 −p0 ) = 0

A bal oldali z´ar´ojelben egy line´aris egyenletrendszert kell megoldani pk+1 -re. Ezt az u ´n. projekci´ os SOR-m´ odszerrel oldjuk meg. Azaz Gauss–Seidel-iter´aci´ot alkalmazunk a megold´asra, ´es minden l´ep´esben vessz¨ uk a pk+1,j+1 ´es p0 vektorok maximum´at addig, m´ıg ”nem konverg´al” az iter´aci´ o. Legyen E − AIE /2 = F − L − R, F a diagon´alis, L ´es R pedig az als´o ´es fels˝o h´aromsz¨ogr´esz -1-szerese. Legyen pk+1,0 = pk , ´es tekints¨ uk az al´abbi iter´aci´ ot: k

pk+1,j+1 = (F−L)−1 (Rpk+1,j +(E+AEE /2)pk +f /2+f pk+1,j+1 = max{pk+1,j+1 , p0 }.

k+1

/2),

A line´aris komplementarit´asi feladat numerikus megold´asa

Addig iter´alunk, m´ıg kpk+1,j+1 − pk+1,j k ≤ ε, ahol ε egy el˝ore r¨ogz´ıtett tolerancia´ert´ek. Ezek ut´an pk+1 := pk+1,j+1 , amivel befejezt¨ uk a (k + 1). l´ep´est. Ekkor nyilv´an a konstrukci´ o miatt pk+1 ≥ p0 is teljes¨ ul, tov´abb´a k

(E − AIE /2)pk+1 − (E + AEE /2)pk − f /2 − f

k+1

/2 ≥ 0,

azaz pk+1 megold´asa lesz a line´aris komplementarit´asi feladatnak.