XProfan

WindowTitle "Demo: Runge-Kutta RK4 für Systeme von Differentialgleichungen höherer Ordnung"
' Angelehnt an: B.Brand: Algorithmen zur praktischen Mathematik, S.240, Oldenbourg
' Testhalber migriert nach XProfan-11 in 2014-09 by P.Specht, Wien.
' Demo für privaten Gebrauch, 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