Matematika | Analízis » Soós Ivett - Közönséges differenciálegyenletek numerikus megoldása

http://www.doksihu Eötvös Loránd Tudományegyetem Természettudományi Kar Közönséges differenciálegyenletek numerikus megoldása Szakdolgozat Soós Ivett Matematika B.Sc, Matematikai elemz® szakirány Témavezet®: Mincsovics Miklós Alkalmazott Analízis és Számításmatematikai Tanszék Budapest 2010 http://www.doksihu Tartalomjegyzék 1. Bevezet® 2 2. Dierenciálegyenletek 3 2.1 A dierenciálegyenletek típusai 2.2 A közönséges dierenciálegyenletek általános alakja . 2.21 A KDE típusai 2.3 A dierenciálegyenletek stabilitása 2.31 Stabilitási alapfogalmak 2.32 Példa 2.4 Dierenciálegyenletek a gyakorlatban 2.41 Példa 1 2.42 Példa 2 2.5 A dierenciálegyenletek megoldhatósága 2.6 A megoldandó probléma: . 3 . . . . . . . . . . 3 4 5 5 6 6 6 7 8 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. Numerikus módszerek 10 3.1 A numerikus módszerek "jósága" 3.11 A numerikus módszerek hibája 3.12 Numerikus módszerek konvergenciája 3.13 A konvergencia rendje 3.14 Numerikus módszerek 0-stabilitása 3.15 Numerikus módszerek konzisztenciája 4. A legegyszer¶bb megoldások . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 11 11 12 12 12 13 4.1 Az explicit Euler módszer 13 1 http://www.doksihu 4.11 Példa 4.12 A hibaegyenlet 4.13 Az explicit Euler módszer konzisztenciája 4.14 Az Euler módszer

lokális hibája 4.15 Az Euler módszer 0-stabilitása 4.16 A "jó" lépésköz megválasztása 4.17 A teszt egyenlet 4.18 Abszolút stabilitás 4.19 A-stabilitás 4.110 Az explicit Euler módszer abszolút stabilitási tartománya 4.111 A módszer hiányosságai 4.2 Az implicit Euler módszer 4.21 Az implicit Euler módszer konzisztenciája 4.22 Az implicit Euler módszer abszolút stabilitási tartománya 4.23 Merev dierenciálegyenletek 4.24 Példa 4.25 Az implicit Euler módszer el®nyei és hátrányai 4.3 Szimmetrikus/ trapéz módszer 4.31 A szimmetrikus módszer konzisztenciája 4.32 A trapéz módszer abszolút stabilitási tartománya . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . 5. Runge-Kutta típusú módszerek 15 16 17 17 18 18 19 19 19 20 21 21 22 22 23 24 25 26 26 26 28 5.1 A kvadratúra formulák 5.11 A javított Euler módszer lokális hibája 5.12 Az RK módszerek 0-stabilitása 5.2 Az implicit Runge-Kutta módszerek 5.21 Az implicit formula el®nyei és hátrányai 5.3 Az explicit Runge-Kutta formulák 5.31 Az els® rend¶ Runge-Kutta formulák 5.32 A megoldhatóság feltételei 5.33 A másodrend¶ Runge-Kutta formulák 5.34 Példák másodrend¶ formulákra 5.35 Harmadrend¶ RK formulák 5.36 A klasszikus negyed rend¶ Runge-Kutta formula 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 29 31 32 32 32 32 33 33 34 34 34 http://www.doksihu 5.37 Az explicit Runge-Kutta módszer abszolút

stabilitási tartománya 35 5.38 RK módszerek összehasonlítása 36 6. Összefoglalás 38 3 http://www.doksihu 1. fejezet Bevezet® "Bár napjaink matematika könyveiben szinte hemzsegnek az absztrakt szimbólumok, ez azonban éppúgy nem jelenti a matematika lényegét, mint ahogy a zene valódi mibenléte sem a hangjegyek jelölésrendszerében keresend®." Keith Devlin A témaválasztásban számomra fontos szerepe volt annak, hogy a matematika olyan területével foglalkozzak, mely közvetlen kapcsolatban áll a gyakorlati problémák megoldásával. A dierenciálegyenletek a tudomány szinte minden területén jelen vannak, ezért a megoldhatóságuk nagyon fontos szerepet tölt be mindennapjainkban. Ám a legtöbb esetben ezeket nem olyan egyszer¶ kiszámítani Ahogy az idézetben is szerepel, a matematika lényege az ehhez hasonló problémák megoldása, ennek ellenére a következ® néhány oldal sem fog sz¶kölködni "absztrakt

szimbólumokban". 4 http://www.doksihu 2. fejezet Dierenciálegyenletek 2.1 A dierenciálegyenletek típusai Egy dierenciálegyenlet egy függvény és annak dierenciáltja közötti kapcsolatot mutatja meg. Ezek az összefüggések sok esetben a tudomány egyéb területein felmerül® problémák matematikai modelljei. Ezek a modellek a gyakorlatban igen hasznosak, ám elég összetettek is ahhoz, hogy a megoldásuk egzakt legyen. Legyen szó zikai, biológiai, vagy közgazdaságtani problémáról, biztosan függenek az id®t®l, vagy esteleg egy másik változó paramétert®l, tehát a dierenciálegyenletek olyan folyamatokat írnak le, melyek nem diszkrét lépésekben zajlanak. A dierenciálegyenleteknek két f® típusa van, nevezetesen 1. Közönséges dierenciálegyenletek 2. Parciális dierenciálegyenletek A f® különbség e két típus között az, hogy a közönséges dierenciálegyenletekben (KDE) az ismeretlen egyváltozós, a parciális

dierenciálegyenletknél pedig többváltozós ismeretlent keresünk. A továbbiakban a közönséges dierenciálegyenletekkel fogunk foglalkozni. 2.2 A közönséges dierenciálegyenletek általános alakja Nézzük meg el®ször az els®rend¶ explicit KDE általános alakját. 5 http://www.doksihu 2.21 Megjegyzés A dierenciálegyenlet rendje a legmagasabb derivált rendje x0 (t) = f (t, x(t)) ahol f : R Rn adott függvény, és az ismeretlenünk pedig x : R Rn . 2.22 Deníció Legyen F : Rn+2 R függvény. Az n-ed rend¶ közönséges dier- enciálegyenlet általános alakja: F (t, x(t), x0 (t) . x(n) (t)) = 0 Legyen f : Rn+1 R adott. Az n-ed rend¶ explicit közönsége dierenciálegyenlet általános alakja x(n) (t) = f (t, x(t), x0 (t) . x(n−1) (t)) 2.21 A KDE típusai A közönséges dierenciálegyenleteket két nagy csoportra oszthatjuk: lineáris és nemlineáris. Ez a megoldhatóság szempontjából igen fontos, hiszen a nemlineáris egyenletek

pontos megoldására nincs bevett módszer, és a közelít® megoldások kiszámítása is komplikáltabb. A lineáris dierenciálegyenletek újabb két nagy csoportra oszthatók. Legyen a lineáris els®rend¶ KDE általalános alakja: y 0 (t) = A(t)y(t) + b(t) Két esetet küönböztetünk meg: 1. ha A(t) nem függ t-t®l, azaz konstans, így az egyenlet állandó együtthatós 2. ha A(t) függ t-t®l, azaz az egyenlet változó együtthatós Az els® esetben létezik egzakt megoldási módszer, de legtöbbször csak nagy nehézségek árán tudjuk meghatározni a pontos megoldást. A második esetre nincs olyan bevett módszer, mellyel kiszámíthatnánk a pontos értékeket. 2.23 Megjegyzés A dierenciálegyenlet típusa, és megoldhatóságának nehézsége, nyílván a közelít® megoldások pontosságát is befolyásolja. 6 http://www.doksihu 2.3 A dierenciálegyenletek stabilitása 2.31 Stabilitási alapfogalmak A stabilitáselméleti alapfogalmak szemléltetésére

tekintsük meg el®ször az alábbi egyszer¶ zikai példát. Képzeljünk el egy golyót az alábi 3 egyensúlyi helyzetben: 1. egy gödör alján 2. egy domb tetején 3. egy vízszintes sík felületen Mindhárom helyzetben egyensúlyban van a golyó (ha nem mozdítjuk meg, helyben marad), azonban ha kicsit elmozdítjuk, majd elengedjük, akkor mindhárom esetben más történik. Az els® esetben a golyó visszagurul a gödör aljára Ezt az egyensúlyi helyzetet nevezzük asszimptotikusan stabilisnak. A második esetben a golyó legurul a domboldalon, egyre jobban eltávolodik az eredeti helyzetét®l. Ezt az egyensúlyi helyzetet nevezzük instabilisnak. A harmadik esetben a golyó ott marad az elmozdítás helyén, azaz nem tér vissza az eredeti helyzetébe, de nem is távolodik el onnan. Ezt az egyensúlyi helyzetet nevezzük stabilisnak Nézzük meg egy példa segítségével, hogy mit is jelentenek ezek a fogalmak a dierenciáálegyenletek megoldásainak körében. 7

http://www.doksihu 2.32 Példa Legyen y 0 (t) = λy , dierenciálegyenlet, y(0) = y0 kezdetiérték feltétetellel. Ha ezt az egyenltet integráljuk, akkor az y(t) = eλt · y0 egyenletet kapjuk megoldásként. Nézzük meg milyen egyensúlyi állapotok állnak fenn λ < 0, λ > 0, illetve λ = 0 esetekben. 1. λ < 0 esetben, ha t +∞, akkor a megoldás határétéke 0, így minden megoldás a 0 egyensúlyi ponthoz közeledik. Ez egy asszimptotikusan stabilis egyensúlyi helyzet. 2. λ > 0 esetben, ha t +∞, akkor a megoldás határértéke ±∞, azaz minden 0-tól különböz® kezdetiérték feltételnek eleget tev® megoldás távolodik 0 egyensúlyi ponttól. Ebben az estben az 0 egyensúlyi pont instabilis 3. λ = 0 esetben a megoldások konstans függvények, azaz nem is közelednek és nem is távolodnak a 0 egyensúlyi ponttól, azaz a megoldás stabilis. 2.4 Dierenciálegyenletek a gyakorlatban A dierenciálegyenletek általában akkor jutnak szerephez,

amikor olyan folyamatot próbálunk modellezni, mely nem diszkrét lépésekben zajlik (mint mondjuk egy sakkjátszma), hanem az id®ben folyamatosan változnak az állapotjelz®k értékei. Ilyen esetekben vagy meggyelések utalnak egy mennyiség és megváltozásának kapcsolatára, vagy feltételeznek egy elméleti relációt a jellemz®k között. Például a természetben populációk növekedésének üteme általában függ magától a populáció nagyságától. Nézzünk egy ilyen példat a gyakorlatban 2.41 Példa 1 Legyen x(t) a populáció mérete t id®pontban. Most azt fogjuk megvizsgálni, hogy egy adott id®intervallumban, milyen mértékkel n® a populáció mérete: x(t + h) − x(t) ≈ h · x(t) · a 8 http://www.doksihu Ebben a modellben nem számolunk a halálozással, csupán a szület® utódok számát vesszük gyelembe. Az eltelt id®t h-val jelöljük, a-val pedig annak a változásnak a mértékét (arányát), mely a növekv® populáció

következtében a szület® utódok számának növekedését mutatja. A fenti egyenletet átrendezve az alábbi dierenciálegyenletet kapjuk x(t + h) − x(t) = a · x(t) h lim x0 (t) = a · x(t) h0 Tehát, ha ezt az egyenletet kiintegráljuk, akkor az alábbi megoldáshoz jutunk: x(t) = eat · x0 Ez az egyenlet azonban csak akkor ad reális eredményt, ha kis id®intervallumokat vizsgálunk Ha hosszú távon szeretnénk populációs modelleket vizsgálni, akkor sajnos ennél egy jóval bonyulultabb egyenletre lesz szükségünk. Például, olyan a-t választunk arányossági tényez®nek, mely függ x(t)-t®l is Például a = r · (K − x(t)), ahol r ismét egy arányossági tényez®, K pedig az eltartó képesség, azaz, hogy egy adott terület egy adott populációnak hány tagját képes "eltartani". Sajnos az ilyen típusú egyenletek már jóval nagyobb m¶veletigénnyel bírnak, mint azt az el®z® egyszer¶ példában láthattuk. 2.42 Példa 2 Ugyancsak a

populáció növekedését leíró folyamat, az úgynevezett "róka-nyúl" modell, ahol már nem csak egy populáció nagyságánának a változását vizsgáljuk, hanem azt, hogy ha egy területen két különböz® faj él, akkor hogyan alakul a populációk mérete. Jelöljük a rókák számát a t id®pntban y(t)-vel, a nyulakét x(t)-vel, a, b, c, d számok pedig pozitív konstansok. x0 (t) = a · x(t) − b · x(t)y(t) y 0 (t) = c · x(t)y(t) − d · y(t) Azaz minden egyes találkozási pontnál, a nyulak száma csökken, illetve a rókák száma "n®", azaz az egyed fejl®dik a tápláléktól. Természetesen itt is gyelembe vehetünk még több paramétert, például, hogy a populációk növekedése milyen tényez®kt®l függ stb, amik tovább nehezítik az egyenlet megoldhatóságát 9 http://www.doksihu 2.5 A dierenciálegyenletek megoldhatósága 2.51 Deníció Kezdetiérték feltétel fogalma: Legyen y0 (t) = f (t, y(t)) egy dier-

enciálegyenlet. Legyen t0 ∈ R, p0 ∈ Rn Egy y megoldás teljesíti az y(t0 ) = p0 kezdeti feltételt, ha átmegy a (t0 , p0 ) ponton. 2.52 Tétel Egy kezdetiérték feladatnak létezik egyértelm¶ megoldása, ha az f : Rn+1 Rn folytonos függvény a 2, 3, . n + 1 -edik változójában Lipschitz folytonos. Tehát a tételb®l biztosan tudjuk, hogy az adott egyenlet megoldható, ám a differenciálegyenleteket kielégít® megoldásfüggvények csak a legegyszer¶bb esetekben fejezhet®k ki zárt alakban. Sok esetben szükségtelen is kiszámolni a konkrét megoldásokat, sokkal többet tudhatunk meg a folyamatokról, ha a megoldások kapcsolatait vizsgáljuk. Más esetben szükséges kiszámítani a megoldás konkrét értékeit Mindkét feladatra számítógépes módszereket használnak, az els® inkább kvalitatív, míg a második kvantitatív eredményt szolgáltat. A dierenciálegyenletek megoldási módszereit három nagy csoportba sorolhatjuk, és ebb®l a három

csoportból mindössze csak egy tartalmazza a pontos megoldások kiszámítását. 1. Analitikus megoldási módszerek Ezek a megoldási módszerek a dierenciálegyenlet pontos megoldásának kiszámítására alkalmasak. A probléma csak az, hogy a dierenciálegyenlteknek csupán egy kis szeletét adják azok a típusok, melyeknek ismerjük a megoldási módszerét, s®t, ha ismerünk is ilyen módszert, a megoldás kiszámítása sok esetben igen költséges. Így más eszközöket kell keresnünk a megoldás kiszámításához, illetve közelítéséhez 2. Kvázianalitikus módszerek Ebbe a kategóriába a Banach xponttételen alapuló iterációs megoldási módszerek tartoznak, ahol y 0 (t) = f (t, y(t)) y(t0 ) = y0 10 http://www.doksihu kezdetiérték problémát a következ® integrálegyenletté írjuk át: Z t y(t) = y(t0 ) + Z 0 t y (s)ds = y0 + t0 f (s, y(s))ds t0 Ekkor a következ® közelít® módszerrel készítünk egy iterációs eljárást, azaz egy y1

, y2 . sorozatot, ahol az n + 1-edik tag a következ®képpen néz ki: Z t f (s, yn (s))ds yn+1 (t) = y0 (t) + t0 Csakhogy a xpont közelít® eljárást minden lépésben közelít® integrálással kell kombinálni, ezért a megoldásunk végül nem egy intervallumon értelmezett függvény lesz, hanem ennek egy diszkretizált alakja. 3. Numerikus módszerek A numerikus módszerek már elég nagy csoportot ölelnek fel. Ezekr®l részletesen olvashatunk a következ® fejezetekben. 2.6 A megoldandó probléma: Legyen f : Rn+1 Rn , olyan folytonos függvény, mely a 2, 3, . n + 1 -edik változójában Lipschitz folytonos. Továbbá legyen u0 ∈ Rn vektor adott Keressük meg azt az u : I Rn dierenciálható függvényt, amelyre igaz, hogy u0 (t) = f (t, u(t)) (2.1) u(0) = u0 (2.2) kezdeiérték feltétel mellett. 11 http://www.doksihu 3. fejezet Numerikus módszerek A numerikus megoldási módszereknek két f®bb osztálya van: az egy- és a többlépéses

módszerek. Az alapvet® különbség az, hogy az egylépéses módszerek csupán egy "lépcs®t" használnak fel a már meglév® n − 1 darabból, azaz yn kiszámításához mindössze az yn−1 közelítést veszi segítségül. Ezzel szemben a többlépéses módszerek az n-edik közelítéshez legalább 2, és legfeljebb n − 1 lépcs®t használnak. Az alábbi fejezetekben az egylépéses numerikus módszerek struktúráját, el®nyeit és hátrányait fogom bemutatni. Egy adott módszer attól lesz diszkrét, hogy a közelít® megoldásokat véges sok pontban keressük. A numerikus módszerek esetében, egy adott [0, b] intervallum t0 < . < tN ekvidisztáns felosztását tekintjük Nézzük meg részletesen, hogyan is épül fel egy numerikus módszer. Legyen u0 (t) = f (t, u(t)) u(0) = u0 kezdeiérték feltétellel. Koordinántánként felírva u0i (t) = fi (t, u1 (t), . un (t)) ui (0) = ui0 Keressük az u(t) = [u1 (t), u2 (t) . un (t)] vekort Tegyük

most fel, hogy n = 1, azaz u0 (t) = f (t, u(t)) f : R2 R, és keressük az u(t) : I R függvényt. Deniáljunk egy rácsháló- függvényt ! Legyen egy tetsz®leges h > 0 esetén ωh := {tn = nh| n = 0, 1 . N } 12 http://www.doksihu Ez lesz az úgy nevezett h lépésköz¶ rácsháló. Deniáljunk az ωh rácshálón egy yh rácshálófüggvényt, és vezessük be a következ® jelölést: yh (tn ) = yn , azaz, yn jelöli azt az értéket, ha a rácshálófüggvényünkbe behelyettesítjük a tn értéket. A cél az, hogy ezt az yn függvényt úgy válasszuk meg, hogy minél közelebb legyen u(tn )-hez, ∀tn ∈ ωh -ra. Ez lesz a numerikus módszerek alapja 3.1 A numerikus módszerek "jósága" A numerikus módszerek pontossága sok paramétert®l függ. Az alábbiakban ezeket a tényez®ket fogjuk vizsgálni 3.11 A numerikus módszerek hibája Legyen u(t) az (2.1)(22) feladat pontos megoldása, ti pontban u(ti ), a továbbiakban jelöljük ui -vel

Hasonlóképpen, legyen yi a ti pontbeli közelít® megoldása az (2.1)(22) feladatnak A két érték közti eltérést, tehát a közelítés hibáját pedig nevezzük di -nek. di := |yi − ui |, feltéve, hogy yi−1 = ui−1 . Ezt nevezzük lokális vagy más néven diszkretizációs hibának A lokális hiba, ahogy azt a neve is jelzi, azt mutatja meg, hogy egy lépés alatt mekkora hiba keletkezik. A globális hibát pedig deniáljuk úgy, hogy ei = |yi − ui |, i = 1, 2.n, tehát itt már nem csak egy lépés alatt keletkez® hibát, hanem n lépés alatt keletkez®t vizsgálunk. 3.12 Numerikus módszerek konvergenciája Azt mondjuk, hogy egy numerikus módszer által el®állított yh (t) rácsháló függvény sorozat, nomodó h lépésközök esetén konvergál az y(t) (2.1)(22) feladat megoldásához, a t∗ ∈ I pontban, ha: • t∗ ∈ ωh , ∀h-ra • limh0 |yh (tn ) − y(t∗ )| = 0, (tn = n · h) 13 http://www.doksihu Általában a numerikus módszert

konvergensnek nevezzük, ha konvergens minden t∗ ∈ I pontban. 3.13 A konvergencia rendje |yh (tn ) − y(t∗ )| = O(hp ), a konvergencia rendjének nevezzük, ezek szerint a numerikus módszer p-ed rendben konvergens. 3.11 Megjegyzés A módszer akkor konvergens, ha limh0 |ei | = 0 3.12 Megjegyzés Egy numerikus módszer csak akkor "jó", ha a módszer konver- gens. A konvegenciát két egyszer¶bben ellen®rizhet® tulajdonsággal lehet garantálni, ezek a 0-stabilitás és a konzisztencia. 3.14 Numerikus módszerek 0-stabilitása Egy numerikus módszer 0-stabil, ha ∃K > 0, melyre |ei | ≤ K(|e0 | + ∀i : 1 ≤ i ≤ N . Pn j=1 |dj |), 3.15 Numerikus módszerek konzisztenciája Azt mondjuk, hogy egy numerikus módszer p-ed rendben konzisztens, ha ∃c > 0 konstans, melyre |di | ≤ c · hp+1 . 3.13 Tétel Ha egy numerikus módszer p-ed rendben konzisztens és 0-stabil, akkor p-ed rendben konvergens. A konvergencia általában nehezen

