2
SUBROUTINE DSTEPS (DF, NEQN, Y, X, H, EPS, WT, START, HOLD, K,
3
+ KOLD, CRASH, PHI, P, YP, PSI, ALPHA, BETA, SIG, V, W, G,
4
+ PHASE1, NS, NORND, KSTEPS, TWOU, FOURU, XOLD, KPREV, IVC, IV,
6
C***BEGIN PROLOGUE DSTEPS
7
C***PURPOSE Integrate a system of first order ordinary differential
9
C***LIBRARY SLATEC (DEPAC)
11
C***TYPE DOUBLE PRECISION (STEPS-S, DSTEPS-D)
12
C***KEYWORDS ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
13
C ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR
14
C***AUTHOR Shampine, L. F., (SNLA)
15
C Gordon, M. K., (SNLA)
16
C MODIFIED BY H.A. WATTS
19
C Written by L. F. Shampine and M. K. Gordon
23
C Subroutine DSTEPS is normally used indirectly through subroutine
24
C DDEABM . Because DDEABM suffices for most problems and is much
25
C easier to use, using it should be considered before using DSTEPS
28
C Subroutine DSTEPS integrates a system of NEQN first order ordinary
29
C differential equations one step, normally from X to X+H, using a
30
C modified divided difference form of the Adams Pece formulas. Local
31
C extrapolation is used to improve absolute stability and accuracy.
32
C The code adjusts its order and step size to control the local error
33
C per unit step in a generalized sense. Special devices are included
34
C to control roundoff error and to detect when the user is requesting
37
C This code is completely explained and documented in the text,
38
C Computer Solution of Ordinary Differential Equations, The Initial
39
C Value Problem by L. F. Shampine and M. K. Gordon.
40
C Further details on use of this code are available in "Solving
41
C Ordinary Differential Equations with ODE, STEP, and INTRP",
42
C by L. F. Shampine and M. K. Gordon, SLA-73-1060.
45
C The parameters represent --
46
C DF -- subroutine to evaluate derivatives
47
C NEQN -- number of equations to be integrated
48
C Y(*) -- solution vector at X
49
C X -- independent variable
50
C H -- appropriate step size for next step. Normally determined by
52
C EPS -- local error tolerance
53
C WT(*) -- vector of weights for error criterion
54
C START -- logical variable set .TRUE. for first step, .FALSE.
56
C HOLD -- step size used for last successful step
57
C K -- appropriate order for next step (determined by code)
58
C KOLD -- order used for last successful step
59
C CRASH -- logical variable set .TRUE. when no step can be taken,
61
C YP(*) -- derivative of solution vector at X after successful
63
C KSTEPS -- counter on attempted steps
64
C TWOU -- 2.*U where U is machine unit roundoff quantity
65
C FOURU -- 4.*U where U is machine unit roundoff quantity
66
C RPAR,IPAR -- parameter arrays which you may choose to use
67
C for communication between your program and subroutine F.
68
C They are not altered or used by DSTEPS.
69
C The variables X,XOLD,KOLD,KGI and IVC and the arrays Y,PHI,ALPHA,G,
70
C W,P,IV and GI are required for the interpolation subroutine SINTRP.
71
C The remaining variables and arrays are included in the call list
72
C only to eliminate local retention of variables between calls.
78
C The user must provide storage in his calling program for all arrays
79
C in the call list, namely
81
C DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12),
82
C 1 ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
87
C The user must also declare START , CRASH , PHASE1 and NORND
88
C logical variables and DF an EXTERNAL subroutine, supply the
89
C subroutine DF(X,Y,YP) to evaluate
90
C DY(I)/DX = YP(I) = DF(X,Y(1),Y(2),...,Y(NEQN))
91
C and initialize only the following parameters.
92
C NEQN -- number of equations to be integrated
93
C Y(*) -- vector of initial values of dependent variables
94
C X -- initial value of the independent variable
95
C H -- nominal step size indicating direction of integration
96
C and maximum size of step. Must be variable
97
C EPS -- local error tolerance per step. Must be variable
98
C WT(*) -- vector of non-zero weights for error criterion
100
C YP(*) -- vector of initial derivative values
101
C KSTEPS -- set KSTEPS to zero
102
C TWOU -- 2.*U where U is machine unit roundoff quantity
103
C FOURU -- 4.*U where U is machine unit roundoff quantity
104
C Define U to be the machine unit roundoff quantity by calling
105
C the function routine D1MACH, U = D1MACH(4), or by
106
C computing U so that U is the smallest positive number such
107
C that 1.0+U .GT. 1.0.
109
C DSTEPS requires that the L2 norm of the vector with components
110
C LOCAL ERROR(L)/WT(L) be less than EPS for a successful step. The
111
C array WT allows the user to specify an error test appropriate
112
C for his problem. For example,
113
C WT(L) = 1.0 specifies absolute error,
114
C = ABS(Y(L)) error relative to the most recent value of the
115
C L-th component of the solution,
116
C = ABS(YP(L)) error relative to the most recent value of
117
C the L-th component of the derivative,
118
C = MAX(WT(L),ABS(Y(L))) error relative to the largest
119
C magnitude of L-th component obtained so far,
120
C = ABS(Y(L))*RELERR/EPS + ABSERR/EPS specifies a mixed
121
C relative-absolute test where RELERR is relative
122
C error, ABSERR is absolute error and EPS =
123
C MAX(RELERR,ABSERR) .
125
C Subsequent calls --
127
C Subroutine DSTEPS is designed so that all information needed to
128
C continue the integration, including the step size H and the order
129
C K , is returned with each step. With the exception of the step
130
C size, the error tolerance, and the weights, none of the parameters
131
C should be altered. The array WT must be updated after each step
132
C to maintain relative error tests like those above. Normally the
133
C integration is continued just beyond the desired endpoint and the
134
C solution interpolated there with subroutine SINTRP . If it is
135
C impossible to integrate beyond the endpoint, the step size may be
136
C reduced to hit the endpoint since the code will not take a step
137
C larger than the H input. Changing the direction of integration,
138
C i.e., the sign of H , requires the user set START = .TRUE. before
139
C calling DSTEPS again. This is the only situation in which START
146
C The subroutine returns after each successful step with START and
147
C CRASH set .FALSE. . X represents the independent variable
148
C advanced one step of length HOLD from its value on input and Y
149
C the solution vector at the new value of X . All other parameters
150
C represent information corresponding to the new X needed to
151
C continue the integration.
153
C Unsuccessful Step --
155
C When the error tolerance is too small for the machine precision,
156
C the subroutine returns without taking a step and CRASH = .TRUE. .
157
C An appropriate step size and error tolerance for continuing are
158
C estimated and all other information is restored as upon input
159
C before returning. To continue with the larger tolerance, the user
160
C just calls the code again. A restart is neither required nor
163
C***REFERENCES L. F. Shampine and M. K. Gordon, Solving ordinary
164
C differential equations with ODE, STEP, and INTRP,
165
C Report SLA-73-1060, Sandia Laboratories, 1973.
166
C***ROUTINES CALLED D1MACH, DHSTRT
167
C***REVISION HISTORY (YYMMDD)
168
C 740101 DATE WRITTEN
169
C 890531 Changed all specific intrinsics to generic. (WRB)
170
C 890831 Modified array declarations. (WRB)
171
C 890831 REVISION DATE from Version 3.2
172
C 891214 Prologue converted to Version 4.0 format. (BAB)
173
C 920501 Reformatted the REFERENCES section. (WRB)
174
C***END PROLOGUE DSTEPS
176
INTEGER I, IFAIL, IM1, IP1, IPAR, IQ, J, K, KM1, KM2, KNEW,
177
1 KOLD, KP1, KP2, KSTEPS, L, LIMIT1, LIMIT2, NEQN, NS, NSM2,
179
DOUBLE PRECISION ABSH, ALPHA, BETA, BIG, D1MACH,
180
1 EPS, ERK, ERKM1, ERKM2, ERKP1, ERR,
181
2 FOURU, G, GI, GSTR, H, HNEW, HOLD, P, P5EPS, PHI, PSI, R,
182
3 REALI, REALNS, RHO, ROUND, RPAR, SIG, TAU, TEMP1,
183
4 TEMP2, TEMP3, TEMP4, TEMP5, TEMP6, TWO, TWOU, U, V, W, WT,
185
LOGICAL START,CRASH,PHASE1,NORND
186
DIMENSION Y(*),WT(*),PHI(NEQN,16),P(*),YP(*),PSI(12),
187
1 ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
189
DIMENSION TWO(13),GSTR(13)
193
DATA TWO(1),TWO(2),TWO(3),TWO(4),TWO(5),TWO(6),TWO(7),TWO(8),
194
1 TWO(9),TWO(10),TWO(11),TWO(12),TWO(13)
195
2 /2.0D0,4.0D0,8.0D0,16.0D0,32.0D0,64.0D0,128.0D0,256.0D0,
196
3 512.0D0,1024.0D0,2048.0D0,4096.0D0,8192.0D0/
197
DATA GSTR(1),GSTR(2),GSTR(3),GSTR(4),GSTR(5),GSTR(6),GSTR(7),
198
1 GSTR(8),GSTR(9),GSTR(10),GSTR(11),GSTR(12),GSTR(13)
199
2 /0.5D0,0.0833D0,0.0417D0,0.0264D0,0.0188D0,0.0143D0,0.0114D0,
200
3 0.00936D0,0.00789D0,0.00679D0,0.00592D0,0.00524D0,0.00468D0/
202
C *** BEGIN BLOCK 0 ***
203
C CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE
204
C PRECISION. IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A
205
C STARTING STEP SIZE.
208
C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE
210
C***FIRST EXECUTABLE STATEMENT DSTEPS
212
IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 5
213
H = SIGN(FOURU*ABS(X),H)
217
C IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE
221
10 ROUND = ROUND + (Y(L)/WT(L))**2
222
ROUND = TWOU*SQRT(ROUND)
223
IF(P5EPS .GE. ROUND) GO TO 15
224
EPS = 2.0D0*ROUND*(1.0D0 + FOURU)
230
IF(.NOT.START) GO TO 99
232
C INITIALIZE. COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP
234
C CALL DF(X,Y,YP,RPAR,IPAR)
239
C20 SUM = SUM + (YP(L)/WT(L))**2
242
C IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM)
243
C H = SIGN(MAX(ABSH,FOURU*ABS(X)),H)
246
BIG = SQRT(D1MACH(2))
247
CALL DHSTRT(DF,NEQN,X,X+H,Y,YP,WT,1,U,BIG,
248
1 PHI(1,3),PHI(1,4),PHI(1,5),PHI(1,6),RPAR,IPAR,H)
257
IF(P5EPS .GT. 100.0D0*ROUND) GO TO 99
262
C *** END BLOCK 0 ***
264
C *** BEGIN BLOCK 1 ***
265
C COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP. AVOID COMPUTING
266
C THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED.
274
C NS IS THE NUMBER OF DSTEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT
275
C ONE. WHEN K.LT.NS, NO COEFFICIENTS CHANGE
277
IF(H .NE. HOLD) NS = 0
278
IF (NS.LE.KOLD) NS = NS+1
280
IF (K .LT. NS) GO TO 199
282
C COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH
287
ALPHA(NS) = 1.0D0/REALNS
290
IF(K .LT. NSP1) GO TO 110
295
BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2
299
105 SIG(I+1) = REALI*ALPHA(I)*SIG(I)
302
C COMPUTE COEFFICIENTS G(*)
304
C INITIALIZE V(*) AND SET W(*).
306
IF(NS .GT. 1) GO TO 120
313
IF (K .EQ. 1) GO TO 140
318
C IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*)
320
120 IF(K .LE. KPREV) GO TO 130
321
IF (IVC .EQ. 0) GO TO 122
329
IF (K .NE. 2) GO TO 123
333
IF(NSM2 .LT. JV) GO TO 130
336
V(I) = V(I) - ALPHA(J+1)*V(I+1)
338
IF (I .NE. 2) GO TO 130
342
C UPDATE V(*) AND SET W(*)
344
130 LIMIT1 = KP1 - NS
347
V(IQ) = V(IQ) - TEMP5*V(IQ+1)
350
IF (LIMIT1 .EQ. 1) GO TO 137
353
137 W(LIMIT1+1) = V(LIMIT1+1)
354
IF (K .GE. KOLD) GO TO 140
358
C COMPUTE THE G(*) IN THE WORK VECTOR W(*)
362
IF(KP1 .LT. NSP2) GO TO 199
367
145 W(IQ) = W(IQ) - TEMP6*W(IQ+1)
370
C *** END BLOCK 1 ***
372
C *** BEGIN BLOCK 2 ***
373
C PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED
374
C SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K,
375
C K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED.
378
C INCREMENT COUNTER ON ATTEMPTED DSTEPS
382
C CHANGE PHI TO PHI STAR
384
IF(K .LT. NSP1) GO TO 215
388
205 PHI(L,I) = TEMP1*PHI(L,I)
391
C PREDICT SOLUTION AND DIFFERENCES
393
215 DO 220 L = 1,NEQN
394
PHI(L,KP2) = PHI(L,KP1)
402
P(L) = P(L) + TEMP2*PHI(L,I)
403
225 PHI(L,I) = PHI(L,I) + PHI(L,IP1)
407
TAU = H*P(L) - PHI(L,15)
409
235 PHI(L,16) = (P(L) - Y(L)) - TAU
411
240 DO 245 L = 1,NEQN
412
245 P(L) = Y(L) + H*P(L)
416
CALL DF(X,P,YP,RPAR,IPAR)
418
C ESTIMATE ERRORS AT ORDERS K,K-1,K-2
425
TEMP4 = YP(L) - PHI(L,1)
427
255 ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2
428
260 ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2
429
265 ERK = ERK + (TEMP4*TEMP3)**2
431
270 ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*SQRT(ERKM2)
432
275 ERKM1 = ABSH*SIG(K)*GSTR(KM1)*SQRT(ERKM1)
433
280 TEMP5 = ABSH*SQRT(ERK)
434
ERR = TEMP5*(G(K)-G(KP1))
435
ERK = TEMP5*SIG(KP1)*GSTR(K)
438
C TEST IF ORDER SHOULD BE LOWERED
441
285 IF(MAX(ERKM1,ERKM2) .LE. ERK) KNEW = KM1
443
290 IF(ERKM1 .LE. 0.5D0*ERK) KNEW = KM1
445
C TEST IF STEP SUCCESSFUL
447
299 IF(ERR .LE. EPS) GO TO 400
448
C *** END BLOCK 2 ***
450
C *** BEGIN BLOCK 3 ***
451
C THE STEP IS UNSUCCESSFUL. RESTORE X, PHI(*,*), PSI(*) .
452
C IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE. IF STEP FAILS MORE
453
C THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE. DOUBLE ERROR
454
C TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE
458
C RESTORE X, PHI(*,*) AND PSI(*)
463
TEMP1 = 1.0D0/BETA(I)
466
305 PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1))
468
IF(K .LT. 2) GO TO 320
470
315 PSI(I-1) = PSI(I) - H
472
C ON THIRD FAILURE, SET ORDER TO ONE. THEREAFTER, USE OPTIMAL STEP
475
320 IFAIL = IFAIL + 1
477
IF(IFAIL - 3) 335,330,325
478
325 IF(P5EPS .LT. 0.25D0*ERK) TEMP2 = SQRT(P5EPS/ERK)
483
IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 340
485
H = SIGN(FOURU*ABS(X),H)
489
C *** END BLOCK 3 ***
491
C *** BEGIN BLOCK 4 ***
492
C THE STEP IS SUCCESSFUL. CORRECT THE PREDICTED SOLUTION, EVALUATE
493
C THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE
494
C DIFFERENCES. DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP.
499
C CORRECT AND EVALUATE
505
RHO = TEMP1*(YP(L) - PHI(L,1)) - PHI(L,16)
507
PHI(L,15) = (Y(L) - P(L)) - RHO
510
410 DO 415 L = 1,NEQN
512
Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1))
514
420 CALL DF(X,Y,YP,RPAR,IPAR)
516
C UPDATE DIFFERENCES FOR NEXT STEP
519
PHI(L,KP1) = YP(L) - PHI(L,1)
520
425 PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2)
523
430 PHI(L,I) = PHI(L,I) + PHI(L,KP1)
526
C ESTIMATE ERROR AT ORDER K+1 UNLESS:
527
C IN FIRST PHASE WHEN ALWAYS RAISE ORDER,
528
C ALREADY DECIDED TO LOWER ORDER,
529
C STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE
532
IF(KNEW .EQ. KM1 .OR. K .EQ. 12) PHASE1 = .FALSE.
534
IF(KNEW .EQ. KM1) GO TO 455
535
IF(KP1 .GT. NS) GO TO 460
537
440 ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2
538
ERKP1 = ABSH*GSTR(KP1)*SQRT(ERKP1)
540
C USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER
543
IF(K .GT. 1) GO TO 445
544
IF(ERKP1 .GE. 0.5D0*ERK) GO TO 460
546
445 IF(ERKM1 .LE. MIN(ERK,ERKP1)) GO TO 455
547
IF(ERKP1 .GE. ERK .OR. K .EQ. 12) GO TO 460
549
C HERE ERKP1 .LT. ERK .LT. MAX(ERKM1,ERKM2) ELSE ORDER WOULD HAVE
550
C BEEN LOWERED IN BLOCK 2. THUS ORDER IS TO BE RAISED
563
C WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP
567
IF(P5EPS .GE. ERK*TWO(K+1)) GO TO 465
569
IF(P5EPS .GE. ERK) GO TO 465
571
R = (P5EPS/ERK)**(1.0D0/TEMP2)
572
HNEW = ABSH*MAX(0.5D0,MIN(0.9D0,R))
573
HNEW = SIGN(MAX(HNEW,FOURU*ABS(X)),H)
576
C *** END BLOCK 4 ***
added
removed
deps/openlibm/slatec/dsteps.f
