p.specht
 | In der Elektrotechnik, aber auch im Maschinenbau (Hydraulik, Thermodynamik) und z.B. in der Regeltechnik, kommen häufig physikalische Systeme vor, die am einfachsten durch lineare Differentialgleichungen 1. Ordnung beschrieben werden können. Eine bekannte, definierte Ausgangssituation vorausgesetzt, kann man dann nach dem konkreten Verhalten bestimmter Systemwerte im Zeitablauf fragen.
Die Methoden zur Nutzung des Computers bei der Beantwortung solcher Fragen sind inzwischen sehr vielfältig und weisen je nach Aufgabenstellung verschiedene Stärken und Schwächen auf. Vier auch heute noch sehr gängige Verfahren sind im nachstehenden Programm anwählbar, sodaß ein Algorithmenvergleich möglich wird.
Die Vorlage zur XProfan-Übersetzung stammt aus 1987 und wurde offenbar zuvor von Fortran77 nach HP-BASIC übersetzt - daher kann man hier mit Fug und Recht von ALTBEWÄHREM reden. Erste Tests hat das Zeug überstanden, Garantie gibts aber wie immer keine! Reine Demo, da Rechtslage unklar.
WindowTitle "Anfangswertprobleme von Differentialgleichungen 1. Ordnung lösen mit verschiedenen Methoden"
Proc FNF : parameters t!,y!
declare DY!
DY! = (T!-Y!)/2
' DY! = 2 * t! * y!
return DY!
endproc
proc s30'SUBROUTINE PRINT FUNCTION
PRINT " Y' = (T-Y)/2"
RETURN
endproc
'{ (TW) 2012-06 Translationware to XProfan 11.2a by P.Specht@gmx.at
' Source: https://read.pudn.com/downloads168/sourcecode/math/775662/BASIC/ALG9%605-8.BAS__.htm
' Keine wie auch immer geartete Gewähr. Rechtsvermerke des Originaltextes:
REM Copyright 1987
REM John H. Mathews
REM Dept. of Mathematics
REM California State University, Fullerton
REM Runge-Kutta-Fehlberg Method, RKF45
REM Adams-Bashforth-Moulton Method
REM Milne-Simpson Method
REM Hamming's Method
Window 0,0 - %maxx,%maxy-44
Font 2:Randomize:set("decimals",15)
Declare METH&,A!,B!,Y0!,M&,STATE&,STAT&,RESP!,DOMO&,I&,J&,K&,ANS$
Declare Valu!,TOL!,H!,TJ!,TK!,YJ!,YK!,K1!,K2!,K3!,K4!,H2!,F0!,F1!,F2!,F3!,F4!,F5!
Declare MEND&,POLD!,YOLD!,PNEW!,PMOD!,P!,COLD!,CNEW!
Declare HMIN!,HMAX!,BR!,Y1!,Y2!,Y3!,Y4!,Y5!,Y6!,K5!,K6!,ERRER!,YNEW!,S!
'}
g100:
REM PROGRAM SOLVE DE'S
var MAXM& = 1000
Declare D![4],T![MAXM&],Y![MAXM&]
METH&=1
A!=0
B!=0
Y0!=0
M&=1
STATE&=1
g118:
REM SUBROUTINE MESSAGE
s1000' Do Subroutine
DOMO&=1
g124:
s2000' Do INPUTS
g128:
s3000' Do END POINTS
IF METH& = 1 :goto "g134" : ELSE : goto "g138": Endif
g134:
s300' Do ADAMS-BASHFORTH-MOULTON
MEND&=M&
GOTO "g176"
g138:
REM ELSE
IF METH& = 2: goto "g142": ELSE : goto "g146" : endif
g142:
s400' Do MILNE
MEND&=M&
GOTO "g176"
g146:
REM ELSE
IF METH& = 3: goto "g150": ELSE : goto "g154" : endif
g150:
s500' Do HAMMING
MEND&=M&
GOTO "g176"
g154:
IF METH& = 4 : goto "g158" : ELSE : goto "g162" : endif
g158:
s600' Do RKF45
MEND&=M&
GOTO "g176"
g162:
REM ELSE
g176:
REM ENDIF
s4000' do RESULTS
PRINT
PRINT " Want to try another initial condition ? <Y/N> ";
INPUT ANS$
IF (ANS$ <> "Y") AND (ANS$ <> "y") : goto "g188" : ELSE : goto "g192" : endif
g188:
STATE&=2
GOTO "g196"
g192:
STATE&=0
g196:
REM ENDIF
IF (STATE& = 0) OR (STATE& = 1) : GOTO "g128" : endif
PRINT " Want to change the differential equation ? <Y/N> ";
INPUT ANS$
IF (ANS$ <> "Y") AND (ANS$ <> "y") :goto "g206" : ELSE : goto "g210":endif
g206:
DOMO&=0
GOTO "g214"
g210:
REM ELSE
STATE&=0
g214:
REM ENDIF
IF DOMO& = 1 :GOTO "g124" : endif
PRINT " Want to try another method of approximation? <Y/N> ";
INPUT ANS$
IF (ANS$ <> "Y") AND (ANS$ <> "y") : goto "g224" : ELSE : goto "g228" : endif
g224:
METH&=0
GOTO "g232"
g228:
REM ELSE
STATE&=0
g232:
REM ENDIF
IF METH& <> 0 :GOTO "g118":endif
GOTO "g5000"
Proc s300
REM SUBROUTINE ADAMS-BASHFORTH-MOULTON
H!=(B!-A!)/M&
T![0]=A!
Y![0]=Y0!
REM Use Runge-Kutta to compute three starting values.
WhileLoop 0,2,1: J& = &Loop
TJ!=T![J&]
YJ!=Y![J&]
K1!=H!*FNF(TJ!,YJ!)
K2!=H!*FNF(TJ!+H!/2,YJ!+0.5*K1!)
K3!=H!*FNF(TJ!+H!/2,YJ!+0.5*K2!)
K4!=H!*FNF(TJ!+H!,YJ!+K3!)
Y![J&+1]=YJ!+(K1!+2*K2!+2*K3!+K4!)/6
T![J&+1]=A!+H!*(J&+1)
Endwhile
REM Start of Procedure Adams-Bashforth-Moulton
F0!=FNF(T![0],Y![0])
F1!=FNF(T![1],Y![1])
F2!=FNF(T![2],Y![2])
F3!=FNF(T![3],Y![3])
H2!=H!/24
WhileLoop 3,m&-1,-1 : K&=&Loop
TK!=T![K&]
YK!=Y![K&]
P!=YK!+H2!*(-9*F0!+37*F1!-59*F2!+55*F3!)
T![K&+1]=A!+H!*(K&+1)
F4!=FNF(T![K&+1],P!)
Y![K&+1]=YK!+H2!*(F1!-5*F2!+19*F3!+9*F4!)
F0!=F1!
g358:
F1!=F2!
F2!=F3!
F3!=FNF(T![K&+1],Y![K&+1])
EndWhile
RETURN
Endproc
proc s400
REM SUBROUTINE MILNE
H!=(B!-A!)/M&
T![0]=A!
Y![0]=Y0!
REM Use Runge-Kutta to compute three starting values.
WhileLoop 0,2,1: J& = &Loop
TJ!=T![J&]
YJ!=Y![J&]
K1!=H!*FNF(TJ!,YJ!)
K2!=H!*FNF(TJ!+H!/2,YJ!+0.5*K1!)
K3!=H!*FNF(TJ!+H!/2,YJ!+0.5*K2!)
K4!=H!*FNF(TJ!+H!,YJ!+K3!)
Y![J&+1]=YJ!+(K1!+2*K2!+2*K3!+K4!)/6
T![J&+1]=A!+H!*(J&+1)
EndWhile
REM Start of Procedure Milne-Simpson
F1!=FNF(T![1],Y![1])
F2!=FNF(T![2],Y![2])
F3!=FNF(T![3],Y![3])
POLD!=0
YOLD!=0
WhileLoop 3,m&-1:K&=&Loop
PNEW!=Y![K&-3]+4*H!*(2*F1!-F2!+2*F3!)/3
PMOD!=PNEW!+28*(YOLD!-POLD!)/29
T![K&+1]=A!+H!*(K&+1)
F4!=FNF(T![K&+1],PMOD!)
Y![K&+1]=Y![K&-1]+H!*(F2!+4*F3!+F4!)/3
POLD!=PNEW!
YOLD!=Y![K&+1]
F1!=F2!
F2!=F3!
F3!=FNF(T![K&+1],Y![K&+1])
Endwhile
RETURN
Endproc
proc s500
REM SUBROUTINE HAMMING
H!=(B!-A!)/M&
T![0]=A!
Y![0]=Y0!
REM Use Runge-Kutta to compute three starting values.
WhileLoop 0,2:J&=&Loop
TJ!=T![J&]
YJ!=Y![J&]
K1!=H!*FNF(TJ!,YJ!)
K2!=H!*FNF(TJ!+H!/2,YJ!+0.5*K1!)
K3!=H!*FNF(TJ!+H!/2,YJ!+0.5*K2!)
K4!=H!*FNF(TJ!+H!,YJ!+K3!)
Y![J&+1]=YJ!+(K1!+2*K2!+2*K3!+K4!)/6
T![J&+1]=A!+H!*(J&+1)
Endwhile
REM Start of Procedure Hamming
F1!=FNF(T![1],Y![1])
F2!=FNF(T![2],Y![2])
F3!=FNF(T![3],Y![3])
POLD!=0
COLD!=0
WhileLoop 3,m&-1 : k&=&Loop
PNEW!=Y![K&-3]+4*H!*(2*F1!-F2!+2*F3!)/3
PMOD!=PNEW!+112*(COLD!-POLD!)/121
T![K&+1]=A!+H!*(K&+1)
F4!=FNF(T![K&+1],PMOD!)
CNEW!=(9*Y![K&]-Y![K&-2]+3*H!*(-F2!+2*F3!+F4!))/8
Y![K&+1]=CNEW!+9*(PNEW!-CNEW!)/121
POLD!=PNEW!
COLD!=CNEW!
F1!=F2!
F2!=F3!
F3!=FNF(T![K&+1],Y![K&+1])
EndWhile
RETURN
endproc
proc s600
REM SUBROUTINE RKF45 Runge Kutta Fehlberg
Declare BIG!,A2!,B2!,A3!,B3!,C3!,A4!,B4!,C4!,D4!,A5!,B5!,C5!,D5!,E5!
Declare A6!,B6!,C6!,D6!,E6!,F6!,R1!,R3!,R4!,R5!,R6!,N1!,N3!,N4!,N5!
BIG!=1*10^13
A2!= 0.25
B2!= 0.25
A3!= 0.375
B3!= 0.09375
C3!= 0.28125
A4!= 0.9230769231
B4!= 0.8793809741
C4!= -3.277196177
D4!= 3.320892126
A5!= 1.0
B5!= 2.032407407
C5!= -8.0
D5!= 7.173489279
E5!= -0.2058966862
A6!= 0.5
B6!= -0.2962962963
C6!= 2.0
D6!= -1.381676413
E6!= 0.4529727096
F6!= -0.275
R1!= 0.002777777778
R3!= -0.02994152047
R4!= -0.02919989367
R5!= 0.02
R6!= 0.03636363636
N1!= 0.1157407407
N3!= 0.5489278752
N4!= 0.535331384
N5!= -0.2
H!=(B!-A!)/M&
HMIN!=H!/1024
HMAX!=H!*64
T![0]=A!
Y![0]=Y0!
T![0]=A!
J&=0
TJ!=A!
BR!=B!-0.0001*ABS(B!)
WHILE TJ! < B!
case (TJ!+H!) > BR! : H!=B!-TJ!
TJ!=T![J&]
YJ!=Y![J&]
Y1!=YJ!
K1!=H!*FNF(TJ!,Y1!)
Y2!=YJ!+B2!*K1!
IF BIG! < ABS(Y2!) : GOTO "g720" : endif
K2!=H!*FNF(TJ!+A2!*H!,Y2!)
Y3!=YJ!+B3!*K1!+C3!*K2!
IF BIG! < ABS(Y3!) : GOTO "g720" : endif
K3!=H!*FNF(TJ!+A3!*H!,Y3!)
Y4!=YJ!+B4!*K1!+C4!*K2!+D4!*K3!
IF BIG! < ABS(Y4!) : GOTO "g720" : endif
K4!=H!*FNF(TJ!+A4!*H!,Y4!)
Y5!=YJ!+B5!*K1!+C5!*K2!+D5!*K3!+E5!*K4!
IF BIG! < ABS(Y5!) : GOTO "g720" : endif
K5!=H!*FNF(TJ!+A5!*H!,Y5!)
Y6!=YJ!+B6!*K1!+C6!*K2!+D6!*K3!+E6!*K4!+F6!*K5!
IF BIG! < ABS(Y6!) : GOTO "g720" : endif
K6!=H!*FNF(TJ!+A6!*H!,Y6!)
ERRER!=ABS(R1!*K1!+R3!*K3!+R4!*K4!+R5!*K5!+R6!*K6!)
YNEW!=YJ!+N1!*K1!+N3!*K3!+N4!*K4!+N5!*K5!
ERRER!=ABS(ERRER!)
IF (ERRER! < TOL!) OR (H! < (2*HMIN!)) :goto "g690" : ELSE : goto "g699" : endif
g690:
Y![J&+1]=YNEW!
IF (TJ!+H!) > BR! :goto "g692" : ELSE : goto "g694" : endif
g692:
T![J&+1]=B!
GOTO "g696"
g694:
REM ELSE
T![J&+1]=TJ!+H!
g696:
REM ENDIF
J&=J&+1
TJ!=T![J&]
g699:
REM CONTINUE
IF ERRER! = 0 :goto "g702" : ELSE : goto "g706" : endif
g702:
S!=0
GOTO "g710"
g706:
REM ELSE
S!=0.84*(TOL!*H!/ERRER!)^0.25
g710:
REM ENDIF
case (S! < 0.75) AND (H! > (2*HMIN!)) : H!=H!/2
case (S! > 1.50) AND ((2*H!) < HMAX!) : H!=H!*2
case (BIG! < ABS(Y![J&])) OR (MAXM& = J&) : GOTO "g720"
Endwhile
g720:
REM CONTINUE
MEND&=J&
IF B! > T![J&] : goto "g726" : ELSE : goto "g730" : endif
g726:
M&=J&+1
GOTO "g734"
g730:
REM ELSE
M&=J&
g734:
REM ENDIF
RETURN
endproc
proc s1000
REM SUBROUTINE MESSAGE
CLS
PRINT
PRINT " Solution of the Differential Equation (DE) Y' = F(T,Y) "
PRINT " with the initial condition Y(A) = Y0 "
PRINT
PRINT " A numerical approximation is computed in Interval [A,B] "
PRINT
PRINT "\n Choose the method of Approximation:"
PRINT
PRINT "\n < 1 > Adams-Bashforth-Moulton method"
PRINT " (Bashford-predicted Adams-Moulton Algorithm)"
PRINT " + Good if many Steps are to be processed, "
PRINT " - Bad if F(t,Y)-Values vary largely. "
PRINT "\n < 2 > Milne-Simpson's method"
PRINT " + Multistep predictor/corrector method using 4 initial values."
PRINT "\n < 3 > Hamming's method"
PRINT " Examines recent behaviour to predict near-term future, then correct it."
PRINT "\n < 4 > Runge-Kutta-Fehlberg"
PRINT " of Order 4, with Fehlberg's Step size adaption for error limitation."
PRINT
PRINT " SELECT < 1 - 4 > ? ";
RESP!=METH&
INPUT RESP!
METH&=INT(RESP!)
case ( (METH& < 1) OR (METH& > 4) ) AND (STATE& <> 0) : METH&=1
RETURN
endproc
proc s2000
REM SUBROUTINE INPUTS
CLS
PRINT " You chose "
IF METH& = 1 : goto "g2040" : ELSE : goto "g2060" : endif
g2040:
PRINT "Adams-Bashforth-Moulton method to solve:"
GOTO "g2250"
g2060:
REM ELSE
IF METH& = 2 : goto "g2080" : ELSE : goto "g2100" : endif
g2080:
PRINT "the Milne-Simpson method to solve"
GOTO "g2250"
g2100:
REM ELSE
IF METH& = 3 :goto "g2120" : ELSE : goto "g2140" : endif
g2120:
PRINT "Hamming's method to solve"
GOTO "g2250"
g2140:
REM ELSE
IF METH& = 4 :goto "g2160" : ELSE : goto "g2180" : endif
g2160:
PRINT "the Runge-Kutta-Fehlberg method to solve"
GOTO "g2250"
g2180:
REM ELSE
g2250:
REM ENDIF
PRINT
s30' Do SUBROUTINE PRINT FUNCTION
cls
PRINT
'WhileLoop 15:I&=&Loop :PRINT:EndWhile
PRINT " You chose ";
IF METH& = 1 :goto "g2340" : ELSE : goto "g2360" : endif
g2340:
PRINT "Adams-Bashforth-Moulton method to solve the D.E."
GOTO "g2550"
g2360:
REM ELSE
IF METH& = 2 :goto "g2380" : ELSE : goto "g2400" : endif
g2380:
PRINT "the Milne-Simpson method to solve the D.E."
GOTO "g2550"
g2400:
REM ELSE
IF METH& = 3 : goto "g2420" : ELSE : goto "g2440" :endif
g2420:
PRINT "Hamming's method to solve the D.E."
GOTO "g2550"
g2440:
REM ELSE
IF METH& = 4 : goto "g2460" : ELSE : goto "g2480" : endif
g2460:
PRINT "the Runge-Kutta-Fehlberg method to solve the D.E."
GOTO "g2550"
g2480:
REM ELSE
g2550:
REM ENDIF
PRINT
s30' do SUBROUTINE PRINT FUNCTION
PRINT
PRINT " with the initial condition Y(A) = Y0."
PRINT
PRINT " A numerical approximation is computed over [A,B]."
PRINT
PRINT " You must enter the endpoints for the interval,"
PRINT
PRINT " the initial condition Y0, and the number of steps M."
PRINT
PRINT
PRINT " Press the <ENTER> key. ";
INPUT ANS$
RETURN
endproc
proc s3000
REM SUBROUTINE END POINTS
g3010:
REM STATUS = (Change,Enter,Done)
STAT&=1 : Case (STATE& = 0) : STAT&=0
WHILE (STAT& = 0) OR (STAT& = 1)
WhileLoop 5 : PRINT : endwhile
s30' Do SUBROUTINE PRINT FUNCTION
PRINT
PRINT
IF STAT& = 1 : goto "g3120" : ELSE : goto "g3180" : endif
g3120:
PRINT " ENTER the left endpoint A = ";
INPUT A!
PRINT
PRINT " ENTER the right endpoint B = ";
INPUT B!
PRINT
g3180:
PRINT " ENTER initial condition Y(A) = ";
INPUT Y0!
PRINT
PRINT " ENTER the number of steps M = ";
VALU!=M&
INPUT VALU!
M&=INT(VALU!)
Case M& < 1 : M&=1
Case (M& > 1000) : M&=1000
PRINT
PRINT " ENTER the tolerance TOL = ";
VALU!=TOL!
INPUT VALU!
TOL!=VALU!
GOTO "g3360"
3280:
REM ELSE
3290:
PRINT " The left endpoint is A = ";A!
PRINT
PRINT " The right endpoint is B = ";B!
PRINT
PRINT " Initial condition is Y(A) = ";Y0!
PRINT
PRINT " The number of steps is M = ";M&
PRINT
PRINT " The tolerance is TOL = ";TOL!
g3360:
REM ENDIF
PRINT
PRINT
PRINT " Do you want to make a change ? <Y/N> ";
INPUT ANS$
IF (ANS$ = "Y") OR (ANS$ = "y"): goto "g3420" : ELSE : goto "g3700" : endif
g3420:
STAT&=0
WhileLoop 15: print: endwhile
PRINT " "
s30' Do SUBROUTINE PRINT FUNCTION
PRINT
PRINT
PRINT " The current left endpoint is A = ";A!
PRINT " ENTER the NEW left endpoint A = ";
INPUT A!
PRINT
PRINT " The current right endpoint is B = ";B!
PRINT " ENTER the NEW right endpoint B = ";
INPUT B!
PRINT
PRINT " The current I. C. is Y(A) = ";Y0!
PRINT " ENTER the NEW I. C. Y(A) = ";
INPUT Y0!
PRINT
PRINT " The current value of M is M = ";M&
PRINT " ENTER the NEW value of M = ";
VALU!=M&
INPUT VALU!
M&=INT(VALU!)
case M& < 1 : M&=1
case M& > 1000 : M&=1000
PRINT
PRINT " The current value of TOL is TOL = ";TOL!
PRINT " ENTER the NEW value of TOL = ";
INPUT VALU!
TOL!=VALU!
GOTO "g3720"
g3700:
REM ELSE
STAT&=2
g3720:
REM ENDIF
EndWhile
RETURN
endproc
proc s4000
REM SUBROUTINE RESULTS
CLS
IF METH& = 1 :goto "g4030" : ELSE : goto "g4050" : endif
g4030:
PRINT "Adams-Bashforth-Moulton method was used to solve the D.E."
GOTO "g4240"
g4050:
REM ELSE
IF METH& = 2 :goto "g4070" : ELSE : goto "g4090" : endif
g4070:
PRINT "The Milne-Simpson method was used to solve the D.E."
GOTO "g4240"
g4090:
REM ELSE
IF METH& = 3 :goto "g4110" : ELSE : goto "g4130" : endif
g4110:
PRINT "Hamming's method was used to solve the D.E."
GOTO "g4240"
g4130:
REM ELSE
IF METH& = 4 : goto "g4150" : ELSE : goto "g4170" : endif
g4150:
PRINT "The Runge-Kutta-Fehlberg method to solve the D.E."
GOTO "g4240"
g4170:
REM ELSE
g4240:
REM ENDIF
s30' Do SUBROUTINE PRINT FUNCTION
PRINT
PRINT "with Y(";T![0];" ) = ";Y![0]
PRINT
PRINT " K"," T(K) "," Y(K) F(T,Y) "
PRINT " -----------------------------------------------------------------"
WhileLoop 0,MEnd&: k&=&Loop
PRINT K&;TAB(15);T![K&];TAB(40);Y![K&]; TAB(65); FNF(T![K&],Y![K&])
if %csrlin>36:waitinput 3000:cls:endif
Endwhile
IF MEND& < M& : goto "g4350" : ELSE : goto "g4370" : endif
g4350:
PRINT
PRINT "The solution points are approaching a pole (POLSTELLE !)"
g4370:
REM CONTINUE
g4380:
RETURN
endproc
g5000:
END
|
|
(939)
Deutsch
English
Français
Español
Italia