megállapítható, hiszen kiszámításához szükséges tudni a globális hibát, és a globális hiba kiszámításához tudni kell a pontos megoldást. A megoldások nyílván nem állnak a rendelkezésünkre, hiszen akkor nem lenne szükség közelít® megoldások kiszámítására. Mivel a konvergenciának kulcsfontosságú szerepe van, hiszen csak akkor "jó" egy adott numerikus módszer, ha az konvergens is, ezért valamilyen módon muszáj információt kapnunk a konvergenciáról. Egy adott módszer 0-stabilitásáról és konzisztenciájáról mindig van információnk, és e két adatból már dönthetünk a konvergenciáról is. 14 http://www.doksihu 4. fejezet A legegyszer¶bb megoldások 4.1 Az explicit Euler módszer Ennek a módszernek az alapötlete igen egyszer¶. Nevezetesen: egy érint® segítségével közelítjük a megoldást a következ®képpen: Húzzuk be az y(0) = y0 kezdetiérték érint®jét, majd kössük össze egy tetsz®leges t1 ponttal,

így megkapjuk az (x1 , y1 ) pontot. Ezek után, ugyanezzel az eljárással, kapjuk meg (t2 , y2 ), . , (tn , yn ) közelít® pontokat Azaz ha ezt felírjuk általánosan, akkor az ábrán látható Φ szög tangensét a következ® képpen kaphatjuk meg: 15 http://www.doksihu 4.1 ábra Az explicit Euler módszer tg(Φ) = yi+1 − yi ti+1 − ti azaz f (ti , yi ) = yi+1 − yi ti+1 − ti yi+1 = yi + f (ti , yi )(ti+1 − ti ) ahol (ti+1 − ti ) = h lépésközzel. Egy másik megközelítséb®l: Legyen a [0, b] interval16 http://www.doksihu lum felosztása a következ®: 0 = t0 < t1 < . < tn−1 < tn = b Legyen hi = ti −ti−1 az i-edik lépésköz. Az egyenlet tn id®pontbeli pontos megoldását jelöljük u(tn )-nel, és a tn id®pontbeli közelít® megoldását pedig yn -nel. A kezdeti érték problémánál, mindig tudjuk, hogy milyen értéket vesz fel a megoldás a t0 id®pontban. Ebben az esetben tegyük fel, hogy ismerjük a tn−1 id®pontbeli

yn−1 közelítést, és ebb®l határozzuk meg a tn -beli yn közelítést. yn+1 − yn = f (tn , yn ) h y(0) = y0 átrendezve, tehát: yn+1 = yn + h · f (tn , yn ). Ezt nevezzük explicit euler módszernek 4.11 Példa Legyen y 0 (t) = 3e−t − 0, 4y(t) y(0) = 5 y(3) =? h=3 Helyettesítsünk be az Explicit Euler módszer egyenletébe: y1 = 5 + f (0, 5) · 3 y1 = 5 + (3e−0 − 0, 4 · 5) · 3 y1 = 5 + (3 − 2) · 3 = 8 Azonban az egyenlet pontos megoldása ebben az esetben 2,763, azaz a numerikus módszer lokiális hibája 5,237, ami igen nagy mérték¶ hibát jelent. Próbáljuk meg a lépésköz csökkentésével, és a lépésszámok növelésével pontosítani a közelít® megoldást. Legyen h = 1, 5 és közelítsük a megoldást két lépésben, el®ször y(0) y(1, 5), majd y(1, 5) y(3)! 1.lépés x0 = 0; y0 = 5; h = 1, 5 17 http://www.doksihu y1 = 5 + f (0, 5) · 1, 5 y1 = 5 + (3e−0 − 0, 4 · 5) · 1, 5 y1 = 5 + 1 · 1, 5 = 6, 5 = y(1, 5) 2. lépés x1

= 1, 5; y1 = 6, 5; h = 1, 5 y2 = 6, 5 + f (1, 5; 6, 5) · 1, 5 y2 = 6, 5 + (3e−1,5 − 0, 4 · 6, 5) · 1, 5 y2 = 6, 5 + (−1, 93061) · 1, 5 = 3, 604 = y(3) Itt a pontos megoldástól való eltérés már mindössze 0,841, és ez az eredmény még tovább csökkenthet®, a lépésköz csökkentésével, illetve a lépésszámok növelésével. 4.12 A hibaegyenlet Mivel yn a pontos megoldás közelítése, ezért felírhatjuk yn = un + zn alakban, azaz pontos megoldás + hiba alakban. Helyettesítsük be ezt a formulát az explicit Euler módszer egyenletébe: un+1 − un zn+1 − zn + − f (tn , un + zn ) = 0 h h un+1 − un zn+1 − zn − f (tn , un ) + f (tn , un ) − f (tn , un + zn ) + =0 h h zn+1 − zn −un+1 − un = + f (tn , un ) − f (tn , un ) + f (tn , un + zn ) h h Ez az úgynevezett hiba egyenlet, ami felírható Ψn,1 + Ψn,2 = 0 alakban, ahol Ψn,1 = un+1 −un + f (tn , un ), a lokális approximációs, vagy más néven reziduális hiba. h Tehát, ha Ψn,1

-be a pontos megoldást helyettesítjük be, akkor nyílvánvalóan 0t kapunk. Így a lokális approximációs hiba azt fejezi ki, hogy az adott numerikus módszer milyen pontosan közelíti a folytonos (2.1)(22) feladatot Ezek alapján pontosítsuk a konzisztencia denícióját: 4.11 Deníció Egy numerikus módszer konzisztens, azaz approximálja a folytonos feladatot, ha limh0 Ψn,1 = 0 Nézzük meg, mit kapunk, ha Taylor sorba fejtjük az u(tn )-et. 18 http://www.doksihu 4.13 Az explicit Euler módszer konzisztenciája un+1 − un + f (tn , un ) = h Tudjuk, hogy a lépésköz hossza h, ezért tn+1 felírható tn + h alakban. Ψn,1 = u(tn ) − u(tn ) + hu0 (tn ) + 21 h2 u00 (tn ) + O(h3 ) + u0 (tn ) = h 1 − hu00 (tn ) + O(h2 ) 2 A konzisztencia rendjére vonatkozó tétel alapján az explicit Euler módszer konzisztens és rendje 1. 4.14 Az Euler módszer lokális hibája Az i-edik lépésben elkövetett lokális hibát jelöljük di -vel,amit nyílván a pontos és a

