|
2329
|
1 SUBROUTINE DDASTP (X, Y, YPRIME, NEQ, RES, JAC, H, WT, JSTART, |
|
|
2 + IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, ALPHA, BETA, GAMMA, |
|
|
3 + PSI, SIGMA, CJ, CJOLD, HOLD, S, HMIN, UROUND, IPHASE, JCALC, |
|
|
4 + K, KOLD, NS, NONNEG, NTEMP) |
|
|
5 C***BEGIN PROLOGUE DDASTP |
|
|
6 C***SUBSIDIARY |
|
|
7 C***PURPOSE Perform one step of the DDASSL integration. |
|
|
8 C***LIBRARY SLATEC (DASSL) |
|
|
9 C***TYPE DOUBLE PRECISION (SDASTP-S, DDASTP-D) |
|
|
10 C***AUTHOR PETZOLD, LINDA R., (LLNL) |
|
|
11 C***DESCRIPTION |
|
|
12 C----------------------------------------------------------------------- |
|
|
13 C DDASTP SOLVES A SYSTEM OF DIFFERENTIAL/ |
|
|
14 C ALGEBRAIC EQUATIONS OF THE FORM |
|
|
15 C G(X,Y,YPRIME) = 0, FOR ONE STEP (NORMALLY |
|
|
16 C FROM X TO X+H). |
|
|
17 C |
|
|
18 C THE METHODS USED ARE MODIFIED DIVIDED |
|
|
19 C DIFFERENCE,FIXED LEADING COEFFICIENT |
|
|
20 C FORMS OF BACKWARD DIFFERENTIATION |
|
|
21 C FORMULAS. THE CODE ADJUSTS THE STEPSIZE |
|
|
22 C AND ORDER TO CONTROL THE LOCAL ERROR PER |
|
|
23 C STEP. |
|
|
24 C |
|
|
25 C |
|
|
26 C THE PARAMETERS REPRESENT |
|
|
27 C X -- INDEPENDENT VARIABLE |
|
|
28 C Y -- SOLUTION VECTOR AT X |
|
|
29 C YPRIME -- DERIVATIVE OF SOLUTION VECTOR |
|
|
30 C AFTER SUCCESSFUL STEP |
|
|
31 C NEQ -- NUMBER OF EQUATIONS TO BE INTEGRATED |
|
|
32 C RES -- EXTERNAL USER-SUPPLIED SUBROUTINE |
|
|
33 C TO EVALUATE THE RESIDUAL. THE CALL IS |
|
|
34 C CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR) |
|
|
35 C X,Y,YPRIME ARE INPUT. DELTA IS OUTPUT. |
|
|
36 C ON INPUT, IRES=0. RES SHOULD ALTER IRES ONLY |
|
|
37 C IF IT ENCOUNTERS AN ILLEGAL VALUE OF Y OR A |
|
|
38 C STOP CONDITION. SET IRES=-1 IF AN INPUT VALUE |
|
|
39 C OF Y IS ILLEGAL, AND DDASTP WILL TRY TO SOLVE |
|
|
40 C THE PROBLEM WITHOUT GETTING IRES = -1. IF |
|
|
41 C IRES=-2, DDASTP RETURNS CONTROL TO THE CALLING |
|
|
42 C PROGRAM WITH IDID = -11. |
|
|
43 C JAC -- EXTERNAL USER-SUPPLIED ROUTINE TO EVALUATE |
|
|
44 C THE ITERATION MATRIX (THIS IS OPTIONAL) |
|
|
45 C THE CALL IS OF THE FORM |
|
|
46 C CALL JAC(X,Y,YPRIME,PD,CJ,RPAR,IPAR) |
|
|
47 C PD IS THE MATRIX OF PARTIAL DERIVATIVES, |
|
|
48 C PD=DG/DY+CJ*DG/DYPRIME |
|
|
49 C H -- APPROPRIATE STEP SIZE FOR NEXT STEP. |
|
|
50 C NORMALLY DETERMINED BY THE CODE |
|
|
51 C WT -- VECTOR OF WEIGHTS FOR ERROR CRITERION. |
|
|
52 C JSTART -- INTEGER VARIABLE SET 0 FOR |
|
|
53 C FIRST STEP, 1 OTHERWISE. |
|
|
54 C IDID -- COMPLETION CODE WITH THE FOLLOWING MEANINGS: |
|
|
55 C IDID= 1 -- THE STEP WAS COMPLETED SUCCESSFULLY |
|
|
56 C IDID=-6 -- THE ERROR TEST FAILED REPEATEDLY |
|
|
57 C IDID=-7 -- THE CORRECTOR COULD NOT CONVERGE |
|
|
58 C IDID=-8 -- THE ITERATION MATRIX IS SINGULAR |
|
|
59 C IDID=-9 -- THE CORRECTOR COULD NOT CONVERGE. |
|
|
60 C THERE WERE REPEATED ERROR TEST |
|
|
61 C FAILURES ON THIS STEP. |
|
|
62 C IDID=-10-- THE CORRECTOR COULD NOT CONVERGE |
|
|
63 C BECAUSE IRES WAS EQUAL TO MINUS ONE |
|
|
64 C IDID=-11-- IRES EQUAL TO -2 WAS ENCOUNTERED, |
|
|
65 C AND CONTROL IS BEING RETURNED TO |
|
|
66 C THE CALLING PROGRAM |
|
|
67 C RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS THAT |
|
|
68 C ARE USED FOR COMMUNICATION BETWEEN THE |
|
|
69 C CALLING PROGRAM AND EXTERNAL USER ROUTINES |
|
|
70 C THEY ARE NOT ALTERED BY DDASTP |
|
|
71 C PHI -- ARRAY OF DIVIDED DIFFERENCES USED BY |
|
|
72 C DDASTP. THE LENGTH IS NEQ*(K+1),WHERE |
|
|
73 C K IS THE MAXIMUM ORDER |
|
|
74 C DELTA,E -- WORK VECTORS FOR DDASTP OF LENGTH NEQ |
|
|
75 C WM,IWM -- REAL AND INTEGER ARRAYS STORING |
|
|
76 C MATRIX INFORMATION SUCH AS THE MATRIX |
|
|
77 C OF PARTIAL DERIVATIVES,PERMUTATION |
|
|
78 C VECTOR,AND VARIOUS OTHER INFORMATION. |
|
|
79 C |
|
|
80 C THE OTHER PARAMETERS ARE INFORMATION |
|
|
81 C WHICH IS NEEDED INTERNALLY BY DDASTP TO |
|
|
82 C CONTINUE FROM STEP TO STEP. |
|
|
83 C |
|
|
84 C----------------------------------------------------------------------- |
|
|
85 C***ROUTINES CALLED DDAJAC, DDANRM, DDASLV, DDATRP |
|
|
86 C***REVISION HISTORY (YYMMDD) |
|
|
87 C 830315 DATE WRITTEN |
|
|
88 C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch) |
|
|
89 C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format. |
|
|
90 C 901026 Added explicit declarations for all variables and minor |
|
|
91 C cosmetic changes to prologue. (FNF) |
|
|
92 C***END PROLOGUE DDASTP |
|
|
93 C |
|
|
94 INTEGER NEQ, JSTART, IDID, IPAR(*), IWM(*), IPHASE, JCALC, K, |
|
|
95 * KOLD, NS, NONNEG, NTEMP |
|
|
96 DOUBLE PRECISION |
|
|
97 * X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*), |
|
|
98 * E(*), WM(*), ALPHA(*), BETA(*), GAMMA(*), PSI(*), SIGMA(*), CJ, |
|
|
99 * CJOLD, HOLD, S, HMIN, UROUND |
|
|
100 EXTERNAL RES, JAC |
|
|
101 C |
|
|
102 EXTERNAL DDAJAC, DDANRM, DDASLV, DDATRP |
|
|
103 DOUBLE PRECISION DDANRM |
|
|
104 C |
|
|
105 INTEGER I, IER, IRES, J, J1, KDIFF, KM1, KNEW, KP1, KP2, LCTF, |
|
|
106 * LETF, LMXORD, LNJE, LNRE, LNST, M, MAXIT, NCF, NEF, NSF, NSP1 |
|
|
107 DOUBLE PRECISION |
|
|
108 * ALPHA0, ALPHAS, CJLAST, CK, DELNRM, ENORM, ERK, ERKM1, |
|
|
109 * ERKM2, ERKP1, ERR, EST, HNEW, OLDNRM, PNORM, R, RATE, TEMP1, |
|
|
110 * TEMP2, TERK, TERKM1, TERKM2, TERKP1, XOLD, XRATE |
|
|
111 LOGICAL CONVGD |
|
|
112 C |
|
|
113 PARAMETER (LMXORD=3) |
|
|
114 PARAMETER (LNST=11) |
|
|
115 PARAMETER (LNRE=12) |
|
|
116 PARAMETER (LNJE=13) |
|
|
117 PARAMETER (LETF=14) |
|
|
118 PARAMETER (LCTF=15) |
|
|
119 C |
|
|
120 DATA MAXIT/4/ |
|
|
121 DATA XRATE/0.25D0/ |
|
|
122 C |
|
|
123 C |
|
|
124 C |
|
|
125 C |
|
|
126 C |
|
|
127 C----------------------------------------------------------------------- |
|
|
128 C BLOCK 1. |
|
|
129 C INITIALIZE. ON THE FIRST CALL,SET |
|
|
130 C THE ORDER TO 1 AND INITIALIZE |
|
|
131 C OTHER VARIABLES. |
|
|
132 C----------------------------------------------------------------------- |
|
|
133 C |
|
|
134 C INITIALIZATIONS FOR ALL CALLS |
|
|
135 C***FIRST EXECUTABLE STATEMENT DDASTP |
|
|
136 IDID=1 |
|
|
137 XOLD=X |
|
|
138 NCF=0 |
|
|
139 NSF=0 |
|
|
140 NEF=0 |
|
|
141 IF(JSTART .NE. 0) GO TO 120 |
|
|
142 C |
|
|
143 C IF THIS IS THE FIRST STEP,PERFORM |
|
|
144 C OTHER INITIALIZATIONS |
|
|
145 IWM(LETF) = 0 |
|
|
146 IWM(LCTF) = 0 |
|
|
147 K=1 |
|
|
148 KOLD=0 |
|
|
149 HOLD=0.0D0 |
|
|
150 JSTART=1 |
|
|
151 PSI(1)=H |
|
|
152 CJOLD = 1.0D0/H |
|
|
153 CJ = CJOLD |
|
|
154 S = 100.D0 |
|
|
155 JCALC = -1 |
|
|
156 DELNRM=1.0D0 |
|
|
157 IPHASE = 0 |
|
|
158 NS=0 |
|
|
159 120 CONTINUE |
|
|
160 C |
|
|
161 C |
|
|
162 C |
|
|
163 C |
|
|
164 C |
|
|
165 C----------------------------------------------------------------------- |
|
|
166 C BLOCK 2 |
|
|
167 C COMPUTE COEFFICIENTS OF FORMULAS FOR |
|
|
168 C THIS STEP. |
|
|
169 C----------------------------------------------------------------------- |
|
|
170 200 CONTINUE |
|
|
171 KP1=K+1 |
|
|
172 KP2=K+2 |
|
|
173 KM1=K-1 |
|
|
174 XOLD=X |
|
|
175 IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0 |
|
|
176 NS=MIN(NS+1,KOLD+2) |
|
|
177 NSP1=NS+1 |
|
|
178 IF(KP1 .LT. NS)GO TO 230 |
|
|
179 C |
|
|
180 BETA(1)=1.0D0 |
|
|
181 ALPHA(1)=1.0D0 |
|
|
182 TEMP1=H |
|
|
183 GAMMA(1)=0.0D0 |
|
|
184 SIGMA(1)=1.0D0 |
|
|
185 DO 210 I=2,KP1 |
|
|
186 TEMP2=PSI(I-1) |
|
|
187 PSI(I-1)=TEMP1 |
|
|
188 BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2 |
|
|
189 TEMP1=TEMP2+H |
|
|
190 ALPHA(I)=H/TEMP1 |
|
|
191 SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I) |
|
|
192 GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H |
|
|
193 210 CONTINUE |
|
|
194 PSI(KP1)=TEMP1 |
|
|
195 230 CONTINUE |
|
|
196 C |
|
|
197 C COMPUTE ALPHAS, ALPHA0 |
|
|
198 ALPHAS = 0.0D0 |
|
|
199 ALPHA0 = 0.0D0 |
|
|
200 DO 240 I = 1,K |
|
|
201 ALPHAS = ALPHAS - 1.0D0/I |
|
|
202 ALPHA0 = ALPHA0 - ALPHA(I) |
|
|
203 240 CONTINUE |
|
|
204 C |
|
|
205 C COMPUTE LEADING COEFFICIENT CJ |
|
|
206 CJLAST = CJ |
|
|
207 CJ = -ALPHAS/H |
|
|
208 C |
|
|
209 C COMPUTE VARIABLE STEPSIZE ERROR COEFFICIENT CK |
|
|
210 CK = ABS(ALPHA(KP1) + ALPHAS - ALPHA0) |
|
|
211 CK = MAX(CK,ALPHA(KP1)) |
|
|
212 C |
|
|
213 C DECIDE WHETHER NEW JACOBIAN IS NEEDED |
|
|
214 TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE) |
|
|
215 TEMP2 = 1.0D0/TEMP1 |
|
|
216 IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1 |
|
|
217 IF (CJ .NE. CJLAST) S = 100.D0 |
|
|
218 C |
|
|
219 C CHANGE PHI TO PHI STAR |
|
|
220 IF(KP1 .LT. NSP1) GO TO 280 |
|
|
221 DO 270 J=NSP1,KP1 |
|
|
222 DO 260 I=1,NEQ |
|
|
223 260 PHI(I,J)=BETA(J)*PHI(I,J) |
|
|
224 270 CONTINUE |
|
|
225 280 CONTINUE |
|
|
226 C |
|
|
227 C UPDATE TIME |
|
|
228 X=X+H |
|
|
229 C |
|
|
230 C |
|
|
231 C |
|
|
232 C |
|
|
233 C |
|
|
234 C----------------------------------------------------------------------- |
|
|
235 C BLOCK 3 |
|
|
236 C PREDICT THE SOLUTION AND DERIVATIVE, |
|
|
237 C AND SOLVE THE CORRECTOR EQUATION |
|
|
238 C----------------------------------------------------------------------- |
|
|
239 C |
|
|
240 C FIRST,PREDICT THE SOLUTION AND DERIVATIVE |
|
|
241 300 CONTINUE |
|
|
242 DO 310 I=1,NEQ |
|
|
243 Y(I)=PHI(I,1) |
|
|
244 310 YPRIME(I)=0.0D0 |
|
|
245 DO 330 J=2,KP1 |
|
|
246 DO 320 I=1,NEQ |
|
|
247 Y(I)=Y(I)+PHI(I,J) |
|
|
248 320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J) |
|
|
249 330 CONTINUE |
|
|
250 PNORM = DDANRM (NEQ,Y,WT,RPAR,IPAR) |
|
|
251 C |
|
|
252 C |
|
|
253 C |
|
|
254 C SOLVE THE CORRECTOR EQUATION USING A |
|
|
255 C MODIFIED NEWTON SCHEME. |
|
|
256 CONVGD= .TRUE. |
|
|
257 M=0 |
|
|
258 IWM(LNRE)=IWM(LNRE)+1 |
|
|
259 IRES = 0 |
|
|
260 CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR) |
|
|
261 IF (IRES .LT. 0) GO TO 380 |
|
|
262 C |
|
|
263 C |
|
|
264 C IF INDICATED,REEVALUATE THE |
|
|
265 C ITERATION MATRIX PD = DG/DY + CJ*DG/DYPRIME |
|
|
266 C (WHERE G(X,Y,YPRIME)=0). SET |
|
|
267 C JCALC TO 0 AS AN INDICATOR THAT |
|
|
268 C THIS HAS BEEN DONE. |
|
|
269 IF(JCALC .NE. -1)GO TO 340 |
|
|
270 IWM(LNJE)=IWM(LNJE)+1 |
|
|
271 JCALC=0 |
|
|
272 CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H, |
|
|
273 * IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR, |
|
|
274 * IPAR,NTEMP) |
|
|
275 CJOLD=CJ |
|
|
276 S = 100.D0 |
|
|
277 IF (IRES .LT. 0) GO TO 380 |
|
|
278 IF(IER .NE. 0)GO TO 380 |
|
|
279 NSF=0 |
|
|
280 C |
|
|
281 C |
|
|
282 C INITIALIZE THE ERROR ACCUMULATION VECTOR E. |
|
|
283 340 CONTINUE |
|
|
284 DO 345 I=1,NEQ |
|
|
285 345 E(I)=0.0D0 |
|
|
286 C |
|
|
287 C |
|
|
288 C CORRECTOR LOOP. |
|
|
289 350 CONTINUE |
|
|
290 C |
|
|
291 C MULTIPLY RESIDUAL BY TEMP1 TO ACCELERATE CONVERGENCE |
|
|
292 TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD) |
|
|
293 DO 355 I = 1,NEQ |
|
|
294 355 DELTA(I) = DELTA(I) * TEMP1 |
|
|
295 C |
|
|
296 C COMPUTE A NEW ITERATE (BACK-SUBSTITUTION). |
|
|
297 C STORE THE CORRECTION IN DELTA. |
|
|
298 CALL DDASLV(NEQ,DELTA,WM,IWM) |
|
|
299 C |
|
|
300 C UPDATE Y,E,AND YPRIME |
|
|
301 DO 360 I=1,NEQ |
|
|
302 Y(I)=Y(I)-DELTA(I) |
|
|
303 E(I)=E(I)-DELTA(I) |
|
|
304 360 YPRIME(I)=YPRIME(I)-CJ*DELTA(I) |
|
|
305 C |
|
|
306 C TEST FOR CONVERGENCE OF THE ITERATION |
|
|
307 DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR) |
|
|
308 IF (DELNRM .LE. 100.D0*UROUND*PNORM) GO TO 375 |
|
|
309 IF (M .GT. 0) GO TO 365 |
|
|
310 OLDNRM = DELNRM |
|
|
311 GO TO 367 |
|
|
312 365 RATE = (DELNRM/OLDNRM)**(1.0D0/M) |
|
|
313 IF (RATE .GT. 0.90D0) GO TO 370 |
|
|
314 S = RATE/(1.0D0 - RATE) |
|
|
315 367 IF (S*DELNRM .LE. 0.33D0) GO TO 375 |
|
|
316 C |
|
|
317 C THE CORRECTOR HAS NOT YET CONVERGED. |
|
|
318 C UPDATE M AND TEST WHETHER THE |
|
|
319 C MAXIMUM NUMBER OF ITERATIONS HAVE |
|
|
320 C BEEN TRIED. |
|
|
321 M=M+1 |
|
|
322 IF(M.GE.MAXIT)GO TO 370 |
|
|
323 C |
|
|
324 C EVALUATE THE RESIDUAL |
|
|
325 C AND GO BACK TO DO ANOTHER ITERATION |
|
|
326 IWM(LNRE)=IWM(LNRE)+1 |
|
|
327 IRES = 0 |
|
|
328 CALL RES(X,Y,YPRIME,DELTA,IRES, |
|
|
329 * RPAR,IPAR) |
|
|
330 IF (IRES .LT. 0) GO TO 380 |
|
|
331 GO TO 350 |
|
|
332 C |
|
|
333 C |
|
|
334 C THE CORRECTOR FAILED TO CONVERGE IN MAXIT |
|
|
335 C ITERATIONS. IF THE ITERATION MATRIX |
|
|
336 C IS NOT CURRENT,RE-DO THE STEP WITH |
|
|
337 C A NEW ITERATION MATRIX. |
|
|
338 370 CONTINUE |
|
|
339 IF(JCALC.EQ.0)GO TO 380 |
|
|
340 JCALC=-1 |
|
|
341 GO TO 300 |
|
|
342 C |
|
|
343 C |
|
|
344 C THE ITERATION HAS CONVERGED. IF NONNEGATIVITY OF SOLUTION IS |
|
|
345 C REQUIRED, SET THE SOLUTION NONNEGATIVE, IF THE PERTURBATION |
|
|
346 C TO DO IT IS SMALL ENOUGH. IF THE CHANGE IS TOO LARGE, THEN |
|
|
347 C CONSIDER THE CORRECTOR ITERATION TO HAVE FAILED. |
|
|
348 375 IF(NONNEG .EQ. 0) GO TO 390 |
|
|
349 DO 377 I = 1,NEQ |
|
|
350 377 DELTA(I) = MIN(Y(I),0.0D0) |
|
|
351 DELNRM = DDANRM(NEQ,DELTA,WT,RPAR,IPAR) |
|
|
352 IF(DELNRM .GT. 0.33D0) GO TO 380 |
|
|
353 DO 378 I = 1,NEQ |
|
|
354 378 E(I) = E(I) - DELTA(I) |
|
|
355 GO TO 390 |
|
|
356 C |
|
|
357 C |
|
|
358 C EXITS FROM BLOCK 3 |
|
|
359 C NO CONVERGENCE WITH CURRENT ITERATION |
|
|
360 C MATRIX,OR SINGULAR ITERATION MATRIX |
|
|
361 380 CONVGD= .FALSE. |
|
|
362 390 JCALC = 1 |
|
|
363 IF(.NOT.CONVGD)GO TO 600 |
|
|
364 C |
|
|
365 C |
|
|
366 C |
|
|
367 C |
|
|
368 C |
|
|
369 C----------------------------------------------------------------------- |
|
|
370 C BLOCK 4 |
|
|
371 C ESTIMATE THE ERRORS AT ORDERS K,K-1,K-2 |
|
|
372 C AS IF CONSTANT STEPSIZE WAS USED. ESTIMATE |
|
|
373 C THE LOCAL ERROR AT ORDER K AND TEST |
|
|
374 C WHETHER THE CURRENT STEP IS SUCCESSFUL. |
|
|
375 C----------------------------------------------------------------------- |
|
|
376 C |
|
|
377 C ESTIMATE ERRORS AT ORDERS K,K-1,K-2 |
|
|
378 ENORM = DDANRM(NEQ,E,WT,RPAR,IPAR) |
|
|
379 ERK = SIGMA(K+1)*ENORM |
|
|
380 TERK = (K+1)*ERK |
|
|
381 EST = ERK |
|
|
382 KNEW=K |
|
|
383 IF(K .EQ. 1)GO TO 430 |
|
|
384 DO 405 I = 1,NEQ |
|
|
385 405 DELTA(I) = PHI(I,KP1) + E(I) |
|
|
386 ERKM1=SIGMA(K)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR) |
|
|
387 TERKM1 = K*ERKM1 |
|
|
388 IF(K .GT. 2)GO TO 410 |
|
|
389 IF(TERKM1 .LE. 0.5D0*TERK)GO TO 420 |
|
|
390 GO TO 430 |
|
|
391 410 CONTINUE |
|
|
392 DO 415 I = 1,NEQ |
|
|
393 415 DELTA(I) = PHI(I,K) + DELTA(I) |
|
|
394 ERKM2=SIGMA(K-1)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR) |
|
|
395 TERKM2 = (K-1)*ERKM2 |
|
|
396 IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430 |
|
|
397 C LOWER THE ORDER |
|
|
398 420 CONTINUE |
|
|
399 KNEW=K-1 |
|
|
400 EST = ERKM1 |
|
|
401 C |
|
|
402 C |
|
|
403 C CALCULATE THE LOCAL ERROR FOR THE CURRENT STEP |
|
|
404 C TO SEE IF THE STEP WAS SUCCESSFUL |
|
|
405 430 CONTINUE |
|
|
406 ERR = CK * ENORM |
|
|
407 IF(ERR .GT. 1.0D0)GO TO 600 |
|
|
408 C |
|
|
409 C |
|
|
410 C |
|
|
411 C |
|
|
412 C |
|
|
413 C----------------------------------------------------------------------- |
|
|
414 C BLOCK 5 |
|
|
415 C THE STEP IS SUCCESSFUL. DETERMINE |
|
|
416 C THE BEST ORDER AND STEPSIZE FOR |
|
|
417 C THE NEXT STEP. UPDATE THE DIFFERENCES |
|
|
418 C FOR THE NEXT STEP. |
|
|
419 C----------------------------------------------------------------------- |
|
|
420 IDID=1 |
|
|
421 IWM(LNST)=IWM(LNST)+1 |
|
|
422 KDIFF=K-KOLD |
|
|
423 KOLD=K |
|
|
424 HOLD=H |
|
|
425 C |
|
|
426 C |
|
|
427 C ESTIMATE THE ERROR AT ORDER K+1 UNLESS: |
|
|
428 C ALREADY DECIDED TO LOWER ORDER, OR |
|
|
429 C ALREADY USING MAXIMUM ORDER, OR |
|
|
430 C STEPSIZE NOT CONSTANT, OR |
|
|
431 C ORDER RAISED IN PREVIOUS STEP |
|
|
432 IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1 |
|
|
433 IF(IPHASE .EQ. 0)GO TO 545 |
|
|
434 IF(KNEW.EQ.KM1)GO TO 540 |
|
|
435 IF(K.EQ.IWM(LMXORD)) GO TO 550 |
|
|
436 IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550 |
|
|
437 DO 510 I=1,NEQ |
|
|
438 510 DELTA(I)=E(I)-PHI(I,KP2) |
|
|
439 ERKP1 = (1.0D0/(K+2))*DDANRM(NEQ,DELTA,WT,RPAR,IPAR) |
|
|
440 TERKP1 = (K+2)*ERKP1 |
|
|
441 IF(K.GT.1)GO TO 520 |
|
|
442 IF(TERKP1.GE.0.5D0*TERK)GO TO 550 |
|
|
443 GO TO 530 |
|
|
444 520 IF(TERKM1.LE.MIN(TERK,TERKP1))GO TO 540 |
|
|
445 IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550 |
|
|
446 C |
|
|
447 C RAISE ORDER |
|
|
448 530 K=KP1 |
|
|
449 EST = ERKP1 |
|
|
450 GO TO 550 |
|
|
451 C |
|
|
452 C LOWER ORDER |
|
|
453 540 K=KM1 |
|
|
454 EST = ERKM1 |
|
|
455 GO TO 550 |
|
|
456 C |
|
|
457 C IF IPHASE = 0, INCREASE ORDER BY ONE AND MULTIPLY STEPSIZE BY |
|
|
458 C FACTOR TWO |
|
|
459 545 K = KP1 |
|
|
460 HNEW = H*2.0D0 |
|
|
461 H = HNEW |
|
|
462 GO TO 575 |
|
|
463 C |
|
|
464 C |
|
|
465 C DETERMINE THE APPROPRIATE STEPSIZE FOR |
|
|
466 C THE NEXT STEP. |
|
|
467 550 HNEW=H |
|
|
468 TEMP2=K+1 |
|
|
469 R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2) |
|
|
470 IF(R .LT. 2.0D0) GO TO 555 |
|
|
471 HNEW = 2.0D0*H |
|
|
472 GO TO 560 |
|
|
473 555 IF(R .GT. 1.0D0) GO TO 560 |
|
|
474 R = MAX(0.5D0,MIN(0.9D0,R)) |
|
|
475 HNEW = H*R |
|
|
476 560 H=HNEW |
|
|
477 C |
|
|
478 C |
|
|
479 C UPDATE DIFFERENCES FOR NEXT STEP |
|
|
480 575 CONTINUE |
|
|
481 IF(KOLD.EQ.IWM(LMXORD))GO TO 585 |
|
|
482 DO 580 I=1,NEQ |
|
|
483 580 PHI(I,KP2)=E(I) |
|
|
484 585 CONTINUE |
|
|
485 DO 590 I=1,NEQ |
|
|
486 590 PHI(I,KP1)=PHI(I,KP1)+E(I) |
|
|
487 DO 595 J1=2,KP1 |
|
|
488 J=KP1-J1+1 |
|
|
489 DO 595 I=1,NEQ |
|
|
490 595 PHI(I,J)=PHI(I,J)+PHI(I,J+1) |
|
|
491 RETURN |
|
|
492 C |
|
|
493 C |
|
|
494 C |
|
|
495 C |
|
|
496 C |
|
|
497 C----------------------------------------------------------------------- |
|
|
498 C BLOCK 6 |
|
|
499 C THE STEP IS UNSUCCESSFUL. RESTORE X,PSI,PHI |
|
|
500 C DETERMINE APPROPRIATE STEPSIZE FOR |
|
|
501 C CONTINUING THE INTEGRATION, OR EXIT WITH |
|
|
502 C AN ERROR FLAG IF THERE HAVE BEEN MANY |
|
|
503 C FAILURES. |
|
|
504 C----------------------------------------------------------------------- |
|
|
505 600 IPHASE = 1 |
|
|
506 C |
|
|
507 C RESTORE X,PHI,PSI |
|
|
508 X=XOLD |
|
|
509 IF(KP1.LT.NSP1)GO TO 630 |
|
|
510 DO 620 J=NSP1,KP1 |
|
|
511 TEMP1=1.0D0/BETA(J) |
|
|
512 DO 610 I=1,NEQ |
|
|
513 610 PHI(I,J)=TEMP1*PHI(I,J) |
|
|
514 620 CONTINUE |
|
|
515 630 CONTINUE |
|
|
516 DO 640 I=2,KP1 |
|
|
517 640 PSI(I-1)=PSI(I)-H |
|
|
518 C |
|
|
519 C |
|
|
520 C TEST WHETHER FAILURE IS DUE TO CORRECTOR ITERATION |
|
|
521 C OR ERROR TEST |
|
|
522 IF(CONVGD)GO TO 660 |
|
|
523 IWM(LCTF)=IWM(LCTF)+1 |
|
|
524 C |
|
|
525 C |
|
|
526 C THE NEWTON ITERATION FAILED TO CONVERGE WITH |
|
|
527 C A CURRENT ITERATION MATRIX. DETERMINE THE CAUSE |
|
|
528 C OF THE FAILURE AND TAKE APPROPRIATE ACTION. |
|
|
529 IF(IER.EQ.0)GO TO 650 |
|
|
530 C |
|
|
531 C THE ITERATION MATRIX IS SINGULAR. REDUCE |
|
|
532 C THE STEPSIZE BY A FACTOR OF 4. IF |
|
|
533 C THIS HAPPENS THREE TIMES IN A ROW ON |
|
|
534 C THE SAME STEP, RETURN WITH AN ERROR FLAG |
|
|
535 NSF=NSF+1 |
|
|
536 R = 0.25D0 |
|
|
537 H=H*R |
|
|
538 IF (NSF .LT. 3 .AND. ABS(H) .GE. HMIN) GO TO 690 |
|
|
539 IDID=-8 |
|
|
540 GO TO 675 |
|
|
541 C |
|
|
542 C |
|
|
543 C THE NEWTON ITERATION FAILED TO CONVERGE FOR A REASON |
|
|
544 C OTHER THAN A SINGULAR ITERATION MATRIX. IF IRES = -2, THEN |
|
|
545 C RETURN. OTHERWISE, REDUCE THE STEPSIZE AND TRY AGAIN, UNLESS |
|
|
546 C TOO MANY FAILURES HAVE OCCURED. |
|
|
547 650 CONTINUE |
|
|
548 IF (IRES .GT. -2) GO TO 655 |
|
|
549 IDID = -11 |
|
|
550 GO TO 675 |
|
|
551 655 NCF = NCF + 1 |
|
|
552 R = 0.25D0 |
|
|
553 H = H*R |
|
|
554 IF (NCF .LT. 10 .AND. ABS(H) .GE. HMIN) GO TO 690 |
|
|
555 IDID = -7 |
|
|
556 IF (IRES .LT. 0) IDID = -10 |
|
|
557 IF (NEF .GE. 3) IDID = -9 |
|
|
558 GO TO 675 |
|
|
559 C |
|
|
560 C |
|
|
561 C THE NEWTON SCHEME CONVERGED,AND THE CAUSE |
|
|
562 C OF THE FAILURE WAS THE ERROR ESTIMATE |
|
|
563 C EXCEEDING THE TOLERANCE. |
|
|
564 660 NEF=NEF+1 |
|
|
565 IWM(LETF)=IWM(LETF)+1 |
|
|
566 IF (NEF .GT. 1) GO TO 665 |
|
|
567 C |
|
|
568 C ON FIRST ERROR TEST FAILURE, KEEP CURRENT ORDER OR LOWER |
|
|
569 C ORDER BY ONE. COMPUTE NEW STEPSIZE BASED ON DIFFERENCES |
|
|
570 C OF THE SOLUTION. |
|
|
571 K = KNEW |
|
|
572 TEMP2 = K + 1 |
|
|
573 R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2) |
|
|
574 R = MAX(0.25D0,MIN(0.9D0,R)) |
|
|
575 H = H*R |
|
|
576 IF (ABS(H) .GE. HMIN) GO TO 690 |
|
|
577 IDID = -6 |
|
|
578 GO TO 675 |
|
|
579 C |
|
|
580 C ON SECOND ERROR TEST FAILURE, USE THE CURRENT ORDER OR |
|
|
581 C DECREASE ORDER BY ONE. REDUCE THE STEPSIZE BY A FACTOR OF |
|
|
582 C FOUR. |
|
|
583 665 IF (NEF .GT. 2) GO TO 670 |
|
|
584 K = KNEW |
|
|
585 H = 0.25D0*H |
|
|
586 IF (ABS(H) .GE. HMIN) GO TO 690 |
|
|
587 IDID = -6 |
|
|
588 GO TO 675 |
|
|
589 C |
|
|
590 C ON THIRD AND SUBSEQUENT ERROR TEST FAILURES, SET THE ORDER TO |
|
|
591 C ONE AND REDUCE THE STEPSIZE BY A FACTOR OF FOUR. |
|
|
592 670 K = 1 |
|
|
593 H = 0.25D0*H |
|
|
594 IF (ABS(H) .GE. HMIN) GO TO 690 |
|
|
595 IDID = -6 |
|
|
596 GO TO 675 |
|
|
597 C |
|
|
598 C |
|
|
599 C |
|
|
600 C |
|
|
601 C FOR ALL CRASHES, RESTORE Y TO ITS LAST VALUE, |
|
|
602 C INTERPOLATE TO FIND YPRIME AT LAST X, AND RETURN |
|
|
603 675 CONTINUE |
|
|
604 CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI) |
|
|
605 RETURN |
|
|
606 C |
|
|
607 C |
|
|
608 C GO BACK AND TRY THIS STEP AGAIN |
|
|
609 690 GO TO 200 |
|
|
610 C |
|
|
611 C------END OF SUBROUTINE DDASTP------ |
|
|
612 END |
