XProfan

WindowTitle "Anfangswertproblem einer Differentalgleichung 2. Ordnung lösen nach Størmer-Verlet"
' ---------------------------------------------------------------------------------------------
' Aus Fortran90 nach XProfan-11 übersetzt 2014-08 von P. Specht, Wien (Österreich)
' Keine wie auch immer geartete Gewähr! Verwendung auf alleinige Gefahr de(s|r) Anwender(s|in)!
' Quelle: https://jean-pierre.moreau.pagesperso-orange.fr/Cplus/stormer.txt
'**********************************************************************************************
'*  Differential equation y"=f(x,y,y') by  Størmer's method  *
' (erweiterbar auf allg. Verlet type, z.B. y'''[...]=f(x,y,y',y''[,y''',...])
'* --------------------------------------------------------- *
'* SAMPLE RUN (BEISPIEL)                                     *
'* Löse  y" = 8*y*y / (1+2*x) from x=0 to x=1                *

proc G :parameters x!,y!,z!'= x, y, y'

    return 8*sqr(y!)/(1+2*x!)

endproc

proc F :parameters x!,y!,z!'= x, y, y'

    return z!

endproc

'* with initial conditions: x(0)=0, y(0)=1 and y'(0)=-2      *
Declare x![4],y![4],z![4],h!' Initial conditions:
x![1]= 0
y![1]= 1
z![1]=-2' dy/dt
h!=0.01' Step
'*  Compare to exact solution   y = 1/(1+2*x)                *

Proc Fx :parameters x!' Exact Solution:

    return 1/(1+2*x!)

endproc

'* --------------------------------------------------------- *
'*   Differential equation y"=f(x,y,y') by Stormer's method  *
'* --------------------------------------------------------- *
'*        X           Y           Y exact         Error      *
'*     0.000       1.000000       1.000000     0.0000000000  *
'*     0.010       0.980392       0.980392     0.0000000001  *
'*     0.020       0.961538       0.961538     0.0000000295  *
'*     0.030       0.943396       0.943396     0.0000000457  *
'*     0.040       0.925926       0.925926     0.0000000974  *
'*     0.050       0.909091       0.909091     0.0000001285  *
'*     ...         ...            ...          ...           *
'*     0.950       0.344866       0.344828     0.0000381695  *
'*     0.960       0.342505       0.342466     0.0000388874  *
'*     0.970       0.340176       0.340136     0.0000396196  *
'*     0.980       0.337878       0.337838     0.0000403406  *
'*     0.990       0.335612       0.335570     0.0000410721  *
'*     1.000       0.333375       0.333333     0.0000418231  *
'*                                                           *
'*   F90 Release By J-P Moreau, Paris (www.jpmoreau.fr).     *
'*************************************************************
WindowStyle 24:font 2:set("decimals",16):Window 0,0-%maxx,%maxy-40
Declare c![4],a1!,a2!,a3!,yex!,k&,er! , x1!,y1!,z1!
a1!= 1.083333333333333333 : a2!=-2*(a1!-1) : a3!= a1!-1
clearclip
print "\n------------------------------------------------------------------------------"
print "       Differentialgleichung y''=f(x,y,y') mit Störmer's Methode lösen "
print "------------------------------------------------------------------------------\n"
print "         X                   Y                  Y exact                Error\n"
putclip " X | Y | Y_exact | Error \n"
yex!=Fx(x![1]):print tab(2);x![1],tab(22);y![1],tab(42);yex!,tab(62);er!
putclip str$(x![1])+"|"+str$(y![1])+"|"+str$(yex!)+"|"+str$(er!)+"\n"
' Runge-Kutta for first 2 steps:

whileloop 1,2:k&=&Loop

    RK4D2 x![k&],y![k&],z![k&],h!
    x![k&+1]=x1!
    y![k&+1]=y1!
    z![k&+1]=z1!
    yex!=Fx(x![k&+1])
    er!=abs(yex!-y![k&+1])
    print tab(2);x![k&+1],tab(22);y![k&+1],tab(42);yex!,tab(62);er!
    putclip str$(x![k&+1])+"|"+str$(y![k&+1])+"|"+str$(yex!)+"|"+str$(er!)+"\n"

endwhile

REPEAT'Main loop G10:

    whileloop 2,4:k&=&Loop

        c![k&]=G( x![5-k&],y![5-k&],z![5-k&] )

    endwhile

    y![4]=2*y![3]-y![2]+sqr(h!)*(a1!*c![2]+a2!*c![3]+a3!*c![4])
    x![4]=x![3]+h!
    yex!=Fx(x![4])
    er!=abs(yex!-y![4])
    print tab(2);x![4],tab(22);y![4],tab(42);yex!,tab(62);er!
    putclip str$(x![4])+"|"+str$(y![4])+"|"+str$(yex!)+"|"+str$(er!)+"\n"

    Whileloop 1,3:k&=&Loop'für nächsten Schritt umkopieren:

        x![k&]=x![k&+1]
        y![k&]=y![k&+1]

    endwhile

    if %csrlin>40:waitinput 4000:cls

        print "\n         X                   Y                  Y-exakt                Error\n"

    endif

UNTIL x![3]>1'end x value = 1

print "\n EOF"
Print " Resultate in Zwischenablage!"
Waitinput
END