közelít® megoldás különbségeként fogunk meghatározni. di = yi − ui yi−1 = h · f (ti−1 , yi−1 ) − u(ti ) u(ti−1 ) = h · f (ti−1 , u(ti−1 )) − u(ti ) u(ti−1 ) = h · f (ti−1 , u(ti−1 )) − u(ti−1 + h) Fejtsük Taylor sorba u(ti−1 + h)-t! u(ti−1 + h) = u(ti−1 ) + h · u0 (ti−1 ) + h2 00 · u (ti−1 ) + . + O(h3 ) 2 Ezt behelyettesítve megkapjuk a módszer lokális hibáját: h2 h2 00 · u (ti−1 ) + o(h3 ) = − · u00 (ti−1 ) + O(h3 ) 2 2 Azaz a lokális hiba másodrend¶, ami a gyakorlatban azt jelenti, hogy a választott lépésköz zsugorításakor annak második hatványával zsugorodik a hibára adott fels® becslés. 19 http://www.doksihu 4.15 Az Euler módszer 0-stabilitása Ha ei (globális hiba) abszolút értékéhez létezik olyan K konstans, hogy a globális P hiba abszolút értéke kisebb legyen, mint K(|e0 | + nj=1 |dj |) , akkor a módszer 0stabil. ei+1 = ei + h(f (ti , y(ti )) − f (ti , yi )) + dj+1 |ei+1 |

= |ei | + hα|ei | + dj+1 = = (1 + hα)|ei | + dj+1 ≤ (1 + hα)(|ei−1 | + hα|ei−1 | + dj ) + dj+1 ei+1 ≤ (1 + hα)|ei | + di+1 ≤ (1 + hα)2 |ei−1 | + (1 + hα)dj + dj+1 ≤ . ≤ ≤ (1 + hα)i+1 |e0 | + X (1 + hα)i−j |dj+1 | (1 + hα)k ≤ ehαK ≤ etk α ≤ eαb ≤ eαb |e0 | + i X eαb |dj+1 | = eαb (|e0 | + j=0 i+1 X |dj |) j=1 Azaz K := eαb választással a módszer 0-stabil. 4.12 Megjegyzés tk -val a k-adik osztópontot jelöltem 4.16 A "jó" lépésköz megválasztása A lépésköz kiválasztása nagy szerepet tölt be a numerikus módszer pontosságában. Az adott egyenlet és a választott numerikus módszer azonban korlátozza a lépéshosszt illet® választási lehet®ségek számát. Így nemcsak megfelel® lépésközr®l, hanem megfelel® numerikus módszerr®l is beszélnünk kell, hiszen a módszert úgy kell megválasztanunk, hogy h lépésközre vonatkozó korlátozások száma és mértéke minimális legyen Ezek a

korlátok általában szoros összefüggésben állnak a numerikus módszer stabilitásával. 20 http://www.doksihu 4.17 A teszt egyenlet Tekintsük azt az esetet, mikor az egyenlet y 0 = Ay alakú, ahol A ∈ Rnxn mátrix. Ha A diagonizálható, akkor ezt az egyenletet felírhatjuk az alábbi módon: w0 = Dw ahol D egy diagonális mátrix. Ha mindezt koordinátákra bontjuk, akkor a következ® n ismeretlenes egyenletrendszert kapjuk: w10 = d1 w1 w20 = d2 w2 . . wn0 = dn wn Itt természetesen di ∀i-re a D diagonális mátrix sajátértékei. Ezt az egyenletrendszert reprezentáljuk a y 0 = λy teszt egyenlettel, ugyanis csak akkor lesz az eredeti egyenlet stabil, ha az egyenletrendszer ∀di sajátértékére stabil a teszt egyenlet. 4.18 Abszolút stabilitás Ha tudjuk, hogy y(0) = c, ahol c > 0 konstans, akkor a teszt egyenlet pontos megoldása y(tn ) = ceλtn lenne. Ezek alapján három esetet különböztethetünk meg: • ha Re(λ) > 0, akkor |y(t)| = ceRe(λ)t

exponenciálisan növekszik t szerint. Ez egy instabilis helyzet. • ha Re(λ) = 0, akkor a megoldás oszcillál. • ha Re(λ) < 0, akkor |y(t)| exponenciálisan csökken, így a megoldások egyre közelednek egymáshoz. Ez egy asszimptotikusan stabil helyzet 4.19 A-stabilitás Vizsgáljuk tovább a tesztegyenletet a harmadik esteben: y(t) = λy(t), ahol λ ∈ C megoldása 0-hoz tart t ∞ határesetben, (Re(λ) < 0). 21 http://www.doksihu 4.13 Deníció Egy numerikus módszer A-stabil, ha a teszt egyenletre alkalmazva rendelkezik a fenti tulajdonsággal, lépéshossztól függetlenül. Másképpen megfogalmazva: Legyen y(1) = 1 a kezdeti feltétel. Ha a tesztegyenletb®l kapott y1 , y2 , . sorozatra limn∞ yn = 0, ∀h-ra, akkor a módszer A-stabil Ez egy nagyon er®s feltétel. Az ilyen módszerek a gyakorlatban nem mutatnak stabilitási problémakat. 4.110 Az explicit Euler módszer abszolút stabilitási tartománya A leírtak alapján már tudjuk, hogy egy

módszernél fontos tulajdonság az Astabilitás. Vizsgáljuk meg az explicit Euler módszert ebb®l a szempontból yn+1 = yn + h · f (tn , yn ), (f (t, y) = λy) yn+1 = yn + h · λyn = yn+1 = (1 + hλ)yn = (1 + hλ)2 yn−1 = . = (1 + hλ)n y0 Tehát, most azt kell megvizsgálnuk, hogy yn sorozat 0-hoz tart-e? yn+1 = (1 + hλ)n y0 0 n ∞ Az alábbi állítás akkor igaz, ha |1 + hλ| < 1, ez például valós λ-ra is megkötés: h < − λ2 . Tehát az explicit Euler módszer nem A-stabil 22 http://www.doksihu 4.2 ábra Az explicit Euler módszer abszolút stabilitási tartománya 4.111 A módszer hiányosságai A módszer nagy el®nye, hogy nagyon egyszer¶, illetve a m¶veletigénye alacsony. Ám sajnos az explicit Euler módszer nyílván nem a "legjobb" közelít® megoldást adja hiszen • csak els®rend¶ • vannak stabilitási problémái, azaz, nem A-stabil 4.2 Az implicit Euler módszer Nézzünk most az (2.1)(22) feladtra egy másik

közelítési módszert Írjuk fel u(tn )et a következ® alakban: u(tn ) = u(tn+1 − h) alakban, és fejtsük Taylor sorba: u(tn+1 − h) = u(tn+1 ) − hu0 (tn+1 ) + O(h2 ) azaz, u(tt+1 ) − u( tn ) = u0 (tn+1 ) + O(h2 ) = f (tn+1 , u(tn+1 )) h yn+1 − yn = f (tn+1 , yn+1 ) h y(0) = y0 23 http://www.doksihu Az egyenletet rendezzük yn+1 -re, így megkapjuk az implicit Euler módszer általános alakját: yn+1 = yn − hf (tn+1 , yn+1 ) Nézzük meg, mit kapunk n = 0 behelyettesítésével: y1 − y0 = f (h, y1 ) h Így láthatjuk, hogy yn kiszámításához, egy általános nemlineáris algebrai egyenlet megoldása szükséges, melyre több módszert is ismerünk (például Newton iteráció). Ez a módszer ugyan költségesebb, mint a fent említett explicit Euler módszer, azonban sok szempontból praktikusabb annál. 4.21 Az implicit Euler módszer konzisztenciája Helyettesítsünk be a lokális approximációs hiba képletébe, majd fejtsük Taylor sorba u(tn+1 −

h)-t! u(tn+1 ) − u(tn+1 − h) + u0 (tn+1 ) h u(tn+1 ) − (u(tn+1 ) − hu0 (tn+1 )) + O(h2 ) + u0 (tn+1 ) =− h = O(h) Ψn,1 = − Az implicit Euler módszer is 1. rendben konzisztens, tehát nem sikerült az explicit formulánál pontosabb közelítést létrehozni, legalábbis a konzisztencia szempontjából Nézzük meg mi a helyzet a stabilitási tartományával 4.22 Az implicit Euler módszer abszolút stabilitási tartománya A teszt egyenletb®l megkapjuk, hogy y 0 (t) = f (t, y(t)) = λy(t) Helyettesítsük be ezt az implicit Euler módszer n + 1-edik egyenletébe: yn+1 = yn + hλyn+1 yn+1 (1 + hλ) = yn yn+1 = 24 1 · yn 1 − hλ http://www.doksihu 1 Legyen z := hλ, így feltételként az 1−z < 1 egyenletet kapjuk. Ebb®l már látható az implicit Euler módszer stabilitási tartománya: ha 1 < |1 − z|, azaz ha h > 0. Tehát az implicit Euler módszer A-stabil. Nézzük meg egy gyakorlati példán keresztül, hogy miért fontos szempont az

A-stabilitás. 4.3 ábra Az implicit Euler módszer stabilitási tartományát a satírozott területen kívüli rész jelöli 4.23 Merev dierenciálegyenletek A merev dierenciálegyenletek és ezek megoldási módszerei egy elég nagy témát ölelnek fel ráadásul nincs általánosan elfogadott egzakt deníció, ezért csak néhány példával érzékeltetném, hogy mely esetekben beszélünk merev derenciálegyenletelr®l: y 0 = Ky + f (t), ahol y ∈ Rm , K ∈ Rmxm és ∃ "nagy" abszolútérték¶ negatív sajátértéke, vagy y 0 = f (t, y) és f Jacobi mátrixának létezik "nagy" abszolútérték¶ negatív sajátértéke. Az ilyen típusú egyenleteknél nagy szerepet játszik, hogy a használni kívánt numerikus módszer A-stabil vagy nem. A továbbiakban nézzünk olyan módszereket, melyek alkalmasak a merev dierenciálegyenletek megoldására. 25 http://www.doksihu 4.24 Példa Az alábbi egyenletet programozzuk le mindkét módszer szerint

Matlab-ban! y 0 (t) = −100(y(t) − sin(t)) 4.4 ábra Explicit Euler 26 http://www.doksihu 4.5 ábra Implicit Euler Láthatjuk, hogy a megadott merev dierenciálegyenletet, hogyan közelíti meg egy olyan módszer, mely A-stabil, és egy olyan, amelyik nem. A különbség szembet¶n® 4.25 Az implicit Euler módszer el®nyei és hátrányai Az implicit módszer tehát stabilitását tekintve jobb az explicitnél, viszont ez még mindig csak egy els®rend¶ közelítés, így a következ® lépésben próbáljunk meg olyan módszert keresni, amely legalább egy másodrend¶ approximáció. 27 http://www.doksihu 4.3 Szimmetrikus/ trapéz módszer Próbáljuk meg a következ®t: vegyük az explicit és az implicit Euler módszer számtani közepét! yn+1 − yn 1 = [f (tn , yn ) + f (tn+1 , yn+1 )] h 2 h yn+1 = yn + [f (tn , yn ) + f (tn+1 , yn+1 )] 2 4.31 A szimmetrikus módszer konzisztenciája u(tn+1 ) − u(tn ) 1 + [f (tn , yn ) + f (tn+1 , yn+1 )] = h 2 h2 00 0 u(tn )

+ hu (tn ) + 2 u (tn ) + O(h3 ) − u(tn ) 1 0 + [u (tn ) + u(tn+1 )] = =− h 2 h 1 = −u0 (tn ) − u00 (tn ) + o(h2 ) + [u0 (tn ) + u0 (tn ) + hu00 (tn ) + O(h2 )] = O(h2 ) 2 2 Ψn,1 = − Tehát sikerült pontosabb módszerhez jutnunk, hiszen a trapéz módszer már 2. rendben konzisztens Vizsgáljuk meg, hogy a stabilitási tartomány hogyan alakul ebben az esetben! 4.32 A trapéz módszer abszolút stabilitási tartománya A szokásos módon helyettesítsünk be a teszt egyenlet alapján a szimmetrikus trapéz formulába. h yn+1 = yn + (λyn + λyn+1 ) 2 1 yn+1 = yn + hλ(yn + yn+1 ) 2 1 1 yn+1 = yn + hλyn + hλyn+1 2 2 1 1 (1 + hλ)yn = yn+1 − (1 − )hλyn+1 2 2 1 + 12 hλ yn+1 = yn 1 − 12 hλ 2+z Azaz a stabilitási tartomány a következ® lesz: | 2−z | < 1 |2 + z| < |2 − z|. Ezt az egyenl®tlenséget próbáljuk tovább alakítani, hogy z -re megoldást kapjunk. Írjuk fel 28 http://www.doksihu z := a + bi alakban, és oldjuk meg az

egyenl®tlenséget: (|2 + a + bi|)2 < (|2 − (a + bi)|)2 (2 + a)2 + b2 < (2 − a)2 + b2 Tehát,a < 0, így Re(z) < 0, azaz a trapéz módszer is A-stabil. 4.6 ábra A szimmetrikus módszer stabilitási tartománya 29 http://www.doksihu 5. fejezet Runge-Kutta típusú módszerek Vizsgáljuk meg újra a fent említett módszereket egy másik megközelítésb®l, ahol a közelít® eljárások alapját az úgynevezett kvadratúra formulák adják, vagyis dolgozzunk numerikus integrálokkal! 5.1 A kvadratúra formulák A Runge-Kutta módszerek ötletének alapját a kvadratúra formulák, azaz a numerikus integrálás szabályai adják. Nézzük meg, mit is jelent ez pontosan! Vegyünk egy [a,b] intervallumot, és készítsünk egy t1 , t2 , . tn ekvidisztáns felosztást Az f függvény az [a, b] intervallumon felírhatjuk az n − 1 darab részintervallum integráljának összegeként is, azaz: Z b f (t) = a n Z X i=1 ti+1 f (t)dt ti Ez akkor hajtható

végre, ha ismerjük a függvény [(t1 , f (t1 )), . , (tn , f (tn ))] alappontjait Ezeket a pontokat felhasználva el®állítunk egy Lagrange interpolációs polinomot, mellyel közelítjük az eredeti függvényünket, majd ezt a polinomot integráljuk: Y t − ti tj − ti Z b ωj = Lj (t)dt Lj (t) = a 30 http://www.doksihu Végül a kapott értékeket szorozzuk f (t1 ), . f (tn ) értékekkel, így megkapjuk f (t) közelít® integrálját. Z b f (t)dt ≈ a n X ωj f (tj ) j=1 n Vegyük az y(tn ) − y(tn−1 ) = ttn−1 y 0 (t)dt egyenletet. A görbe alatti területet approximálhatjuk bal (Explicit -Euler) és jobb (Implicit Euler) közelít® integrálokkal Ezek az els®rend¶ módszerek. Próbáljuk meg az explicit Euler módszert "nomítani"! Eddig yn kiszámításához csak az yn−1 értéket használtuk fel. Azonban ha felvennénk egy köztes pontot, például y(n−1)/2 -t, akkor elméletileg pontosabb megoldást kapnánk. Vizsgáljuk meg ezt az

