📄 simulation.f
字号:
1 HOLD,K,KOLD,CRASH,PHI,P,YP,PSI) STEP 2
C STEP 3
C SUBROUTINE STEP INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY STEP 4
C DIFFERENTIAL EQUATIONS ONE STEP, NORMALLY FROM X TO X+H, USING A STEP 5
C MODIFIED DIVIDED DIFFERENCE FORM OF THE ADAMS PECE FORMULAS. LOCAL STEP 6
C EXTRAPOLATION IS USED TO IMPROVE ABSOLUTE STABILITY AND ACCURACY. STEP 7
C THE CODE ADJUSTS ITS ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR STEP 8
c per unit step in a generalized sense. special devices are included step 9
C TO CONTROL ROUNDOFF ERROR AND TO DETECT WHEN THE USER IS REQUESTING STEP 10
C TOO MUCH ACCURACY. STEP 11
C STEP 12
C THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, STEP 13
C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS: THE INITIAL STEP 14
C VALUE PROBLEM BY L. F. SHAMPINE AND M. K. GORDON. STEP 15
C STEP 16
C STEP 17
C THE PARAMETERS REPRESENT: STEP 18
C X -- INDEPENDENT VARIABLE STEP 19
C Y(*) -- SOLUTION VECTOR AT X STEP 20
C YP(*) -- DERIVATIVE OF SOLUTION VECTOR AT X AFTER SUCCESSFUL STEP 21
C STEP STEP 22
C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED STEP 23
C H -- APPROPRIATE STEP SIZE FOR NEXT STEP. NORMALLY DETERMINED BYSTEP 24
C CODE STEP 25
C EPS -- LOCAL ERROR TOLERANCE. MUST BE VARIABLE STEP 26
C WT(*) -- VECTOR OF WEIGHTS FOR ERROR CRITERION STEP 27
C START -- LOGICAL VARIABLE SET .TRUE. FOR FIRST STEP, .FALSE. STEP 28
C OTHERWISE STEP 29
C HOLD -- STEP SIZE USED FOR LAST SUCCESSFUL STEP STEP 30
C K -- APPROPRIATE ORDER FOR NEXT STEP (DETERMINED BY CODE) STEP 31
C KOLD -- ORDER USED FOR LAST SUCCESSFUL STEP STEP 32
C CRASH -- LOGICAL VARIABLE SET .TRUE. WHEN NO STEP CAN BE TAKEN, STEP 33
C .FALSE. OTHERWISE. STEP 34
C THE ARRAYS PHI, PSI ARE REQUIRED FOR THE INTERPOLATION SUBROUTINE STEP 35
C INTRP . THE ARRAY P IS INTERNAL TO THE CODE. STEP 36
C STEP 37
C INPUT TO STEP STEP 38
C STEP 39
C FIRST CALL -- STEP 40
C STEP 41
C THE USER MUST PROVIDE STORAGE IN HIS DRIVER PROGRAM FOR ALL ARRAYS STEP 42
C IN THE CALL LIST, NAMELY STEP 43
C STEP 44
C DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12) STEP 45
C STEP 46
C THE USER MUST ALSO DECLARE START AND CRASH LOGICAL VARIABLES STEP 47
C AND F AN EXTERNAL SUBROUTINE, SUPPLY THE SUBROUTINE F(X,Y,YP) STEP 48
C TO EVALUATE STEP 49
C DY(I)/DX = YP(I) = F(X,Y(1),Y(2),...,Y(NEQN)) STEP 50
C AND INITIALIZE ONLY THE FOLLOWING PARAMETERS: STEP 51
C X -- INITIAL VALUE OF THE INDEPENDENT VARIABLE STEP 52
C Y(*) -- VECTOR OF INITIAL VALUES OF DEPENDENT VARIABLES STEP 53
C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED STEP 54
C H -- NOMINAL STEP SIZE INDICATING DIRECTION OF INTEGRATION STEP 55
C AND MAXIMUM SIZE OF STEP. MUST BE VARIABLE STEP 56
C EPS -- LOCAL ERROR TOLERANCE PER STEP. MUST BE VARIABLE STEP 57
C WT(*) -- VECTOR OF NON-ZERO WEIGHTS FOR ERROR CRITERION STEP 58
C START -- .TRUE. STEP 59
C STEP 60
C STEP REQUIRES THE L2 NORM OF THE VECTOR WITH COMPONENTS STEP 61
C LOCAL ERROR(L)/WT(L) BE LESS THAN EPS FOR A SUCCESSFUL STEP. THESTEP 62
C ARRAY WT ALLOWS THE USER TO SPECIFY AN ERROR TEST APPROPRIATE STEP 63
C FOR HIS PROBLEM. FOR EXAMPLE, STEP 64
C WT(L) = 1.0 SPECIFIES ABSOLUTE ERROR, STEP 65
C = ABS(Y(L)) ERROR RELATIVE TO THE MOST RECENT VALUE OF THESTEP 66
C L-TH COMPONENT OF THE SOLUTION, STEP 67
C = ABS(YP(L)) ERROR RELATIVE TO THE MOST RECENT VALUE OF STEP 68
C THE L-TH COMPONENT OF THE DERIVATIVE, STEP 69
C = AMAX1(WT(L),ABS(Y(L))) ERROR RELATIVE TO THE LARGEST STEP 70
C MAGNITUDE OF L-TH COMPONENT OBTAINED SO FAR, STEP 71
C = ABS(Y(L))*RELERR/EPS + ABSERR/EPS SPECIFIES A MIXED STEP 72
C RELATIVE-ABSOLUTE TEST WHERE RELERR IS RELATIVE STEP 73
C ERROR, ABSERR IS ABSOLUTE ERROR AND EPS = STEP 74
C AMAX1(RELERR,ABSERR) . STEP 75
C STEP 76
C SUBSEQUENT CALLS -- STEP 77
C STEP 78
C SUBROUTINE STEP IS DESIGNED SO THAT ALL INFORMATION NEEDED TO STEP 79
C CONTINUE THE INTEGRATION, INCLUDING THE STEP SIZE H AND THE ORDER STEP 80
C K , IS RETURNED WITH EACH STEP. WITH THE EXCEPTION OF THE STEP STEP 81
C SIZE, THE ERROR TOLERANCE, AND THE WEIGHTS, NONE OF THE PARAMETERS STEP 82
C SHOULD BE ALTERED. THE ARRAY WT MUST BE UPDATED AFTER EACH STEP STEP 83
C TO MAINTAIN RELATIVE ERROR TESTS LIKE THOSE ABOVE. NORMALLY THE STEP 84
C INTEGRATION IS CONTINUED JUST BEYOND THE DESIRED ENDPOINT AND THE STEP 85
C SOLUTION INTERPOLATED THERE WITH SUBROUTINE INTRP . IF IT IS STEP 86
C IMPOSSIBLE TO INTEGRATE BEYOND THE ENDPOINT, THE STEP SIZE MAY BE STEP 87
C REDUCED TO HIT THE ENDPOINT SINCE THE CODE WILL NOT TAKE A STEP STEP 88
C LARGER THAN THE H INPUT. CHANGING THE DIRECTION OF INTEGRATION, STEP 89
C I.E., THE SIGN OF H , REQUIRES THE USER SET START = .TRUE. BEFORE STEP 90
C CALLING STEP AGAIN. THIS IS THE ONLY SITUATION IN WHICH START STEP 91
C SHOULD BE ALTERED. STEP 92
C STEP 93
C OUTPUT FROM STEP STEP 94
C STEP 95
C SUCCESSFUL STEP -- STEP 96
C STEP 97
C THE SUBROUTINE RETURNS AFTER EACH SUCCESSFUL STEP WITH START AND STEP 98
C CRASH SET .FALSE. . X REPRESENTS THE INDEPENDENT VARIABLE STEP 99
C ADVANCED ONE STEP OF LENGTH HOLD FROM ITS VALUE ON INPUT AND Y STEP 100
C THE SOLUTION VECTOR AT THE NEW VALUE OF X . ALL OTHER PARAMETERS STEP 101
C REPRESENT INFORMATION CORRESPONDING TO THE NEW X NEEDED TO STEP 102
C CONTINUE THE INTEGRATION. STEP 103
C STEP 104
C UNSUCCESSFUL STEP -- STEP 105
C STEP 106
C WHEN THE ERROR TOLERANCE IS TOO SMALL FOR THE MACHINE PRECISION, STEP 107
C THE SUBROUTINE RETURNS WITHOUT TAKING A STEP AND CRASH = .TRUE. . STEP 108
C AN APPROPRIATE STEP SIZE AND ERROR TOLERANCE FOR CONTINUING ARE STEP 109
C ESTIMATED AND ALL OTHER INFORMATION IS RESTORED AS UPON INPUT STEP 110
C BEFORE RETURNING. TO CONTINUE WITH THE LARGER TOLERANCE, THE USER STEP 111
C JUST CALLS THE CODE AGAIN. A RESTART IS NEITHER REQUIRED NOR STEP 112
C DESIRABLE. STEP 113
C STEP 114
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
LOGICAL START,CRASH,PHASE1,NORND STEP 115
DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12) STEP 116
DIMENSION ALPHA(12),BETA(12),SIG(13),W(12),V(12),G(13), STEP 117
1 GSTR(13),TWO(13) STEP 118
COMMON/IPRDCT/IPRDCT
EXTERNAL F STEP 119
SAVE
C***********************************************************************STEP 120
C* THE ONLY MACHINE DEPENDENT CONSTANTS ARE BASED ON THE MACHINE UNIT *STEP 121
C* ROUNDOFF ERROR U WHICH IS THE SMALLEST POSITIVE NUMBER SUCH THAT *STEP 122
C* 1.0+U .GT. 1.0 . THE USER MUST CALCULATE U AND INSERT *STEP 123
C* TWOU=2.0*U AND FOURU=4.0*U IN THE DATA STATEMENT BEFORE CALLING *STEP 124
C* THE CODE. THE ROUTINE MACHIN CALCULATES U . *STEP 125
DATA TWOU,FOURU/2.8D-14,5.68D-14/ STEP 126
C***********************************************************************STEP 127
DATA TWO/2D0,4D0,8D0,16D0,32D0,64D0,128D0,256D0,512D0,1024D0, STEP 128
1 2048D0,4096D0,8192D0/ STEP 129
DATA GSTR/.5D0,.0833D0,.0417D0,.0264D0,.0188D0,.0143D0,.0114D0, STEP 130
1 .936D-2,.789D-2,.679D-2,.592D-2,.524D-2,.468D-2/ STEP 131
DATA G(1),G(2)/1.0D0,0.5D0/,SIG(1)/1.0D0/ STEP 132
C STEP 133
C STEP 134
C *** BEGIN BLOCK 0 *** STEP 135
C CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE STEP 136
C PRECISION. IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A STEP 137
C STARTING STEP SIZE. STEP 138
C *** STEP 139
C STEP 140
C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE STEP 141
C STEP 142
CRASH = .TRUE. STEP 143
IF(DABS(H) .GE. FOURU*DABS(X)) GO TO 5 STEP 144
H =DSIGN(FOURU*DABS(X),H) STEP 145
RETURN STEP 146
5 P5EPS = 0.5D0*EPS STEP 147
C STEP 148
C IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE STEP 149
C STEP 150
ROUND = 0.0D0 STEP 151
DO 10 L = 1,NEQN STEP 152
10 ROUND = ROUND + (Y(L)/WT(L))**2 STEP 153
ROUND =TWOU*DSQRT(ROUND) STEP 154
IF(P5EPS .GE. ROUND) GO TO 15 STEP 155
EPS = 2D0*ROUND*(1D0 + FOURU) STEP 156
RETURN STEP 157
15 CRASH = .FALSE. STEP 158
IF(.NOT.START) GO TO 99 STEP 159
C STEP 160
C INITIALIZE. COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP STEP 161
C STEP 162
IPRDCT=0
CALL F(X,Y,YP) STEP 163
SUM = 0.0D0 STEP 164
DO 20 L = 1,NEQN STEP 165
PHI(L,1) = YP(L) STEP 166
PHI(L,2) = 0.0D0 STEP 167
20 SUM = SUM + (YP(L)/WT(L))**2 STEP 168
SUM =DSQRT(SUM) STEP 169
ABSH =DABS(H) STEP 170
IF(EPS .LT. 16D0*SUM*H*H) ABSH = 0.25D0*DSQRT(EPS/SUM) STEP 171
H =DSIGN(DMAX1(ABSH,FOURU*DABS(X)),H) STEP 172
HOLD = 0.0D0 STEP 173
K = 1 STEP 174
KOLD = 0 STEP 175
START = .FALSE. STEP 176
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -