Numerikus Matematika a gyakorlatban Dr. Vertse Tam´ as Jegyzet Szeptember 12. K¨ ovetelm´ enyek: jelenl´ et aj´ anlott, saj´ at TIGRIS account kell, ´ erdemjegyet a programra kapnak Bevezet´ es a VMS-be : file szerkezet, DEL, PU utasit´ asok.
Bevezet´ es a FORTRAN nyelvbe FOR LINK COP.COM F.COM szeptemeber 19 JEGYZET.txt
F¨ uggv´ enyk¨ ozelit´ esek, inter- ´ es extrapol´ aci´ o Lagrange interpol´ aci´ o LAG.FOR Spline interpol´ aci´ o @COP TRIDAG.FOR @COP TRIDAG.DAT ED TRIDAG.RES TRIDAG.FOR VSW.FOR Febru´ ar 15. Az ED editor haszn´ alata Egyszer˝ u Fortran program bemutat´ asa: LAG.FOR LAG.COM ASS DEASS LAGRANGE.FOR Febru´ ar 22. A Lagrange interpol´ aci´ oval nem kapunk becsl´ est az interpol´ aci´ o hib´ aj´ ara. Ebben seg´ıt a Nevill f´ ele algoritmus. SUBROUTINE POLINT(XA,YA,N,X,Y,DY) @COP POLINT.FOR @COP MAINPOLINT.FOR Interpol´ aci´ os polinom egy¨ utthat´ oinak sz´ am´ıt´ asa Vandermonde matrix egyenletrendszer´ enek megold´ asa MAINPOLCOF.FOR POLCOF.FOR FOR MAINPOLCOF,POLCOF LIN MAINPOLCOG,POLCOF 1 undefined symbol: POLINT LIN MAINPOLCOF,POLCOF,POLINT
R MAINPOLCOF FORTRAN PAUSE Hibakeres´ es DEBUGER-rel ED FDBG.COM *C $FOR/NOOP/DEB ’P1 $LIN/DEB ’P1 CNRL Z *EX @FDBG MAINPOLCOF,POLCOF,POLINT R MAINPOLCOF implicit real*8 (a-h,o-z) besz´ ur´ asa a MAINPOLINT f˝ oprogramba ´ es a POLCOF subroutine-ba @F JO erredm´ eny: p(x) = −6.05 + 0.44167x + 0.7x2 − 0.09167x3
M´ arcius 1. Inverzinterpol´ aci´ o @COP INVERZPOL.FOR
A legkisebb n´ egyzetek m´ odszere Line´ aris eset: norm´ alegyenletek 1. p´ elda illeszt´ es x=0-ra centr´ alt polinommal φj (x) = xj
@COP NORMAL.FOR @COP NORMAL.DAT FOR/D LINES NORMAL ED fd.com *C Polinom helyettes´ıt´ esi ´ ert´ ek´ enek kisz´ am´ıt´ asa Horner elrendez´ essel FUNCTION HORN(N,A,X) Visszat´ er´ es a legkisebb n´ egyzetek m´ odszer´ ehez Tov´ abbra is line´ aris eset. 2. p´ elda at a-ra M´ odositsuk a NORMAL.FOR programot arra az esetre, ha a φj (x) = (x − a)j teh´
centr´ alt polinom read(5,*) a do i=1,n x(i)=x(i)-a enddo Megold´ as:
NORMAL1.FOR NORMAL1.DAT 3. p´ elda exponenci´ alis f¨ uggv´ enyek line´ arkombin´ aci´ oja: φj (x) = e−λj x
a λj boml´ asi ´ alland´ ok ´ ert´ ekei adottak. @COP NORMAL2.FOR @COP NORMAL2.DAT @F NORMAL2 R NORMAL2 fortran stop ED NORMAL2.RES c(0)=5.0503 c(1)=2.9718 y(x) = 5.0503e−x + 2.9718e−0.5x
Ha a λ-kat is ismeretleneknek tekinten´ enk, akkor ez a probl´ ema a nemline´ aris egyenletrendszerek megold´ as´ at k¨ oveteli. aci´ oval A norm´ alegyenletek line´ aris egyenletrendszerek. Ezek megold´ as´ at Gauss Jordan elimin´ (kik¨ usz¨ ob¨ ol´ essel) v´ egezt¨ uk. Ezt a feladatot a SUBROUTINE GAUSSJ(A,N,NP,B,M, MP) v´ egezte el. A param´ eterek jelent´ ese. Egy¨ utthat´ ok m´ atrixa: A(NP,NP) ismeretlenek sz´ ama: N¡=NP T¨ obbsz¨ or¨ os jobboldal: B(NP,MP) A sz´ amol´ as v´ egen ez fogja tartalmazni a megold´ ast. Line´ aris egyenletrendszerek megold´ asa
Line´ aris egyenletrendszerek megold´ asa Direkt m´ odszerek: Gauss kik¨ usz¨ ob¨ ol´ eses m´ odszer: ismeretlenek fokozatos kik¨ usz¨ ob¨ ol´ ese, amig fels˝ o h´ aromsz¨ ogm´ atrixot kapunk. Ezut´ an a kapott egyenletrendszert visszahelyettesit´ essel oldjuk meg. Gauss m´ odszer determin´ ans kisz´ amit´ as´ ara: f˝ oa ´tl´ obeli elemek szorzata (−1)k , ahol k a sorcser´ ek sz´ ama a h´ aromsz¨ ogm´ atrix´ a alakit´ as sor´ an. M´ atrixinverzi´ o Gauss m´ odszerrel: Inverzm´ atrix oszlopvektorai az ismeretlenek, jobboldalon az egys´ egvektorok. Gauss-Jordan kik¨ usz¨ ob¨ ol´ es (elimin´ aci´ o). Hasonl´ o a Gauss m´ odszerhez, de enn´ el az xk ismeretleneket valamennyi egyenletb˝ ol kik¨ usz¨ ob¨ olj¨ uk (nemcsak az adott egyenlet alattiakb´ ol). A k=n-1 edik lesz´ armaztatott egyenletrendszer m´ atrixa diagon´ alis lesz, igy a f˝ oa ´tl´ obeli elemekkel beosztva egys´ egmatrixot kapunk. A jobboldal fogja az inverzm´ atrixot tartalmazni. F˝ oelem kiv´ alaszt´ as(pivoting) r´ eszleges teljes Line´ aris egyenletrendszer megold´ asa Gauss-Jordan elimin´ aci´ oval @COP MAINGAUSJ.FOR
@COP GAUSJ.DAT @COP GAUSJ.COM M´ atrixinverzi´ o Gauss-Jordan elimin´ aci´ oval @COP GAUSJINV.FOR @COP GAUSJINV.DAT @F GAUSJINV R GAUSJINV Line´ aris egyenletrendszer megold´ asa h´ aromsz¨ ogm´ atrixokra bont´ assal A=L . U Ax = b A.x = L.(Ux) = L.y = b, ahol U.x = y L.y = b el˝ orehelyettesit´ essel oldhat´ o meg, y ismeret´ eben az U.x = y visszahelyettesit´ essel. A = L.U felbont´ as a Crout f´ ele algoritmussal t¨ ort´ enik. T´ arol´ as az A-m´ atrix hely´ en (ez´ ert
el˝ osz¨ or A-t c´ elszer˝ u kimenteni). @COP LU.FOR @COP LU.DAT @F LU R LU ED LU.RES Determin´ ans sz´ am´ıt´ asa hasonl´ oan t¨ ort´ enhet, mint a Gauss modszern´ el. M´ atrixinvert´ al´ as LU felbont´ assal A.X = E egyenletrendszerek megold´ asa r´ ev´ en
@COP LUINVERZ.FOR @COP LUINVERZ.DAT @F LUINVERZ R LUINVERZ ED LUINVERZ.RES V´ ege a direkt m´ odszereknek M´ arcius 8. Iter´ aci´ os m´ odszerek line´ aris egyenletrendszerek megold´ as´ ara. Ezeket a direkt m´ odszerrel kapott k¨ ozelit˝ o megold´ as pontos´ıt´ as´ ara haszn´ aljuk. K¨ ozel´ıt˝ o megold´ as:x0 Eredeti egyenlet: A.x = b A.x0 = b0 = b + c0 = A.x − (b − b0 )
Lin. egyenletrendszer az x − x0 korrekci´ ora: A.(x − x0 ) = b − b0 = c
Mivel ennek megold´ asa is numerikus, ez´ ert ezt is csak k¨ ozelit˝ oen kapjuk. Iter´ aci´ o a korrekci´ okra @COP MPROVE.FOR @COP MPROVE.DAT @F MPROVE ED MPROVE.RES El˝ osz¨ or k¨ ozel´ıt˝ o megold´ ast gener´ alunk, u ´ gy, hogy az LU felbont´ assal kapott eredm´ enyt kiss´ e elrontjuk. Ez lesz az x0 kezd˝ o megold´ asunk. Az iteraci´ ot addig folytatjuk, amig k´ et egymasut´ ani megold´ as elt´ er´ ese adott pontoss´ agn´ al kisebb lesz.
Nemline´ aris egyenletek megold´ asa. Ha f(x) folytonos ´ es sign(f (a)) 6= sign(f (b)) akkor l´ etezik gy¨ ok [a,b]-ben. Intervallum felez´ esi elj´ ar´ assal (bisection) Ez biztos, de lass´ u m´ odszer. P´ elda: f (x) = −1 − x ∗ (1 − x2 ) [a, b] = [1, 1.5]
acc=1.E-7
@cop rtbis.for @f rtbis r rtbis 1.,1.5,1.e-7
root is in the interval [ 1.00000, 1.50000] the root is: x= 1.32472+-0.1E-06 M´ asik p´ elda:f (x, c) = e−c∗x @cop rtbis1.for @cop khi2.dat @f rtbis1 r rtbis1 0.5,2.5,1.e-7 root is in the interval [ 0.50000, 2.50000] the root is: x= 0.95099,+-0.1E-06 H´ urm´ odszer, szel˝ o m´ odszer: line´ aris interpol´ aci´ o ´ Erint˝ o m´ odszer, Newton-Raphson m´ odszer: xk+1 = xk − f (xk )/f ′ (xk )
kvadratikusan konverg´ al, ha j´ o kezd˝ o´ erteket adunk. Nem mindig konverg´ al. Gener´ aljunk kezd˝ o´ ert´ eket az intervallum felez´ esi m´ odszerrel f (x, c) = e−c∗x
@cop nlinkhi2.for @f nlinkhi2 r nlinkhi2 the root is: 0.95099+-0.1E-5 FORTRAN STOP Nemline´ aris egyenletrendszerek megold´ asa Feladat: megoldand´ o a k¨ ovetkez˝ o k´ etismeretlenes nemline´ aris egyenletrendszer F1 (x1 , x2 ) = x21 + 3x1 − x2 − 1 = 0 F2 (x1 , x2 ) = x1 x2 + 3x2 + 9 = 0 β1 = −F1 β2 = −F2
α1,1 = ∂F1 ∂x1 = 2x1 + 3
α1,2 = ∂F1 ∂x2 = −1
α2,1 = ∂F2 ∂x1 = x2
α2,2 = ∂F2 ∂xx = x1 + 3
tehat a Jacobi matrix: α1,1 α1,2 J(x) = α2,1 α2,2 (0)
(0)
kezdoertekek: x1 = −4x2 = 6
Megoldas: x1 = −4.5251023x2 = 5.9012437
A t¨ obbparam´ eteres legkisebb n´ egyzetes illeszt´ esi feladat egyenletrendszerre vezet. Nemline´ aris egyenletrendszerre nincs biztosan m˝ uk¨ od˝ o m´ odszer.
Ha van elfogadhat´ o
kezd˝ o´ ert´ ek sorozatunk, akkor az ´ altal´ anos´ıtott Newton-Raphson m´ odszerrel pr´ ob´ alkozhatunk. Feladat: f (x, c(0 : 1)) = c(0) ∗ e−c(1)∗x illeszt´ ese az xi , yi alappontokra.
Ez k´ etismeretlenes nemline´ aris egyenletrendszerre vezet. @cop mnewt1.for @cop mnewt1.dat @f mnewt1 r mnewt1 FORTRAN STOP ed mnewt1.res *c M´ arcius 22 K¨ ozel´ıt˝ o (numerikus) differenci´ al´ as Differenciah´ anyadossal val´ o k¨ ozel´ıt´ es.
el˝ ore(progressz´ıv), h´ atra(retrogr´ ad) ´ es centr´ alis differenci´ ak. T´ abl´ azatosan (alappontokban) adott f¨ uggv´ eny deriv´ al´ asa: Interpol´ aci´ os polinommal, majd a polinom deriv´ al´ asa. Spline interpol´ aci´ oval pl. Ekkor a deriv´ alt m´ asodfok´ u polinomja x-nek, ez´ ert a magasabbrend˝ u deriv´ altak sz´ amol´ asa probl´ em´ as. K¨ ozel´ıt´ es Csebishev polinomokkal: SUBROUTINE CHEBFT(A,B,C,N,FUNC) A FUNC(X) f¨ uggv´ enyt k¨ ozel´ıti N-edfok´ u Csebisev polinommal, teh´ at kisz´ am´ıtja a C[1:N] sorfejt´ esi egy¨ utthat´ okat. Ezut´ an megvizsg´ aljuk az egy¨ utthat´ okat ´ es esetleg lev´ agjuk a sorfejt´ est. K¨ ozel´ıt˝ o f¨ uggv´ eny sz´ amol´ asa Clenshow rekurzi´ oval. FUNCTION CHEBEV(A,B,C,M,X) A deriv´ alt polinom egy¨ utthat´ oinak sz´ amol´ asa. Mi a teend˝ o, ha a f¨ uggv´ eny csak t´ abl´ azatosan adott? i) interpol´ aci´ o ii) interpol´ aci´ os polinom Csebisev sorfejt´ esi egy¨ utthat´ oinak meghat´ aroz´ asa. iii) a deriv´ alt f¨ uggv´ eny Csebisev sorfejt´ esi egy¨ utthat´ oinak sz´ amol´ asa.
SUBROUTINE CHDER(A,B,C,CDER,N) a deriv´ alt egy¨ utthat´ oi a CDER t¨ ombben @COP CHEBDIF.FOR @COP WERNER.NUM input file @F CHEBDIF R CHEBDIF ed chebdif.res *c deriv´ alt f¨ uggv´ eny ki´ır´ asa m´ eg egyszer a WERNER.DEN1 fileba is. L´ asd k´ es˝ obb. ´ Aprilis 12 Ha a Csebisev sorfejt´ esi egy¨ utthat´ okat ismerj¨ uk, akkor sz´ amolhatjuk az integr´ al sorfejt´ esi egy¨ utthat´ oit is. @COP CHEBINT.FOR @COP WERNER.DEN input file R CHEBINT ED CHEBINT.RES *C ED WERNER.NUM1 *C Ez k¨ ozel´ıt˝ oleg megegyezik a WERNER.NUM-mal, aminek deriv´ altja a werner.den1, amit el˝ oz˝ oleg meg˝ or´ızt¨ unk. Numerikus integr´ al´ as (Kvadrat´ ura) A Lagrange interpol´ aci´ os polinom kiintegr´ al´ asa Newton-Cotes formul´ ara vezet. Egyenl˝ ok¨ oz˝ u beoszt´ asok: trap´ ez-formula, hibatag f (x) = ex
@COP TRAGYAK.FOR a Simpson formula, mint a trap´ ez-formula javit´ oformul´ aja. @COP SIMP.FOR f (x) = sin(x)
Jav´ıt´ oformula a Simpson m´ odszerre. Feladat: ´ırjunk programot az f (x) = x3 f¨ uggv´ eny ekvidiszt´ ans pontokban vett Simpson k¨ ozel´ıt´ essel sz´ amolt integr´ alj´ anak ´ es becs¨ ult hib´ aj´ anak a kisz´ am´ıt´ as´ ara. megold´ as: SIMPGYAK.FOR Romberg integr´ al´ as @COP MAINROM.FOR @COP ROMB.DAT @F MAINROM R MAINROM ED ROMB.RES *C Gauss kvadrat´ ur´ ak:
V´ eges [a, b] intervallumon a Gauss-Legendre kvadrat´ ur´ at haszn´ aljuk. n-pontos kvadrat´ ur´ an´ al az oszt´ aspontok az n-edfok´ u Legendre polinom z´ erushelyei. Az alappontok ´ es a s´ ulyok megtal´ alhat´ ok pl. Abramowitz-Stegun, Handbook of Mathematical Functions, Table 25.4 916.oldal Pl. 10 pontos Gauss-Legendre kvadrat´ ura: @COP QGAUS.FOR Gauss kvadrat´ ur´ ak m´ as orthogon´ alis polinomok segits´ eg´ evel. Gauss Csebisev kvadrat´ ura, Gauss Hermite, Gauss Laguerre kvadrat´ ur´ ak. Hogyan sz´ am´ıtsuk az alappontokat ´ es a s´ ulyokat Gauss-Legendre kvadrat´ ur´ ara? Mintaprogram: @COP GAULEG.FOR @COP GAULEG.DAT input file a,b,n @F GAULEG R GAULEG ED GAULEG.RES *C GL.FOR fileb´ ol beszerkeszthetj¨ uk az alappontok ´ es a s´ ulyok pontos ´ ert´ ekeit a programba. ´ Aprilis 19 T¨ obbsz¨ or¨ os integr´ alok kisz´ am´ıt´ asa K´ etszeres integr´ alok I=
Z Z
f (x, y)dxdy, T
ahol T k´ etdimenzi´ os z´ art halmaz. Jel¨ olj¨ uk a = min(x,y)∈T {y} ´ es b = max(x,y)∈T {y}. Tegy¨ uk
fel, hogy b´ armely y0 ∈ [a, b]-re az y = y0 egyenes maximum k´ et ponban metszi T hat´ ar´ at.
Ha ez nem teljes¨ ul, akkor darabol´ assal alak´ıtjuk T -t olyan T1 , ..., Tk tartom´ anyok ¨ osszeg´ ere,
melyekre ez a felt´ etel m´ ar teljes¨ ul. Jel¨ olje az y = y0 egyenes balodali metsz´ espontj´ at a ϕ(y0 ), a jobboldali metsz´ espontj´ at ψ(y0 ). Ekkor I=
Z
b
F (y)dy,
F (y) =
a
Z
ψ(y)
f (x, y)dx ϕ(y)
Ezzel a k´ etszeres integr´ al kisz´ am´ıt´ as´ at egyszeres integr´ alok kisz´ am´ıt´ as´ ara vezett¨ uk vissza: I=
n X
(n)
F (yk )Ak
k=1
F (yk ) =
mk X
(k)
(mk )
f (xi , yk )Ai
,
i=1
ahol az mk ,
k = 1, ..., n alappont a [ϕ(yk ), ψ(yk )] intervallumon. Ezekkel kifejezve I=
mk n X X
k=1 i=1
(k)
(n)
(mk )
f (xi , yk )Ak Ai
.
A legegyszer¨ ubb esetben, amikor T t´ eglalap, akkor az x beli alappontrendszer f¨ uggetlen yk -t´ ol: mk = m
(k)
xi
= xi ,
i = 1, .., m I=
m n X X
(n)
(m)
f (xi , yk )Ak Ai
.
k=1 i=1
P´ eldak´ eppen a k´ tszeres integr´ al kisz´ am´ıt´ as´ ara, sz´ am´ıtsuk ki a f´ elg¨ omb t´ erfogat´ at numerikus integr´ al´ assal. V´ egezz¨ unk 10 pontos Gauss-Legendre kvadrat´ ur´ at mindk´ et integr´ alra: @COP QGAUS2D.FOR @F QGAUS2D R QGAUS2D L´ athat´ o, hogy a 10 pontos Gauss kvadrat´ ura nem el´ eg pontos.
Ha az integr´ al ´ ert´ ek´ et
pontos´ıtani akarjuk, akkor t¨ obbpontos kvadrat´ ur´ ara van sz¨ uks´ eg¨ unk. 16 pontos Gauss kvadrat´ ur´ at haszn´ al a QGAUS2D1.FOR 32 pontos Gauss kvadrat´ ur´ at haszn´ al a QGAUS2D2.FOR Integr´ al´ as t¨ obb (d) dimenzi´ oban. Monte-Carlo integr´ al´ as T¨ obbsz¨ or¨ os integr´ al´ as ´ altal´ anos´ıthat´ o d > 2 esetre is, de a sz¨ uks´ eges id˝ o rohamosan n˝ o a d dimenzi´ o n¨ oveked´ es´ evel. A hiba becsl´ ese is egyre bonyolultabb. Ez´ ert ´ erdemes lesz a sztochasztikus szimul´ aci´ ot haszn´ al´ o u.n. Monte-Carlo m´ odszert haszn´ alni. Ebben az oszt´ aspontokat v´ eletlenszer˝ uen v´ alasztjuk. Ha az ri sz´ amok a [0,1] intervallumon egyenletes eloszl´ as´ u v´ eletlen sz´ amok, akkor az I=
Z
1
... 0
Z
1
F (x1 , x2 , ..., xn )dx1 dx2 ... dxn 0
integr´ al k¨ ozel´ıt´ ese az F (r1 , r2 , ..., rn ) f¨ uggv´ eny´ ert´ ekekb˝ ol k´ epezett ¨ osszeg. El˝ onye m´ eg, hogy egyszer˝ uen alkalmazhat´ o t¨ obbdimenzi´ os komplik´ alt integr´ aci´ os tartom´ anyra is.
Fel kell
venni azt a legkisebb V hipert´ eglalapot, ami a T = W integr´ aci´ os tartom´ anyt k¨ or¨ ul¨ oleli. Az integr´ al´ ast erre a V hipert´ eglalapra kell elv´ egezni, u ´ gy, hogy a W -n k´ıv¨ uli pontokat nem adjuk hozz´ a az ¨ osszeghez. Ez ekvivalens egy olyan f (xi ) f¨ uggv´ eny defini´ al´ as´ aval, ami azonosan zerus, hacsak xi ∈ W nem teljes¨ ul, mikoris f (xi ) = F (xi ) ´ all fenn. Az I integr´ al Monte-Carlo k¨ ozel´ıt´ esben N esem´ eny gener´ al´ asa mellett: Z V p < f 2 > − < f >2 , f (x)dV = V < f > ± √ I= N V
ahol dV = dx1 ...dxn , ´ es < f >=
< f 2 >=
Mivel a hiba
√1 -nel N
N 1 X f (xi ) N i=1 N 1 X f (xi )2 . N i=1
ar´ anyos, aj´ anlatos a V hipert´ eglalapot a lehet˝ o legkisebbre venni, mert
f zerus ´ ert´ ekei n¨ ovelik a hibatagot, u.i. a W -n k´ıv¨ uli pontoknak nincs inform´ aci´ o tartalma, 1 ez´ ert az esem´ enyek hat´ asos Nef f sz´ ama N -n´ el kisebb, ´ es a hiba √ -vel ar´ anyos. Nef f
V´ eletlensz´ amok gener´ al´ asi m´ odszerei szam´ıt´ og´ eppel:
A sz´ am´ıt´ og´ ep nem k´ epes igazi v´ eletlen sz´ amokat gener´ alni, csak ehhez k¨ ozel´ıt˝ o u.n. pseudo v´ eletlensz´ amokat. Igen fontos, hogy ezek eloszl´ asa a lehet˝ o legjobban k¨ ozel´ıtse az egyenletes eloszl´ ast. Neumann J´ anos f´ ele k¨ ozbens˝ o-n´ egyzet m´ odszer: r-jegy˝ u sz´ am k¨ oz´ eps˝ o r/2 jegye D.H. Lehmer szorzat-marad´ ek m´ odszer: Adva van egy a szorz´ o ´ es egy m oszt´ o (eg´ eszek), valamint egy ugyancsak eg´ esz r0 kezd˝ o´ ert´ ek. A tov´ abbi elemeket az ri = a × ri−1 (mod m)
rekurzi´ oval ´ all´ıtjuk el˝ o. M´ odos´ıtott form´ aja, ha van m´ eg egy b konstans is, ´ es a rekurzi´ o
ri = (a × ri−1 + b)(mod m) alak´ u, ahol m = 2t t az integer sz´ o bitjeinek sz´ ama. Szorz´ as 2 × t
jegy˝ u, a hozz´ aadott b t jegy˝ u, a k¨ ovetkez˝ o v´ eletlen sz´ am az eredm´ eny kev´ esb´ e szignifik´ ans
t jegye, amit m-mel val´ o oszt´ assal alak´ıtunk [0,1] k¨ oz´ e es˝ o lebeg˝ opontos sz´ amm´ a.
Ahrens-Dieter m´ odszer: ri = (a × ri−1 + b × ri−2 )(mod m) NAG Lib. SUBROUTINE G05CAF m = 259
a = 1313
r0 = 123456789 × (232 + 1)
Az eddigi m´ odszerek k¨ oz¨ os h´ atr´ anya az, hogy korrel´ aci´ o van az egym´ ast k¨ ovet˝ o sz´ amok k¨ oz¨ ott. Jav´ıt´ as m´ odja: el˝ osz¨ or gener´ alunk egy sz´ amsorozatot, majd ebb˝ ol v´ alasztunk egy elemet v´ eletlen¨ ul, majd ennek a hely´ et felt¨ oltj¨ uk egy u ´ jonnan gener´ alt v´ eletlen sz´ ammal. ´Igy m˝ uk¨ odik a RAN2 function. FUNCTION RAN2(IDUM) c c [0,1] intervallumon egyenletes eloszlasu veletlen szamokat generalo rutin c PARAMETER (M=714025,IA=1366,IC=150889,RM=1.4005112E-6) DIMENSION IR(97) DATA IFF /0/ IF(IDUM.LT.0.OR.IFF.EQ.0)THEN IFF=1 IDUM=MOD(IC-IDUM,M) DO 11 J=1,97 IDUM=MOD(IA*IDUM+IC,M) IR(J)=IDUM 11
CONTINUE IDUM=MOD(IA*IDUM+IC,M) IY=IDUM ENDIF J=1+(97*IY)/M IF(J.GT.97.OR.J.LT.1)PAUSE IY=IR(J) RAN2=IY*RM IDUM=MOD(IA*IDUM+IC,M) IR(J)=IDUM RETURN END
A Monte-Carlo integr´ al´ as sokdimenzi´ os integr´ alok kisz´ am´ıt´ as´ ara gazdas´ agos.
P´ elda Monte-Carlo integr´ al´ asra: t´ oruszdarab t´ erfogat´ anak sz´ am´ıt´ asa: Hat´ arozzuk meg egy 4 egys´ eg k¨ uls˝ o ´ es 2 egys´ eg bels˝ o sugar´ u t´ orusz ´ es egy x ≥ 1, y ≥ −3
felt´ etelekkel adott doboz k¨ oz¨ os r´ esz´ et k´ epez˝ o t´ argy s´ ulypontj´ at ´ es t¨ omeg´ et.
A t´ orusz
belsej´ ebe es˝ o pontra teljes˝ ul, hogy z2 + (
p
x2 + y 2 − 3)2 ≤ 1.
Tegy¨ uk fel, hogy a t´ argy ρ s˝ ur˝ us´ ege ´ alland´ o.
A feladat megold´ as´ ahoz a k¨ ovetkez˝ o in-
tegr´ alokat kell a test belsej´ ere kisz´ am´ıtani: Z Z Z M= ρdV, Mx = xρdV, My = W
W
yρdV, Mz = W
Z
zρdV, W
ahol dV = dx dy dz . A s´ ulypont koordin´ at´ ai a fenti integr´ alokkal: xs =
My Mz Mx , ys = , zs = . M M M
A testet mag´ abafoglal´ o V -vel jel¨ olt t´ eglatest: x ∈ [1, 4], y ∈ [−3, 4], z ∈ [−1, 1], V = 3 × 7 × 2 = 42.
@COP MONTC.FOR @F MONTC R MONTC adja meg az esem´ enyek sz´ am´ at 100000 n= 100000 a test t¨ omege= 22.1579 hiba= 0.0663 a s´ ulypont x,y,z koordin´ at´ ai: x= 2.4076 hiba= 0.0079 y= 0.1629 hiba= 0.0080 z= 0.0002 hiba= 0.0022 FORTRAN STOP Szimmetria okokb´ ol a s´ ulypont z koordin´ at´ aja 0, s l´ athat´ o, hogy a k¨ ozel´ıt˝ o ´ ert´ ek a hibahat´ aron bel¨ ul ezt teljes´ıti. H´ azi feladatk´ ent hat´ arozzuk meg Monte-Carlo m´ odszerrel az egys´ egnyi sugar´ u f´ elg¨ omb t´ erfogat´ at. @cop halfsphmc @f halfsphmc r halfsphmc adja meg az esemenyek szamat 100000 n= 100000 a test tomege= 2.0982 hiba= 0.0063 a sulypont x,y,z koordinatai: x= 0.0020 hiba= 0.0019 y= -0.0010 hiba= 0.0019
z= 0.3753 hiba= 0.0016 felgomb terfogata= 2.09440 FORTRAN STOP egzakt ´ erte´ ekek:
2π 3
= 2.09440, valamint xs = ys = 0. L´ athat´ o, hogy a k¨ ozel´ıt˝ o ´ ert´ ekek a
hibahat´ arokon bel¨ ul egyeznek az egzakt ´ ert´ ekekkel.
Differenci´ alegyenletek numerikus megold´ asa Differenci´ alegyenletek oszt´ alyoz´ asa: - k¨ oz¨ ons´ eges differenci´ alegyenletek (egyv´ altoz´ os f¨ uggv´ enyekre: y(x)), - parci´ alis differenci´ alegyenlet(t¨ obbv´ altoz´ os f¨ uggv´ enyekre(pl. u(x, y)) El˝ osz¨ or a k¨ oz¨ ons´ eges differenci´ alegyenletekkel foglalkozunk. Alapfeladatok k¨ oz¨ ons´ eges differenci´ alegyenletekre: - Kezdeti´ ert´ ek-probl´ ema: a differenci´ alegyenlet megold´ asa (´ es a deriv´ altjai) adott(ak) az x0 = a kezd˝ opontban, meghat´ arozand´ o bizonyos diszkr´ et xi pontokban az eg´ esz xi ∈ [a, b] intervallumon.
- Perem´ ert´ ek-probl´ ema: a differenci´ alegyenlet megold´ asa adott az [a, b] intervallum a kezd˝ o´ es b v´ egpontj´ aban ( m´ asodrend˝ u, vagy magasabbrend˝ u differenci´ alegyenletekre). Pl. Ly(x) = y ′′ (x) + p(x)y ′ (x) + q(x)y(x) = 0,
ahol az y(a) = A ´ es az y(b) = B perem´ ert´ ekek adottak, az y , (a) = α azonban ismeretlen. A perem´ ert´ ek-probl´ ema visszavezethet˝ o kezdeti´ ert´ ek-probl´ em´ ara. - Saj´ at´ ert´ ek-probl´ ema: Ismeretlen param´ eter van a differenci´ alegyenletben. Saj´ at´ ert´ ek: azon param´ eter ´ ert´ ek, amin´ el a perem´ ert´ ek-probl´ em´ anak van megold´ asa.
Ez a
megold´ as a saj´ atf¨ uggv´ eny. Pl. Ly − λy = 0,
ahol λ ismeretlen param´ eter. Az y(a) = A ´ es az y(b) = B perem´ ert´ ekek adottak. Ez is visszavezethet˝ o a kezdeti´ ert´ ek-probl´ em´ ak, vagy matrix saj´ at´ ert´ ek probl´ em´ aj´ anak megold´ as´ ara.
Kezdeti´ ert´ ek probl´ em´ ak numerikus megold´ asi m´ odszerei
A legt¨ obb differenci´ alegyenletb˝ ol kifejezhet˝ ok a deriv´ altak, ´ es a differenci´ alegyenlet fel´ırhat´ o u.n. standard alakba, ami a k¨ ovetkez˝ o: y′ = f (x, y), ahol y ´ es f N -dimenzi´ os vektorok. A magasabb rend˝ u differenci´ alegyenletek visszavezethet˝ ok a standard alakra, azaz els˝ orend˝ u differenci´ alegyenletek rendszer´ ere. Pl. az y ′′ + py ′ + qy = 0
m´ asodrend˝ u differenci´ alegyenlet eset´ en vezess¨ uk be: az y1 = y ,´ es az y1′ = f1 = y2 = y ′ jel¨ ol´ eseket.
Ezek felhaszn´ al´ as´ aval y2′ = f2 = y ′′ = −py2 − qy1 , azaz a m´ asodrend˝ u
differenci´ alegyenlet a standard alakot veszi fel N = 2-vel , ugyanis y′ = (y1′ , y2′ )T = (f1 , f2 )T ,
ahol f1 = y2 ´ es f2 = −py2 − qy1 , ´ es a T a transzpon´ altat jel¨ oli.
A numerikus megold´ asi m´ odszerek oszt´ alyoz´ asa: A differenci´ alegyenletek megold´ asi m´ odszereit szok´ as felosztani egyl´ ep´ eses ´ es t¨ obbl´ ep´ eses m´ odszerekre. Vegy¨ unk egyenl˝ o, h = (b−a)/n l´ ep´ esk¨ oz˝ u beoszt´ as´ at az [a, b] intervallumnak, azaz legyen xi = a + ih, i = 0, 1, ...n.
- t¨ obbl´ ep´ eses m´ odszerr˝ ol akkor besz´ el¨ unk, ha az yi+1 kisz´ am´ıt´ as´ ahoz egyn´ el t¨ obb oszt´ aspontbeli megold´ ast haszn´ alunk fel, azaz k > 0 mellett felhaszn´ aljuk az yi , yi−1 , ..., yi−k ´ ert´ ekeket. - egyl´ ep´ eses a m´ odszer, ha k=0, teh´ at csak az yi ´ ert´ eket haszn´ aljuk fel. K¨ oz¨ ons´ eges differenci´ alegyenletek megold´ asa egyl´ ep´ eses k¨ ozel´ıt˝ o m´ odszerrel: az Euler m´ odszer. ′
Tekints¨ uk az y = f (x, y) differenci´ alegyenlet-rendszer x0 pontbeli kezdeti´ ert´ ek-probl´ em´ aj´ at, ahol adottak az y(x0 ) = y0 ´ es y′ (x0 ) = f (x0 , y0 ) kezdeti´ ert´ ek vektorok ´ ert´ ekei. Sz´ am´ıtsuk ki a differenci´ alegyenlet megold´ as´ at az x1 = x0 + h pontban. A legdurv´ abb k¨ ozel´ıt´ es az els˝ orend˝ u u.n. Euler m´ odszer, amiben az [x0 , x1 ] intervallumban a keresett megold´ ast egyenesekkel k¨ ozel´ıtj¨ uk. Ha N = 1, azaz csak egyetlen els˝ orend˝ u differenci´ alegyenlettel van dolgunk, akkor y1 = y0 + hy ′ (x0 ) = y0 + hf (x0 , y0 )
vagy ezt ´ atrendezve,
∆y0 y1 − y0 = = y ′ (x0 ) = f (x0 , y0 ). h h
Vagyis a differenciah´ anyadost az intervallum kezd˝ opontjabeli deriv´ alttal helyettes´ıtett¨ uk. Ez a k¨ ozel´ıt´ es nyilv´ anval´ oan ann´ al rosszabb, min´ el ink´ abb elt´ er az y(x) megold´ as az x ∈ [x0 , x1 ] intervallumban az egyenest˝ ol, vagyis, min´ el nagyobb a h l´ ep´ esk¨ oz. Az Euler
m´ odszer teh´ at val´ odi egyl´ ep´ eses m´ odszer, hiszen az y1 kisz´ am´ıt´ as´ ahoz csak az x0 pontbeli
mennyis´ egeket y0 -at ´ es f (x0 , y0 )-at haszn´ alja fel. Ha a m´ odszert az x0 helyett az xi kezd˝ opontb´ ol ind´ıtjuk, ´ es t¨ obb egyenlet¨ unk van (N > 1), akkor az Euler m´ odszer term´ eszetesen yi+1 = yi + hf (xi , yi )
alak´ u lesz. P´ elda az Euler m´ odsszer haszn´ alat´ ara: x0 = 0.d0, y0 = 0.d0, b = 1.d0 ´ es h = 0.1
@COP EULER.FOR @COP EULER.DAT @F EULER R EULER FORTRAN STOP ED EULER.RES *C y(x=1)=0.29254 , ha h=0.1 y(x=1)=0.34402 , ha h=0.01 y(x=1)=0.34961 , ha h=0.001
N = 1 y ′ = f (x, y) = x2 + y 2 adottak:
y(x=1)=0.35040 , ha h=0.0001 Ahhoz, hogy az els˝ orend˝ u Euler m´ odszern´ el pontosabb eredm´ enyt sz´ am´ıthassunk, k¨ ozbens˝ o pontokat kell felvenn¨ unk az [xi , xi+1 ] intervallumban. Jel¨ olj¨ uk ezeket a k¨ ozbens˝ o pontokat x ˜j -vel. Vegy¨ unk t darab k¨ ozbens˝ o pontot, jel¨ olj¨ uk az ezekben a pontokban vett ˜ j -vel ´ k¨ ozel´ıt˝ o megold´ asokat y es az itteni k¨ ozel´ıt˝ o deriv´ altakat kj -vel: ˜ j ). kj = f (˜ xj , y
Term´ eszetesen j=1,..,t. Az
∆yi yi+1 − yi = h h
differenciah´ anyadost pontosabban k¨ ozel´ıthetj¨ uk a k¨ ozel´ıt˝ o deriv´ altak s´ ulyozott ¨ osszeg´ evel, ha az αj s´ ulyokat megfelel˝ oen v´ alasztjuk: t
X yi+1 − yi ˜ j ) ≡ g(xi , yi , h). αj f (˜ xj , y = h j=1
A most k¨ ovetkez˝ o Runge-Kutta tipus´ u m´ odszerek u ´ gy hat´ arozz´ ak meg az αj s´ ulyokat, valamint az x ˜j y˜j ´ ert´ ekeket, hogy az adott t-hez tartoz´ o m´ odszer H hibavektora a h-ban p-edrend˝ u legyen: H(xi , h) ≡ y(xi ) + hg(xi , y(xi ), h) − y(xi + h) = O(hp+1 )1,
ahol 1 egy olyan N -dimenzi´ os vektor, melynek minden komponense 1.
Runge-Kutta tipus´ u m´ odszerek. V´ alasszunk egyenl˝ ok¨ oz˝ u beoszt´ ast az [x0 , b] intervallumon: xi = x0 + ih,
ahol i = 0, 1, ..., n.
Oldjuk meg az y′ = f (x, y) differenci´ alegyenlet-rendszer x0 pontbeli
kezdeti´ ert´ ek-probl´ emaj´ at, ahol adottak az y(x0 ) = y0 ´ es y, (x0 ) = f (x0 , y0 ) ´ ert´ ekei. A megold´ as meghat´ aroz´ as´ at x0 -b´ ol kiind´ ulva, h hossz´ us´ ag´ u l´ ep´ esenk´ ent v´ egezz¨ uk.
Tegy¨ uk fel, hogy
m´ ar megtett¨ unk k l´ ep´ est ´ es meghat´ aroztuk az xk pontbeli yk k¨ ozel´ıt˝ o megold´ ast. Ekkor a soronk¨ ovetkez˝ o feladat, ennek felhaszn´ al´ as´ aval kisz´ am´ıtani az xk+1 = xk + h pontbeli yk+1 k¨ ozel´ıt˝ o megold´ ast. Vegy¨ unk fel ehhez t darab x ˜ j = x k + µj h
k¨ ozbens˝ opontot az [xk , xk+1 ] intervallumban, azaz az egyel˝ ore ismeretlen µj param´ eterekre k¨ oss¨ uk ki, hogy µj ∈ [0, 1]. Vezess¨ unk be tov´ abbi ismeretlen λji param´ etereket, melyekkel meghat´ arozzuk a k¨ ozbens˝ opontbeli (k¨ ozel´ıt˝ o) megold´ asokat ´ es deriv´ altakat ˜ j = yk + h y
t X
λji ki
i=1
˜ j ) = f (xk + µj h, yk + h kj = f (˜ xj , y
t X i=1
λji ki ).
Ha a λji param´ etereket u ´ gy v´ alasztjuk meg, hogy i ≥ j eset´ en λji = 0 legyen, akkor a j -t 1-t˝ ol t-ig n¨ ovelve, a
kj = f (˜ xj , yk + h
j−1 X
λji ki )
i=1
kifejez´ es mindig csak a m´ ar kisz´ am´ıtott kj -ket fogja tartalmazni, teh´ at ebb˝ ol a k1 ,k2 ,...,kt ´ ert´ ekeket egym´ asut´ an ki tudjuk sz´ am´ıtani.
A xk+1 pontbeli k¨ ozel´ıt˝ o megold´ as ezek is-
meret´ eben a k¨ ovetkez˝ ok´ eppen sz´ am´ıthat´ o: yk+1 = yk + h
t X
αj k j .
j=1
Id˝ ohi´ any miatt itt nem foglalkozunk azzal, hogy milyen egyenletekb˝ ol lehet meghat´ arozni a Runge-Kutta m´ odszerekben szerepl˝ o param´ eterek ´ ert´ ekeit, csup´ an n´ eh´ any fontosabb speci´ alis esettel foglalkozunk.
A speci´ alis eseteket a benn¨ uk szerepl˝ o k¨ ozbens˝ o pontok
sz´ ama: t, vagy a k¨ ozel´ıt´ es p rendje szerint csoportos´ıthatjuk. Ha t = 1, akkor csak egyetlen pont van, ami az intervallum kezd˝ opontja: x ˜ 1 = xk , λ11 = 0, ´ıgy µ1 = 0 ´ es k1 = f (xk , yk ) miatt, valamint α1 = 1 felhaszn´ al´ as´ aval az Euler m´ odszert
kapjuk vissza: yk+1 = yk + hf (xk , yk ).
Az Euler m´ odszer els˝ orend˝ u Runge-Kutta m´ odszer (p = 1), melynek hibatagja O(h2 ) (teh´ at h-ban m´ asodrend˝ u).
Ez igen pontatlan m´ odszer, ez´ ert a gyakorlatban nem is c´ elszer˝ u
haszn´ alni. Gyakorlatban haszn´ alatos m´ odszer a negyedrend˝ u (p = 4) Runge-Kutta m´ odszer. Itt t=4´ es a k¨ ozbens˝ oponti k¨ ozel´ıt˝ o deriv´ atakat a k¨ ovetkez˝ ok´ eppen sz´ am´ıtjuk: k1 = f (xk , yk ) h , yk + 2 h k3 = f (xk + , yk + 2 k2 = f (xk +
h k1 ) 2 h k2 ) 2
k4 = f (xk + h, yk + hk3 ).
Ezek seg´ıts´ eg´ evel az xk+1 pontbeli megold´ as: yk+1 = yk +
h (k1 + 2k2 + 2k3 + k4 ). 6
A negyedrend˝ u Runge-Kutta m´ odszer hibatagja: O(h5 ) , teh´ at egzakt maximum negyedfok´ u polinomra. Alkalmazzuk a negyedrend˝ u Runge-Kutta m´ odszert az k¨ ovetkez˝ o feladat megold´ as´ ara: N =1 y′ = 2 ∗
p
1 − y 2 y(0) = sin(0.25) egzakt megold´ asa a k.´ e.p.-nak: y = sin(2x + 0.25)
@COP RK4.FOR @COP RK4.DAT @F RK4 R RK4
FORTRAN STOP
ED RK4.RES *C ED RK4.FOR *C Az f (x, y) jobboldalvektor komponenseit a SUBROUTINE DERIVS sz´ amolja.
M´ asik
feladatra ezt kell ´ at´ırni. H´ azi feladat: oldjuk meg az y ′ = x2 + y 2 kezdeti´ ert´ ek-probl´ ema-t y(0) = 0-ra az RK4 felhaszn´ al´ as´ aval. a megold´ as: RK41.FOR ED RK41.res *C x approximate y(x) 0.000000000 0.000000000 0.100000000 0.000333335 0.200000000 0.002666875 0.300000000 0.009003498 0.400000000 0.021359447 0.500000000 0.041791288 0.600000000 0.072448125 0.700000000 0.115660305 0.800000000 0.174081004 0.900000000 0.250907869 1.000000000 0.350233742 A negyedrend˝ u Runge-Kutta m´ odszer alkalmaz´ asa m´ asodrend˝ u differenci´ alegyenlet-re: y ′′ = −y ∗ K 2 ,
y(0) = 1 , y ′ (0) = 0 kezdeti´ ert´ ek-probl´ ema megold´ asa:
Els˝ o l´ ep´ es a standard alakra hoz´ as ´ es az ennek megfelel˝ o DERIVS subroutine elk´ esz´ıt´ ese. @COP RK24.FOR c main program for calling rk4 c integration y”=-xk**2*y c with y(0)=1 and y’(0)=0 c c exact solution: y(x)=cos(xk*x) implicit real*8 (a-h,o-z) parameter (nmax=1000) dimension yy(nmax),xx(nmax),f(nmax),dev(nmax), 1 y(2),yout(2),dydx(2) common/par/xk external derivs open(unit=5,file=’rk24.dat’,status=’old’) open(unit=6,file=’rk24.res’,status=’unknown’) read(5,*) n,xmax,xk c number of equations nn=2 c converting second order diff. eqn into the system of first order diff.eqns
c y(1)=y(x), y(2)=y’(x)=dydx(1) y”(x)=dydx(2)=-xk**2*y(1) c initial values: at x=0 y(1)=1 y(2)=0 c c setting initial condition x=0.d0 y(1)=1.d0 y(2)=0.d0 c h=xmax/n xx(1)=x yy(1)=y(1) c exact solution in x=h f(1)=cos(xk*x) c step by step integration from x=0 to x=xmax in step size h do i=1,n ip=i+1 c define k1 = f (x, y) call derivs(x,y,dydx) call rk4(y,dydx,nn,x,h,yout,derivs) x=x+h xx(ip)=x c store approximate solution yy(ip)=yout(1) y(1)=yout(1) y(2)=yout(2) c store exact solution and deviation f(ip)=cos(xk*x) dev(ip)=abs((yy(ip)-f(ip))/f(ip)) enddo write(6,100) write(6,101)(xx(i),yy(i),f(i),dev(i),i=1,n+1) 100 format(7x,’x’,10x,’approximate y(x)’,3x,’exact y(x)’,4x, 1 ’relative error’) 101 format(1x,3f15.9,f19.16) stop end subroutine derivs(x,y,dydx) implicit real*8 (a-h,o-z) dimension y(2),dydx(2) common/par/xk c
yi′ (x)
= dydxi = fi (x, y1 , y2 )i = 1, 2
dydx(1)=y(2) dydx(2)=-xk**2*y(1) return
end @COP RK24.DAT TYPE RK24.DAT 10,1.8,2.d0 @F RK24 R RK24 TYPE RK24.RES x approximate y(x) exact y(x) relative error 0.000000000 1.000000000 1.000000000 0.0000000000000000 0.180000000 0.935899840 0.935896824 0.0000032229215751 0.360000000 0.751846764 0.751805729 0.0000545821793318 0.540000000 0.471434530 0.471328364 0.0002252473826228 0.720000000 0.130606600 0.130423709 0.0014022840479742 0.900000000 -0.226951116 -0.227202095 0.0011046501167430 1.080000000 -0.555409741 -0.555699146 0.0005207944531408 1.260000000 -0.812671410 -0.812952037 0.0003451947985342 1.440000000 -0.965764864 -0.965979312 0.0002220007929906 1.620000000 -0.995071125 -0.995161903 0.0000912195880846 1.800000000 -0.896837674 -0.896758416 0.0000883827341111 A negyedrend˝ u Runge-Kutta m´ odszer jav´ıt´ o formul´ aja. Jel¨ olje yh a h l´ eo´ esk¨ oz hszn´ alat´ aval, y2h pedig a 2h l´ ep´ esk¨ oz haszn´ alata mellett kapott k¨ oze´ıt˝ o megold´ ast az x = xk pontban.
Kimutathat´ o, hogy ugyanitt az y=
1 16 yh − y2h 15 16
jav´ıt´ o formula jobb k¨ ozel´ıt´ ese az y(xk ) egzakt megold´ asnak, mint az yh . Ezt a jav´ıt´ o formul´ at teh´ at felhszn´ alhatjuk a h l´ ep´ esk¨ ozzel kapott k¨ ozel´ıt´ es tov´ abbjav´ıt´ as´ ara.
Ezt azonban
ritk´ an alkalmazzuk. Ink´ abb haszn´ aljuk a jav´ıt´ oformul´ ab´ ol sz´ armaztathat´ o hibabecsl´ est: a h l´ ep´ sk¨ ozzel kapott yh k¨ ozel´ıt´ es hib´ aja eszerint δ =
1 15 (yh
− y2h ). Ennek seg´ıts´ eg´ evel a h
integr´ al´ asi l´ ep´ esk¨ ozt be´ allithatjuk egy olyan ´ ert´ ekre, hogy a k¨ ozel´ıt´ es megadott pontoss´ ag´ u legyen. Automatikus l´ ep´ esk¨ oz-szab´ alyoz´ as adott pontoss´ ag´ u k¨ ozel´ıt´ eshez
B´ ar az automatikus l´ ep´ esk¨ oz-szab´ alyoz´ ast az el˝ obbi hibabecsl´ es seg´ıts´ eg´ evel is elv´ egezhetn´ enk, ¨ tekints¨ unk itt egy magasabbrend˝ u k¨ ozel´ıt´ est. Ot¨ odrend˝ u (p = 5), vagy be´ agyazott RungeKutta formula: Ez a k¨ ozel´ıt´ es hat k¨ ozbens˝ o pontot haszn´ al (t = 6). Jel¨ olj¨ uk az ezekben a pontokban vett k¨ ozel´ıt˝ o deriv´ altakat ki -vel (i = 1, ..., 6). Mivel ez a k¨ ozel´ıt´ es ¨ ot¨ odrend˝ u, a m´ odszer hibatagja: O(h6 ), felt´ eve, hogy megfelel˝ ok´ eppen (optim´ alisan) v´ alasztjuk az αi s´ ulyokat: yk+1 = yk + h
6 X
αi ki + O(h6 ).
i=1
Ha az optim´ alis αi s´ ulyok helyett kev´ esb´ e optim´ alis αi∗ s´ ulyokat haszn´ alunk, akkor egy kev´ esb´ e pontos, negyedrend˝ u k¨ ozel´ıt´ est kapunk: ∗ yk+1 = yk + h
6 X i=1
αi∗ ki + O(h5 ).
A k´ et k¨ ozel´ıt´ es k¨ ul¨ onbs´ eg´ et k´ epezve hibabecsl´ est kaphatunk yk+1 -re : ∆ = yk+1 −
∗ yk+1
=
6 X i=1
ki (αi − αi∗ ) = O(h5 ).
Itt term´ eszetesen az αi valamint az αi∗ konstansok megadott ´ ert´ ekek. A hibabecsl´ es lehet˝ ov´ e teszi, hogy adott pontoss´ ag´ u k¨ ozel´ıt´ est ´ all´ıtsunk el˝ o a h integr´ al´ asi l´ ep´ esk¨ oz automatikus szab´ alyoz´ as´ aval. Ehhez kihaszn´ aljuk, hogy a −6
tekinthet˝ o. Tegy¨ uk fel, hogy ǫ = 10
∆ h5
j´ o k¨ ozel´ıt´ essel h-t´ ol f¨ uggetlen ´ alland´ onak
pontoss´ ag´ u k¨ ozel´ıt˝ o megold´ ast akarunk gener´ alni.
Ekkor, ha egy h1 l´ ep´ esk¨ ozzel sz´ amolva a kapott hiba ´ ert´ eke ∆1 , akkor az ǫ = ∆o hib´ at h o = 1
ep´ esk¨ ozzel sz´ amolva ´ erhetj¨ uk el. Ha teh´ at |∆1 | > ǫ, akkor a sz´ amol´ ast ho l´ ep´ esk¨ ozzel h1 | ∆ǫ1 | 5 l´ meg kell ism´ eteln¨ unk. Ellenkez˝ o esetben nincs sz¨ uks´ eg a sz´ amol´ as megism´ etl´ es´ ere, mert az nem baj, ha a k¨ ozel´ıt´ es¨ unk pontosabb a megk´ıv´ antn´ al. A m´ odszer alkalmaz´ as´ ara tekints¨ uk a kor´ abbi p´ eld´ at: y ′′ = −y ∗ K 2 ,
y(0) = 1,
y ′ (0) = 0
@COP ADAPT2.FOR @COP ADAPT2.DAT @f ADAPT2 R ADAPT2 FORTRAN STOP ED ADAPT2.RES *C
x1 x2 approximate y(x) exact y(x) 0.000000000 1.800000000 -0.896758490 -0.896758416 L´ athat´ o, hogy a k¨ ozel´ıt˝ o megold´ as val´ oban ǫ = 10−8 pontoss´ aggal egyezik az egzakt eredm´ ennyel. Az ADAPT2 programban a SUBROUTINE RKQS v´ egzi az automatikus l´ ep´ esk¨ oz be´ all´ıt´ ast a kiv´ ant eps pontoss´ ag el´ er´ es´ ehez. A komponensek pontoss´ ag´ at YSCAL sk´ al´ azza. RKQSnek meg kell adni aktu´ alis param´ eterk´ ent azt az egyl´ ep´ eses numerikus integr´ al´ asi rutin nev´ et, amit haszn´ alni akarunk. Jelen esetben ez az RKCK t=5 R-K m´ odszer. Az x-ben val´ o el˝ orehalad´ ast az intervallum v´ egpontj´ aig az ODEINT subroutine szervezi. Az egyl´ ep´ eses m´ odszerekhez sorolhatjuk m´ eg az u.n. m´ odos´ıtott k¨ ozbens˝ oponti m´ odszert (MMID). Ennek l´ enyege az, hogy a H hossz´ us´ ag´ u egyetlen l´ ep´ est h =
H M
hossz´ us´ ag´ u elemi
l´ ep´ esekre bontjuk, s ´ıgy M +1 darab k¨ ozbens˝ opontot vezet¨ unk be. Ennek tov´ abbfejleszt´ ese a Bulirsch-Stoer m´ odszer, ami Richardson extrapol´ aci´ ot haszn´ al (mint a Romberg integr´ al´ as). A h=0-ra val´ o extrapol´ al´ ast a SUBROUTINE RZEXTR v´ egzi racion´ alis extrapol´ aci´ oval. Ennek a m´ odszernek a m˝ uk¨ od´ es´ et az al´ abbi p´ eldaprogram mutatja. @COP BESTADAPT.FOR adatfile a m´ ar ´ atm´ asolt ADAPT2.DAT A feladat is ugyanaz, mint az ADAPT2.FOR-ban specifik´ alt feladat, teh´ at y ′′ (x) = −k 2 ∗ y(x) , y(0) = 1, y ′ (0) = 0 egzakt m.o.: y = cos(kx)
A standard alak ism´ et: y1 = y
y2 = y ′ = y1′ = f1
ez´ ert a DERIVS ugyanaz, mint kor´ abban.
y2′ = y ′′ = −k 2 ∗ y1 = f2 ,
F˝ oprogram: I/O fileok megnyit´ asa, beolvas´ as, kezdeti´ ert´ ekek megad´ asa: ystart(1) = 1
ystart(2) = 0
A m´ odos´ıtott k¨ ozbens˝ oponti m´ odszer: A megold´ as adott X-ben:
Y(1:N), ´ es keresett X+H-ban, amit majd YOUT(1:N)-ben
t´ arolunk. A k¨ ozbens˝ o pontok sz´ ama: NSTEP, s ezek egyenl˝ ok¨ ozzel vannak elhelyezve: h = H/N ST EP l´ ep´ esk¨ ozzel. A k¨ ozel´ıt˝ o megold´ asok a k¨ ozbens˝ o pontokban a k¨ ovetkez˝ oek: z0 = y(X) z1 = z0 + h ∗ f (X, z0 ) z2 = z0 + 2h ∗ f (X + h, z1 ) z3 = z1 + 2h ∗ f (X + 2h, z1 ) .. zm+1 = zm−1 + 2h ∗ f (X + mh, zm ),
m = 1, ...n − 1
Az els˝ o´ es az utols´ o pontt´ ol eltekintve ez a k¨ ozbens˝ oponti vagy centr´ alis differencia m´ odszer a z(x) k¨ ozel´ıt˝ o megold´ asra. y(m)′ =
zm+1 − zm−1 2h
Az els˝ o´ es az utols´ o pontban azonban a szimmetria megt¨ orik ´ es n´ emileg m´ odosul a k´ eplet, ez indokolja a m´ odszer nev´ et. A pontos megold´ as X + H pontban: y(X + H). A k¨ ozel´ıt˝ o megold´ as pedig ugyanitt: y(n) = 0.5 ∗ [zn + zn−1 + h ∗ f (X + H, zn )]
A m´ odos´ıtott k¨ ozbens˝ oponti m´ odszer m´ asodrend˝ u m´ odszer, hibatagja O(h3 ) rend˝ u, ´ ert´ eke: P∞ 2i y(n) − y(X + H) = i=1 αi ∗ h , teh´ at csak p´ aros hatv´ anyokat tartalmaz
Jav´ıt´ o formul´ aja:
y = 4/3 ∗ yn − 1/3 ∗ yn/2 ,
ahol yn h, yn/2 pedig 2h l´ ep´ esk¨ ozzel volt sz´ amolva. y(X + H) = yn + O(h5 )
A Bulirsch-Stoer m´ odszer h´ıvja a SUBROUTINE MMID-et NSTEP=2,4,6,8,12,14,24 beoszt´ asokkal, mindaddig, m´ıg a Richardsson extrapol´ aci´ oval h=0 -ra kapott ´ ert´ ek hibabecsl´ ese kisebb nem lesz a megk´ıv´ ant ǫ (relat´ıv) pontoss´ agn´ al. A Richardsson extrapol´ aci´ ot a SUBROUTINE RZEXTR v´ egzi polinom helyett pn (x)/qm (x) racion´ alis t¨ ortkifejez´ essel. A SUBROUTINE BSSTEP(Y,DYDX,NV,X,HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS)-be val´ o bel´ ep´ eskor az Y t¨ omb tartalma az X-beli megold´ as, kil´ ep´ eskor pedig az X+H pontbeli megold´ as. A megold´ as propag´ al´ as´ at x1 -bol x2 -be az ODEINT v´ egzi. Ez gondoskodik arr´ ol is, hogy a h l´ ep´ esk¨ oz ne legyen hmin -n´ el kisebb. Optimaliz´ alja a sz¨ uks´ eges sz´ amol´ ast H n¨ ovel´ es´ evel ill. cs¨ okkent´ es´ evel.
@F BESTADAPT R BESTADAPT FORTRAN STOP ED ADAPT2.RES *C x1 x2 approximate y(x) exact y(x) 0.000000000 1.800000000 -0.896758414 -0.896758416
Line´ aris t¨ obbl´ ep´ eses m´ odszerek
Az egyszer˝ us´ eg kedv´ e´ ert tekints¨ uk az x0 pontbeli kezdeti´ ert´ ek-probl´ em´ at N = 1 esetben. Enn´ el y ′ = f (x, y),
ahol y(x0 ) = y0 adott.
Vegy¨ unk h l´ ep´ eses egyenl˝ ok¨ oz˝ u beoszt´ ast:
xi = x0 + ih, ahol
i = 1, 2, .., m. L´ attuk, hogy egy k -l´ ep´ eses m´ odszer k darab f¨ uggv´ eny´ ert´ eket haszn´ al fel az yn+1
kisz´ am´ıt´ as´ ahoz, ezek a k¨ ovetkez˝ ok: yn ,yn−1 ,...,yn−k+1 . Teh´ at a k l´ ep´ eses m´ odszer ´ altal´ anos alakja: yn+1 +
k X
αi yn+1−i = h
i=1
k X
βi fn+1−i ,
i=0
ahol az αi ´ es βi mennyis´ egek megadott s´ ulyok az i = 0, 1, ..., k ´ ert´ ekekre. Itt bevezett¨ uk az ft = f (xt , yt ) jel¨ ol´ est. Ha nem haszn´ aln´ ank ezt a t¨ om¨ or jel¨ ol´ est, akkor jobban l´ atszana, hogy
a jobboldalon ´ all´ o kifejez´ es is tartalmazza a meghat´ arozand´ o yn+1 mennyis´ eget, kiv´ eve, ha β0 = 0. Ha β0 = 0, akkor explicit egyenlet¨ unk van az ismeretlen yn+1 meghat´ aroz´ as´ ara, s
ezt az egyenletet extrapol´ aci´ os formul´ anak nevezz¨ uk. Ellenkez˝ o esetben, teh´ at, ha β0 6= 0, implicit egyenlet¨ unk van ´ es ezt interpol´ aci´ os formul´ anak nevezz¨ uk.
Itt sem foglalkozunk azzal, hogy hogyan lehet adott rend˝ u k¨ ozel´ıt´ eshez meghat´ arozni az αi , βi sz´ am´ ert´ ekeit, hanem ezeket ismerteknek tekintj¨ uk.
Az Adams f´ ele m´ odszerek
azok, amik r¨ ogz´ıtett k mellett a legmagasabbrend˝ u k¨ ozel´ıt´ eseket szolg´ altatj´ ak. Ez´ ert a tov´ abbiakban csak ezekkel foglalkozunk. Kimutathat´ o, hogy a k -adrend˝ u Adams m´ odszern´ el az extrapol´ aci´ os formula hibatagja O(hk+1 ) rend˝ u, m´ıg a pontosabb interpol´ aci´ os formula hib´ aja O(hk+2 ) rend˝ u. P´ eldak´ eppen fel´ırjuk a k = 3 esetre vonatkoz´ o Adams formul´ akat. Kimutathat´ o, hogy k = 3-ra az extrapol´ aci´ os formula: yn+1 = yn +
h (23fn − 16fn−1 + 5fn−2 ) + O(h4 ). 12
M´ıg az interpol´ aci´ os formula: yn+1 = yn +
h (9fn+1 + 19fn − 5fn−1 + fn−2 ) + O(h5 ). 24
L´ athat´ o, hogy az interpol´ aci´ os formula alkalmaz´ as´ ahoz ismerni kell egy k¨ ozel´ıt˝ o ´ ert´ ek´ et a meghat´ arozand´ o mennyis´ egnek.
K´ ezenfekv˝ o, hogy ezt az extrapol´ aci´ os formula fel-
haszn´ al´ as´ aval sz´ am´ıtsuk ki. Ezt az elj´ ar´ ast prediktor-korrektor m´ odszernek h´ıvj´ ak . A prediktor-korrektor m´ odszer a k¨ ovetkez˝ o h´ arom l´ ep´ esb˝ ol ´ all: 1.) Az yn+1 meghat´ aroz´ asa exrapol´ aci´ os formul´ aval (prediktor l´ ep´ es, jele: P) 2.) Az fn+1 = f (xn+1 , yn+1 ) kisz´ am´ıt´ asa az el˝ oz˝ o l´ ep´ sben kapott yn+1 seg´ıts´ eg´ evel (jele: E) 3.) Az yn+1 meghat´ aroz´ asa interpol´ aci´ os formul´ aval, az el˝ oz˝ o l´ ep´ esben meghat´ arozott fn+1 felhaszn´ al´ as´ aval (korrektor l´ ep´ es, jele: C). Elvileg lehets´ eges az interpol´ aci´ os m´ odszer ism´ etelt felhaszn´ al´ as´ aval pontosabb k¨ ozel´ıt˝ o megold´ ast sz´ amolni yn+1 -re. Azaz a PECECEC.... iter´ aci´ ot alkalmazni a korrektor l´ ep´ es implicit egyenlet´ enek megold´ as´ ara. Ez azonban nem szok´ asos. A PEC l´ ep´ esek elv´ egz´ ese ut´ an tov´ abbl´ ep¨ unk ´ es az yn+2 -t hat´ arozzuk meg. Ennek az az oka, hogy nincs sok ´ ertelme pontosabban meghat´ arozni egy am´ ugy is csak k¨ ozel´ıt˝ o mennyis´ eget. A most k¨ ovetkez˝ o prop ′ 2 gram a k = 1, ..4-es az Adams formul´ ak haszn´ alat´ at mutatja az y = 2∗ 1 − y , y(0) = sin(0.25)
kezdeti´ ert´ ekprobl´ em´ ara. Mivel az indit´ ashoz k kezd˝ opontbeli f¨ uggv´ eny´ ert´ ek sz¨ uks´ eges, erre
felhaszn´ aljuk , hogy ismerj¨ uk az egzakt megold´ ast: y = sin(2 ∗ x + 0.25)
@COP MULTISTEP.FOR @COP MULTISTEP.DAT @F MULTISTEP R MULTISTEP FORTRAN STOP ED MULTISTEP.RES *C Feladat: Irjuk ´ at a programot az y ′′ = −y ∗ K ∗ ∗2
y(0) = 1,
y ′ (0) = 0 m´ asodrend˝ u d.e.
megold´ as´ ara. egzakt megold´ asa: y = cos(x ∗ K)
megold´ as: MULTI2.FOR adatfile: MULTI2.DAT A k¨ ozel´ıt˝ o megold´ as teljes hib´ aja. Az eddigiekben csak egyetlen h l´ epes megt´ etel´ en´ el vizsg´ altuk az elk¨ ovetett hib´ at a dif-
ferenci´ alegyenlet numerikus integr´ al´ asa sor´ an . Ez a hiba alapvet˝ oen az alkalmazott k¨ ozel´ıt˝ o m´ odszer p rendj´ et˝ ol, valamint a h l´ ep´ eshosszt´ ol f¨ ugg¨ ott. Ezt a hib´ at k´ eplethib´ anak nevezz¨ uk, ´ es ez csak egy r´ esze a megold´ as teljes hib´ aj´ anak. A kezdeti´ ert´ ek-probl´ ema megold´ as´ anak hib´ aj´ at az is befoly´ asolja, hogy milyen hib´ at k¨ ovet¨ unk el a kezdeti´ ert´ ek(ek) megad´ as´ an´ al az x0 pontban. Ahogy t´ avolodunk a kezd˝ opontt´ ol a k´ etf´ ele hiba tovaterjed ´ es kieg´ esz¨ ul min-
den l´ ep´ esn´ el az abban elv´ egzett aritmetikai m˝ uveletek kerek´ıt´ esi hib´ aival. Teh´ at, amikor eljutunk a b v´ egpontba, az itt sz´ am´ıtott megold´ as hib´ aj´ aban mindh´ arom fajta hiba szerpet j´ atszik. A k´ eplethiba ann´ al kisebb, min´ el r¨ ovidebb a h l´ ep´ esk¨ oz. Ennek ellen´ ere nem v´ alaszthatjuk h ´ ert´ ek´ et nagyon kicsinek, mert akkor a kerek´ıt´ esi hib´ ak fognak feln˝ oni a ´ sok integr´ aci´ os l´ ep´ es miatt, ami a b v´ egpont el´ er´ es´ ehez sz¨ uks´ eges. Altal´ anos javaslat, hogy becs¨ ulj¨ uk meg azt a minim´ alis pontoss´ agot, amire a megold´ as meghat´ aroz´ as´ an´ al sz¨ uks´ eg¨ unk van, ´ es haszn´ aljuk az automatikus l´ ep´ esk¨ oz szab´ alyoz´ asi lehet˝ os´ eget.
Perem´ ert´ ek-feladatok megold´ asa A gyakorlatban leggyakrabban felmer¨ ul˝ o perem´ ert´ ek-probl´ em´ ak line´ arisak, s ezeket egyszer˝ uen visszavezethetj¨ uk kezdeti´ ert´ ek-probl´ em´ ak megold´ asaira. Ennek m´ odj´ at a k¨ ovetkez˝ o egyszer˝ u p´ eld´ an mutatjuk be. Tekints¨ uk az L oper´ atorral defini´ alt Ly(x) = y ′′ (x) + p(x)y ′ (x) + q(x)y(x) = r(x),
m´ asodrend˝ u, line´ aris, inhomog´ en differenci´ alegyenletet, ami az y(a) = A ´ es az y(b) = B hat´ arfelt´ etelekkel perem´ ert´ ek-probl´ emat hat´ aroz meg.
A perem´ ert´ ek-probl´ ema y(x)
megold´ as´ at az al´ abbi k´ et kezdeti´ ert´ ek-probl´ ema megold´ as´ anak line´ arkombin´ aci´ ojak´ ent ´ all´ıtjuk el˝ o: 1 Az u(x) f¨ uggv´ eny jel¨ olje az Lu = 0 homog´ en differenci´ alegyenletmegold´ as´ at u(a) = 0 ´ es u′ (a) = 1 kezdeti´ ert´ ekek mellett. Ez teh´ at egy kezdeti´ ert´ ek-probl´ ema,
melynek
megold´ as´ ara sz´ amos m´ odszert ismert¨ unk meg az el˝ oz˝ oekben. 2 A z(x) f¨ uggv´ eny jel¨ olje az Lz = r inhomog´ en differenci´ alegyenlet megold´ as´ at z(a) = A ´ es z ′ (a) = 0 kezdeti´ ert´ ekek mellett. Ez egy m´ asik kezdeti´ ert´ ek-probl´ ema. K¨ onny˝ u bel´ atni, hogy az y(x) = c1 u(x) + c2 z(x)
line´ arkombin´ aci´ o c1 =
B−z(b) u(b)
´ e s c2 = 1 ´ ert´ ekek mellett a perem´ ert´ ek-probl´ ema megold´ as´ at
szolg´ altatja.
Nemline´ aris differenci´ alegyenlettel kapcsolatos perem´ ert´ ek-probl´ ema eset´ eben pr´ ob´ alkozhatunk az u.n. bel¨ ov´ eses m´ odszerrel (shooting method). Noha ennek a m´ odszernek a konvergenci´ aj´ ara ´ altal´ aban nincs garancia, ha azonban konverg´ al, akkor ´ altal´ aban gyorsabb a konkurens m´ odszerekn´ el. A bel¨ ov´ eses m´ oszer nev´ et arr´ ol a technik´ ar´ ol kapta, amit r´ egen a t˝ uz´ erek haszn´ altak, amikor ´ agy´ uval akartak egy t´ avoli c´ elpontot eltal´ alni.
Az ´ agy´ u
v´ızszinteshez m´ ert sz¨ og´ enak v´ altoztat´ as´ aval l˝ ott´ ek be a c´ elt, azaz kerestek olyan α1 ´ e s α2 sz¨ ogeket, hogy az egyik sz¨ ogn´ el a l¨ oved´ ek a c´ elpont el˝ ott, a m´ asikn´ al a c´ elpont m¨ og¨ ott csap´ odott be. M´ asodrend˝ u differenci´ alegyenletre vonatkoz´ o perem´ ert´ ek-probl´ ema eset´ en ennek az felel meg, hogy az a kezd˝ opontban az y ′ (a) mennyis´ egre tal´ alunk olyan α1 ´ e s α2 ´ ert´ ekeket, hogy az ezekkel ind´ıtott kezdeti´ ert´ ek-probl´ em´ ak megold´ asai x = b-ben k¨ ozrefogj´ ak a b pontbeli y(b) = B peremfelt´ etelt. Ekkor a perem´ ert´ ek-probl´ ema megold´ as´ ahoz tartoz´ o igazi y ′ (a) ´ ert´ ekre fenn´ all, hogy y ′ (a) ∈ [α1 , α2 ]
Az igazi y ′ (a) ´ ert´ ek meghat´ aroz´ asa ezut´ an a nemline´ aris egyenletek megold´ as´ ara vonatkoz´ o m´ odszerek (intervallum felez´ es, h´ urm´ odszer, stb) valamelyik´ evel t¨ ort´ enhet. V´ eges differenci´ ak m´ odszere
Oldjuk meg az el˝ oz˝ oekben vizsg´ alt perem´ ert´ ek feladatot a v´ eges differenci´ ak m´ odszer´ evel! Ekkor teh´ at az y ′′ (x) + p(x)y ′ (x) + q(x)y(x) = r(x),
differenci´ alegyenletet, kell megoldanunk az y(a) = A ´ es az y(b) = B hat´ arfelt´ etelek mellett. Ehhez vegy¨ uk fel az [a, b] intervallumnak a h =
b−a n+1
(egyenl˝ o) l´ ep´ esk¨ oz˝ u beoszt´ as´ at, melyre
xi = a + ih valamennyi i = 0, 1, ..., n + 1 ´ ert´ ekre. Alkalmazzuk az yi = y(xi ) jel¨ ol´ est, amivel y0 = A ´ es yn+1 = B fenn kell ´ alljon. Az xk oszt´ aspontokban a megold´ as els˝ o, illetve m´ asodik
deriv´ altjait a k¨ ovetkez˝ o v´ eges differencia h´ anyadosokkal k¨ ozel´ıthetj¨ uk: y ′ (xk ) =
valamint y ′′ (xk ) =
yk+1 − yk−1 + O(h), 2h
yk+1 − 2yk + yk−1 + O(h2 ). h2
Jel¨ olj¨ uk az oszt´ aspontokban a megfelel˝ o f¨ uggv´ eny´ ert´ ekeket:
pk = p(xk ), qk = q(xk ) ´ es
rk = r(xk )-val! Ezeket a differenci´ alegyenletbe helyettes´ıtve a k¨ ovetkez˝ o line´ aris egyenle-
trendszert kaphatjuk az ismeretlen y1 , y2 , ..., yn mennyis´ egekre: yk+1 (2 + hpk ) + yk (2h2 qk − 4) + yk−1 (2 − hpk ) = 2h2 rk
valamennyi k = 1, ..., n ´ ert´ ekre. Ennek a line´ aris egyenletrendszernek a megold´ asa szolg´ altatja a perem´ ert´ ek-probl´ ema megold´ as´ at. Megjegyezz¨ uk, hogy ennek a line´ aris egyenletrendszernek a matrixa u.n. h´ arom´ atl´ os matrix, melyn´ el csak a f˝ oa ´tl´ oban ´ es az alatta, illetve felette elhelyezked˝ o mell´ ek´ atl´ okban lev˝ o elemek k¨ ul¨ onb¨ oznek null´ at´ ol. A m´ odszert az al´ abbi p´ eldaprogrammal mutatjuk be.
Feladat: oldjuk meg az y ′′ − y ∗ (1 + x2 ) = 1 diff. e. y(0) = 1 ´ es y(1) = 0 ´ ert´ ekekkel adott
prem´ ert´ ekfeladat´ at a v´ eges differenci´ ak m´ odszer´ evel: @COP PEREMERT.FOR @COP PEREMERT.DAT @f PEREMERT R PEREMERT FORTRAN STOP A megold´ as:: ED PEREMERT.RES *C 1.0000
1.0746
1.1701
1.2877
1.4294
1.5976
1.7958
2.0285
2.3013
2.6219
3.0000
Saj´ at´ ert´ ek-probl´ em´ ak megold´ asa v´ eges differenci´ ak m´ odszer´ evel.
Tekints¨ uk az y ′′ (x) + p(x)y ′ (x) + (q(x) − λ)y(x) = 0,
m´ asodrend˝ u differenci´ alegyenletnek az y(a) = 0 ´ es az y(b) = 0 hat´ arfelt´ etelekkel defini´ alt saj´ at´ ert´ ek-feladat´ at! Nyilv´ anval´ o, hogy a differenci´ alegyenlet megold´ asa f¨ ugg az ismeretlen λ param´ eter ´ ert´ ek´ et˝ ol, teh´ at y = y(x, λ).
perem´ ert´ ek-probl´ ema kiel´ eg¨ ulj¨ on.
A λ ´ ert´ ek´ et u ´ gy kell meghat´ arozni, hogy a
Jel¨ olj¨ uk λ param´ eternek ezeket az ´ ert´ ekeit, teh´ at a
saj´ at´ ert´ ekeket λn -nel. A λ = λn saj´ at´ ert´ ekek mellett kapott y(x, λn ) megold´ asok a probl´ ema saj´ atf¨ uggv´ enyei.
Vegy¨ uk az el˝ oz˝ oekben haszn´ alt egyenl˝ ok¨ oz˝ u beoszt´ ast ´ es alkalmazzuk
a deriv´ altakra a v´ eges differenciah´ anyadosokkal val´ o k¨ ozel´ıt´ est!
Irjuk be ezt a differ-
enci´ alegyenletbe , ekkor a k¨ ovetkez˝ o line´ aris egyenletrendszert kapjuk: yk+1 (
pk 1 pk 2 1 + ) + yk [(qk − 2 ) − λ] + yk−1 ( 2 − )=0 h2 2h h h 2h
k = 1, ..., n-re. Ez ism´ et egy h´ arom´ atl´ os matrix´ u, homog´ en, line´ aris egyenletrendszer az y1 , y2 , ..., yn ismeretlenekre. A λ-t´ ol f¨ ugg˝ o tagot az egyenlet jobboldal´ ara rendezve Ay = λy
alak´ u egyenletrendszert kapunk, ami nem m´ as, mint az A matrix saj´ at´ ert´ ek-probl´ em´ aja. Ezzel a differenci´ alegyenlet saj´ at´ ert´ ek-probl´ em´ aj´ at matrix saj´ at´ ert´ ek-probl´ em´ aj´ anak megold´ as´ ara vezett¨ uk vissza, amire standard numerikus m´ odszerek ´ allnak rendelkez´ esre. Meg kell jegyezni, hogy a matrix saj´ at´ ert´ ek-probl´ em´ aj´ anak a megold´ asa a differenci´ alegyenlet saj´ at´ ert´ ekprobl´ em´ aja megold´ as´ anak csup´ an k¨ ozel´ıt´ ese, aminek pontoss´ aga f¨ ugg p´ eld´ aul a h beoszt´ as nagys´ ag´ at´ ol. (Hiszen a deriv´ altakat a h l´ ep´ esk¨ ozh¨ oz tartoz´ o v´ eges differenci´ ak h´ anyadosaival k¨ ozel´ıtett¨ uk.) Bizonyos hat´ arig pontos´ıthatjuk a differenci´ alegyenlet saj´ at´ ert´ ek-probl´ ema megold´ as´ at a h cs¨ okkent´ es´ evel, ennek ´ ara azonban az, hogy egyre nagyobb m´ eret˝ u A matrix saj´ at´ ert´ ek-probl´ em´ aj´ at kell megoldanunk. Bizonyos hat´ aron t´ ul azonban a m´ odszer pontoss´ aga nem jav´ıthat´ o, mert a kerek´ıt´ esi hib´ ak halmoz´ od´ asa a k¨ ozel´ıt´ es roml´ as´ at okozza.
Saj´ at´ ert´ ek-probl´ em´ ak megold´ asa iter´ aci´ oval.
Egy m´ asik m´ odszer az el˝ oz˝ o saj´ at´ ert´ ek-probl´ emamegold´ as´ ara az, hogy v´ alasztunk egy a < xm < b illeszt´ esi pontot, ´ es megoldjuk a k¨ ovetkez˝ o k´ et kezdeti´ ert´ ek-probl´ em´ at:
1.) Kisz´ am´ıtunk egy baloldali megold´ ast az [a, xm ] intervallumon, az yl (a, λ) = 0 ´ es yl′ (a, λ) = 1 kezdeti felt´ etelekkel, majd
2.) kisz´ am´ıtunk egy jobboldali megold´ ast az [xm , b] intervallumon, az yr (b, λ) = 0 ´ es yr′ (b, λ) = 1 kezdeti felt´ etelekkel.
Ezut´ an k´ epezz¨ uk az u.n. logaritmikus deriv´ altak k¨ ul¨ onbs´ eg´ et: f (λ) =
yl′ (xm , λ) yr′ (xm , λ) − . yl (xm , λ) yr (xm , λ)
Az f (λ) nemline´ aris f¨ uggv´ eny, ´ es a λn saj´ at´ ert´ ekek f (λn ) = 0 egyenlet gy¨ okei.
Ugyanis
ezek mellett mindk´ et megold´ as az y(x, λn ) saj´ atf¨ uggv´ ennyel ar´ anyos, ´ıgy logaritmikus deriv´ altjaik egyenl˝ oek. Tegy¨ uk fel, hogy tudunk olyan λ(0) ´ ert´ eket adni, mely k¨ ozel´ ebe esik valamelyik λn saj´ at´ ert´ eknek. Ekkor a λ(0) kezd˝ o´ ert´ ekkel elind´ıtott Newton iter´ aci´ o´ altal´ aban a λn saj´ at´ ert´ ekhez fog konverg´ alni. Term´ eszetesen az f (λ) = 0 egyenlet gy¨ okeit m´ as, a nemline´ aris egyenlet megold´ as´ ara vonatkoz´ o numerikus m´ odszerrel is megkereshetj¨ uk. A saj´ at´ ert´ ek meghat´ aroz´ asa ut´ an a saj´ atf¨ uggv´ enyt u ´ gy kaphatjuk meg, hogy a jobboldali megold´ ast az xm illeszt´ esi pontban csatlakoztatjuk a baloldali megold´ ashoz, majd ha sz¨ uks´ eges, akkor az ´ıgy kapott folytonos f¨ uggv´ enyt 1-re norm´ aljuk.
M´ atrixok saj´ at´ ert´ ekprobl´ em´ aja M´ atrixokr´ ol ´ altal´ aban: transzpon´ al´ as, n´ egyzetes m´ atrixok: egys´ eg, diagon´ alis, tridiagon´ alis, kontinu´ ans m´ atrixok fels˝ o h´ aromsz¨ og ´ es Hessenberg m´ atrixok als´ o””” Szimmetrikus, ¨ onadjung´ alt, ortogon´ alis, unit´ er, norm´ alis m´ atrixok Val´ os elem˝ u m´ atrixokra Hermitikus=szimmetrikus ortogon´ alis=unit´ er=norm´ alis M´ atrixok saj´ at´ ert´ ek feladat´ anak numerikus megold´ asa jobbodali es baloldali saj´ atvektorok Szimmetrikus A matrix eset´ en baloldali saj´ atvektor a jobboldali transzpon´ altja Csak jobboldali (oszlop) saj´ atvektorokkal foglalkozunk. Ax = y trafo a saj´ atvektorok ir´ any´ at v´ altozatlanul hagyja. Ax = λx(A − λE)x = 0Det(A − λE) = 0φ(λ) = 0
Karakterisztikus polinom (n-edfok´ u). maximum n db. k¨ ul¨ onb¨ oz˝ o gy¨ oke lehet a karakterisztikus egyenletnek. Ha A diagonalis vagy U , vagy L, akkor determin´ ansa l´ enyeg´ eben a f˝ oa ´tl´ obeli elemeinek szorzata, teh´ at a saj´ at´ ert´ ekek a f˝ oa ´tl´ obeli elemek.
Saj´ atf¨ uggv´ enyek, ha a saj´ at´ ert´ ekek
k¨ ul¨ onb¨ oz˝ oek (nem degener´ alt eset), akkor az ei egys´ egvektorok. Ha T regul´ aris (Det(T ) 6= 0), akkor B = T × A × T −1 hasonl´ os´ agi trafo v´ altozatlanul hagyja
A saj´ at´ ert´ ekeit. Ha xi A saj´ atvektora, akkor B -nek ugyanahhoz a saj´ at´ ert´ ekhez tartoz´ o
saj´ atvektora: T × xi Teh´ at, ha siker¨ ul az A m´ atrixot hasonl´ os´ agi traf´ oval diagonaliz´ alni,
akkor megkapjuk a saj´ at´ ert´ ekeket a f˝ oa ´tl´ oban, ´ es a saj´ atvektor is rekonstru´ alhat´ o a hasonl´ os´ agi traf´ ob´ ol. A hasonl´ os´ agi traf´ o oszlopvektorai a saj´ atvektorok. Val´ os, szimmetrikus m´ atrixra a saj´ atvektorok val´ osak ´ es ortogon´ alisak, ´ıgy a hasonl´ os´ agi traf´ o is ortogon´ alis traf´ o. A = Z T .A.Z
Val´ os, nem szimmetrikus m´ atrixok is diagonaliz´ alhat´ ok a saj´ atvektorokb´ ol ¨ ossze´ all´ o m´ atrixal, de a saj´ atvektorok ´ altal´ aban komplexek.
Val´ os hasonl´ os´ agi traf´ oval kv´ azi-diagonaliss´ a
tehet˝ ok (2 × 2-es al-m´ atrixokt´ ol eltekintve).
A numerikus megold´ asn´ al az A m´ atrixot
P1 , P2 , ... hasonl´ os´ agi traf´ ok sorozat´ aval a diagon´ alis alak fel´ e l¨ okd¨ oss¨ uk.
Ha ez sikeres,
akkor XR = P1 .P2 .P3 ... a saj´ atvektorokb´ ol alkotott m´ atrix lesz. Numerikus m´ odszerek fajt´ ai: I./ Kik¨ usz¨ ob¨ ol´ esi m´ odszerek: (Jacobi, Gauss, Householder) egyl´ ep´ esben kinull´ aznak egy, vagy t´ obb elemet, ´ es ezeket a l´ ep´ eseket iter´ aljuk. II./ Faktoriz´ aci´ os m´ odszerek: k¨ oz¨ os, hogy A = FL .FR szorzatra bontjuk, s FR .F= FL−1 .A.FL hasonl´ os´ agi traf´ ot kapunk. Ide tartozik a QR ´ es a QL m´ odszer. I m´ odszerrel j´ o kezdeti ´ allapotot ´ all´ıtunk el˝ o, majd folytatjuk a II m´ odszerrel, ami ekkor gyorsabban konverg´ al. Householder m´ odszer (szimmetrikus m´ atrixokra) n − 2 db ortogon´ alis traf´ oval tridiagon´ alis
alakra hoz.
Ezt valos´ıtja meg a SUBROUTINE TRED2(A,N,NP,D,E) rutin A diagonaliz´ alni kiv´ ant m´ atrixot az A N × N -es t¨ ombbe helyezz¨ uk,
N <= N P
Az eredm´ eny: f˝ oa ´tl´ obeli elemek a D(NP) val´ os vektorba, D(1),...,D(N) lesznek elt´ arolva, a f˝ oa ´tl´ o alatti elemek E(2),...,E(N)-ben, E(1)=0 mindig. Az A t¨ omb fel¨ ul´ır´ odik. A kapott tridiagon´ alis m´ atrix saj´ at´ ert´ ekeinek ( ´ es esetleg saj´ at vektorainak) a meghat´ aroz´ asa: fel´ırjuk a φ(λ) karakterisztikus polinomot, majd a φ(x) = 0 egyenlet gy¨ okeit intervallumfelez´ essel, vagy NR m´ odszerrel meghat´ arozzuk.
A hozz´ ajuk tartoz´ o saj´ atvektorokat in-
verziter´ aci´ oval hat´ arozhatjuk meg. Ha nincs sz¨ uks´ eg valamennyi saj´ at´ ert´ ekre ´ es vektorra, akkor hat´ ekonyabb egy faktoriz´ aci´ os elj´ ar´ as alkalmaz´ asa. A QR ´ es a QL algoritmus B´ armely val´ os A m´ atrix fel´ırhat´ o A = Q.R alakba, ahol Q ortogon´ alis m´ atrix, R pedig fels˝ o h´ aromsz¨ og m´ atrix. Irjuk most ˝ oket
ford´ıtott sorrendben, azaz sz´ am´ıtsuk ki az A′ = R.Q m´ atrixot. Mivel Q ortogon´ alis, ez´ ert Q−1 = QT , s ´ıgy QT .A = Q−1 .Q.R = R, teh´ at A′ = QT .A.Q A-nak ortogon´ alis transzform´ altja.
Bizony´ıthat´ o, hogy a QR traf´ o meg˝ orzi az A m´ atrix szimmetri´ aj´ at, teh´ at tridiagonalit´ as´ at ´ es Hessenberg alakj´ at (ha volt neki). Hasonl´ ok´ eppen kimutathat´ o, hogy b´ armely A m´ atrix fel´ırhat´ o A = Q.L alakba is, ahol L als´ o h´ aromsz¨ ogm´ atrix ´ es Q ortogon´ alis m´ atrix.
A
tov´ abbiakban a QL algoritmussal foglalkozunk, ami ortogon´ alis traf´ ok al´ abbi sorozata: az s-edik l´ ep´ esben: As = Qs .Ls a traf´ o: As+1 = Ls .Qs = QTs .As .Qs
Ha s v´ egtelenhez tart, akkor As als´ o h´ aromsz¨ ogm´ atrixhoz tart, ha A saj´ at´ ert´ ekei nem degener´ altak. Ellenkez˝ o esetben, ha λi p > 1 multiplicit´ as´ u, akkor As egy p-edrend˝ u m´ atrixtol (melynek saj´ at´ ert´ ekei a λi ) eltekintve tart als´ o h´ aromsz¨ ogm´ atrixhoz. M˝ uvelet´ıg´ eny: O(n3 ) iter´ aci´ onk´ ent ´ altal´ anos m´ atrixra, Hesenberg alakra O(n2 ), tridiagon´ alis m´ atrixra O(n) iter´ aci´ onk´ ent. Al´ abb a kerek´ıt´ esi hiba szempontj´ ab´ ol el˝ ony¨ osebb u.n. implicit eltol´ asos QL algoritmust alkalmaz´ o szubrutint adjuk: SUBROUTINE TQLI(D,E,N,NP,Z) Bemen˝ o param´ eterek: D(NP),E(NP) , ami a tridiagon´ alis m´ atrix diagon´ alis ´ es f˝ oa ´tl´ o alatti elemeit tartalmaz´ o vektorok. N ×N -es eredeti A m´ atrix eset´ en D(1),..,D(N) ´ es E(2),...,E(N)
van felt¨ oltve. Ez az eredm´ enye a TRED2 szubrutinnak. Z egy N P × N P m´ eret˝ u val´ os t¨ omb,
amit egys´ egm´ atrix-szal kell felt¨ olteni, ha a saj´ atvektorokat is meg akarjuk hat´ arozni. Mivel TRED2 az A m´ atrixot ´ıgy t¨ olti fel, c´ elszer˝ u ezt megadni TQLI-ben Z-nek. @COP SYMSAJ.FOR ld. Szidarovszki 7.9 288 oldal @COP SYMSAJ.DAT mivel A szimmetrikus, ez´ ert el´ eg az als´ o h´ aromsz¨ oget megadni, teh´ at @F SYMSAJ 1 1 0.5 R SYMSAJ A= 1 1 0.25 FORTRAN STOP 0.5 0.25 2 eset´ en 1,1,1,0.5,0.25,2 -t ED SYMSAJ.DAT 3,1
1.d0,1.d0,1.d0,0.5d0,0.25d0,2.d0, ED SYMSAJ.RES eigenvalue= -0.0166472836 eigenvector: 0.72
-0.69
-0.09
eigenvalue= 1.4801214232 eigenvector: -0.44
-0.56
0.70
eigenvalue= 2.5365258604 eigenvector: 0.53
0.46
0.71
Nemszimmetrikus val´ os m´ atrixok saj´ at´ ert´ ekeinek meghat´ aroz´ asa 3 l´ ep´ esben:
1./
el˝ osz¨ or c´ elszer˝ u kiegyens´ ulyozni az A m´ atrixot, hogy sorainak ´ es oszlopainak norm´ aja k¨ ozel´ıt˝ oleg azonos legyen, 2./ ezut´ an Hessenberg alakra hozzuk (fels˝ o Hessenberg alak: als´ o mell´ ek´ atl´ o alatt valamennyi elem zerus) 3./ veg¨ ul alkalmazzuk a QR algoritmust a val´ os Hessenberg m´ atrixra. b˝ ovebben: 1./ Hasonl´ os´ agi transzform´ aci´ oval v´ egezz¨ uk a SUBROUTINE BALANC(A,N,NP) seg´ıts´ eg´ evel. Param´ eterek: A(NP,NP) val´ os t¨ omb, bemenetkor az egyens´ ulyozand´ o, kimenetkor a kiegyens´ ulyozott A m´ atrixot tartalmazza. A maxim´ alis m´ erete.
N az A m´ atrix aktu´ alis m´ erete, N P az
2./ A Hessenberg alakra hoz´ ast vagy f˝ oelemkiv´ alaszt´ asos Gauss
m´ odszerrel vagy Householder redukci´ okkal ´ erhetj¨ uk el. Mivel a Gauss elimin´ aci´ o nem hasonl´ os´ agi traf´ oja az A m´ atrixnak, ez´ ert azt kicsit m´ odos´ıtani kell. Az r-edik l´ ep´ es el˝ ott az eredeti A = A1 m´ atrix Ar m´ atrix´ a v´ altozott, ´ es ez egy fels˝ o Hessenberg m´ atrix az els˝ o r−1 sor´ aban ´ es oszlop´ aban. Az r-edik l´ ep´ es a k¨ ovetkez˝ o m˝ uveleteket tartalmazza: a./ megker-
ess¨ uk az r-edik oszlopnak a f˝ oa ´tl´ o alatti maxim´ alis abszol´ ut´ ert´ ek˝ u elem´ et. Ha ez zerus, akkor atugorjuk a tov´ ´ abbi k´ et l´ ep´ est. Egy´ ebk´ ent jel¨ olj¨ uk ezen elem sor´ at r′ -vel. b./ Felcser´ elj¨ uk az r′ -edik sort az r + 1-edikkel. (Ez a f˝ oelemkiv´ alaszt´ as). Az´ ert, hogy hasonl´ os´ agi traf´ o legyen az r′ oszlopot is felcser´ elj¨ uk az r + 1-edikkel. c./ i = r + 2, r + 3, ..., N -re kisz´ am´ıtjuk az ni,r+1 = ai,r /ar+1,r szorz´ ofaktort ´ es ezzel megszorozva az r + 1-edik sort, azt kivonjuk az iedik sorb´ ol. Az´ ert, hogy hasonl´ os´ agi traf´ o legyen az i-edik oszlop ni,r+1 szeres´ et hozz´ aadjuk az r + 1-edik oszlophoz. Az al´ abbi rutin ezeket a l´ ep´ eseket tartalmazza ´ es mivel csak a saj´ at´ ert´ ekeket hat´ arozza meg, nem k¨ oveti nyomon a hasonl´ os´ agi traf´ ok sorrendj´ et. SUBROUTINE ELMHES(A,N,NP) param´ eterek: mint kor´ abban, az A t¨ omb bemenetkor a m´ ar feltehet˝ oleg kiegyens´ ulyozott nem szimmetrikus val´ os m´ atrixot, kimenetkor pedig annak fels˝ o Hessenberg alakj´ at tartalmazza, teh´ at az A(i, j) m´ atrix i > j + 1 elemei gyakorlatilag 0-nak tekinthet˝ ok. A legkisebb lebeg˝ opontos sz´ am nagysagrendj´ ebe es˝ o v´ eletlen sz´ amok.
3./ A QR algoritmust az eloz˝ oekben kapott fels˝ o Hessenberg alak´ u val´ os A m´ atrixra a SUBROUTINE HQR(A,N,NP,WR,WI) rutine-nal v´ egezhetj¨ uk. Itt az els˝ o h´ arom param´ eter ugyanaz, mint kor´ abban, a WR(1),...,WR(N) val´ os vektor elemei kimeneten a komplex sajat´ ert´ ekek val´ os r´ esz´ et, a WI(1),..,WI(N) elemek pedig a saj´ at´ ert´ ekek k´ epzetes r´ eszeit tartalmazz´ ak. Az ´ıgy kapott saj´ at´ ert´ ekek pontoss´ aga tov´ abb n¨ ovelhet˝ o, ´ es a saj´ atvektorok is meghat´ arozhat´ ok inverz iter´ aci´ oval. P´ elda: Hat´ arozzuk meg a k¨ ovetkez˝ o 1.5
0.1
4.5
− 1.5
−22.5
3.5
12.5
−2.5
0.3
4.5
−2.5
0.1
4.5
− 2.5 − 2.5 2.5
nemszimmetrikus val´ os m´ atrix saj´ at´ ert´ ekeit. Ezt a feladatot a fent v´ azolt m´ odszerrel oldjuk meg. @COP NONSY.FOR @COP NONSY.DAT @F NONSY R NONSY FORTRAN STOP ED NONSY.RES *C sajatertek: 3.0000 +i* -4.0000 sajatertek: 3.0000 +i* 4.0000 sajatertek: 2.0000 +i* 0.0000 sajatertek: 4.0000 +i* 0.0000 L´ athat´ o, hogy az utols´ o k´ et saj´ at´ ert´ ek val´ os, m´ıg az els˝ o kett˝ o komplex ´ es konjug´ altjai egym´ asnak. A saj´ at´ ert´ ekek finom´ıt´ asa ´ es a saj´ atvektorok kisz´ am´ıt´ asa inverz iter´ aci´ oval t¨ otrt´ enhet. A megoldand´ o saj´ at´ ert´ ekprobl´ ema: (A − λ.E).x = 0
Az inverz iter´ aci´ o k -adik l´ ep´ es´ eben a k¨ ovetkez˝ o egyenletet kell megoldanunk: (A − τk .E).y = bk Itt a τk a k -adik l´ ep´ esbeli k¨ ozel´ıt´ ese az egyik saj´ at´ ert´ eknek, bk pedig egyre
cs¨ okken˝ o elemekkel rendelkez˝ o marad´ ekvektor (ez tart zerushoz, ha a pontos mego Miut´ an
az inhomog´ en line´ aris egyenletrendszert az y vektor komponenseire megoldottuk, a k + 1edik k¨ ozel´ıt´ es´ et az x saj´ atvektornak 1-re norm´ al´ assal kapjuk: y/|y|. A saj´ at´ ert´ ek k¨ ozel´ıt´ es´ et pedig τk+1 = τk +1/(bk .y) szolg´ altatja. L´ enyeges, hogy a kiindul´ o saj´ at´ ert´ ek¨ unk el´ eg j´ o legyen, m´ ask¨ ul¨ onben a m´ odszer t´ ul lassan, vagy egy´ altal´ an nem konverg´ al.
Irodalom
Szidarovszky Ferenc, Bevezet´ es a numerikus m´ odszerekbe, K¨ ozgazdas´ agi ´ es Jogi K¨ onyvkiad´ o, Budapest, 1974 William H. Press, Saul A. Teukolsky, William T. Vettering, Brian P Flannery, Numerical Recipes in FORTRAN, University Press, Cambridge, 1992. elemei szint´ en ismertek.
ld´ ast k¨ ozel´ıtj¨ uk), melynek