esetet. R h ỹ(n−1)/2 = yn−1 + f (tn−1 , yn−1 ) 2 yn = yn−1 + hf (t(n−1)/2 , ỹ(n−1)/2 ) Ezt a módszert javított Euler módszernek nevezzük. Nézzük meg, hogy mennyivel közelíti jobban a pontos megoldást! 5.11 A javított Euler módszer lokális hibája y(tn ) − y(tn−1 ) h − f (tn−1 ), y(tn−1 + f (tn−1 , y(tn−1 ))) h 2 2 h h h h2 = y 0 + y 00 + y 000 − (f + (ft + fy f ) + (ftt + 2fty f + fyy f 2 )) + O(h3 ) 2 6 2 8 Ψn,1 = Tehát a módszer 2. rendben konzisztens, így valóban sikerült növelni a pontosságot Ennek alapján próbáljunk meg minél magasabb rend¶ módszereket létrehozni Ehhez vezessük be az alábbi jelöléseket. k1 :=f (tn , yn ) k2 :=f (tn + 0, 5h, yn + 0, 5h · k1 ) yn+1 =yn + hk2 31 http://www.doksihu A javított Euler módszer példájára, általánosítsuk ezt jelölést m lépésre. k1 = f (tn , yn ) k2 = f (tn + a2 h, yn + hb21 k1 ) k3 = f (tn + a3 h, yn + h(b31 k1 + b32 k2 )) . . km = f (tn + am h, yn +

h(bm1 k1 + . + bm,m−1 km−1 )) yn+1 = yn + h(c1 k1 + c2 k2 + . + cm km ) Ez az úgynevezett általános m lépéses (explicit) Runge-Kutta módszer, melyet 1900 körül dolgozotak ki Karl Runge és Martin Kutta német matematikusok. A javított Euler módszernél láthattuk, hogy a "javítással" sikerült növelni a módszer rendjét. Itt az a2 , , am ; bik ; és c1 , , cm tetsz®leges paraméterek, illetve (bik )-k által el®állított B mátrixot a Runge-Kutta módszer mátixának, ak -kat a módusainak, P ci -ket pedig a mátrix súlyainak nevezzük. (Éppen ezért m i=1 ci = 1.) A Runge-Kutta formula lényege, hogy minél magasabb rend¶ módszereket tudjunk létrehozni. Az adott módszerhez kapott konstansokat egy úgynevezett Butcher-tábla foglalja össze. (Butcher tableau, John C. Butcher neve után) a2 b21 0 a3 b31 a4 am . . . b32 0 0 b41 b42 b43 . 0 0 0 bm1 bm2 bm3 . bm,m−1 c1 c2 c3 . cm . . . . Az egyszer¶ség

kedvéért vezessük be az alábbi jelöléseket: yn+1 = yn + h · φ(xn , yn , h) m X φ(x, y, h) = = ci ki ki = f (tn + ai h, yn + i=1 m X bij kj ) (5.1) (5.2) (5.3) i=1 5.11 Megjegyzés Ahhoz, hogy p-ed rend¶ RK módszert kapjunk a a2 , . , a m ; bik ; és c1 , . , cm együtthatókat úgy kell megválasztani, hogy a módszer lokális hibája 32 http://www.doksihu p + 1 rend¶ legyen. Azaz, egy p-ed rend¶ a következ® egyenl®ségnek kell teljesülnie: d(t, y, h) = u(t + h) − u(t) − hφ(t, y, h) = O(hp+1 ) Ahol d(t, y, h) a lokális hibát, u(t) pedig az egyenlet pontos megoldását jelöli. Els® lépésként fejtsük Taylor sorba d(t, y, h)-t, u(t + h) − u(t)-t, és φ(t, y, h)-t h szerint h = 0 pontban.(Ahhoz, hogy az egyenl®ség teljesüljön, az kell, hogy a d lokális hiba Taylor-sorának els® p tagja t¶njön el.) d(t, y, h) = ∞ X ∂ id i=1 ∂hi (5.4) · hi z (p) (t) u00 (t) + . + hp + O(hp+1 ) 2! p! h2 φ(t, y, h) = φ(t, y, 0) + hφ0