PROC RK4D2 :parameters x!,y!,z!,h!',x1!,y1!,z1!

    declare c1!,d1!,c2!,d2!,c3!,d3!,c4!,d4!
    c1!=F(x!,y!,z!)
    d1!=G(x!,y!,z!)
    c2!=F(x!+h!/2, y!+h!/2*c1!, z!+h!/2*d1!)
    d2!=G(x!+h!/2, y!+h!/2*c1!, z!+h!/2*d1!)
    c3!=F(x!+h!/2, y!+h!/2*c2!, z!+h!/2*d2!)
    d3!=G(x!+h!/2, y!+h!/2*c2!, z!+h!/2*d2!)
    c4!=F(x!+h!,   y!+h!*c3!,   z!+h!*d3!  )
    d4!=G(x!+h!,   y!+h!*c3!,   z!+h!*d3!  )
    x1!=x!+h!
    y1!=y!+h!*(c1!+2*c2!+2*c3!+c4!)/6
    z1!=z!+h!*(d1!+2*d2!+2*d3!+d4!)/6

endproc

PROGEND
'{=============================================================================================
'    Stormer, aka Newton-Størmer-Verlet leapfrog method
' =============================================================================================
' Namensgebung: Der Algorithmus selbst geht zurück auf Delambre (1791), oftmals wiederentdeckt,
' u.a. 1907 durch Carl Størmer, 1909 durch Cowell und Cromlin, 1967 populär durch Loup Verlet.
' Verbessert die explizite Euler-Methode ohne allzu hohen Mehraufwand. Oft in Physik-Engines
' zu finden. Details siehe https://en.wikipedia.org/wiki/Verlet_integration
' =============================================================================================
' Methode STORMER (aka Newton-Störmer-Verlet-leapfrog method):
' löst Y"=f(x,y,y') mit gegebenen Anfangsbedingungen
' ----------------------------------------------------------------------------------
' Eine Differentialgleichung 2. Ordnung (Y") kann ersetzt werden durch ein
' Differentialgleichungs-SYSTEM bestehend aus 2 Gleichungen 1. Ordnung!
' Wenn also y"=f(x,y,y') gegeben ist mit y(a=Anfang) und y'(a=Anfang), dann kann man
' diese Formel durch Substitution u = y' umformen zu:
'
'                    | u'=f(x,y,u)
'                    | y'=u
'                    | mit den neuen Anfangsbedingungen y(a), u(a)
'
' Beginnen wir die Betrachtung mit jener speziellen Form der Taylor-Entwicklung,
' bei der der Rest als Integral dargestellt wird (Wiki: Taylorentwicklung)
'                                       x+h
'           y(x+h) = y(x) + h * y'(x) + INTEGRAL (x+h-t)*y"(t).dt
'                                        x
' Durch Verbringung der Nichtintegral-Terme auf die linke Seite und Aufspaltung
' des 'Restfehler-Integrals' in zwei Teile kann diese Form auch folgender-
' maßen dargestellt werden:
'                                        x+h              x-h
'           y(x+h) - 2y(x) + y(x-h) = INTEGRAL ...  +  INTEGRAL ...
'                                         x                x
' Für das zweite Integral wird die Variable x nun substituiert durch:  u = 2*x-t
' Dann gilt auf Grund der Kettenregel:
'           x-h                              x+h
'         INTEGRAL (x-h-t) * y"(t).dt = - INTEGRAL (x+h-u)*y"(u).du
'            x                                x
' und die Ableitung .du rück-substituiert (zu .-dt) liefert dann:
'                                    x+h
'           du=-dt ==>           = INTEGRAL (x+h-t)*y"(2x-t).dt
'                                     x
' Somit kann man schreiben:              x+h
'           y(x+h) - 2*y(x) + y(x-h) = INTEGRAL (x+h-t)*[y"(t)+y"(2x-t)].dt
'                                         x
' Nun verwenden wir unter Hinweis, daß ja y"(x)=f(x,y,y') gilt, folgenden Ausdruck
' als Interpolationspolynom:
'                             x_n+1
'     y_n+1 - y_n + y_n-1 = INTEGRAL (x_n+1 - t)*[P(t)+P(2*x_n - 1)].dt  = ...
'                              x_n
'               x_n+1
'    ...= h^2*INTEGRAL (x_n+1 -t)*[O_00*Div(f_n)+O_11*Div(f_n)+O_22*Div(f_n)].dt ,
'                x_n
'                         1           |(-s) (m)|
' wobei O_mm  = -1^m * INTEGRAL (1-s)*|    +   | .ds
'                         0           |( m) (s)|
'                                               (vgl. auch 'Adams-Bashford-Methode')
' Letzteres bringt uns schließlich zu STORMER's Formeln:
' ==================================================================================
' Explizit formuliert:  y_n+1 -2*y_n + y_n-1 = h^2/12 * [13 * f_n - 2*f_n-1 + f_n-2]
' Das interessierende   y_n+1 = h^2/12 * [13 * f_n - 2*f_n-1 + f_n-2]+ 2*y_n - y_n-1
' (f_n+1 muss nicht bekannt sein: Prognoserechnung bei aktuell ergänzten Zeitreihen)
' ----------------------------------------------------------------------------------
' oder implizit:        y_n+1 -2*y_n + y_n-1 = h^2/12 * [f_n+1 + 10 * f_n + f_n-1 ]
'}==================================================================================