Abt. Runge-Kutta RK4 für Nichtlineare Differentalgleichungssysteme
============================================
Wie schon erwähnt, läßt sich das Verfahren von Runge und Kutta auch auf DGS höherer Ordnung anwenden. Allerdings nimmt die Präzision (trotz des im linearen Bereich exzellenten Verhaltens, immerhin Fehlerordnung 5!) mit der Anzahl R der Dimensionen rapide ab - genauer gesagt mit der R-ten Wurzel aus 'Schrittweite^5'. Man sollte diesfalls also nicht allzu hohe Erwartungen hegen .
Gruss
P.S.: Zur Erhaltung maßvoller technisch Genauigkeitsanforderungen wäre bereits erheblicher :such03: Zusatzaufwand erforderlich. Besser rückt man solchen Problemen dann mit sog. 'Monte Carlo-Methoden' zu Leibe. Literatur dazu: P. Zinterhof: Integrale und mehrdimensionale Funktionen, Aufsatz in OCG-Schriftenreihe Band 12 'Zahlentheoretische Methoden in der Numerischen Mathematik', S.132ff
WindowTitle "Runge-Kutta RK4 für Systeme von Differentialgleichungen höherer Ordnung"
' BASIC-Original: B.Brand: Algorithmen zur praktischen Mathematik, S.240, Oldenbourg Vlg.
' Migriert nach XProfan-11 in 2014-09 by P.Specht, Wien.
' Keine wie auch immer geartete Gewähr!
WindowStyle 24:Window 0,0-%maxx,%maxy
declare k1!,k2!,k3!,k4!,l1!,l2!,l3!,l4!
':: Beispiel: 'Nichtlinares Differentialgleichungssystem'
':: Zwei Testfunktionen, die eine gekrümmte Fläche bilden (X: Abszisse, Y: Ordinate, Z: Höhe)
declare x!,y!,z!,f!,g!
declare i&,h!,x![],y![],z![],xend!
proc F :parameters x!,y!,z!
return (y!-z!)/x!
endproc
proc G :parameters x!,y!,z!
return (y!+z!)/x!
endproc
' Hinweis: Zeilen, die mit :: beginnen, stammen vom Übersetzter nach XProfan-11
print "\n Im Beispiel wird ein System aus '2 nichtlinearen Differentialgleichungen "
print " in 3 Variablen' gelöst. Die Gleichungen lauten F=(y+z)/x und G=(y-z)/x. "
print " Die Anfangswerte für x, y und z sowie die Schrittweite h sind vorzugeben. "
print " Ferner wird ein Abbruchkriterium, z.B. der Maximale x-Endwert eingegeben. \n"
print " Aus Gründen meiner Bequemlichkeit werden bei Eingabe 'x=0' die Beispiels- "
print " werte der Originalvorlage genommen (Ergebnisse einsehbar im Programmtext).\n"
' Die Input-Werte des Original-Beispiels:
' x!=1 : y!=1 : z!=0 : h!=0.2 : xend!=3
' führen zu folgender Ausgabe:
'Schritt X Y Z
' 0 1 1 0
' 1 1.2 1.18011937557392 0.217584940312213
' 2 1.4 1.32150759540037 0.462240502741668
' 3 1.6 1.42651937035306 0.724648391007437
' 4 1.8 1.49791768072086 0.998168848790414
' 5 2 1.53848713605723 1.27796059664382
' 6 2.2 1.55087497139788 1.56043022121442
' 7 2.4 1.5375315540692 1.84287395296026
' 8 2.6 1.50069575905948 2.12323677414932
' 9 2.8 1.44240050556552 2.39994581694282
' 10 3 1.36448684585809 2.67179234356735
' ---
':: ACHTUNG! Auch das Endkriterium ist jeweils an die Gleichungsstruktur anzupassen,
':: sonst rechnet das Ding, bis der Speicher voll ist!!!
print " X = ";:input x! : if x!=0:x!=1:y!=1:z!=0:h!=0.2:xend!=3.0:goto "skip":endif
print " Y = ";:input y!
print " Z = ";:input z!
print " Schrittweite h = ";:input h!
print " Endwert (hier: X_end) = ";:input xend!
skip:
set("decimals",17)
clear x![],y![],z![]:clearclip
x![0]=x!:y![0]=y!:z![0]=z!
ausgabe:
print "\n Schritt";Tab(25);" X";Tab(50);" Y";Tab(75);" Z":print mkstr$("-",90)
Repeat
print " ";i&,Tab(25);format$("%g",X![i&]);Tab(50);format$("%g",Y![i&]);Tab(75);format$("%g",Z![i&])
::putclip str$(i&)+" "+format$("%g",X![i&])+" "+format$("%g",Y![i&])+" "+format$("%g",Z![i&])+"\n"
k1!=h!*F(x![i&],y![i&],z![i&])
l1!=h!*G(x![i&],y![i&],z![i&])
k2!=h!*F(x![i&]+h!/2,y![i&]+k1!/2,z![i&]+l1!/2)
l2!=h!*G(x![i&]+h!/2,y![i&]+k1!/2,z![i&]+l1!/2)
k3!=h!*F(x![i&]+h!/2,y![i&]+k2!/2,z![i&]+l2!/2)
l3!=h!*G(x![i&]+h!/2,y![i&]+k2!/2,z![i&]+l2!/2)
k4!=h!*F(x![i&]+h!, y![i&]+k3!,z![i&]+l3!)
l4!=h!*G(x![i&]+h!, y![i&]+k3!,z![i&]+l3!)
y![i&+1]=y![i&]+1/6*(k1!+2*k2!+2*k3!+k4!)
z![i&+1]=z![i&]+1/6*(l1!+2*l2!+2*l3!+l4!)
x![i&+1]=x![i&]+h!
inc i&
::if %csrlin>55:waitinput ' 55 und 85 an aktuelle Schirmdimensionen anpassen!
::cls:print:print " Schritt";Tab(25);" X";Tab(50);" Y";Tab(75);" Z":print mkstr$("-",90)
::endif
until x![i&]>(xend!*1.000000000000001) '(wegen höherer Genauigkeit von Intel-FPUs)
::print "OK.":beep:print " Ausgabe steht auch in Zwischenablage!":waitinput
End
Alles anzeigen