GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
stode.f
Go to the documentation of this file.
1  SUBROUTINE stode (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
2  1 wm, iwm, f, jac, pjac, slvs, ierr)
3 CLLL. OPTIMIZE
4  EXTERNAL f, jac, pjac, slvs
5  INTEGER NEQ, NYH, IWM
6  INTEGER ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
7  1 mxstep, mxhnil, nhnil, ntrep, nslast, cnyh,
8  2 ialth, ipup, lmax, meo, nqnyh, nslp
9  INTEGER ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
10  1 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
11  INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
12  DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM
13  DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
14  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
15  DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
16  1 r, rh, rhdn, rhsm, rhup, told, vnorm
17  dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
18  1 acor(*), wm(*), iwm(*)
19  COMMON /ls0001/ conit, crate, el(13), elco(13,12),
20  1 hold, rmax, tesco(3,12),
21  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
22  2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
23  3 mxstep, mxhnil, nhnil, ntrep, nslast, cnyh,
24  3 ialth, ipup, lmax, meo, nqnyh, nslp,
25  4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
26  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
27 C-----------------------------------------------------------------------
28 C STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
29 C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
30 C NOTE.. STODE IS INDEPENDENT OF THE VALUE OF THE ITERATION METHOD
31 C INDICATOR MITER, WHEN THIS IS .NE. 0, AND HENCE IS INDEPENDENT
32 C OF THE TYPE OF CHORD METHOD USED, OR THE JACOBIAN STRUCTURE.
33 C COMMUNICATION WITH STODE IS DONE WITH THE FOLLOWING VARIABLES..
34 C
35 C NEQ = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND
36 C PASSED AS THE NEQ ARGUMENT IN ALL CALLS TO F AND JAC.
37 C Y = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN
38 C ALL CALLS TO F AND JAC.
39 C YH = AN NYH BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
40 C AND THEIR APPROXIMATE SCALED DERIVATIVES, WHERE
41 C LMAX = MAXORD + 1. YH(I,J+1) CONTAINS THE APPROXIMATE
42 C J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J)
43 C (J = 0,1,...,NQ). ON ENTRY FOR THE FIRST STEP, THE FIRST
44 C TWO COLUMNS OF YH MUST BE SET FROM THE INITIAL VALUES.
45 C NYH = A CONSTANT INTEGER .GE. N, THE FIRST DIMENSION OF YH.
46 C YH1 = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH.
47 C EWT = AN ARRAY OF LENGTH N CONTAINING MULTIPLICATIVE WEIGHTS
48 C FOR LOCAL ERROR MEASUREMENTS. LOCAL ERRORS IN Y(I) ARE
49 C COMPARED TO 1.0/EWT(I) IN VARIOUS ERROR TESTS.
50 C SAVF = AN ARRAY OF WORKING STORAGE, OF LENGTH N.
51 C ALSO USED FOR INPUT OF YH(*,MAXORD+2) WHEN JSTART = -1
52 C AND MAXORD .LT. THE CURRENT ORDER NQ.
53 C ACOR = A WORK ARRAY OF LENGTH N, USED FOR THE ACCUMULATED
54 C CORRECTIONS. ON A SUCCESSFUL RETURN, ACOR(I) CONTAINS
55 C THE ESTIMATED ONE-STEP LOCAL ERROR IN Y(I).
56 C WM,IWM = REAL AND INTEGER WORK ARRAYS ASSOCIATED WITH MATRIX
57 C OPERATIONS IN CHORD ITERATION (MITER .NE. 0).
58 C PJAC = NAME OF ROUTINE TO EVALUATE AND PREPROCESS JACOBIAN MATRIX
59 C AND P = I - H*EL0*JAC, IF A CHORD METHOD IS BEING USED.
60 C SLVS = NAME OF ROUTINE TO SOLVE LINEAR SYSTEM IN CHORD ITERATION.
61 C CCMAX = MAXIMUM RELATIVE CHANGE IN H*EL0 BEFORE PJAC IS CALLED.
62 C H = THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
63 C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
64 C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
65 C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
66 C HMIN = THE MINIMUM ABSOLUTE VALUE OF THE STEP SIZE H TO BE USED.
67 C HMXI = INVERSE OF THE MAXIMUM ABSOLUTE VALUE OF H TO BE USED.
68 C HMXI = 0.0 IS ALLOWED AND CORRESPONDS TO AN INFINITE HMAX.
69 C HMIN AND HMXI MAY BE CHANGED AT ANY TIME, BUT WILL NOT
70 C TAKE EFFECT UNTIL THE NEXT CHANGE OF H IS CONSIDERED.
71 C TN = THE INDEPENDENT VARIABLE. TN IS UPDATED ON EACH STEP TAKEN.
72 C JSTART = AN INTEGER USED FOR INPUT ONLY, WITH THE FOLLOWING
73 C VALUES AND MEANINGS..
74 C 0 PERFORM THE FIRST STEP.
75 C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST.
76 C -1 TAKE THE NEXT STEP WITH A NEW VALUE OF H, MAXORD,
77 C N, METH, MITER, AND/OR MATRIX PARAMETERS.
78 C -2 TAKE THE NEXT STEP WITH A NEW VALUE OF H,
79 C BUT WITH OTHER INPUTS UNCHANGED.
80 C ON RETURN, JSTART IS SET TO 1 TO FACILITATE CONTINUATION.
81 C KFLAG = A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
82 C 0 THE STEP WAS SUCCESFUL.
83 C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED.
84 C -2 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED.
85 C -3 FATAL ERROR IN PJAC OR SLVS.
86 C A RETURN WITH KFLAG = -1 OR -2 MEANS EITHER
87 C ABS(H) = HMIN OR 10 CONSECUTIVE FAILURES OCCURRED.
88 C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF TN AND
89 C THE YH ARRAY ARE AS OF THE BEGINNING OF THE LAST
90 C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
91 C MAXORD = THE MAXIMUM ORDER OF INTEGRATION METHOD TO BE ALLOWED.
92 C MAXCOR = THE MAXIMUM NUMBER OF CORRECTOR ITERATIONS ALLOWED.
93 C MSBP = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0).
94 C MXNCF = MAXIMUM NUMBER OF CONVERGENCE FAILURES ALLOWED.
95 C METH/MITER = THE METHOD FLAGS. SEE DESCRIPTION IN DRIVER.
96 C N = THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
97 C IERR = ERROR FLAG FROM USER-SUPPLIED FUNCTION
98 C-----------------------------------------------------------------------
99  kflag = 0
100  told = tn
101  ncf = 0
102  ierpj = 0
103  iersl = 0
104  jcur = 0
105  icf = 0
106  delp = 0.0d0
107  IF (jstart .GT. 0) go to 200
108  IF (jstart .EQ. -1) go to 100
109  IF (jstart .EQ. -2) go to 160
110 C-----------------------------------------------------------------------
111 C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE
112 C INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
113 C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
114 C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE
115 C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
116 C FOR THE NEXT INCREASE.
117 C-----------------------------------------------------------------------
118  lmax = maxord + 1
119  nq = 1
120  l = 2
121  ialth = 2
122  rmax = 10000.0d0
123  rc = 0.0d0
124  el0 = 1.0d0
125  crate = 0.7d0
126  hold = h
127  meo = meth
128  nslp = 0
129  ipup = miter
130  iret = 3
131  go to 140
132 C-----------------------------------------------------------------------
133 C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1.
134 C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE.
135 C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1),
136 C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP.
137 C IF THE CALLER HAS CHANGED METH, CFODE IS CALLED TO RESET
138 C THE COEFFICIENTS OF THE METHOD.
139 C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
140 C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY.
141 C IF H IS TO BE CHANGED, YH MUST BE RESCALED.
142 C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1
143 C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS.
144 C-----------------------------------------------------------------------
145  100 ipup = miter
146  lmax = maxord + 1
147  IF (ialth .EQ. 1) ialth = 2
148  IF (meth .EQ. meo) go to 110
149  CALL cfode(meth, elco, tesco)
150  meo = meth
151  IF (nq .GT. maxord) go to 120
152  ialth = l
153  iret = 1
154  go to 150
155  110 IF (nq .LE. maxord) go to 160
156  120 nq = maxord
157  l = lmax
158  DO 125 i = 1,l
159  125 el(i) = elco(i,nq)
160  nqnyh = nq*nyh
161  rc = rc*el(1)/el0
162  el0 = el(1)
163  conit = 0.5d0/dble(nq+2)
164  ddn = vnorm(n, savf, ewt)/tesco(1,l)
165  exdn = 1.0d0/dble(l)
166  rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
167  rh = dmin1(rhdn,1.0d0)
168  iredo = 3
169  IF (h .EQ. hold) go to 170
170  rh = dmin1(rh,dabs(h/hold))
171  h = hold
172  go to 175
173 C-----------------------------------------------------------------------
174 C CFODE IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE
175 C CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET
176 C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM.
177 C-----------------------------------------------------------------------
178  140 CALL cfode(meth, elco, tesco)
179  150 DO 155 i = 1,l
180  155 el(i) = elco(i,nq)
181  nqnyh = nq*nyh
182  rc = rc*el(1)/el0
183  el0 = el(1)
184  conit = 0.5d0/dble(nq+2)
185  go to(160, 170, 200), iret
186 C-----------------------------------------------------------------------
187 C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
188 C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO
189 C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS
190 C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE.
191 C-----------------------------------------------------------------------
192  160 IF (h .EQ. hold) go to 200
193  rh = h/hold
194  h = hold
195  iredo = 3
196  go to 175
197  170 rh = dmax1(rh,hmin/dabs(h))
198  175 rh = dmin1(rh,rmax)
199  rh = rh/dmax1(1.0d0,dabs(h)*hmxi*rh)
200  r = 1.0d0
201  DO 180 j = 2,l
202  r = r*rh
203  DO 180 i = 1,n
204  180 yh(i,j) = yh(i,j)*r
205  h = h*rh
206  rc = rc*rh
207  ialth = l
208  IF (iredo .EQ. 0) go to 690
209 C-----------------------------------------------------------------------
210 C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
211 C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX.
212 C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1).
213 C WHEN RC DIFFERS FROM 1 BY MORE THAN CCMAX, IPUP IS SET TO MITER
214 C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED.
215 C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS.
216 C-----------------------------------------------------------------------
217  200 IF (dabs(rc-1.0d0) .GT. ccmax) ipup = miter
218  IF (nst .GE. nslp+msbp) ipup = miter
219  tn = tn + h
220  i1 = nqnyh + 1
221  DO 215 jb = 1,nq
222  i1 = i1 - nyh
223 CDIR$ IVDEP
224  DO 210 i = i1,nqnyh
225  210 yh1(i) = yh1(i) + yh1(i+nyh)
226  215 CONTINUE
227 C-----------------------------------------------------------------------
228 C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS
229 C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR
230 C WEIGHT VECTOR EWT. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE
231 C VECTOR ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP.
232 C-----------------------------------------------------------------------
233  220 m = 0
234  DO 230 i = 1,n
235  230 y(i) = yh(i,1)
236  ierr = 0
237  CALL f(neq, tn, y, savf, ierr)
238  IF (ierr .LT. 0) RETURN
239  nfe = nfe + 1
240  IF (ipup .LE. 0) go to 250
241 C-----------------------------------------------------------------------
242 C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND
243 C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET
244 C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE.
245 C-----------------------------------------------------------------------
246  ierr = 0
247  CALL pjac(neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac,
248  1 ierr)
249  IF (ierr .LT. 0) RETURN
250  ipup = 0
251  rc = 1.0d0
252  nslp = nst
253  crate = 0.7d0
254  IF (ierpj .NE. 0) go to 430
255  250 DO 260 i = 1,n
256  260 acor(i) = 0.0d0
257  270 IF (miter .NE. 0) go to 350
258 C-----------------------------------------------------------------------
259 C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM
260 C THE RESULT OF THE LAST FUNCTION EVALUATION.
261 C-----------------------------------------------------------------------
262  DO 290 i = 1,n
263  savf(i) = h*savf(i) - yh(i,2)
264  290 y(i) = savf(i) - acor(i)
265  del = vnorm(n, y, ewt)
266  DO 300 i = 1,n
267  y(i) = yh(i,1) + el(1)*savf(i)
268  300 acor(i) = savf(i)
269  go to 400
270 C-----------------------------------------------------------------------
271 C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
272 C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND
273 C P AS COEFFICIENT MATRIX.
274 C-----------------------------------------------------------------------
275  350 DO 360 i = 1,n
276  360 y(i) = h*savf(i) - (yh(i,2) + acor(i))
277  CALL slvs(wm, iwm, y, savf)
278  IF (iersl .LT. 0) go to 430
279  IF (iersl .GT. 0) go to 410
280  del = vnorm(n, y, ewt)
281  DO 380 i = 1,n
282  acor(i) = acor(i) + y(i)
283  380 y(i) = yh(i,1) + el(1)*acor(i)
284 C-----------------------------------------------------------------------
285 C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
286 C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
287 C-----------------------------------------------------------------------
288  400 IF (m .NE. 0) crate = dmax1(0.2d0*crate,del/delp)
289  dcon = del*dmin1(1.0d0,1.5d0*crate)/(tesco(2,nq)*conit)
290  IF (dcon .LE. 1.0d0) go to 450
291  m = m + 1
292  IF (m .EQ. maxcor) go to 410
293  IF (m .GE. 2 .AND. del .GT. 2.0d0*delp) go to 410
294  delp = del
295  ierr = 0
296  CALL f(neq, tn, y, savf, ierr)
297  IF (ierr .LT. 0) RETURN
298  nfe = nfe + 1
299  go to 270
300 C-----------------------------------------------------------------------
301 C THE CORRECTOR ITERATION FAILED TO CONVERGE.
302 C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR
303 C THE NEXT TRY. OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES
304 C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF H CANNOT BE
305 C REDUCED OR MXNCF FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2.
306 C-----------------------------------------------------------------------
307  410 IF (miter .EQ. 0 .OR. jcur .EQ. 1) go to 430
308  icf = 1
309  ipup = miter
310  go to 220
311  430 icf = 2
312  ncf = ncf + 1
313  rmax = 2.0d0
314  tn = told
315  i1 = nqnyh + 1
316  DO 445 jb = 1,nq
317  i1 = i1 - nyh
318 CDIR$ IVDEP
319  DO 440 i = i1,nqnyh
320  440 yh1(i) = yh1(i) - yh1(i+nyh)
321  445 CONTINUE
322  IF (ierpj .LT. 0 .OR. iersl .LT. 0) go to 680
323  IF (dabs(h) .LE. hmin*1.00001d0) go to 670
324  IF (ncf .EQ. mxncf) go to 670
325  rh = 0.25d0
326  ipup = miter
327  iredo = 1
328  go to 170
329 C-----------------------------------------------------------------------
330 C THE CORRECTOR HAS CONVERGED. JCUR IS SET TO 0
331 C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER.
332 C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500
333 C IF IT FAILS.
334 C-----------------------------------------------------------------------
335  450 jcur = 0
336  IF (m .EQ. 0) dsm = del/tesco(2,nq)
337  IF (m .GT. 0) dsm = vnorm(n, acor, ewt)/tesco(2,nq)
338  IF (dsm .GT. 1.0d0) go to 500
339 C-----------------------------------------------------------------------
340 C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY.
341 C CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1.
342 C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR
343 C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
344 C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
345 C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A
346 C FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT
347 C TESTING FOR THAT MANY STEPS.
348 C-----------------------------------------------------------------------
349  kflag = 0
350  iredo = 0
351  nst = nst + 1
352  hu = h
353  nqu = nq
354  DO 470 j = 1,l
355  DO 470 i = 1,n
356  470 yh(i,j) = yh(i,j) + el(j)*acor(i)
357  ialth = ialth - 1
358  IF (ialth .EQ. 0) go to 520
359  IF (ialth .GT. 1) go to 700
360  IF (l .EQ. lmax) go to 700
361  DO 490 i = 1,n
362  490 yh(i,lmax) = acor(i)
363  go to 700
364 C-----------------------------------------------------------------------
365 C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
366 C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
367 C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
368 C ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE
369 C BY A FACTOR OF 0.2 OR LESS.
370 C-----------------------------------------------------------------------
371  500 kflag = kflag - 1
372  tn = told
373  i1 = nqnyh + 1
374  DO 515 jb = 1,nq
375  i1 = i1 - nyh
376 CDIR$ IVDEP
377  DO 510 i = i1,nqnyh
378  510 yh1(i) = yh1(i) - yh1(i+nyh)
379  515 CONTINUE
380  rmax = 2.0d0
381  IF (dabs(h) .LE. hmin*1.00001d0) go to 660
382  IF (kflag .LE. -3) go to 640
383  iredo = 2
384  rhup = 0.0d0
385  go to 540
386 C-----------------------------------------------------------------------
387 C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
388 C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED
389 C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
390 C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE.
391 C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
392 C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
393 C ADDITIONAL SCALED DERIVATIVE.
394 C-----------------------------------------------------------------------
395  520 rhup = 0.0d0
396  IF (l .EQ. lmax) go to 540
397  DO 530 i = 1,n
398  530 savf(i) = acor(i) - yh(i,lmax)
399  dup = vnorm(n, savf, ewt)/tesco(3,nq)
400  exup = 1.0d0/dble(l+1)
401  rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0)
402  540 exsm = 1.0d0/dble(l)
403  rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
404  rhdn = 0.0d0
405  IF (nq .EQ. 1) go to 560
406  ddn = vnorm(n, yh(1,l), ewt)/tesco(1,nq)
407  exdn = 1.0d0/dble(nq)
408  rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
409  560 IF (rhsm .GE. rhup) go to 570
410  IF (rhup .GT. rhdn) go to 590
411  go to 580
412  570 IF (rhsm .LT. rhdn) go to 580
413  newq = nq
414  rh = rhsm
415  go to 620
416  580 newq = nq - 1
417  rh = rhdn
418  IF (kflag .LT. 0 .AND. rh .GT. 1.0d0) rh = 1.0d0
419  go to 620
420  590 newq = l
421  rh = rhup
422  IF (rh .LT. 1.1d0) go to 610
423  r = el(l)/dble(l)
424  DO 600 i = 1,n
425  600 yh(i,newq+1) = acor(i)*r
426  go to 630
427  610 ialth = 3
428  go to 700
429  620 IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1d0)) go to 610
430  IF (kflag .LE. -2) rh = dmin1(rh,0.2d0)
431 C-----------------------------------------------------------------------
432 C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
433 C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED.
434 C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
435 C-----------------------------------------------------------------------
436  IF (newq .EQ. nq) go to 170
437  630 nq = newq
438  l = nq + 1
439  iret = 2
440  go to 150
441 C-----------------------------------------------------------------------
442 C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURRED.
443 C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1.
444 C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
445 C YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
446 C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN
447 C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED,
448 C UNTIL IT SUCCEEDS OR H REACHES HMIN.
449 C-----------------------------------------------------------------------
450  640 IF (kflag .EQ. -10) go to 660
451  rh = 0.1d0
452  rh = dmax1(hmin/dabs(h),rh)
453  h = h*rh
454  DO 645 i = 1,n
455  645 y(i) = yh(i,1)
456  ierr = 0
457  CALL f(neq, tn, y, savf, ierr)
458  IF (ierr .LT. 0) RETURN
459  nfe = nfe + 1
460  DO 650 i = 1,n
461  650 yh(i,2) = h*savf(i)
462  ipup = miter
463  ialth = 5
464  IF (nq .EQ. 1) go to 200
465  nq = 1
466  l = 2
467  iret = 3
468  go to 150
469 C-----------------------------------------------------------------------
470 C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD
471 C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
472 C-----------------------------------------------------------------------
473  660 kflag = -1
474  go to 720
475  670 kflag = -2
476  go to 720
477  680 kflag = -3
478  go to 720
479  690 rmax = 10.0d0
480  700 r = 1.0d0/tesco(2,nqu)
481  DO 710 i = 1,n
482  710 acor(i) = acor(i)*r
483  720 hold = h
484  jstart = 1
485  RETURN
486 C----------------------- END OF SUBROUTINE STODE -----------------------
487  END
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
subroutine stode(NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, F, JAC, PJAC, SLVS, IERR)
Definition: stode.f:1
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
subroutine cfode(METH, ELCO, TESCO)
Definition: cfode.f:1
double precision function vnorm(N, V, W)
Definition: vnorm.f:1