(t, y, h) + φ00 (t, y, h) + . 2! hp−1 p−1 + φ (t, y, h) + O(hp ) (p − 1)! P ∂id Ahhoz, hogy a módszer p-ed rend¶ legyen az szükséges, hogy pi=1 ∂h i = 0 u(t + h) = u(t) + hu0 (t) + h2 (5.5) (5.6) (5.7) feltétel teljesüljön. Második lépésként szorozzuk be mindkét oldalt h-val u00 (t) z (p) (t) + . + hp + O(hp+1 )− 2! p! hp − φh(t, y, 0) − h2 φ0 (t, y, h) − . − φ(p−1) (t, y, h) − O(hp ) (p − 1)! 1 d(x, y, 0) = h(u0 (t) − φ(t, y, 0)) + h2 ( u00 (t)− 2 0 p 1 (p) − φ (t, y, 0)) + . + h ( z − φ(p−1) (t, y, 0) + O(hp+1 ) p! d(t, y, h) = hu0 (t) + h2 Így egy p egyenletb®l álló egyenletrendszert kapunk, mely megoldható, ha ai = Pp j=1 bij . 5.12 Az RK módszerek 0-stabilitása A Runge-Kutta formulát az 5.1 egyenletben, olyan alakra hoztuk, mely nagyon hasonlít az explicit Euler módszerhez. Ennek alapján, hasonló levezetéssel megkaphatjuk, hogy a Runge-Kutta módszer 0-stabil Természetesen most is beszélhetünk

explicit és implicit módszerekr®l. Nézzünk most mindkét esetre néhány példát! 33 http://www.doksihu 5.2 Az implicit Runge-Kutta módszerek Az implicit és explicit módszerek közti legszembet¶n®bb különbség, a Butcher tábla, hiszen implicit esetben a B mátrix nem feltétlenül alsó háromszög mátrix. Éppen ezért az implicit formulákat sokszor nevezik általános RK módszereknek. 5.21 Az implicit formula el®nyei és hátrányai Az implicit módszer nagy hátránya, hogy a (4.3) egyenletrendszer egy nemlineáris rendszer lesz ki -kre, amit minden egyes lépésben iterációval kell megoldani, így a m¶veletigénye igen nagy. El®nyei közül a két legjelent®sebbet említeném meg: • stabilitás tulajdonságuk • 5.21 Tétel Minden p ≥ 1 számhoz létezik pontosan egy 2p-ed rend¶ implicit Runge-Kutta módszer. 5.3 Az explicit Runge-Kutta formulák 5.31 Deníció Egy Runge-Kutta módszer explicit, ha bik = 0; i < k Ekkor a ki -k explicit módon

számolhatóak. 5.31 Az els® rend¶ Runge-Kutta formulák Az els®rend¶ RK módszert fel tudom írni a következ® alakban: φ(x, y, h) = c1 k1 = f (t, y), illetve az el®z®ek alapján u(t + h) − u(t) = hu0 (t) + O(h2 ). Ezek alapján írjuk fel az els®rend¶ módszerek lokális hibáját. d(t, y, h) = u(t + h) − u(t) − hφ(t, y, h) = hf (t, y) − hc1 + f (t, y) + O(h2 ) = h(1 − c1 )f (t, y) + O(h2 ) Az egyenletet megoldva azt kapjuk, hogy c1 = 1 választással kapunk els® rend¶ módszert, és ezt visszahelyettesítve a yn+1 = yn +hf (tn , yn ) formulát, azaz az explicit Euler módszert kapjuk. Ez ez egyetlen els® rend¶ (explicit) RK módszer 34 http://www.doksihu 5.32 A megoldhatóság feltételei A ai = pj=1 bij feltétel mellett, minden p-ed rend¶ formulának vannak csak rá vonatkozó feltételei is. Az els®rend¶ módszer pontosan akkor oldható meg, ha teljesül a P ce = 1 (5.8) feltétel, ahol e az (1 . 1)T vektort jelenti, c pedig a Butcher tábla

ci elemib®l álló vektor. 5.33 A másodrend¶ Runge-Kutta formulák A másodrend¶ Runge-Kutta formulák pontosan akkor oldhatók meg, ha teljesül az 5.7 feltétel, és c·a= 1 2 (5.9) egyenl®ség is. (Az aj elemekb®l álló vektort jelöltük a-val) Hasonlóképpen mint az els®rend¶nél φ(t, y, h) = c1 k1 + c2 k2 = c1 f (t, y) + c2 f (t + a2 h, yn + hb21 f (t, y)) Ha Taylor sorba fejtjük k2 -t, és d(t, y, h)-t, akkor megkapjuk h0 , h1 , h2 együtthatóit. (A Taylor sorfejtés ebben az esetben igen hosszadalmas, ezért csak a megoldásokat írom le). b1 = 0 1 − c1 − c2 f = 0 1 1 ( − c2 a2 )∂t f + ( − c2 b21 )∂y f = 0 2 2 Így a másodrend¶ RK módszerre a következ® együtthatókat kapjuk. a1 = 0 a2 = b21 1 a2 = 2c2 c1 + c2 = 1 Látható, hogy az egyenletünk határozatlan, c1 , c2 szabad paraméterek, így igen sok másodrend¶ formula létrehozható. 35 http://www.doksihu 5.34 Példák másodrend¶ formulákra 0 0 0 1 1 0 2 2 0 1 0 0 0 2 2 0 3 3 0

0 0 1 1 0 1 4 1 2 3 4 1 2 5.35 Harmadrend¶ RK formulák A harmadrend¶ RK formulák megoldhatóságának feltételei, az 5.7 és az 58 feltételek, tehát minden, ami az alacsonyabb rend¶ formuláknál kellett és teljesülnie kell az 1 3 1 cBa = 6 (5.10) ca2 = (5.11) egyenleteknek is, ahol B a Butcher tábla bij elemeib®l álló mátrixa. 0 0 0 0 0 0 0 0 1 2 1 2 0 0 0 0 3 3 Példák harmadrend¶ formulákra: 2 2 2 1 -1 2 0 0 32 0 3 1 6 2 3 1 6 1 4 3 8 3 8 5.36 A klasszikus negyed rend¶ Runge-Kutta formula A megoldhatóság feltételei természetesen az összes eddigi feltétel és az alábbiak együttes teljesülése: 1 4 1 cBa2 = 24 1 cABc = 8 ca3 = 36 (5.12) (5.13) (5.14) http://www.doksihu A-val az ai értékekb®l álló mxm-es diagonális mátrixot jelöltem. k1 = f (tn , yn ) h h k2 = f (tn + , yn + k1 ) 2 2 h h k3 = f (tn + , yn + k2 ) 2 2 k4 = f (tn + h, yn + hk3 ) h yn+1 = yn + [k1 + 2k2 + 2k3 + k4 ] 6 0 0 0 1 1 0 2 2 1 0 12 2 1 0 0 0 0 0 1 0 0 0 0

1 6 1 3 1 6 1 3 Ha jobban megnézzük ezt a módszert, akkor láthatjuk, hogy az y(tn ) − y(tn−1 ) közelít® integrál a Simpson formulával felírva közel hasonló eredményt ad: y(tn ) − y(tn−1 ) ≈ h 0 (y (tn−1 ) + 4y 0 (t(n−1)/2 ) + y 0 (tn )) 6 Ha összevetjük az els®-, a másod-, a harmad-, és a negyed rend¶ formulákat , akkor azt a hasonlóságot fedezhetjük fel, hogy p = 1, 2, 3, 4 rend¶ formuláknál a módszerek lépésszáma is rendre m = 1, 2, 3, 4. Valójában a lépések számából nem következik a módszer rendje, amit már m = 5 lépés esetén is láthatunk. lépések száma 1 2 3 4 5 6 7 8 9 10 módszer rendje 1 2 3 4 4 5 6 6 7 7 5.37 Az explicit Runge-Kutta módszer abszolút stabilitási tartománya A Runge-Kutta módszerek családja széleskörben elterjedt közelítési eljárás, annak ellenére, hogy ezeknek a módszereknek is vannak stabilitási problémáik. Abszolút stabilitási tartományuk korlátos, így sajnos nem

A-stabilak 37 http://www.doksihu 5.38 RK módszerek összehasonlítása Nézzük meg az alábbi példán keresztül, hogy milyen pontossággal közelít egy els®-, egy másod-, illetve egy negyed rend¶ Runge-Kutta formula: y 0 (t) = −5t(y(t))2 + 5/t − 1/(t2 ), y(1) = 1 5.1 ábra Els®rend¶ RK formula 5.2 ábra Másodrend¶ RK formula 38 http://www.doksihu 5.3 ábra Negyedrend¶ RK formula Els®rend¶ RK Másodrend¶ RK Negyedrend¶ RK maximális hiba 0,03333 0,56108 0,00533 39 http://www.doksihu 6. fejezet Összefoglalás Az említett módszerek egyt®l egyig alkalmasak közönséges dierenciálegyenletek közelít® megoldásainak kiszámításához. Természetesen ezek a módszerek nem minden esetben "egyformán jók", de a példák alapján tehetünk néhány általános megállapítást: • Egy módszer minél magasabb rend¶, annál pontosabban közelíti a megoldást kivéve, • ha nem megfelel® módszert használunk, azaz implicit helyett

explicitet, vagy fordítva. • Merev dierenciálegyenletek közelítésénél sokkal fontosabb a numerikus mód- szer abszolút stabilitási tartománya, a módszer rendjénél. • Az implicit módszereknél gyelembe kell venni a nagy m¶veletigényüket. A téma kimeríthetetlen, hiszen az említett módszereket csupán az egylépéses közelítések családjából válogattam, és csupán az els®rend¶ egyváltozós kezdetiérték problémát vizsgáltuk meg. A többlépéses módszerekr®l, a merev dierenciálegyenletekr®l, vagy akár a peremérték feladatok megoldásáról még sok-sok oldalt lehetne megtölteni. Irodalomjegyzék: Simon L. Péter, Tóth János: Dierenciálegyenletek Uri M Ascher, Linda R. Petzold: Computer Methods of Ordinary Dierental Equations and Dierential Algebric Equations 40