1
subroutine rkf45(fydot,neqn,y,t,tout,itol,rerr,aerr,
2
1 itask,iflag,iopt,work,lrw,iwork,liw,bjac,mf)
4
c fehlberg fourth-fifth order runge-kutta method
6
c written by h.a.watts and l.f.shampine
8
c albuquerque,new mexico
10
c rkf45 is primarily designed to solve non-stiff and mildly stiff
11
c differential equations when derivative evaluations are inexpensive.
12
c rkf45 should generally not be used when the user is demanding
17
c subroutine rkf45 integrates a system of neqn first order
18
c ordinary differential equations of the form
19
c dy(i)/dt = fydot(t,y(1),y(2),...,y(neqn))
20
c where the y(i) are given at t .
21
c typically the subroutine is used to integrate from t to tout but it
22
c can be used as a one-step integrator to advance the solution a
23
c single step in the direction of tout. on return the parameters in
24
c the call list are set for continuing the integration. the user has
25
c only to call rkf45 again (and perhaps define a new value for tout).
26
c actually, rkf45 is an interfacing routine which calls subroutine
27
c rkfs for the solution. rkfs in turn calls subroutine fehl which
28
c computes an approximate solution over one step.
30
c rkf45 uses the runge-kutta-fehlberg (4,5) method described
32
c e.fehlberg , low-order classical runge-kutta formulas with stepsize
33
c control , nasa tr r-315
35
c the performance of rkf45 is illustrated in the reference
36
c l.f.shampine,h.a.watts,s.davenport, solving non-stiff ordinary
37
c differential equations-the state of the art ,
38
c sandia laboratories report sand75-0182 ,
39
c to appear in siam review.
42
c the parameters represent-
43
c fydot -- subroutine fydot(neqn,t,y,yp) to evaluate derivatives
45
c neqn -- number of equations to be integrated
46
c y(*) -- solution vector at t
47
c t -- independent variable
48
c tout -- output point at which solution is desired
49
c rerr,aerr -- relative and absolute error tolerances for local
50
c error test. at each step the code requires that
51
c abs(local error) .le. rerr*abs(y) + aerr
52
c for each component of the local error and solution vectors
53
c iflag -- indicator for status of integration
54
c work(*) -- array to hold information internal to rkf45 which is
55
c necessary for subsequent calls. must be dimensioned
57
c iwork(*) -- integer array used to hold information internal to
58
c rkf45 which is necessary for subsequent calls. must be
59
c dimensioned at least 5
64
c the user must provide storage in his calling program for the arrays
65
c in the call list - y(neqn) , work(3+6*neqn) , iwork(5) ,
66
c declare fydot in an external statement, supply
67
c subroutine fydot(neqn,t,y,yp)
68
c and initialize the following parameters-
70
c neqn -- number of equations to be integrated. (neqn .ge. 1)
71
c y(*) -- vector of initial conditions
72
c t -- starting point of integration , must be a variable
73
c tout -- output point at which solution is desired.
74
c t=tout is allowed on the first call only, in which case
75
c rkf45 returns with iflag=2 if continuation is possible.
76
c rerr,aerr -- relative and absolute local error tolerances
77
c which must be non-negative. rerr must be a variable while
78
c aerr may be a constant. the code should normally not be
79
c used with relative error control smaller than about 1.e-8 .
80
c to avoid limiting precision difficulties the code requires
81
c rerr to be larger than an internally computed relative
82
c error parameter which is machine dependent. in particular,
83
c pure absolute error is not permitted. if a smaller than
84
c allowable value of rerr is attempted, rkf45 increases
85
c rerr appropriately and returns control to the user before
86
c continuing the integration.
87
c iflag -- +1,-1 indicator to initialize the code for each new
88
c problem. normal input is +1. the user should set iflag=-1
89
c only when one-step integrator control is essential. in this
90
c case, rkf45 attempts to advance the solution a single step
91
c in the direction of tout each time it is called. since this
92
c mode of operation results in extra computing overhead, it
93
c should be avoided unless needed.
98
c y(*) -- solution at t
99
c t -- last point reached in integration.
100
c iflag = 2 -- integration reached tout. indicates successful retur
101
c and is the normal mode for continuing integration.
102
c =-2 -- a single successful step in the direction of tout
103
c has been taken. normal mode for continuing
104
c integration one step at a time.
105
c = 3 -- integration was not completed because relative error
106
c tolerance was too small. rerr has been increased
107
c appropriately for continuing.
108
c = 4 -- integration was not completed because more than
109
c 3000 derivative evaluations were needed. this
110
c is approximately 500 steps.
111
c = 5 -- integration was not completed because solution
112
c vanished making a pure relative error test
113
c impossible. must use non-zero aerr to continue.
114
c using the one-step integration mode for one step
115
c is a good way to proceed.
116
c = 6 -- integration was not completed because requested
117
c accuracy could not be achieved using smallest
118
c allowable stepsize. user must increase the error
119
c tolerance before continued integration can be
121
c = 7 -- it is likely that rkf45 is inefficient for solving
122
c this problem. too much output is restricting the
123
c natural stepsize choice. use the one-step integrator
125
c = 8 -- invalid input parameters
126
c this indicator occurs if any of the following is
127
c satisfied - neqn .le. 0
128
c t=tout and iflag .ne. +1 or -1
129
c rerr or aerr .lt. 0.
130
c iflag .eq. 0 or .lt. -2 or .gt. 8
131
c work(*),iwork(*) -- information which is usually of no interest
132
c to the user but necessary for subsequent calls.
133
c work(1),...,work(neqn) contain the first derivatives
134
c of the solution vector y at t. work(neqn+1) contains
135
c the stepsize h to be attempted on the next step.
136
c iwork(1) contains the derivative evaluation counter.
139
c subsequent calls to rkf45
141
c subroutine rkf45 returns with all information needed to continue
142
c the integration. if the integration reached tout, the user need onl
143
c define a new tout and call rkf45 again. in the one-step integrator
144
c mode (iflag=-2) the user must keep in mind that each step taken is
145
c in the direction of the current tout. upon reaching tout (indicated
146
c by changing iflag to 2),the user must then define a new tout and
147
c reset iflag to -2 to continue in the one-step integrator mode.
149
c if the integration was not completed but the user still wants to
150
c continue (iflag=3,4 cases), he just calls rkf45 again. with iflag=3
151
c the rerr parameter has been adjusted appropriately for continuing
152
c the integration. in the case of iflag=4 the function counter will
153
c be reset to 0 and another 3000 function evaluations are allowed.
155
c however,in the case iflag=5, the user must first alter the error
156
c criterion to use a positive value of aerr before integration can
157
c proceed. if he does not,execution is terminated.
159
c also,in the case iflag=6, it is necessary for the user to reset
160
c iflag to 2 (or -2 when the one-step integration mode is being used)
161
c as well as increasing either aerr,rerr or both before the
162
c integration can be continued. if this is not done, execution will
163
c be terminated. the occurrence of iflag=6 indicates a trouble spot
164
c (solution is changing rapidly,singularity may be present) and it
165
c often is inadvisable to continue.
167
c if iflag=7 is encountered, the user should use the one-step
168
c integration mode with the stepsize determined by the code or
169
c consider switching to the adams codes de/step,intrp. if the user
170
c insists upon continuing the integration with rkf45, he must reset
171
c iflag to 2 before calling rkf45 again. otherwise,execution will be
174
c if iflag=8 is obtained, integration can not be continued unless
175
c the invalid input parameters are corrected.
177
c it should be noted that the arrays work,iwork contain information
178
c required for subsequent integration. accordingly, work and iwork
179
c should not be altered.
182
integer neqn,iflag,iwork(5)
183
double precision y(neqn),t,tout,rerr,aerr,work(1)
188
integer k1,k2,k3,k4,k5,k6,k1m
189
if(itask.eq.1) then iflag=1
190
if(itask.eq.2) then iflag=-1
193
c compute indices for the splitting of the work array
204
c this interfacing routine merely relieves the user of a long
205
c calling list via the splitting apart of two working storage
206
c arrays. if this is not compatible with the users compiler,
207
c he must use rkfs directly.
209
call rkfs(fydot,neqn,y,t,tout,rerr,aerr,iflag,work(1),work(k1m),
210
1 work(k1),work(k2),work(k3),work(k4),work(k5),work(k6),
211
2 work(k6+1),work(k7),
212
3 iwork(1),iwork(2),iwork(3),iwork(4),iwork(5))
216
subroutine rkfs(fydot,neqn,y,t,tout,rerr,aerr,iflag,yp,h,f1,f2,f3,
217
1 f4,f5,savre,savae,savey,nfe,kop,init,jflag,kflag)
219
c fehlberg fourth-fifth order runge-kutta method
222
c rkfs integrates a system of first order ordinary differential
223
c equations as described in the comments for rkf45 .
224
c the arrays yp,f1,f2,f3,f4,and f5 (of dimension at least neqn) and
225
c the variables h,savre,savae,nfe,kop,init,jflag,and kflag are used
226
c internally by the code and appear in the call list to eliminate
227
c local retention of variables between calls. accordingly, they
228
c should not be altered. items of possible interest are
229
c yp - derivative of solution vector at t
230
c h - an appropriate stepsize to be used for the next step
231
c nfe- counter on the number of derivative function evaluations
234
logical hfaild,output
236
integer neqn,iflag,nfe,kop,init,jflag,kflag
237
double precision y(neqn),t,tout,rerr,aerr,h,yp(neqn),
238
1 f1(neqn),f2(neqn),f3(neqn),f4(neqn),f5(neqn),savre,
244
double precision a,ae,dt,ee,eeoet,esttol,et,hmin,remin,rer,s,
245
1 scale,tol,toln,twoeps,u26,ypk
247
integer k,maxnfe,mflag
249
double precision dabs,dmax1,dmin1,dsign,dlamch
251
c remin is the minimum acceptable value of rerr. attempts
252
c to obtain higher accuracy with this subroutine are usually
253
c very expensive and often unsuccessful.
258
c the expense is controlled by restricting the number
259
c of function evaluations to be approximately maxnfe.
260
c as set, this corresponds to about 500 steps.
264
c here two constants emboding the machine epsilon is present
265
c twoesp is set to twice the machine epsilon while u26 is set
266
c to 26 times the machine epsilon
268
c data twoeps, u26/4.4d-16, 5.72d-15/
269
twoeps = 2.*dlamch('p')
273
c check input parameters
276
if (neqn .lt. 1) go to 10
277
if ((rerr .lt. 0.0d0) .or. (aerr .lt. 0.0d0)) go to 10
279
if ((mflag .ge. 1) .and. (mflag .le. 8)) go to 20
285
c is this the first call
286
20 if (mflag .eq. 1) go to 50
288
c check continuation possibilities
290
if ((t .eq. tout) .and. (kflag .ne. 3)) go to 10
291
if (mflag .ne. 2) go to 25
294
if (kflag .eq. 3) go to 45
295
if (init .eq. 0) go to 45
296
if (kflag .eq. 4) go to 40
297
if ((kflag .eq. 5) .and. (aerr .eq. 0.0d0)) go to 30
298
if ((kflag .eq. 6) .and. (rerr .le. savre) .and.
299
1 (aerr .le. savae)) go to 30
302
c iflag = 3,4,5,6,7 or 8
303
25 if (iflag .eq. 3) go to 45
304
if (iflag .eq. 4) go to 40
305
if ((iflag .eq. 5) .and. (aerr .gt. 0.0d0)) go to 45
307
c integration cannot be continued since user did not respond to
308
c the instructions pertaining to iflag=5,6,7 or 8
311
c reset function evaluation counter
313
if (mflag .eq. 2) go to 50
315
c reset flag value from previous call
317
if (kflag .eq. 3) mflag=iabs(iflag)
319
c save input iflag and set continuation flag value for subsequent
324
c save rerr and aerr for checking input on subsequent calls
328
c restrict relative error tolerance to be at least as large as
329
c 2*eps+remin to avoid limiting precision difficulties arising
330
c from impossible accuracy requests
333
if (rerr .ge. rer) go to 55
335
c relative error tolerance too small
343
if (mflag .eq. 1) go to 60
344
if (init .eq. 0) go to 65
348
c set initialization completion indicator,init
349
c set indicator for too many output points,kop
350
c evaluate initial derivatives
351
c set counter for function evaluations,nfe
352
c evaluate initial derivatives
353
c set counter for function evaluations,nfe
354
c estimate starting stepsize
360
call fydot(neqn,a,y,yp)
363
if (t .ne. tout) go to 65
372
tol=rerr*dabs(y(k))+aerr
373
if (tol .le. 0.) go to 70
376
if (ypk*h**5 .gt. tol) h=(tol/ypk)**0.2d0
378
if (toln .le. 0.0d0) h=0.0d0
379
h=dmax1(h,u26*dmax1(dabs(t),dabs(dt)))
383
c set stepsize for integration in the direction from t to tout
387
c test to see if rkf45 is being severely impacted by too many
390
if (dabs(h) .ge. 2.0d0*dabs(dt)) kop=kop+1
391
if (kop .ne. 100) go to 85
393
c unnecessary frequency of output
398
85 if (dabs(dt) .gt. u26*dabs(t)) go to 95
400
c if too close to output point,extrapolate and return
403
90 y(k)=y(k)+dt*yp(k)
405
call fydot(neqn,a,y,yp)
411
c initialize output point indicator
415
c to avoid premature underflow in the error tolerance function,
416
c scale the error tolerances
422
c step by step integration
426
c set smallest allowable stepsize
430
c adjust stepsize if necessary to hit the output point.
431
c look ahead two steps to avoid drastic changes in the stepsize and
432
c thus lessen the impact of output points on the code.
435
if (dabs(dt) .ge. 2.0d0*dabs(h)) go to 200
436
if (dabs(dt) .gt. dabs(h)) go to 150
438
c the next successful step will complete the integration to the
449
c core integrator for taking a single step
451
c the tolerances have been scaled to avoid premature underflow in
452
c computing the error tolerance function et.
453
c to avoid problems with zero crossings,relative error is measured
454
c using the average of the magnitudes of the solution at the
455
c beginning and end of a step.
456
c the error estimate formula has been grouped to control loss of
458
c to distinguish the various arguments, h is not permitted
459
c to become smaller than 26 units of roundoff in t.
460
c practical limits on the change in the stepsize are enforced to
461
c smooth the stepsize selection process and to avoid excessive
462
c chattering on problems having discontinuities.
463
c to prevent unnecessary failures, the code uses 9/10 the stepsize
464
c it estimates will succeed.
465
c after a step failure, the stepsize is not allowed to increase for
466
c the next attempted step. this makes the code more efficient on
467
c problems having discontinuities and more effective in general
468
c since local extrapolation is being used and extra caution seems
472
c test number of derivative function evaluations.
473
c if okay,try to advance the integration from t to t+h
475
200 if (nfe .le. maxnfe) go to 220
482
c advance an approximate solution over one step of length h
487
call fehl(fydot,neqn,y,t,h,yp,f1,f2,f3,f4,f5,f1,savey)
492
c compute and test allowable tolerances versus local error estimates
493
c and remove scaling of tolerances. note that relative error is
494
c measured with respect to the average of the magnitudes of the
495
c solution at the beginning and end of the step.
499
et=dabs(savey(k))+dabs(f1(k))+ae
500
if (et .gt. 0.0d0) go to 240
502
c inappropriate error tolerance
506
240 ee=dabs((-2090.0d0*yp(k)+(21970.0d0*f3(k)-15048.0d0*f4(k)))+
507
1 (22528.0d0*f2(k)-27360.0d0*f5(k)))
508
250 eeoet=dmax1(eeoet,ee/et)
510
esttol=dabs(h)*eeoet*scale/752400.0d0
512
if (esttol .le. 1.0d0) go to 260
516
c reduce the stepsize , try again
517
c the decrease is limited to a factor of 1/10
522
if (esttol .lt. 59049.0d0) s=0.9d0/esttol**0.2d0
524
if (dabs(h) .gt. hmin) go to 200
526
c requested error unattainable at smallest allowable stepsize
533
c store solution at t+h
534
c and evaluate derivatives there
540
call fydot(neqn,a,y,yp)
545
c choose next stepsize
546
c the increase is limited to a factor of 5
547
c if step failure has just occurred, next
548
c stepsize is not allowed to increase
551
if (esttol .gt. 1.889568d-4) s=0.9d0/esttol**0.2d0
552
if (hfaild) s=dmin1(s,1.0d0)
553
h=dsign(dmax1(s*dabs(h),hmin),h)
555
c end of core integrator
558
c should we take another step
560
if (output) go to 300
561
if (iflag .gt. 0) go to 100
564
c integration successfully completed
576
subroutine fehl(fydot,neqn,y,t,h,yp,f1,f2,f3,f4,f5,s,savey)
577
c subroutine fehl(fydot,neqn,y,t,h,yp,f1,f2,f3,f4,f5,s)
579
c fehlberg fourth-fifth order runge-kutta method
581
c fehl integrates a system of neqn first order
582
c ordinary differential equations of the form
583
c dy(i)/dt=fydot(t,y(1),---,y(neqn))
584
c where the initial values y(i) and the initial derivatives
585
c yp(i) are specified at the starting point t. fehl advances
586
c the solution over the fixed step h and returns
587
c the fifth order (sixth order accurate locally) solution
588
c approximation at t+h in array s(i).
589
c f1,---,f5 are arrays of dimension neqn which are needed
590
c for internal storage.
591
c the formulas have been grouped to control loss of significance.
592
c fehl should be called with an h not smaller than 13 units of
593
c roundoff in t so that the various independent arguments can be
598
double precision y(neqn),t,h,yp(neqn),f1(neqn),f2(neqn),
599
1 f3(neqn),f4(neqn),f5(neqn),s(neqn),savey(neqn)
608
221 y(k)=savey(k)+ch*yp(k)
609
call fydot(neqn,t+ch,y,f1)
614
222 y(k)=savey(k)+ch*(yp(k)+3.0d0*f1(k))
615
call fydot(neqn,t+3.0d0*h/8.0d0,y,f2)
620
223 y(k)=savey(k)+ch*(1932.0d0*yp(k)+
621
1 (7296.0d0*f2(k)-7200.0d0*f1(k)))
622
call fydot(neqn,t+12.0d0*h/13.0d0,y,f3)
627
224 y(k)=savey(k)+ch*((8341.0d0*yp(k)-845.0d0*f3(k))+
628
1 (29440.0d0*f2(k)-32832.0d0*f1(k)))
629
call fydot(neqn,t+h,y,f4)
634
225 y(k)=savey(k)+ch*((-6080.0d0*yp(k)+(9295.0d0*f3(k)-
635
1 5643.0d0*f4(k)))+(41040.0d0*f1(k)-28352.0d0*f2(k)))
636
call fydot(neqn,t+h/2.0d0,y,f5)
639
c compute approximate solution at t+h
643
230 s(k)=savey(k)+ch*((902880.0d0*yp(k)+(3855735.0d0*f3(k)-
644
1 1371249.0d0*f4(k)))+(3953664.0d0*f2(k)+
added
removed
routines/integ/rkf45.f
