p.specht
 | ... kann man z.B. mit dem Runge Kutta Fehlberg -Algorithmus (RKF45). Das Verfahren nutzt Schrittweitensteuerung nach dem Deutschen Mathematiker Erwin Fehlberg [Fehlberg 1970], und schafft durch Nutzung von XProfan's Double-Precision Arithmetik eine Genauigkeit von immerhin fünf Kommastellen, was bei solchen "expliziten" Verfahren absolut nicht selbstverständlich ist!
Hinweis: Bei sog. "steifen" (= widerspenstigen) Problemen wird die Sache allerdings sehr sehr langsam, dazu gibt es dann bessere Verfahren, etwa nach dem Prädiktor-Korrektor-Prinzip. Dazu an anderer Stelle mehr ...
Ich darf hier ein von Algol nach Fortran über QBasic nach C++ nach Pascal nach XProfan 11.2a übersetztes Standardstück der Numerischen Mathematik als DEMO zum privaten experimentieren OHNE JEDE GEWÄHR präsentieren, mit dem z.B. ein in den ersten 20 errechneten Sekunden realistisches Doppelpendel-Verhalten berechnet werden könnte (Mit Profan leider nicht in Echtzeit, aber nach 20 Sekunden schlägt die Chaostheorie mit ihrem Schmetterlingseffekt ohnehin voll zu!). Berechnet werden drei Test-Differentialgleichungen unterschiedlicher Komplexität. Die korrekten Ausgabewerte sind im Programmtext zu finden.
'RUNGE KUTTA der Ordnung 4 mit Schrittweitensteuerung der Ordnung 5 nach FEHLBERG"
'***********************************************************************
'* (DEMO ONLY!) 2014-05 Transcribed to XProfan 11 by P.Specht, Vienna *
'* ALL COPYRIGHTS BY THE RESPECTIVE OWNERS! RECHTE DRITTER BEACHTEN! *
'***********************************************************************
'* Integrate a System of Ordinary Differential Equations By the *
'* Runge-Kutta-Fehlberg method (double precision) *
'* ------------------------------------------------------------------- *
'* REFERENCE: H A Watts and L F Shampine, *
'* Sandia Laboratories, *
'* Albuquerque, New Mexico. *
'* ------------------------------------------------------------------- *
'* SAMPLE RUN: *
'* *
'* PROGRAM TRKF45 *
'* Demonstrate the RKF45 ODE integrator. *
'* *
'* TEST01 *
'* Solve a scalar equation: *
'* *
'* Y' = 0.25 * Y * ( 1 - Y / 20 ) *
'* *
'* T Y Y_Exact Error *
'* *
'* 0.00000 1.00000 1.00000 0.0000000 *
'* 4.00000 2.50321 2.50322 -0.0000087 *
'* 8.00000 5.60007 5.60009 -0.0000193 *
'* 12.00000 10.27774 10.27773 -0.0000069 *
'* 16.00000 14.83682 14.83683 -0.0000038 *
'* 20.00000 17.73017 17.73017 -0.0000084 *
'* *
'* TEST02 *
'* Solve a vector equation: *
'* *
'* Y'(1) = Y(2) *
'* Y'(2) = - Y(1) *
'* *
'* T Y1 Y2 *
'* *
'* 0.00000 1.00000 0.00000 *
'* 0.52360 0.86603 -0.50000 *
'* 1.04720 0.50000 -0.86603 *
'* 1.57080 0.00000 -1.00000 *
'* 2.09440 -0.50000 -0.86603 *
'* 2.61799 -0.86603 -0.50000 *
'* 3.14159 -1.00000 -0.00000 *
'* 3.66519 -0.86603 0.50000 *
'* 4.18879 -0.50000 0.86603 *
'* 4.71239 -0.00000 1.00001 *
'* 5.23599 0.50000 0.86604 *
'* 5.75959 0.86604 0.50001 *
'* 6.28319 1.00002 0.00000 *
'* *
'* TEST03 *
'* Solve a vector equation: *
'* *
'* Y'(1) = Y(2) *
'* Y'(2) = Y(3) *
'* Y'(3) = Y(4) *
'* Y'(4) = Y(5) *
'* Y'(5) = (45 * Y(3) * Y(4) * Y(5) - 40 * Y(4)^3) / (9 * Y(3)^2) *
'* *
'* T Y1 Y2 Y3 Y4 Y5 *
'* *
'* 0.00000 1.00000 1.00000 1.00000 1.00000 1.00000 *
'* 0.13636 1.14610 1.14609 1.14587 1.14068 1.05604 *
'* 0.27273 1.31354 1.31340 1.31128 1.28532 1.05248 *
'* 0.40909 1.50538 1.50460 1.49612 1.42333 0.95209 *
'* 0.54545 1.72508 1.72223 1.69844 1.53859 0.71111 *
'* 0.68182 1.97638 1.96840 1.91370 1.60897 0.28809 *
'* 0.81818 2.26328 2.24438 2.13400 1.60781 -0.33918 *
'* 0.95455 2.58984 2.55011 2.34770 1.50801 -1.15027 *
'* 1.09091 2.96003 2.88369 2.53985 1.28946 -2.06094 *
'* 1.22727 3.37739 3.24105 2.69376 0.94820 -2.92046 *
'* 1.36364 3.84475 3.61589 2.79372 0.50379 -3.54302 *
'* 1.50000 4.36396 4.00000 2.82843 -0.00000 -3.77124 *
'* ------------------------------------------------------------------- *
'* *
'* Basic Release 1.1 by J-P Moreau, Paris. *
'* (www.jpmoreau.fr) *
'* Release 1.1: added test #3. *
'***********************************************************************
'* Source: Homepg. of J-P. Moreau: DIFFERENTIAL EQUATIONS IN BASIC.htm *
'* ------------------------------------------------------------------- *
'* All Copyrights (C) 2014- are attributed to their respective Owners *
'***********************************************************************
Windowstyle 24:randomize:font 2
Window 0,0-%maxx,%maxy-40:cls co()
var xx&=width(%hwnd)\2:var yy&=height(%hwnd)\2
set("decimals",17):set("numwidth",24)
proc co
return rgb(200+rnd(56),200+rnd(56),200+rnd(56))
endproc
Declare R$,num&,neqn&,F$,G$,yexact!,istep&,t!,remin!,maxnfe&,eps!
declare mflag&,jflag&,kflag&,savr!,save!,rer!,dt!,init&,kop&,a!,b!
declare h!,i&,k&,nfe&,toln!,tol!,ypk!,tmp!,ep!,IA&,IB&,ISign&,sign!
declare kop!,IOUTPUT&,scale!,ae!,IHFAILD&,hmin!,ch!
declare eedet!,eeoet!,et!,ee!,eecet!,xmax!,esttol!,s!
var abserr! = 0.000001
var relerr! = 0.000001
var iflag& = 1
var tstart! = 0.0
var tstop! = 20.0
var nstep& = 5
var tout! = 0
var tt! = tout!
var NEQ& = 5' Maximum number of equations
def !Pi Pi()
var EPSILON! = val("2.22E-16")' Small number
Declare y![NEQ&], yp![NEQ&]'auxiliary variables used by 400
Declare yy![NEQ&], yyp![NEQ&]'work space used by successive calls to 1000 and 2000
Declare f1![NEQ&], f2![NEQ&], f3![NEQ&], f4![NEQ&], f5![NEQ&]
F$ = "%g"'"####0.0########"
G$ = "%g"'"####0.0##############"
CLS co()
PRINT "\n PROGRAM RKF45 AND TRKF45 FOR TESTING RKF45 \n\n"
PRINT " Demonstrates the RKF45 ODE integrator\n"
'call test01
S600
Print "\n Press any key to continue... ",:beep:Waitinput:cls co()
'call test02
S700
Print "\n Press any key to continue... ",:beep:Waitinput:cls co()
'call test03
S800
Print "\n Press ENTER to end the PGM. ",:beep:Input R$
END'of Main Pgm
S400:
'{ User defined system of Differential Equations
' Subroutine S400f :parameters tt!, yy!, yyp!
'--------------------------------------------------------------------
' Fn evaluates the derivative for the ODE (TESTS #1 and #2).
'--------------------------------------------------------------------
IF num& = 1' Test-No.
yyp![1] = 0.25 * yy![1] * (1.0 - yy![1]/20.0)
ELSEIF num& = 2
yyp![1] = yy![2]
yyp![2] = -1*yy![1]
ELSE
yyp![1] = yy![2]
yyp![2] = yy![3]
yyp![3] = yy![4]
yyp![4] = yy![5]
yyp![5] = (45*yy![3]*yy![4]*yy![5] - 40*yy![4]*sqr(yy![4])) / (9*sqr(yy![3]))
ENDIF
RETURN
S500:
'--------------------------------------------------------------------
' function yexact(tt)
' YEXACT evaluates the exact solution of the ODE (For TEST #1).
'--------------------------------------------------------------------
yexact! = 20/(1+19*EXP(-0.25*tt!))
RETURN
'}
proc S600
'--------------------------------------------------------------------
' TEST01 solves a scalar ODE in double precision.
'--------------------------------------------------------------------
num& = 1'example #1
neqn& = 1'one equation
PRINT " TEST01 \n\n"
PRINT " Solve a scalar equation:"
PRINT
PRINT " Y' = 0.25 * Y * ( 1 - Y / 20 )\n\n"
PRINT
abserr! = .000001
relerr! = .000001
iflag& = 1
tstart! = 0.0
tstop! = 20.0
nstep& = 5
tout! = 0
y![1] = 1
PRINT " T Y Y_Exact Error "
PRINT
tt! = tout!: GOSUB "S500"'calculate yexact(tout)
PRINT " ";format$(F$, tout!),
PRINT tab(26);format$(F$,y![1]),
PRINT tab(52);format$(F$,yexact!),
PRINT tab(77);format$(G$,y![1] - yexact!)
Whileloop nstep& : istep&=&Loop
t! = ((nstep& - istep& + 1) * tstart! + (istep& - 1) * tstop!) / nstep&
tout! = ((nstep& - istep&) * tstart! + (istep&) * tstop!) / nstep&
'Bracket error in French source corrected!
'call rkfs (neqn,y,t,tout,relerr,abserr,iflag,yp,h,f1..f5,savr,save,nfe,kop,init,jflag,kflag)
GOSUB "S2000"
tt! = tout!
GOSUB "S500"' calculate yexact(tout)
PRINT " ";format$(F$, tout!),
PRINT tab(26);format$(F$, y![1]),
PRINT tab(52);format$(F$, yexact!),
PRINT tab(77);format$(G$, y![1] - yexact!)
endwhile
endproc
Proc S700
'--------------------------------------------------------------------
' TEST02 solves a vector ODE (Ordinary Differential Equation)
'--------------------------------------------------------------------
num& = 2'Example #2
neqn& = 2'2 equations
PRINT
PRINT " TEST02\n"
PRINT " Solve a vector equation:"
PRINT
PRINT " Y'(1) = Y(2)"
PRINT " Y'(2) = - Y(1)\n\n"
abserr! = 0.000001
relerr! = 0.000001
iflag& = 1
tstart! = 0
tstop! = 2*PI()
nstep& = 12
tout! = 0
y![1] = 1
y![2] = 0
PRINT
PRINT " T Y1 Y2 "
PRINT
PRINT " ";format$(F$, tout!),
PRINT tab(26);format$(F$, y![1]),
PRINT tab(52);format$(F$, y![2])
Whileloop nstep&
istep&=&Loop
t! = ((nstep& - istep& + 1) * tstart! + (istep& - 1) * tstop!) / nstep&
tout! = ((nstep& - istep&) * tstart! + (istep& * tstop!)) / nstep&')
GOSUB "S2000"'call rkfs (neqn,y,t,tout,relerr,abserr,iflag,yp,h,f1..f5,savr,save,
'nfe,kop,init,jflag,kflag)
PRINT " ";format$(F$, tout!),
PRINT tab(26);format$(F$,y![1]),
PRINT tab(52);format$(F$,y![2])
endwhile
EndProc
Proc S800
'--------------------------------------------------------------------
' TEST03 solves a vector ODE.
'--------------------------------------------------------------------
num& = 3'Example #3
neqn& = 5'5 equations
PRINT "\n TEST03 \n"
PRINT " Solve a vector equation: \n"
PRINT " Y'(1) = Y(2)"
PRINT " Y'(2) = Y(3)"
PRINT " Y'(3) = Y(4)"
PRINT " Y'(4) = Y(5)"
PRINT " Y'(5) = (45 * Y(3) * Y(4) * Y(5) - 40 * Y(4)^3) / (9 * Y(3)^2) \n"
abserr! = .000001
relerr! = .000001
iflag& = 1
tstart! = 0
tstop! = 1.5
nstep& = 11
tout! = 0
y![1] = 1
y![2] = 1
y![3] = 1
y![4] = 1
y![5] = 1
PRINT "\n T Y1 Y2 Y3 Y4 Y5 \n"
PRINT " ";format$(F$,tout!);
PRINT tab(26);format$(F$,y![1]);
PRINT tab(52);format$(F$,y![2]);
PRINT tab(77);format$(F$,y![3]);
PRINT tab(104);format$(F$,y![4]);
PRINT tab(130);format$(F$,y![5])
Whileloop nstep&
istep&=&Loop
t! = ((nstep& - istep& + 1) * tstart! + (istep& - 1) * tstop!) / nstep&
tout! = ((nstep& - istep&) * tstart! + (istep& * tstop!) / nstep&)
GOSUB "S2000"'call rkfs (neqn,y,t,tout,relerr,abserr,iflag,yp,h,f1..f5,savr,save,nfe,kop,init,jflag,kflag)
PRINT " ";format$(F$,tout!);
PRINT tab(26);format$(F$,y![1]);
PRINT tab(52);format$(F$,y![2]);
PRINT tab(77);format$(F$,y![3]);
PRINT tab(104);format$(F$,y![4]);
PRINT tab(130);format$(F$,y![5])
Endwhile
EndProc
'{* Differential Equations By the Runge-Kutta-Fehlberg method (double precision)
S1000:
'Subroutine fehl(neqn, y, t, h, yp, f1, f2, f3, f4, f5, s)
ch! = h! / 4
whileloop neqn&:i&=&Loop
f5![i&] = y![i&] + ch! * yp![i&]
endwhile
tt! = t! + ch!
whileloop neqn&:i&=&Loop
yy![i&] = f5![i&]
endwhile
GOSUB "S400"'call f(t+ch, f5, f1)
whileloop neqn&:i&=&Loop
f1![i&] = yyp![i&]
endwhile
ch! = 3*h!/32
whileloop neqn&:i&=&Loop
f5![i&] = y![i&] + ch! * (yp![i&] + 3 * f1![i&])
endwhile
tt! = t!+3*h!/8
whileloop neqn&:i&=&Loop
yy![i&] = f5![i&]
endwhile
GOSUB "S400"' call f(t+3*h/8, f5, f2)
whileloop neqn&:i&=&Loop
f2![i&] = yyp![i&]
endwhile
ch! = h!/2197
whileloop neqn&:i&=&Loop
f5![i&] = y![i&] + ch! * (1932 * yp![i&] + (7296 * f2![i&] - 7200 * f1![i&]))
endwhile
tt! = t!+12*h!/13
whileloop neqn&:i&=&Loop
yy![i&] = f5![i&]
endwhile
GOSUB "S400"'call f(t+12#*h/13#, f5, f3)
whileloop neqn&:i&=&Loop
f3![i&] = yyp![i&]
endwhile
ch! = h!/4104
whileloop neqn&:i&=&Loop
f5![i&] = y![i&] + ch! * ((8341 * yp![i&] - 845 * f3![i&]) + (29440 * f2![i&] - 32832 * f1![i&]))
endwhile
tt! = t! + h!
whileloop neqn&:i&=&Loop
yy![i&] = f5![i&]
endwhile
GOSUB "S400"'call f(t+h, f5, f4)
whileloop neqn&:i&=&Loop
f4![i&] = yyp![i&]
endwhile
ch! = h!/20520
whileloop neqn&:i&=&Loop
tmp! = (41040 * f1![i&] - 28352 * f2![i&])
f1![i&] = y![i&] + ch! * ((-6080 * yp![i&] + (9295 * f3![i&] - 5643 * f4![i&])) + tmp!)
endwhile
tt! = t!+h!/2
whileloop neqn&:i&=&Loop
yy![i&] = f1![i&]
endwhile
GOSUB "S400"'call f(t+h/2#, f1, f5)
whileloop neqn&:i&=&Loop
f5![i&] = yyp![i&]
endwhile
ch! = h! / 7618050
whileloop neqn&:i&=&Loop
tmp! = (3953664 * f2![i&] + 277020 * f5![i&])
f1![i&] = y![i&] + ch! * ((902880 * yp![i&] + (3855735 * f3![i&] - 1371249 * f4![i&])) + tmp!)
endwhile
RETURN
S1200:
IF ib& < 0 : ISign& = -ABS(ia&)
ELSE : ISign& = ABS(ia&)
ENDIF
RETURN
S1210:
IF b!<0 :Sign! = -1*ABS(a!)
ELSE :Sign! = ABS(a!)
ENDIF
RETURN
S1300:
IF a!>=b!:XMax! = a!
ELSE :XMax! = b!
ENDIF
RETURN
S1310:
IF a!<=b!:XMIN! = a!
ELSE :XMIN! = b!
ENDIF
RETURN
'}
S2000:
'{Subroutine rkfs(neqn, y, t, tout, relerr, abserr, iflag, yp, h, f1, f2, f3, f4,f5,
' savr, save, nfe, kop, init, jflag, kflag)
'***************************************************************************************
' RKFS implements the Runge-Kutta-Fehlberg method (double precision).
'***************************************************************************************
' Labels: 25,40,45,50,60,65,80,100,200,260
remin! = .000000000001
maxnfe& = 3000
eps! = EPSILON!
IF neqn& < 1:iflag& = 8:RETURN
ENDIF
IF relerr! < 0:iflag& = 8:RETURN
ENDIF
IF abserr! < 0:iflag& = 8:RETURN
ENDIF
mflag& = ABS(iflag&)
IF (ABS(iflag&) < 1) OR (ABS(iflag&) > 8)
iflag& = 8
RETURN
ENDIF
case mflag& = 1 : GOTO "G50"
IF (t! = tout!) AND (kflag& <> 3)
iflag& = 8
RETURN
ENDIF
Case mflag& <> 2 : GOTO "G25"
Case kflag& = 3 : GOTO "G45"
Case init& = 0 : GOTO "G45"
Case kflag& = 4 : GOTO "G40"
Case (kflag& = 5) AND (abserr! = 0): END
Case (kflag& = 6) AND (relerr! <= savr!) AND (abserr! <= save!): END
GOTO "G50"
' iflag = 3,4,5,6,7 or 8
G25:
Case iflag& = 3 : GOTO "G45"
Case iflag& = 4 : GOTO "G40"
Case (iflag& = 5) AND (abserr > 0) : GOTO "G45"
END
G40:
nfe& = 0
Case mflag& = 2 : GOTO "G50"
G45:
iflag& = jflag&
Case kflag& = 3 : mflag& = ABS(iflag&)
G50:
jflag& = iflag&
kflag& = 0
savr! = relerr!
save! = abserr!
rer! = 2 * EPSILON! + remin!
IF relerr! < rer!
relerr! = rer!
iflag& = 3
kflag& = 3
RETURN
ENDIF
dt! = tout! - t!
Case mflag& = 1 : GOTO "G60"
Case init& = 0 : GOTO "G65"
GOTO "G80"
G60:
init& = 0
kop& = 0
a! = t!
whileloop neqn&:i&=&Loop
yy![i&] = y![i&]
endwhile
GOSUB "S400"'call f(a, y, yp)
whileloop neqn&:i&=&Loop
yp![i&] = yyp![i&]
endwhile
nfe& = 1
IF t! = tout! : iflag& = 2 : RETURN
ENDIF
G65:
init& = 1
h! = ABS(dt!)
toln! = 0
whileloop neqn&:k&=&Loop
tol! = relerr! * ABS(y![k&]) + abserr!
IF tol! > 0
toln! = tol!
ypk! = ABS(yp![k&])
IF (ypk!*h!^5) > tol!
h! = (tol! / ypk!)^0.2
ENDIF
ENDIF
endwhile
case toln! <= 0 : h! = 0
IF ABS(t!) > ABS(dt!) : tmp! = ABS(t!)
ELSE : tmp! = ABS(dt!)
ENDIF
IF h! < (26 * ep! * tmp!) : h! = 26 * eps! * tmp!
ENDIF
ia& = 2: ib& = iflag&: GOSUB "S1200"
jflag& = ISign&
G80:
a! = h!: b! = dt!: GOSUB "S1210"
h! = Sign!
Case ABS(h!) >= (2*ABS(dt!)): kop& = kop& + 1
IF kop! = 100
kop! = 0
iflag& = 7
RETURN
ENDIF
IF ABS(dt!) <= (26 * eps! * ABS(t!))
WhileLoop neqn&:i&=&Loop
y![i&] = y![i&] + dt! * yp![i&]
Endwhile
a! = tout!
WhileLoop neqn&:i&=&Loop
yy(i) = y(i)
endwhile
GOSUB "S400"'call f(a, y, yp)
WhileLoop neqn&:i&=&Loop
yp![i&] = yyp![i&]
endwhile
nfe& = nfe& + 1
t! = tout!
iflag! = 2
RETURN
ENDIF
ioutput& = 0
scale! = 2 / relerr!
ae! = scale! * abserr!
G100:
ihfaild& = 0
hmin! = 26 * eps! * ABS(t!)
dt! = tout! - t!
Case ABS(dt!) >= (2 * ABS(h!)) : GOTO "G200"
IF ABS(dt!) <= ABS(h!)
ioutput& = 1
h! = dt!
GOTO "G200"
ENDIF
h! = 0.5 * dt!'reduce step
G200:
IF nfe& > maxnfe&
iflag& = 4
kflag& = 4
RETURN
ENDIF
GOSUB "S1000"'call fehl(neqn, y, t, h, yp, f1, f2, f3, f4, f5, f1)
nfe& = nfe& + 5
eeoet! = 0
Whileloop neqn&:k&=&Loop
et! = ABS(y![k&]) + ABS(f1![k&]) + ae!
IF et! <= 0
iflag& = 5
RETURN
ENDIF
tmp! = (22528 * f2![k&] - 27360 * f5![k&])
ee! = ABS((-2090 * yp![k&] + (21970 * f3![k&] - 15048 * f4![k&])) + tmp!)
a! = eecet!: b! = ee! / et!: GOSUB "S1300"
eecet! = XMax!
endwhile
esttol! = ABS(h!) * eeoet! * scale! / 752400
Case esttol! <= 1 : GOTO "G260"
ihfaild& = 1
ioutput& = 0
IF esttol! < 59049
s! = 0.9 / (esttol!^0.2)
ELSE
s! = 0.1
ENDIF
h! = s! * h!
IF ABS(h!) < hmin!
iflag& = 6
kflag& = 6
RETURN
ELSE
GOTO "G200"
ENDIF
G260:
t! = t! + h!
whileloop neqn&:i&=&Loop
y![i&] = f1![i&]
endwhile
a! = t!
whileloop neqn&:i&=&Loop
yy![i&] = y![i&]
endwhile
GOSUB "S400"' call f(a, y, yp)
whileloop neqn&:i&=&Loop
yp![i&] = yyp![i&]
endwhile
nfe& = nfe& + 1
IF esttol! > .0001889568
s! = 0.9/esttol^0.2
ELSE
s! = 5
ENDIF
IF ihfaild& <> 0
a! = s!: b! = 1: GOSUB "S1310": s! = XMIN!
ENDIF
a! = s! * ABS(h!): b! = hmin!: GOSUB "S1300"
a! = XMax!: b! = h!: GOSUB "S1210": h! = Sign!
IF ioutput& <> 0
t! = tout!
iflag& = 2
ENDIF
Case iflag& > 0 : GOTO "G100"
iflag& = -2
RETURN
'}
|
|
(851)
Deutsch
English
Français
Español
Italia