📄 simulation.f
字号:
C RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES DE 50
C IFLAG -- +1,-1. INDICATOR TO INITIALIZE THE CODE. NORMAL INPUT DE 51
C IS +1. THE USER SHOULD SET IFLAG=-1 ONLY IF IT IS DE 52
C IMPOSSIBLE TO CONTINUE THE INTEGRATION BEYOND TOUT . DE 53
C ALL PARAMETERS EXCEPT F , NEQN AND TOUT MAY BE ALTERED BY THE DE 54
C CODE ON OUTPUT SO MUST BE VARIABLES IN THE CALLING PROGRAM. DE 55
C DE 56
C OUTPUT FROM DE -- DE 57
C DE 58
C NEQN -- UNCHANGED DE 59
C Y(*) -- SOLUTION AT T DE 60
C T -- LAST POINT REACHED IN INTEGRATION. NORMAL RETURN HAS DE 61
C T = TOUT . DE 62
C TOUT -- UNCHANGED DE 63
C RELERR,ABSERR -- NORMAL RETURN HAS TOLERANCES UNCHANGED. IFLAG=3DE 64
C SIGNALS TOLERANCES INCREASED DE 65
C IFLAG = 2 -- NORMAL RETURN. INTEGRATION REACHED TOUT DE 66
C = 3 -- INTEGRATION DID NOT REACH TOUT BECAUSE ERROR DE 67
C TOLERANCES TOO SMALL. RELERR , ABSERR INCREASED DE 68
C APPROPRIATELY FOR CONTINUING DE 69
C = 4 -- INTEGRATION DID NOT REACH TOUT BECAUSE MORE THAN DE 70
C MAXNUM STEPS NEEDED DE 71
C = 5 -- INTEGRATION DID NOT REACH TOUT BECAUSE EQUATIONS DE 72
C APPEAR TO BE STIFF DE 73
C = 6 -- INVALID INPUT PARAMETERS (FATAL ERROR) DE 74
C THE VALUE OF IFLAG IS RETURNED NEGATIVE WHEN THE INPUT DE 75
C VALUE IS NEGATIVE AND THE INTEGRATION DOES NOT REACH TOUT ,DE 76
C I.E., -3, -4, -5. DE 77
C DE 78
C SUBSEQUENT CALLS TO DE -- DE 79
C DE 80
C SUBROUTINE DE RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE DE 81
C THE INTEGRATION. IF THE INTEGRATION REACHED TOUT , THE USER NEED DE 82
C ONLY DEFINE A NEW TOUT AND CALL AGAIN. IF THE INTEGRATION DID NOTDE 83
C REACH TOUT AND THE USER WANTS TO CONTINUE, HE JUST CALLS AGAIN. DE 84
C THE OUTPUT VALUE OF IFLAG IS THE APPROPRIATE INPUT VALUE FOR DE 85
C SUBSEQUENT CALLS. THE ONLY SITUATION IN WHICH IT SHOULD BE ALTERED DE 86
C IS TO STOP THE INTEGRATION INTERNALLY AT THE NEW TOUT , I.E., DE 87
C CHANGE OUTPUT IFLAG=2 TO INPUT IFLAG=-2 . ERROR TOLERANCES MAY DE 88
C BE CHANGED BY THE USER BEFORE CONTINUING. ALL OTHER PARAMETERS MUSTDE 89
C REMAIN UNCHANGED. DE 90
C DE 91
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
LOGICAL START,CRASH,STIFF DE 92
DIMENSION Y(300),PSI(12) DE 93
DIMENSION YY(300),WT(300),PHI(300,16),P(300),YP(300),YPOUT(300) DE 94
EXTERNAL F DE 95
SAVE
COMMON/IPRDCT/IPRDCT
C DE 96
C***********************************************************************DE 97
C* THE ONLY MACHINE DEPENDENT CONSTANT IS BASED ON THE MACHINE UNIT *DE 98
C* ROUNDOFF ERROR U WHICH IS THE SMALLEST POSITIVE NUMBER SUCH THAT *DE 99
C* 1.0+U .GT. 1.0 . U MUST BE CALCULATED AND FOURU=4.0*U INSERTED *DE 100
C* IN THE FOLLOWING DATA STATEMENT BEFORE USING DE . THE ROUTINE *DE 101
C* MACHIN CALCULATES U . FOURU AND TWOU=2.0*U MUST ALSO BE *DE 102
C* INSERTED IN SUBROUTINE STEP BEFORE CALLING DE . *DE 103
DATA FOURU/8.88D-16/ DE 104
C***********************************************************************DE 105
C DE 106
C THE CONSTANT MAXNUM IS THE MAXIMUM NUMBER OF STEPS ALLOWED IN ONE DE 107
C CALL TO DE . THE USER MAY CHANGE THIS LIMIT BY ALTERING THE DE 108
C FOLLOWING STATEMENT DE 109
DATA MAXNUM/500/ DE 110
C THIS VERSION OF DE WILL HANDLE UP TO 20 EQUATIONS. TO CHANGE THISDE 111
C NUMBER, ONLY THE NUMBER 20 IN THE DIMENSION STATEMENT ABOVE AND THE DE 112
C NUMBER 20 IN THE FOLLOWING STATEMENT NEED BE CHANGED DE 113
C DE 114
C *** *** *** DE 115
C TEST FOR IMPROPER PARAMETERS DE 116
C DE 117
IF(NEQN .LT. 1 .OR. NEQN .GT. 300) GO TO 10 DE 118
IF(T .EQ. TOUT) GO TO 10 DE 119
IF(RELERR .LT. 0.D0 .OR. ABSERR .LT. 0.D0)GO TO 10 DE 120
EPS = DMAX1(RELERR,ABSERR) DE 121
IF(EPS .LE. 0.D0)GO TO 10 DE 122
IF(IFLAG .EQ. 0) GO TO 10 DE 123
ISN = ISIGN(1,IFLAG) DE 124
IFLAG = IABS(IFLAG) DE 125
IF(IFLAG .EQ. 1) GO TO 20 DE 126
IF(T .NE. TOLD) GO TO 10 DE 127
IF(IFLAG .GE. 2 .AND. IFLAG .LE. 5) GO TO 20 DE 128
10 IFLAG = 6 DE 129
RETURN DE 130
C DE 131
C ON EACH CALL SET INTERVAL OF INTEGRATION AND COUNTER FOR NUMBER OF DE 132
C STEPS. ADJUST INPUT ERROR TOLERANCES TO DEFINE WEIGHT VECTOR FOR DE 133
C SUBROUTINE STEP DE 134
C DE 135
20 DEL = TOUT - T DE 136
ABSDEL =DABS(DEL) DE 137
TEND = T + 10.0D0*DEL DE 138
IF(ISN .LT. 0) TEND = TOUT DE 139
NOSTEP = 0 DE 140
KLE4 = 0 DE 141
STIFF = .FALSE. DE 142
RELEPS = RELERR/EPS DE 143
ABSEPS = ABSERR/EPS DE 144
IF(IFLAG .EQ. 1) GO TO 30 DE 145
IF(ISNOLD .LT. 0) GO TO 30 DE 146
IF(DELSGN*DEL .GT. 0.D0)GO TO 50 DE 147
C DE 148
C ON START AND RESTART ALSO SET WORK VARIABLES X AND YY(*), STORE THE DE 149
C DIRECTION OF INTEGRATION AND INITIALIZE THE STEP SIZE DE 150
C DE 151
30 START = .TRUE. DE 152
X = T DE 153
DO 40 L = 1,NEQN DE 154
40 YY(L) = Y(L) DE 155
DELSGN =DSIGN(1.D0,DEL) DE 156
H=DSIGN(DMAX1(DABS(TOUT-X),FOURU*DABS(X)),TOUT-X) DE 157
C DE 158
C IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN DE 159
C DE 160
50 IF(DABS(X-T).LT. ABSDEL) GO TO 60 DE 161
CALL INTRP2(X,YY,TOUT,Y,YPOUT,NEQN,KOLD,PHI,PSI) DE 162
IFLAG = 2 DE 163
T = TOUT DE 164
TOLD = T DE 165
ISNOLD = ISN DE 166
RETURN DE 167
C DE 168
C IF CANNOT GO PAST OUTPUT POINT AND SUFFICIENTLY CLOSE, DE 169
C EXTRAPOLATE AND RETURN DE 170
C DE 171
60 IF(ISN .GT. 0 .OR. DABS(TOUT-X) .GE. FOURU*DABS(X))GO TO 80 DE 172
H = TOUT - X DE 173
IPRDCT=0
CALL F(X,YY,YP) DE 174
DO 70 L = 1,NEQN DE 175
70 Y(L) = YY(L) + H*YP(L) DE 176
IFLAG = 2 DE 177
T = TOUT DE 178
TOLD = T DE 179
ISNOLD = ISN DE 180
RETURN DE 181
C DE 182
C TEST FOR TOO MUCH WORK DE 183
C DE 184
80 IF(NOSTEP .LT. MAXNUM) GO TO 100 DE 185
IFLAG = ISN*4 DE 186
IF(STIFF) IFLAG = ISN*5 DE 187
DO 90 L = 1,NEQN DE 188
90 Y(L) = YY(L) DE 189
T = X DE 190
TOLD = T DE 191
ISNOLD = 1 DE 192
RETURN DE 193
C DE 194
C LIMIT STEP SIZE, SET WEIGHT VECTOR AND TAKE A STEP DE 195
C DE 196
100 H =DSIGN(DMIN1(DABS(H),DABS(TEND-X)),H) DE 197
DO 110 L = 1,NEQN DE 198
110 WT(L) =RELEPS*DABS(YY(L)) + ABSEPS DE 199
CALL STEPOCP(X,YY,F,NEQN,H,EPS,WT,START, DE 200
1 HOLD,K,KOLD,CRASH,PHI,P,YP,PSI) DE 201
C DE 202
C TEST FOR TOLERANCES TOO SMALL DE 203
C DE 204
IF(.NOT.CRASH) GO TO 130 DE 205
IFLAG = ISN*3 DE 206
RELERR = EPS*RELEPS DE 207
ABSERR = EPS*ABSEPS DE 208
DO 120 L = 1,NEQN DE 209
120 Y(L) = YY(L) DE 210
T = X DE 211
TOLD = T DE 212
ISNOLD = 1 DE 213
RETURN DE 214
C DE 215
C AUGMENT COUNTER ON WORK AND TEST FOR STIFFNESS DE 216
C DE 217
130 NOSTEP = NOSTEP + 1 DE 218
KLE4 = KLE4 + 1 DE 219
IF(KOLD .GT. 4) KLE4 = 0 DE 220
IF(KLE4 .GE. 50) STIFF = .TRUE. DE 221
GO TO 50 DE 222
END DE 223
SUBROUTINE STEPOCP(X,Y,F,NEQN,H,EPS,WT,START, STEP 1
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -