Beentjes, P.A.

Dekker, K.

NARCIS 

Een 5e orde 6-punts Runge-Kutta formule met optimale stabiliteitsgrens

CWI's Institutional Repository

REKENAFDE LI NG stichting mathematisch centrum P.A. BEENTJES en K. DEKKER EEN 5e ORDE 6-PUNTS RUNGE-KUTTA FORMULE MET OPTlMALE STABILITEITSGRENS NR 27/72 RA ~ MC JUNE 2e boerhaavestraat 49 amsterdam MATHEMATl~l'i AMSTEFlDAfli CENU.UM PJUn:ted a.:t the Ma.:the.ma:ti.c.ai. Cen:tJL.e, 49, 2e Bavr.haa.ve6:tJr.aa.:t 49, Am6teJtda.m. The Ma.:the.ma;U.c.ai. Cen.t.tz.e, 6au.nded the 11-th a 6 Febtw.aJLy 1946, ,i.l, a. nan.-PM 6U .ln-6:tLtuti.an cumlng a.:t the pJtama.ti..an a 6 pUll.e malhema;U.e6 a.nd .lt6 a.ppUc.a..t..lan-6. It ,i.l, -6pan6aJted by the Ne.t.hV'l.i.a.nd6 GaveJLnmen:t th/tough the Ne.t.helti.a.nd6 01tga.vu.za.t.lan 601t the Adva.nc.e.men:t 06 Pu.Jte Re6ea.Jtc.h (Z.W.O.), by the Mu.n.lclpai..l:ty 06 Am6te1tda.m, by the Un.lveM.l.t.y 06 Am6teJtda.m, by the F1tee Un.lveM.l.t.y a.:t Am6te1tda.m, a.nd by .lndU6:t.lUe6. Inhoud 1. Inleiding e 2. 5 orde Runge-Kutta methoden met een schatting voor de discrepantie 3. De twee-parameter farnilie van oplossingen van Zonneveld 4. Het stabiliteitspolynoom van een RK(5,6) methode 5. Een RK(5,6) methode met maximale reele stabiliteitsgrens 6. Stabiliteitsgebieden bij RK(5,6) methoden 7. Testvoorbeeld blz. 2 6 8 10 12 14 1. Inleiding In van der Houwen [2] ziJn een aantal expliciete gestabiliseerde Runge-Kutta formules, met een orde van ten hoogste 4, gegeven. Als aanvulling daarop geeft dit rapport een 5e orde gestabiliseerde Runge-Kutta formule die voor vergelijkingen van het type ( 1. 1) y' = f(x,y), y(O) = y0, een grotere integratie-stap toelaat dan de bestaande 5e orde 6-punts Runge-Kutta (RK (5,6)) methoden. De constructie van dit Runge-Kutta schema is gebaseerd op de manier waarop Zonneveld [5] zijn RK(5,6) formule afleidde. 2 2. 5e orde Runge-Kutta methoden met een schatti~g voor de discrepantie Een m-punts Runge-Kutta methode voor vergelijkingen van het type (1.1) wordt gedefinieerd door m-1 k(i) Yn+1 = y + I a. ' n i=O i n n = 0, 1 ,2, ... , k(O) = h f(x ,Y ) , n n n n ( 2. 1) i-1 ( ) = h f(x + M.h, y + ' L .. k j ) n n in n l J.J n ' j=O ]. = 1 , .•. ,m-1 , h = xn+1 - xn, n i-1 M. = I L .•• ]. j=O J.J e We noemen methode (2.1) p orde consistent indien geldt (2.2) n = 0, 1 ,2, •.• , waarin y (x) de locaal analytische oplossing is dat wil zeggen de oplos-n sing van het beginwaarde probleem y' = f(x,y), n n = 0, 1 ,2,. . . . Door ontwikkeling van het linkerlid van (2.2) in een Taylor-reeks in h n vinden we voor 5e orde consistentie de volgende betrekkingen tussen de parameters a., M. en L.. (i = O, •.• ,m-1; j = 0, •.. ,i-1) ]. ]. J.J (2.3a) (2.4a) (2.5a) (2.6a) (2.8a) (2. 1 Oa) (2.11a) (2. 12a) ( 2. 13a) ( 2. 14a) ( 2. 15a) m-1 I a. ]. i=O m-1 I i=1 a.M. ]. ]. m-1 I a.if. . 1 ]. ]. i= m-1 I i=2 m-1 I i=1 m-1 I i=2 3 a.M. ]. ]. i-1 2 a. l Li.lMl 1 1=1 m-1 i-1 1-1 l a. l Lil l Ll M i=3 1 1=2 p=1 pp m-1 I i=1 4 a.M. ]. ]. 3 = 1 ; 1 =-12 4 (2. 16a) m-1 i-1 1-1 L M2 I a. I Lil I i=3 ]. 1=2 p=l lp p (2.17a) m-1 i-1 1-1 I a. I LilMl I L1 M i=3 ]. 1=2 p=1 pp m-1 i-1 1-1 (2. 18a) I a.M. I Lil I L1 M i=3 ]. ]. 1=2 p=1 pp (2.19a) m-1 i-1 1-1 p-1 l ai l Lil l L1 2 i=4 1=3 p=2 p q=l 1 = 60 ; 1 = 40 ; 1 =- ; 30 1 L M = -- ; pq q 120 Analoog aan Zonneveld [5] proberen we, simu.ltaan met het oplossen van deze consistentie-voorwaarden, te voldoen aan de relatie (2.20) m-1 ( ) h5 5 l b.k i = 2f L y (x) i=0 in 5, dx5 n . . . ukk" . k(i) . . Deze lineaire uitdr ing in de , i = 0, ••. ,m-1, kan dienen als maat n voor de discrepantie van de numerieke oplossing. Relatie (2.20) leidt tot het volgende additionele stelsel vergelij-kingen m-1 (2.3b) I b. i=0 ]. = 0 m-1 (2.4b) I b.M. i=l ]. ]. = 0 ( 2. 5b) m-1 2 I b.M. i=1 ]. ]. = 0 m-1 i-1 (2.6b) I b. I LilMl i=2 ]. 1=1 = 0 (2.7b) m-1 3 I b.M. i=1 ]. ]. = 0 (2.8b) (2.9b) (2. 1 Ob) (2.11b) (2. 12b) (2.13b) ( 2. 14b) 2. 15b) ( 2. 16b) ( 2. 17b) ( 2. 18b) (2. 19b) m-1 i-1 2 b. 2 L-1Mi . 1 1 1 1 1=2 = m~1 I i=1 4 b.M. 1 1 m-1 i-1 2 2 b.M. l L.1Mi i=2 1 1 1=1 1 L M lp p = 0 = 0 ; = 0 ; 1 = -5 1 = - ; 20 1 =-10 1 = - ; 20 5 6 3. De twee-parameter familie van oplossingen van Zonneveld Voor het geval m = 6 heeft Zonneveld [5] de stelsels (2.3a) - ( 2. 19a), (2.3b) - (2.19b) vereenvoudigd tot een stelsel van 33 vergelijkingen met 35 onbekenden, te weten 6 6 ( 3. 1 ) I a. = 1 ; ( 3. 2) I b. = 0 ; i=O 1 i=O 1 6 1 6 ( 3. 3) l a.M. = - (3.4) I b.M. = 0 i=1 1 1 2 i=1 1 1 6 2 1 6 2 (3.5) I a.M. = - ; (3.6) I b.M. = 0 ; i=1 1 1 3 i=1 1 1 6 3 1 6 3 ( 3. 7) J, a.M. =4 ( 3 .8) I b.M. = 0 1 1 i=1 1 1 6 4 1 6 4 1 (3.9) I a.M. = - (3.10) I b.M. i=1 1 1 5 i=1 1 1 5 6 6 (3.11) I a.L. 1 = 0 ; (3. 12) l b.L.1 = 0 ; 1 1 . 2 1 1 i=2 1= 6 6 (3.13) I a.L. 2 = a (1-M) ; (3.14) l b.L.2 = -b2(1-M2)/4 ; i=3 1 1 2 2 i=3 1 1 6 6 (3.15) I a.L. 3 = a (1-M) ; (3.16) l b.L.3 = -b (1-M )/4 i=4 1 1 3 3 i=4 1 1 3 3 6 6 (3.17) I a.L. 4 = a4(1-M4 ) . (3.18) I b.L. 4 = -b4(1-M4 )/4 1 1 , 1 l i=5 i=5 6 6 (3.19) I a.L. 1M. = 0 ; (3.20) l b.L. ,M. = 0 i=2 1 1 1 i=2 1 1 1 1 = M2/2 1 L2.~ = M3/3 (3.21) I L2 .M. (3.22) I ; 1 1 2 1 1 2 i=1 i=1 7 2 = M2/2 2 L3.~ = ~/3 (3.23) I L3 .M. (3.24) I ; i=1 l l 3 i=1 l l 3 3 2 3 - 3 (3.25) I L4.M. = M4/2 ; (3.26) l L4-~ - M4/3 i=1 l l . 1 l l 1= (3.27) a1 = 0 ; (3.28) b1 = 0 (3.29) = 1 (3.30) = 1 (3.31) ( 3. 32) = 0 (3.33) = 0 Iedere oplossing van dit stelsel geeft volgens (3.33) een 6-punts in-tegratie formule terwijl met een extra punt een schatting gemaakt kan worden van de discrepantie. De oplossing die door Zonneveld wordt aangegeven wordt bepaald door 1 4 M3 = 2 en M4 = 5 . De waarden van de parameters zijn dan (3.34) a= --12.. a -a -a -o· a - 162. a - 125. a5 = 33146,· bo= -¥-r.-11,. b1=b5=0,· o 336' 1- 3 6- ' 2 - 336 ' 4 - 336' l'f -162 224 -125 42 2 1 1 b2= 14; b3= ~; b4= 14; b6= "'f4; L10= 9; L20= 12; L21= 4; L = l. L =0 · L = l. L = --21. L - .=.Ll.2.. L - 126 · L - 2£._, 30 8' 31 ' 32 8' 40 125' 41- 125' 42- 125' 43- 125' -63 189 -36 -112 50 133 L50- 28; L51- 28; L52- 28' L53- 28' L54- 28' L60- 168' _ -378. _ 276. _ 112. _ 25. 2 1 L61- 168' L62- 168' L63- 168' L64- 168' L65=0; M,= 9; M2= 3; 1 4 M3= 2; M4= 5; M5=M6=1· 8 4. Het stabiliteitspolynoom van een RK(5,6) methode Beschouw het volgende lineaire beginwaarde-probleem y' = oy, ( 4. 1 ) Passen we hierop een RK(5,6) schema toe dan vinden we met 5 zi 6 ( 4 .2) P6(z) = I -=t + (3z , i=O 1.. waarin 5 ( 4. 3) (3 = a II L. . 1. 5 i=1 1 1-We noemen P6(z) het stabiliteitspolynoom. We definieren het stabiliteitsgebied S van P6(z) als volgt Verder definieren we de reele stabiliteitsgrens a van P6(z) door middel van a= inf{x I x > O, IP6(-x)I > 1}. Het is duidelijk data een functie is van (3. Van de ons bekende RK(5,6) formules had de methode van Lawson [4] de grootste reele stabiliteitsgrens, namelijk a = 5,6 L Met behulp van een methode, gegeven in van der Houwen [1] berekenden we ec°'ter dat (4.4) en a. max = max a.(S) = a 6.26, In paragraaf 5 geven we een oplossing van de vergelijkingen (3.1) -(3.33) en (4.3), waarbij we in de laatste vergelijking S = S t stellen. op Het is duidelijk dat deze oplossing correspondeert met een RK(5,6) schema waarvan de reele stabiliteitsgrens gelijk is aan a. • max 9 10 5. Een RK(5 26) methode met maximale reele stabiliteitsgrens Zonneveld [5] heeft het stelsel vergelijkingen (3.1) - (3.33) opge-lost, met M3 en M4 als vrije parameters. Daarbij geeft hij de uitdrukkingen (5.4) Uit (3.17) en (3.33) volgt en uit (3.1), (3.3), (3.5), (3,7), (3.9), (3.27) en (3.33) vinden we Substitutie van (5.2), ... ,(5.7) in (4.3) gee~t de vergelijking Uit (5.1) en (5.8) leiden we een vierkantsvergelijking in M2 af, De grootste wortel van M2 leidt tot minimale waarden van de overige para-meters. De nog overgebleven vrijheidsgraad hebben we benut door a3 = 0 te kiezen. Hierdoor wordt M4 eenduidig bepaald: 11 (5.10) De waarden van de parameters a., b. en M. geven we in tabel 5.1, terwijl 1 1 1 in tabel 5.2 de waarden van L .. opgenomen zijn. 1J Tabel 5.1 De waarden van a., b. en M .• 1 1 1 1 a. b. 1 1 0 .1013 8388 8447 4274 +.1687 9574 4506 3191 10+1 M. 1 0 1 0 0 .2166 3751 5122 2449 2 .4710 9636 5451 7556 3 0 4 ,3760 3358 8853 7316 5 .5148 6157 2470 8536 10-1 6 0 -.1422 2671 9546 3261 +.1797 2905 6842 5553 -.8017 3755 1797 5110 0 +.2579 1843 4328 8994 Tabel 5 .2 De waarden van L ... 1J 10+2 .3249 5627 2683 3674 10+2 .4641 0728 0027 7517 10+1 .7856 4291 2011 1007 1 10+1 1 L10 = +.2166 3751 5122 2449. ; 120 = +.8123 9068 1708 4184 10-1 L21 = +.2437 1720 4512 5255 130 = +.1088 9359 0760 4054 ; L31 = +,7137 3905 6569 5119 10-1 ; 132 = +.2838 3978 3610 3951 L40 = +.4136 4798 7348 0195 ; 141 = -.9615 3115 2649 3416 ; L42 = +,7328 5885 8261 3591 143 = +.6006 6721 9051 0636 ; L50 = -,1795 2996 1930 4468 10+1 L51 = +.4792 6226 0139 7445 10+1 ; 152 = +.8241 2636 9753 6218 153 = -.4387 0248 2693 7314 10+1 154 = +.1565 5754 7509 0715 10+1 L60 = +.8113 7444 5235 0849 161 = -.2142 3218 4008 3255 10+1 L62 = +.1230 7807 2667 0698 10+1 ; L63 = +,9335 8481 1221 5743 ; L64 = +.1665 8185 6955 8982 ; 165 = 0 & 12 6. Stabiliteitsgebieden bij RK(5,6) methoden In deze paragraaf geven we het stabiliteitsgebied in een omgeving van de oorsprong van enige RK(5,6) methoden. In verband met de symmetrie van het stabiliteitsgebied is slechts het gedeelte met positieve y-coordinaat getekend. In figuur 1 geven we het stabiliteitsgebied behorende bij de methode van Zonneveld; het stabiliteitspolynoom hiervan is P ( ) 1 + l z2 + l z3 1 4 + 1 5 1 6 6 z = + z 2 6 + 2.4 z 120 z + 1 440 z ' Het stabiliteitsgebied van de methode van Lawson (hierbij is SL= 1/1280) is in figuur 2 getekend, terwijl figuur 3 het stabiliteitsgebied van de in dit rapport beschreven methode laat zien (S t = ,725591 -3). op 0 Figuur 1 Stabiliteitsgebied van Zonneveld's methode -1 -6 -s Figuur 2 Stabiliteitsgebied van Lawson's methode -7 -6 -s -a Figuur 3 Stabiliteits ebied van de nieuwe methode -a -2 -1 1_3 +2 +1 Re(z) +O +1 +2 +I Re(z) +O +1 +2 +1 Re(z) +O +1 14 7. Testvoorbeeld Als testvoorbeeld kozen we de stijve differentiaal-vergelijking (zie Lapidus en Seinfeld [3]) y'(x) = -200(y(x)-F(x)) + F'(x), F(x) = 10 - (10+x)e-x, y(O) = 10. De exacte oplossing van (7.1) luidt y(x) = F(x) + 10e-2oox. De jacobiaan van deze differentiaal-vergelijking is de scalar -200. Dit houdt in dat RK(5,6) methoden stabiel zijn bij integratie met staplengten h, Cl. h<h =--- max 200 De RK(5,6) methode van Zonneveld heeft een maximale integratiestap h = .022; terwijl de in dit rapport voorgestelde RK(5,6) methode stappen max toelaat, die 40% grater zijn. Dit blijkt duidelijk uit tabel 7.1, waarin we de relatieve fout in x = .4 en x = 10 geven bij integratie met constante stap. 15 Tabel J. 1 Relatieve fout in x = .4 en x = 10 bij integratie van (7.1) stap Zonneveld's metbode Nieuwe methode X = .4 X = 10 X = .4 X = 10 .020 1.9 10-3 7.7 10-10 4.5 10-5 5.7 10-10 .021 1.0 10-1 4.7 10-10 1.4 10-3 3.3 10-10 .022 1 .6 5.8 10-1 1.6 10-2 1.3 10-10 .023 7.3 10+29 5,3 10-2 1. 1 10-10 .024 1.8 10+3 10+52 2.8 10-1 2.2 10-10 .025 4.4 10+2 10+69 3 .1 2.0 10-10 .026 4.7 10+3 10+78 3 .1 10-1 3.2 10-4 .027 9.9 10+2 10+83 6.6 10-1 3.0 10-10 .028 1. 8 10+3 10+81 7.7 10-3 1.8 10-9 .029 2.9 10+1 10+71 1.0 10-4 1.7 10-9 .030 8.2 10+48 1.3 10-5 5,5 1C-9 .031 6.o 10-2 2.4 10-8 .032 312 10+70 16 Referenties [1] Houwen, P.J. van der, One step methods for linear initial value problems I, Polynomial methods, TW report 119, Mathematisch Centrum, .Amsterdam, (1970). [2] Houwen, P.J. van der, Stabilized Runge-Kutta methods with limited storage requirements, TW report 124, Mathematisch Centrum, .Amsterdam, (1971). [3] Lapidus L., J.H. Seinfeld, Numerical solution of ordinary differential equations, Acadamic Press, New York, (1971). [4] Lawson, J.D., An order five Runge-Kutta process with extended region of stability, SIAM J. Numer. Anal. 3, 593, (1966). [5] Zonneveld, J.A., Automatic numerical integration, MC Tract 8, Mathematisch Centrum, .Amsterdam, (1964). 

Een 5e orde 6-punts Runge-Kutta formule met optimale stabiliteitsgrens

Abstract

Similar works

Full text

Available Versions

NARCIS

CWI's Institutional Repository