{VERSION 3 0 "IBM RISC UNIX" "3.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 } {PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 8 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Title" 0 18 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 }{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 8 8 0 0 0 0 0 0 -1 0 }} {SECT 0 {PARA 18 "" 0 "" {TEXT 256 30 "Ein eingebettetes RK-Verfahren " }{TEXT -1 0 "" }}{PARA 19 "" 0 "" {TEXT -1 11 "Simon Giehl" }}{PARA 0 "" 0 "" {TEXT -1 44 "Seminar: Software fuer Diffentialgleichungen" } }{PARA 0 "" 0 "" {TEXT -1 29 "Thema: Schrittweitensteuerung" }}{PARA 0 "" 0 "" {TEXT -1 152 "Worksheet: Eine Implementierung eines eingebet tetes RK-Verfahren mit DP5(4)-Koeffinzientensatz und Anwendung auf das restringierte Drei-Koerper-Problem." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 43 " Der Dormand-Prince-5(4)-Koeffizientensatz." }}{PARA 0 "" 0 "" {TEXT -1 60 "* eingebettetes RK-Schema\n* Ordnung:\n p1 = 4\n p2 = 5 " }}{PARA 0 "" 0 "" {TEXT -1 24 "* Stufenzahl\n s1 = 6" }}{PARA 0 "" 0 "" {TEXT -1 306 " s2 = 7\n* Verwendung des Fehlber-Tricks. So mit hat das Verfahren nur effektiv\n s = 6\n Stufen.\n* Zusaetz lich zu den Konvergenzbedingungen wurde der Koeffizientensatz nach Kri terien ausgewaehlt, die die Schrittweitensteuerung betreffen. So wurde auf die Qualitaet des Fehlerschaetzers Wert gelegt." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 443 "A := matrix( [[0,0,0,0,0,0,0], [1/5,0,0, 0,0,0,0], [3/40,9/40,0,0,0,0,0], [44/45,-56/15,32/9,0,0,0,0], [19372/6 561,-25360/2187,64448/6561,-212/729,0,0,0], [9017/3168,-355/33,46732/5 247,49/176,-5103/18656,0,0], [35/384,0,500/1113,125/192,-2187/6784,11/ 84,0]] ):\nc := vector( [0,1/5,3/10,4/5,8/9,1,1] ):\nb1 := vector( [35 /384,0,500/1113,125/192,-2187/6784,11/84,0] ):\nb2 := vector( [5179/57 600,0,7571/16695,393/640,-92097/339200,187/2100,1/40] ):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "print(c);\nprint(A);\nprint(b1);\np rint(b2);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 17 " Implementierung." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 16 "GRUNDALGORIHMUS." }} {PARA 0 "" 0 "" {TEXT -1 54 "In jedem Schritt werden folgende Schritte durchlaufen:" }}{PARA 0 "" 0 "" {TEXT -1 39 "1) Berechnung des naec hsten Schrittes" }}{PARA 0 "" 0 "" {TEXT -1 21 "2) Fehlerschaetzung " }}{PARA 0 "" 0 "" {TEXT -1 186 "3) Schrittweitenvorschlag\n4a) Ann ahme des Schrittes: Naechster Schritt mit vorgeschlagener Schrittweite \n4b) Ablehnung des Schrittes: Wiederhole Schritt mit vorgeschlagener \+ Schrittweite" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 10 "NAME: DP54" }}{PARA 0 "" 0 "" {TEXT -1 110 "BESCHREIBUNG: Einge bettetes RK-Verfahren nach DP5(4)-Koeffizientensatz fuer mehrdimension ale GDGen 1. Ordnung." }}{PARA 0 "" 0 "" {TEXT -1 8 "EINGABE:" }} {PARA 0 "" 0 "" {TEXT -1 54 " F - Prozedur zur Auswertung der \+ rechten Seite" }}{PARA 0 "" 0 "" {TEXT -1 246 " t0 - Startzeitp unkt\n tend - Endzeitpunkt\n X0 - Startvektor\n hstart - Sta rtschrittweite\n TOL - Tolerance: erwuenschte lokale Genauigkeit\n \+ AUSGABE - der name der table-Datenstruktur, in der die Ausgabe gespe ichert werden soll" }}{PARA 0 "" 0 "" {TEXT -1 253 "AUSGABE: Eine tabl e-Datenstruktur, die folgende Daten enthaelt:\n* zu jedem Schritt den \+ Zeitpunkt t und den Vektor X \n* die Anzahl der akzeptierten und die A nzahl der zurueckgewiesenen Schritte\n* der Aufwand, gemessen durch di e Anzahl der f-Auswertungen\n" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1645 "DP54 :=\nproc( F::procedure, t0::numeric, tend::numeric, X0::v ector, hstart::numeric, TOL::numeric, AUSGABE )\nlocal t, T, h, hnext, normEST, \n rho, hMAX, hMIN, q, beta,\n ENDE, ACCEPT, REJEC T, d, i, \n X, EST, TEMP, k1, k2, k3, k4, k5, k6, k7;\nwith(linal g,norm);\n\nrho := 0.9;\nbeta := 1./(4+1);\n\nd := op(2,op(2,eval(X0)) );\n\nk1 := vector(d);\nk2 := vector(d);\nk3 := vector(d);\nk4 := vect or(d);\nk5 := vector(d);\nk6 := vector(d);\nk7 := vector(d);\nEST := v ector(d);\nTEMP := vector(d);\n\nX := map(evalf,X0);\nt := evalf(t0); \nT := evalf(tend);\nh := evalf(hstart);\n\nACCEPT := 1; REJECT := 0; \nAUSGABE[ACCEPT] := [t,seq(X[i],i=1..d)];\nF(t,X,k1);\n\nwhile t " 0 "" {MPLTEXT 1 0 505 "GetSteps := proc( Ausgabe::table ) \n local Steps, N, i, t, h;\n N := Ausgabe[Acc];\n Steps := NULL;\n for i from 1 to N-1 do\n t := Ausgabe[i][1];\n h := Ausgabe[i+ 1][1]-t;\n Steps := Steps, [t,h];\n od;\n RETURN( [Steps] );\nend : GetGraph := proc( Ausgabe::table, liste::list )\n local Graph, N, i , j, v;\n Graph := NULL;\n N := Ausgabe[Acc];\n for i from 1 to N d o\n v := NULL;\n for j in liste do\n v := v, Ausgabe[i][j+1 ];\n od;\n Graph := Graph, [v];\n od;\n RETURN( [Graph] );\nen d: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 75 " Anwendung auf periodische Orbits des restringi erten Drei-Koerper-Problems." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 95 "Die Prozedur F implementiert die rechte Seite d er DGL des restringierten Drei-Koerper-Problems." }}{PARA 0 "" 0 "" {TEXT -1 64 "Sie wird dem implementierten Verfahren als Parameter uebe rgeben." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 375 "F := proc( t::nu meric, X::vector, K::vector )\n local mu1, mu2, D1, D2;\n mu1 : = 0.012277471;\n mu2 := 1-mu1;\n D1 := ((X[1]+mu1)^2 + X[3]^2)^( 3/2);\n D2 := ((X[1]-mu2)^2 + X[3]^2)^(3/2);\n K[1] := X[2];\n \+ K[2] := X[1] + 2*X[4] - mu2*(X[1]+mu1)/D1 - mu1*(X[1]-mu2)/D2;\n \+ K[3] := X[4];\n K[4] := X[3] - 2*X[2] - mu2*X[3]/D1 - mu1*X[3]/D2; \n RETURN();\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 18 "Es kann \+ los gehen!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 99 "X0 := vector( 4,[0.994,0,0,-2.00158510637908252240537862224]);\nT := 17.065216560157 9625588917206249;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "DP54(F ,0,T,X0,0.001,0.00001, G);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 10 "Sta tistik:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "G[Acc];\nG[Rej]; \nG[Aufwand];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 30 "Die Regelung der Schrittweite:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "plot( Get Steps( G ) );" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 83 "Graphische Darst ellung der numerischen Loesung: Extraktion der x- und y-Koordinaten" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "plot( GetGraph( G, [1,3] ) , title=\"Arenstorf-Orbit\");" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 32 " Ein weiterer periodischer Orbit:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "X0 := vector(4,[1.2,0,0,-1.04936]);\nT := 6.19217;" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "DP54(F,0,T,X0,0.001,0.0001 , H);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "H[Acc];\nH[Rej];\n H[Aufwand];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "plot( GetSte ps( H ));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "plot( GetGraph ( H, [1,3] ));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}} {MARK "10" 0 }{VIEWOPTS 1 1 0 1 1 1803 }