📄 simulation.f
字号:
PHASE1 = .TRUE. STEP 177
NORND = .TRUE. STEP 178
IF(P5EPS .GT. 100D0*ROUND) GO TO 99 STEP 179
NORND = .FALSE. STEP 180
DO 25 L = 1,NEQN STEP 181
25 PHI(L,15) = 0.0D0 STEP 182
99 IFAIL = 0 STEP 183
C *** END BLOCK 0 *** STEP 184
C STEP 185
C *** BEGIN BLOCK 1 *** STEP 186
C COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP. AVOID COMPUTING STEP 187
C THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED. STEP 188
C *** STEP 189
C STEP 190
100 KP1 = K+1 STEP 191
KP2 = K+2 STEP 192
KM1 = K-1 STEP 193
KM2 = K-2 STEP 194
C STEP 195
C NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT STEP 196
C ONE. WHEN K.LT.NS, NO COEFFICIENTS CHANGE STEP 197
C STEP 198
IF(H .NE. HOLD) NS = 0 STEP 199
IF(NS .LE. KOLD) NS=NS+1 STEP 200
NSP1 = NS+1 STEP 201
IF (K .LT. NS) GO TO 199 STEP 202
C STEP 203
C COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH STEP 204
C ARE CHANGED STEP 205
C STEP 206
BETA(NS) = 1.0D0 STEP 207
REALNS = NS STEP 208
ALPHA(NS) =1.D0/REALNS STEP 209
TEMP1 = H*REALNS STEP 210
SIG(NSP1) = 1.0D0 STEP 211
IF(K .LT. NSP1) GO TO 110 STEP 212
DO 105 I = NSP1,K STEP 213
IM1 = I-1 STEP 214
TEMP2 = PSI(IM1) STEP 215
PSI(IM1) = TEMP1 STEP 216
BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2 STEP 217
TEMP1 = TEMP2 + H STEP 218
ALPHA(I) = H/TEMP1 STEP 219
REALI = I STEP 220
105 SIG(I+1) = REALI*ALPHA(I)*SIG(I) STEP 221
110 PSI(K) = TEMP1 STEP 222
C STEP 223
C COMPUTE COEFFICIENTS G(*) STEP 224
C STEP 225
C INITIALIZE V(*) AND SET W(*). G(2) IS SET IN DATA STATEMENT STEP 226
C STEP 227
IF(NS .GT. 1) GO TO 120 STEP 228
DO 115 IQ = 1,K STEP 229
TEMP3 = IQ*(IQ+1) STEP 230
V(IQ) =1.D0/TEMP3 STEP 231
115 W(IQ) = V(IQ) STEP 232
GO TO 140 STEP 233
C STEP 234
C IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*) STEP 235
C STEP 236
120 IF(K .LE. KOLD) GO TO 130 STEP 237
TEMP4 = K*KP1 STEP 238
V(K) =1.D0/TEMP4 STEP 239
NSM2 = NS-2 STEP 240
IF(NSM2 .LT. 1) GO TO 130 STEP 241
DO 125 J = 1,NSM2 STEP 242
I = K-J STEP 243
125 V(I) = V(I) - ALPHA(J+1)*V(I+1) STEP 244
C STEP 245
C UPDATE V(*) AND SET W(*) STEP 246
C STEP 247
130 LIMIT1 = KP1 - NS STEP 248
TEMP5 = ALPHA(NS) STEP 249
DO 135 IQ = 1,LIMIT1 STEP 250
V(IQ) = V(IQ) - TEMP5*V(IQ+1) STEP 251
135 W(IQ) = V(IQ) STEP 252
G(NSP1) = W(1) STEP 253
C STEP 254
C COMPUTE THE G(*) IN THE WORK VECTOR W(*) STEP 255
C STEP 256
140 NSP2 = NS + 2 STEP 257
IF(KP1 .LT. NSP2) GO TO 199 STEP 258
DO 150 I = NSP2,KP1 STEP 259
LIMIT2 = KP2 - I STEP 260
TEMP6 = ALPHA(I-1) STEP 261
DO 145 IQ = 1,LIMIT2 STEP 262
145 W(IQ) = W(IQ) - TEMP6*W(IQ+1) STEP 263
150 G(I) = W(1) STEP 264
199 CONTINUE STEP 265
C *** END BLOCK 1 *** STEP 266
C STEP 267
C *** BEGIN BLOCK 2 *** STEP 268
C PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED STEP 269
C SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K, STEP 270
C K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED. STEP 271
C *** STEP 272
C STEP 273
C CHANGE PHI TO PHI STAR STEP 274
C STEP 275
IF(K .LT. NSP1) GO TO 215 STEP 276
DO 210 I = NSP1,K STEP 277
TEMP1 = BETA(I) STEP 278
DO 205 L = 1,NEQN STEP 279
205 PHI(L,I) = TEMP1*PHI(L,I) STEP 280
210 CONTINUE STEP 281
C STEP 282
C PREDICT SOLUTION AND DIFFERENCES STEP 283
C STEP 284
215 DO 220 L = 1,NEQN STEP 285
PHI(L,KP2) = PHI(L,KP1) STEP 286
PHI(L,KP1) = 0.0D0 STEP 287
220 P(L) = 0.0D0 STEP 288
DO 230 J = 1,K STEP 289
I = KP1 - J STEP 290
IP1 = I+1 STEP 291
TEMP2 = G(I) STEP 292
DO 225 L = 1,NEQN STEP 293
P(L) = P(L) + TEMP2*PHI(L,I) STEP 294
225 PHI(L,I) = PHI(L,I) + PHI(L,IP1) STEP 295
230 CONTINUE STEP 296
IF(NORND) GO TO 240 STEP 297
DO 235 L = 1,NEQN STEP 298
TAU = H*P(L) - PHI(L,15) STEP 299
P(L) = Y(L) + TAU STEP 300
235 PHI(L,16) = (P(L) - Y(L)) - TAU STEP 301
GO TO 250 STEP 302
240 DO 245 L = 1,NEQN STEP 303
245 P(L) = Y(L) + H*P(L) STEP 304
250 XOLD = X STEP 305
X = X + H STEP 306
ABSH =DABS(H) STEP 307
IPRDCT=IPRDCT+1
CALL F(X,P,YP) STEP 308
C STEP 309
C ESTIMATE ERRORS AT ORDERS K,K-1,K-2 STEP 310
C STEP 311
ERKM2 = 0.0D0 STEP 312
ERKM1 = 0.0D0 STEP 313
ERK = 0.0D0 STEP 314
DO 265 L = 1,NEQN STEP 315
TEMP3 =1.D0/WT(L) STEP 316
TEMP4 = YP(L) - PHI(L,1) STEP 317
IF(KM2)265,260,255 STEP 318
255 ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2 STEP 319
260 ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2 STEP 320
265 ERK = ERK + (TEMP4*TEMP3)**2 STEP 321
IF(KM2)280,275,270 STEP 322
270 ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*DSQRT(ERKM2) STEP 323
275 ERKM1 = ABSH*SIG(K)*GSTR(KM1)*DSQRT(ERKM1) STEP 324
280 TEMP5 = ABSH*DSQRT(ERK) STEP 325
ERR = TEMP5*(G(K)-G(KP1)) STEP 326
ERK = TEMP5*SIG(KP1)*GSTR(K) STEP 327
KNEW = K STEP 328
C STEP 329
C TEST IF ORDER SHOULD BE LOWERED STEP 330
C STEP 331
IF(KM2)299,290,285 STEP 332
285 IF(DMAX1(ERKM1,ERKM2) .LE. ERK) KNEW = KM1 STEP 333
GO TO 299 STEP 334
290 IF(ERKM1 .LE. .5D0*ERK)KNEW = KM1 STEP 335
C STEP 336
C TEST IF STEP SUCCESSFUL STEP 337
C STEP 338
299 IF(ERR .LE. EPS) GO TO 400 STEP 339
C *** END BLOCK 2 *** STEP 340
C STEP 341
C *** BEGIN BLOCK 3 *** STEP 342
C THE STEP IS UNSUCCESSFUL. RESTORE X, PHI(*,*), PSI(*) . STEP 343
C IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE. IF STEP FAILS MORE STEP 344
C THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE. DOUBLE ERROR STEP 345
C TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINESTEP 346
C PRECISION. STEP 347
C *** STEP 348
C STEP 349
C RESTORE X, PHI(*,*) AND PSI(*) STEP 350
C STEP 351
PHASE1 = .FALSE. STEP 352
X = XOLD STEP 353
DO 310 I = 1,K STEP 354
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -