c APPLIED NUMERICAL METHODS, EXAMPLE 6.3
c ELECTRICAL TRANSIENTS USING THE RUNGE-KUTTA METHOD
IMPLICIT REAL*8( A-H, O-Z)
INTEGER RUNGE
DIMMENSION F(2), V(2), IMAGE(1500)
c …..READ AND PRINT DATA
1 READ (5,100) R, HL, C, VZERO, H, TMAX, IFREQ, IFPLOT
WRITE (6,200) R, HL, C, VZERO, H, TMAX, IFREQ, IFPLOT
c …..INITIALIZE T, V(1), V(2) AND ICOUNT
T = 0.0
V(1) = VZERO
V(2) = 0.0
ICOUNT = 0
c …..COMPUTE ALPHSQ AND ALPHA, PRINT HEADINGS
ALPHSQ = 1./(C*HL) – R* R /(4.0*HL*HL)
ALPHA = DSQRT (DABX(ALPHSQ))
RO2L = R/(2.*HL)
IF ( IFPLOT .NE. 1 ) GO TO 3
CALL PLOT1( 0, J, 11, 6, 19 )
CALL PLOT2( IMAGE, TMAX, 0., DABS(VZERO), -DABS(VZERO) )
3 WRITE (6,201)
c …..IS CIRCUIT OVER-, CRITICALLY-, OR UNDER-DAMPED
c …..COMPUTE ANALYTICAL SOLUTION, PRINT AND PLOT
4 IF (ALPHSQ) 5, 6, 7
5 TRUV =VZERO*DEXP(-RO2L*T)*((1.+RO2L/ALPHA)*
1 DEXP(ALPHA*T)+(1.-RO2L/ALPFA)*DEXP(-ALPHA*T))/2.0
GO TO 8
6 TRUV =VZERO*DEXP(-RO2L*T)*(1.+RO2L*T)
GO TO 8
7 TRUV=VZERO*DEXP(-RO2L*T)*DCOS(ALPHA*T-DATAN(RO2L/ALP
1 HA))/(ALPHA*DSQRT(C*HL))
8 IF ( IFPLOT .NE. 1 ) GO TO 10
CALL PLOT3( 1H*, T, V(1), 1, 8)
10 WRITE (6,202) T, V(1), TRUV, V(2)
c …..CALL ON THE FOURTH-ORDER RUNGE-KUTTA FUNCTION
11 K = RUNGE(2,V, F, T, H )
c …..WHENEVER K=1, COMPUTE DERIVATIVE VALUES
IF ( K .NE. 1 ) GO TO 13
F(1)=V(2)
F(2)=-R*V(2)/HL –V(1)/(HL*C)
GO TO 11
c …..IF T EXCEEDS TMAX, TERMINATE INTEGRATION
13 IF ( T .LE. TMAX ) GO TO 16
IF ( IFPLOT .NE. 1 ) GO TO 1
WRITE ( 6,203)
CALL PLOT4( 7, 7HVOLTAGE )
WRITE ( 6,204)
GO TO 1
16 ICOUNT=ICOUNT+1
c …..PRINT RESULTS OR CALL DIRECTLY ON RUNGE
IF ( ICOUNT .NE. IFREQ ) GO TO 11
ICOUNT = 0
GO TO 4
c …..FORMATS FOR INPUT AND OUTPUT STATEMENTS
100 FORMAT ( …………..)
200 FORMAT (……………)
201 FORMAT (…全略……)
202 FORMAT (……………)
203 FORMAT (……………)
204 FORMAT (……………)
END
FUNCTION RUNGE( N, Y, F, X, H )
IMPLICIT REAL*8(A-H, O-Z)
REAL*8 Y, F, X, H
INTEGER RUNGE
DIMENSION PHI(50), SAVEY(50), Y(N), F(N)
DATA M/0/
M=M+1
GO TO (1,2,3,4,5), M
c …..PASS 1
1 RUNGE=1
RETURN
c …..PASS 2
2 DO 22 J = 1 , N
SAVEY(J)=Y(J)
PHI(J)=F(J)
22 Y(J)=SAVEY(J)+0.5*H*F(J)
X=X+0.5*H
RUNGE=1
RETURN
c …..PASS 3
3 DO 33 J = 1, N
PHI(J)=PHI(J)+2.0*F(J)
33 Y(J)=SAVEY(J)+0.5*H*F(J)
RUNGE=1
RETURN
c …..PASS 4
4 DO 44 J=1, N
PHI(J)=PHI(J)+2.0*F(J)
44 Y(J)=SAVEY(J)+H*F(J)
X=X+0.5*H
RUNGE=1
RUTURN
c …..PASS 5
5 DO 55 J = 1, N
55 Y(J)=SAVEY(J)+(PHI(J)+F(J))*H/6.0
M=0
RUNGE=0
RETURN
END
改寫完成的 ABC FORTH 完整程式及十群範例輸入數據如下:
\ Runge-Kutta Method for ODE
INTEGER M INTEGER J INTEGER RUNGE
INTEGER N
10 ARRAY Y 10 ARRAY F
REAL X REAL H
50 ARRAY PHI 50 ARRAY SAVEY
: F(RK4) BASIC
10 REM ***F(RUNGE) = FUNCTION RUNGE( N, Y, F, X, H )
20 LET M = M + 1
30 IF M = 1 THEN 60
40 GOTO 80
50 REM ***PASS 1
60 LET RUNGE = 1
70 GOTO 9999
80 IF M = 2 THEN 110
90 GOTO 180
100 REM ***PASS 2
110 FOR J = 1 TO N
120 LET { SAVEY ( J ) = Y ( J ) }
125 LET { PHI ( J ) = F ( J ) }
130 LET { Y ( J ) = SAVEY ( J ) + 0.5 * H * F ( J ) }
140 NEXT J
150 LET { X = X + 0.5 * H }
160 LET RUNGE = 1
170 GOTO 9999
180 IF M = 3 THEN 210
190 GOTO 270
200 REM ***PASS 3
210 FOR J = 1 TO N
220 LET { PHI ( J ) = PHI ( J ) + 2.0 * F ( J ) }
230 LET { Y ( J ) = SAVEY ( J ) + 0.5 * H * F ( J ) }
240 NEXT J
250 LET RUNGE = 1
260 GOTO 9999
270 IF M = 4 THEN 300
280 GOTO 370
290 REM ***PASS 4
300 FOR J = 1 TO N
310 LET { PHI ( J ) = PHI ( J ) + 2.0 * F ( J ) }
320 LET { Y ( J ) = SAVEY ( J ) + H * F ( J ) }
330 NEXT J
340 LET { X = X + 0.5 * H }
350 LET RUNGE = 1
360 GOTO 9999
370 IF M = 5 THEN 400
380 GOTO 430
390 REM ***PASS 5
400 FOR J = 1 TO N
410 LET { Y ( J ) = SAVEY ( J ) + ( PHI ( J ) + F ( J ) ) * H / 6.0 }
420 NEXT J
430 LET M = 0
440 LET RUNGE = 0
9999 END ;
INTEGER K INTEGER ICOUNT常微方程式之 Runge-Kutta 解法
INTEGER IFREQ INTEGER IFPLOT
REAL RT1 REAL RT2 REAL RT3
REAL R REAL HL REAL C REAL VZERO
REAL RO2L REAL ALPHA REAL ALPHSQ
REAL TRUY REAL TMAX
: GROUP1
{{ R = 100.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.00 }}
{{ H = 0.00001 }} {{ TMAX = 0.025 }}
[[ IFREQ = 25 ]] [[ IFPLOT = 0 ]]
;
: GROUP2
{{ R = 0.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.0 }}
{{ H = 0.00010 }} {{ TMAX = 0.020 }}
[[ IFREQ = 10 ]] [[ IFPLOT = 0 ]]
;
: GROUP3
{{ R = 1000.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.0 }}
{{ H = 0.00010 }} {{ TMAX = 0.020 }}
[[ IFREQ = 10 ]] [[ IFPLOT = 0 ]]
;
: GROUP4
{{ R = 1500.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.0 }}
{{ H = 0.00010 }} {{ TMAX = 0.020 }}
[[ IFREQ = 10 ]] [[ IFPLOT = 0 ]]
;
: GROUP5
{{ R = 100.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.0 }}
{{ H = 0.00001 }} {{ TMAX = 0.020 }}
[[ IFREQ = 100 ]] [[ IFPLOT = 0 ]]
;
: GROUP6
{{ R = 100.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.0 }}
{{ H = 0.00010 }} {{ TMAX = 0.020 }}
[[ IFREQ = 10 ]] [[ IFPLOT = 0 ]]
;
: GROUP7
{{ R = 100.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.0 }}
{{ H = 0.00100 }} {{ TMAX = 0.020 }}
[[ IFREQ = 1 ]] [[ IFPLOT = 0 ]]
;
: GROUP8
{{ R = 100.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.0 }}
{{ H = 0.00200 }} {{ TMAX = 0.020 }}
[[ IFREQ = 1 ]] [[ IFPLOT = 0 ]]
;
: GROUP9
{{ R = 100.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.0 }}
{{ H = 0.00500 }} {{ TMAX = 0.020 }}
[[ IFREQ = 1 ]] [[ IFPLOT = 0 ]]
;
: GROUP10
{{ R = 100.00 }} {{ HL = 0.50 }}
{{ C = 2.0E-6 }} {{ VZERO = 10.0 }}
{{ H = 0.01000 }} {{ TMAX = 0.020 }}
[[ IFREQ = 1 ]] [[ IFPLOT = 0 ]]
;
: PRINT-INPUT-DATA BASIC
10 REM ***PRINT DATA
20 PRINT " R = " ; { R }
30 PRINT " HL = " ; { HL }
40 PRINT " C = " ; { C }
50 PRINT " VZERO = " ; { VZERO }
60 PRINT " H = " ; { H }
70 PRINT " TMAX = " ; { TMAX }
80 PRINT " IFREQ = " , IFREQ
90 PRINT " IFPLOT = " , IFPLOT
100 END ;
: INIT BASIC
10 REM ***INITIALIZE X, Y(1), Y(2), AND ICOUNT
20 LET { X = 0 }
30 LET { Y ( 1 ) = VZERO }
40 LET { Y ( 2 ) = 0.0 }
50 LET ICOUNT = 0
60 LET N = 2
70 END ;
: WRITE(6,201)
\ PRINT HEADINGS
CR 10 SPACES ." X " 10 SPACES ." CALC. Y "
10 SPACES ." TRUE Y " 10 SPACES ." CALC. Y' " CR
;
: COMPUTE-CONSTANT BASIC
10 REM ***COMPUTE ALPHSQ AND ALPHA, PRINT HEADINGS
20 LET { ALPHSQ = 1. / ( C * HL ) - R * R / ( 4.0 * HL * HL ) }
30 LET { ALPHA = SQRT ( ABS ( ALPHSQ ) ) }
40 LET { RO2L = R / ( 2. * HL ) }
50 IF IFPLOT <> 1 THEN 80
60 REM ***CALL PLOT1
70 REM ***CALL PLOT2
80 RUN WRITE(6,201)
90 END ;
: WRITE(6,202) BASIC
10 REM ***PRINT SOLUTION
20 PRINT { X , Y ( 1 ) , TRUY , Y ( 2 ) }
30 END ;
: DAMPED BASIC
10 REM ***IS CIRCUIT OVER-, CRITICALLY-, OR UNDER-DAMPED
20 REM ***COMPUTE ANALYTICAL SOLUTION, PRINT OR PLOT
30 IF { ALPHSQ < 0 } THEN 50
40 GOTO 100
50 LET { RT1 = ( 1. + RO2L / ALPHA ) * EXP ( ALPHA * X ) }
60 LET { RT2 = ( 1. - RO2L / ALPHA ) * EXP ( NEGATE ( ALPHA ) * X ) }
70 LET { RT3 = ( RT1 + RT2 ) / 2.0 }
80 LET { TRUY = VZERO * EXP ( NEGATE ( RO2L ) * X ) * RT3 }
90 GOTO 170
100 IF { ALPHSQ = 0 } THEN 120
110 GOTO 140
120 LET { TRUY = VZERO * EXP ( NEGATE ( RO2L ) * X ) * ( 1. + RO2L * X ) }
130 GOTO 170
140 LET { RT1 = COS ( ALPHA * X - ATAN ( RO2L / ALPHA ) ) }
150 LET { RT2 = ALPHA * SQRT ( C * HL ) }
160 LET { TRUY = VZERO * EXP ( NEGATE ( RO2L ) * X ) * RT1 / RT2 }
170 IF IFPLOT <> 1 THEN 190
180 REM ***CALL PLOT3
190 RUN WRITE(6,202)
200 END ;
: CALL-FUNCTION-F(RK4)
F(RK4)
;
: MAIN BASIC \ 執行Runge-Kutta的主要指令
10 RUN GROUP2 \ (1)輸入數據建立了十組,於此處自行選擇更換
20 RUN PRINT-INPUT-DATA
30 RUN INIT
40 RUN COMPUTE-CONSTANT
50 RUN DAMPED
60 REM ***CALL ON THE FOURTH-ORDER RUNGE-KUTTA FUNCTION
70 LET M = 0
80 RUN CALL-FUNCTION-F(RK4)
90 LET K = RUNGE
100 REM ***WHENEVER K=1, COMPUTE DERIVATIVE VALUES
110 IF RUNGE <> 1 THEN 150
120 LET { F ( 1 ) = Y ( 2 ) }
130 LET { F ( 2 ) = NEGATE ( R ) * Y ( 2 ) / HL - Y ( 1 ) / ( HL * C ) }
140 GOTO -80
150 REM *** IF X EXCEEDS TMAX, TEERMINATE INTEGRATION
160 IF { X <= TMAX } THEN 210
170 GOTO 260
\ 180 IF IFPLOT <> 1 THEN -10
\ 190 REM ***WRITE(6,203),CALL PLOT4,WRITE(6,204)
\ 200 GOTO -10
210 LET ICOUNT = ICOUNT + 1
220 REM ***PRINT RESULTS OR CALL DIRECTLY ON RUNGE
230 IF ICOUNT <> IFREQ THEN -70
240 LET ICOUNT = 0
250 GOTO -50
260 PRINT " ***Run next GROUPn by change ( 10 RUN GROUPn )*** "
270 END ;
直接操作ABC FORTH系統,執行MAIN指令,所得結果顯示如下:
MAIN
R = 0.000000000E-1
HL = .5000000000
C = 2.000000000E-6
VZERO = 10.00000000
H = .0001000000
TMAX = .0200000000
IFREQ = 10
IFPLOT = 0
X CALC. Y TRUE Y CALC. Y'
0.000000000E-1 10.00000000 10.00000000 0.000000000E-1
.0010000000 5.403029671 5.403023059 -8414.704778
.0020000000 -4.161452687 -4.161468365 -9092.979918
.0030000000 -9.899919391 -9.899924966 -1411.224448
.0040000000 -6.536459532 -6.536436209 7568.001143
.0050000000 2.836581058 2.836621855 9589.251198
.0060000000 9.601684950 9.601702867 2794.201656
.0070000000 7.539057071 7.539022543 -6569.818977
.0080000000 -1.454933809 -1.455000338 -9893.586642
.0090000000 -9.111266133 -9.111302619 -4121.250371
.0100000000 -8.390754644 -8.390715291 5440.137662
.0110000000 .0441656076 .0442569799 9999.894840
.0120000000 8.438479098 8.438539587 5365.808798
.0130000000 9.074504988 9.074467815 -4201.568624
.0140000000 1.367486013 1.367372182 -9906.048042
.0150000000 -7.596790229 -7.596879129 -6502.966258
.0160000000 -9.576622425 -9.576594803 2878.902739
.0170000000 -2.751765848 -2.751633381 9613.924740
.0180000000 6.603046592 6.603167082 7509.961785
.0190000000 9.887056797 9.887046182 -1498.614135
***Run next GROUPn by change ( 10 RUN GROUPn )*** ok
執行後所得結果,與原始資料內容比較,毫不遜色。