Code: Alles auswählen
; Autor Funktionen: alter Mann , adaptiert von Engeln-Müllges/Uhlig "Numerical Algorithms with C"
; Autor Beispiel : Froggerprogger
; Date: 23. 10. 2008
; OS: Windows
; Demo: Ja
;
#DBL_EPSILON =2.2204460492503131e-016 ; kleinster Wert, so dass 1.0+DBL_EPSILON != 1.0
EnableExplicit
Procedure.d MathMin (Wert1.d,Wert2.d)
If Wert1 < Wert2
ProcedureReturn Wert1
EndIf
ProcedureReturn Wert2
EndProcedure
Procedure.d MathMax (Wert1.d,Wert2.d)
If Wert1 > Wert2
ProcedureReturn Wert1
EndIf
ProcedureReturn Wert2
EndProcedure
Prototype DglSysFkt ( x.d, y.d(1), f.d(1) )
;************************************************************************
;* *
;* Loesung eines gewoehnlichen Differentialgleichungssystems 1. Ordnung *
;* -------------------------------------------------------------------- *
;* nach der Runge-Kutta-Fehlberg-Methode [O(h^6)] mit Schaetzung des *
;* ----------------------------------------------------------------- *
;* lokalen Fehlers und Schrittweitensteuerung *
;* ------------------------------------------ *
;************************************************************************
#MAX_FCT = 10000 ; maximale Zahl an Auswertungen der rechten Seite des Differentialgleichungssystems
;************************************************************************
;* Ein System gewoehnlicher Differentialgleichungen 1. Ordnung wird *
;* nach der Runge-Kutta-Fehlberg-Methode [O(h^6)] mit Schaetzung des *
;* lokalen Fehlers und Schrittweitensteuerung integriert. *
;* *
;* Eingabeparameter: *
;* ================= *
;* x Anfangspunkt des Integrationsintervalls *
;* xend rechter Rand des Integrationsintervalls (xend < x erlaubt) *
;* n Anzahl der Differentialgleichungen *
;* y Loesung des Differentialgleichungssystems bei x *
;* dgl Zeiger auf eine Funktion, die die rechte Seite des Diffen- *
;* tialgleichungssystems y' = f(x,y) auswertet *
;* h Anfangsschrittweite; falls h unrealistisch vorgegeben wird, *
;* wird h intern veraendert; h kann negativ sein, wenn *
;* xend < x. *
;* hmax obere Grenze fuer den Betrag der waehrend der Rechnung be- *
;* nutzten Schrittweiten (hmax > 0.0) *
;* epsabs \ Grenzen fuer den zulaessigen lokalen Fehler, relativ zur *
;* epsrel / aktuellen Schrittweite. Gilt fuer jede Komponente der er- *
;* rechneten Loesung y[i] *
;* |Schaetzung des lokalen Fehlers| <= *
;* |h| * (epsrel * |y[i]| + epsabs), *
;* dann wird die Loesung im momentanen Schritt akzeptiert. *
;* epsabs = 0.0 => nur Test fuer relativen Fehler *
;* epsrel = 0.0 => nur Test fuer absoluten Fehler *
;* *
;* Ausgabeparameter: *
;* ================= *
;* x letzte Stelle, an der eine Loesung erfolgreich berechnet *
;* wurde. Normalerweise ist dann x = xend. *
;* y Naeherungswert fuer die Loesung des DGL-Systems an der neuen *
;* Stelle x *
;* h optimale, im letzten Schritt benutzte Schrittweite *
;* *
;* Funktionswert: *
;* ============== *
;* = 0: alles in Ordnung (Loesung errechnet an der Stelle xend) *
;* = 1: Nach MAX_FCT Aufrufen von dgl() wird abgebrochen, ohne dass *
;* i. a. der Endpunkt xend erreicht wurde. Soll weitergerechnet *
;* werden, muss man rk_fehl() mit unveraenderten Parametern *
;* erneut aufrufen. *
;* = 2: falsche Eingabedaten, d. h. *
;* epsabs < 0 oder *
;* epsrel < 0 oder *
;* epsabs + epsrel = 0 oder *
;* hmax <= 0 *
;* = 3: Die zur Weiterrechnung optimale Schrittweite kann im Rechner *
;* nicht dargestellt werden. *
;* = 4: Der Test auf den lokalen Fehler kann nicht ausgefuehrt werden, *
;* werden, weil epsabs und eine Komponente der berechneten *
;* Naeherungsloesung gleichzeitig Null sind. Daher ist ein reiner *
;* Test auf den relativen Fehler fuer diese Komponente unmoeglich. *
;***********************************************************************/
Procedure.l Num_RungKuttaFehlberg (*x.d,xend.d,n.l,y.d(1),*dgl.DglSysFkt,*h.d,hmax.d,epsabs.d,epsrel.d)
Protected Dim yt.d(n) ; Hilfsvektor fuer den Runge-Kutta-Schritt
Protected Dim t.d(n) ; Hilfsvektor fuer den letzten Teil des Runge-Kutta-Schritts
Protected Dim r.d(n) ; Vektor mit den Schaetzwerten fuer den lokalen Fehler
Protected Dim k1.d(n) ; \
Protected Dim k2.d(n) ; \
Protected Dim k3.d(n) ; \ Hilfsvektoren fuer den
Protected Dim k4.d(n) ; / Runge-Kutta-Schritt
Protected Dim k5.d(n) ; /
Protected Dim k6.d(n) ; /
Protected da.d ; Laenge des Integrationsintervalls
Protected xt.d ; Hilfsvariable fuer den eigentlichen Runge-Kutta-Schritt
Protected hmx.d ; maximale Schrittweite, hoechstens jedoch Laenge des Integrationsintervalls
Protected hf.d = 0.0 ; Hilfsvariable zur Aufbewahrung der Schrittweite
Protected quot.d ; Mass fuer die Genauigkeit des Integrationsschritts
Protected et.d ; epsrel * |yt[i]| + epsabs
Protected tr.d ; Hilfsvariable zur Berechnung von quot
Protected i.l ; Laufvariable
Protected aufrufe.l ; Anzahl benoetigter Aufrufe von dgl()
Protected Rep.l ; Flagge, die anzeigt, ob ein Schritt mit der gleichen Schrittweite wiederholt werden soll
Protected ende.l ; Flagge, die anzeigt, dass das Verfahrenbeendet werden kann
Protected ax.d = PeekD(*x)
Protected ah.d = PeekD(*h)
aufrufe = 0
Rep = #False
ende = #False
da = xend - ax
If epsrel < 0.0 Or epsabs < 0.0 Or epsrel + epsabs = 0.0 Or hmax <= 0.0 ; Plausibilitaetskontrollen
ProcedureReturn 2
EndIf
If Abs(da) <= 13.0 * #DBL_EPSILON * MathMax(Abs(ax), Abs(xend))
ProcedureReturn 3
EndIf
hmx = MathMin(hmax, Abs(da))
If Abs(ah) <= 13.0 * #DBL_EPSILON * Abs(ax)
ah = hmx
EndIf
Repeat
If Rep = #False ; h auf hmx begrenzen
ah = MathMin(Abs(ah), hmx) ; und so waehlen, dass der
If da < 0.0 ; Endpunkt xend erreicht wird
ah = -ah
EndIf
If Abs(xend - ax) <= 1.25 * Abs(ah) ; falls moeglich
hf = ah ; Falls ende gesetzt und h = xend - x
ende = #True ; zulaessig, wird nach dem naechsten
ah = xend - ax ; Integrationsschritt abgebrochen.
EndIf
*dgl(ax, y(), k1()) ; einen Integrationsschritt ausfuehren
aufrufe + 1
EndIf
xt = ah * 0.25
For i = 0 To n-1 Step 1
yt(i) = y(i) + xt * k1(i)
Next i
xt + ax
*dgl(xt, yt(), k2())
For i = 0 To n-1 Step 1
yt(i) = y(i) + ah * (k1(i) * 3.0/32.0 + k2(i) * 9.0/32.0)
Next i
xt = ax + ah * 0.375
*dgl(xt, yt(), k3())
For i = 0 To n-1 Step 1
yt(i) = y(i) + ah * (k1(i) * 1932.0 / 2197.0 - k2(i) * 7200.0 / 2197.0 + k3(i) * 7296.0 / 2197.0 )
Next i
xt = ax + ah * 12.0 / 13.0
*dgl(xt, yt(), k4())
For i = 0 To n-1 Step 1
yt(i) = y(i) + ah * (k1(i) * 439.0 / 216.0 - k2(i) * 8.0 + k3(i) * 3680.0 / 513.0 - k4(i) * 845.0 / 4104.0 )
Next i
xt = ax + ah
*dgl(xt, yt(), k5())
For i = 0 To n-1 Step 1
yt(i) = y(i) + ah * (-k1(i) * 8.0 / 27.0 + k2(i) * 2.0 - k3(i) * 3544.0 / 2565.0 + k4(i) * 1859.0 / 4104.0 - k5(i) * 11.0 / 40.0)
Next i
xt = ax + 0.5 * ah
*dgl(xt, yt(), k6())
For i = 0 To n-1 Step 1
t(i) = k1(i) * 25.0 / 216.0 + k3(i) * 1408.0 / 2565.0 + k4(i) * 2197.0 / 4104.0 - k5(i) * 0.2
yt(i) = y(i) + ah * t(i)
Next i
; yt ist jetzt das vorlaeufige Ergebnis des Schrittes. Nun wird
; r, die Schaetzung des lokalen Fehlers, relativ zur aktuellen
; Schrittweite berechnet.
For i = 0 To n-1 Step 1
r(i) = k1(i) / 360.0 - k3(i) * 128.0 / 4275.0 - k4(i) * 2197.0 / 75240.0 + k5(i) / 50.0 + k6(i) * 2.0 / 55.0
Next i
quot = 0.0
For i = 0 To n-1 Step 1 ; Genauigkeitstest
et = epsrel * Abs(yt(i)) + epsabs
If et = 0.0 ;Wenn die Naeherung verschwindet, ist eine reine Abfrage auf den relativen Fehler ungeeignet.
ProcedureReturn 5
EndIf
tr = Abs(r(i)) / et
quot = MathMax(quot, tr)
Next i
If quot <= 1.0 ;Ergebnis akzeptabel?
For i = 0 To n-1 Step 1
y(i) = yt(i)
Next i
ax + ah
If ende = #True ; Wenn xend erreicht wurde, wird das Verfahren beendet.
ah = hf
Break
EndIf
quot = MathMax(quot,0.00065336) ; naechsten Schritt vorbereiten
EndIf
quot = MathMin(quot,4096.0) ;h maximal um den Faktor 5 vergroessern bzw. um den Faktor 10 verkleinern
ah = 0.8 * ah / Sqr(Sqr(quot))
If Abs(ah) <= 13.0 * #DBL_EPSILON * Abs(ax)
ProcedureReturn 3
EndIf
aufrufe + 5
If aufrufe >= #MAX_FCT
ProcedureReturn 1
EndIf
If quot > 1.0 ;Schritt mit kleinerem h wiederholen?
Rep = #True
ende = #False
Else ;Schritt akzeptabel?
Rep = #False
EndIf
ForEver
PokeD (*x,ax)
PokeD (*h,ah)
ProcedureReturn 0
EndProcedure
;-
;- example
;-
Define i.l
; ### 1. Formulierung des mathematischen Problems
#E = 2.718281828459
#N = 3
Define x0.d = 0.0 ; x-Wert der Anfangswerte
Dim AW.d(#N) : AW(0)=1 : AW(1)=-1 : AW(2)=3 ; beliebige Anfangswerte
Global Dim a.d(#N) : a(0)=1.0 : a(1)=2.0 : a(2)=4.0 ; die a's in unseren DGLn y'=a*y
;Prototype DglSysFkt ( x.d, y.d(1), f.d(1) )
Procedure ExpFuncDgl(x.d, y.d(1), f.d(1))
; y' = a*y
; definiert implizit y als Exponentialfunktion: y(x)=C*e^(a*x)
; weil: y'/y=a => ln(y)'=a => ln(y)(x)=a*x+c => y(x)=e^(a*x+c) mit C := e^c
Protected i.l
For i=0 To #N-1
f(i) = a(i)*y(i)
Next
EndProcedure
; ### 2. Funktionsaufruf
Define h.d = 0.01 ; Anfangsschrittweite
Define x.d = x0 ; Kopie des Anfangs-x-Wert (wird vom Algorithmus als Ausgabe verwendet)
Dim y.d(#N) ; Kopien der Anfangs-y-Werte (werden vom Algorithmus als Ausgabe verwendet)
For i=0 To #N-1
y(i)=AW(i)
Next
;Procedure.l Num_RungKuttaFehlberg (*x.d,xend.d,n.l,y.d(1),*dgl.DglSysFkt,*h.d,hmax.d,epsabs.d,epsrel.d)
Define res.l = Num_RungKuttaFehlberg(@x, 1.0, #N, y(), @ExpFuncDgl(), @h, 10, 0.0, 0.0001)
; ### 3. Ergebnisausgabe
If res <> 0
Debug "Fehlercode " + Str(res)
Debug "xCheck = " + StrD(x)
Debug "hCheck = " + StrD(h)
Else
Debug "maximale Schrittweite bei Berechnung: " + StrD(h)
Debug "Auswertung an x = " + StrD(x)
For i=0 To #N-1
Debug "#### " + Str(i+1) + ". Gleichung: y'=" + StrD(a(i), 2) + "*y mit AW y(" + StrD(x0, 2) + ")=" + StrD(AW(i), 2)
Debug StrD(y(i)) + " (Runge-Kutta-Fehlberg)"
Debug StrD(AW(i) * Pow(#E, a(i)*1.0)) + " (direkt per Exponentialfunktion)"
Next
EndIf