subroutine rkf45(f,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
c
c fehlberg fourth-fifth order runge-kutta method
c
c written by h.a.watts and l.f.shampine
c sandia laboratories
c albuquerque,new mexico
c
c rkf45 is primarily designed to solve non-stiff and mildly stiff
c differential equations when derivative evaluations are inexpensive.
c rkf45 should generally not be used when the user is demanding
c high accuracy.
c
c abstract
c
c subroutine rkf45 integrates a system of neqn first order
c ordinary differential equations of the form
c dy(i)/dt = f(t,y(1),y(2),...,y(neqn))
c where the y(i) are given at t .
c typically the subroutine is used to integrate from t to tout but it
c can be used as a one-step integrator to advance the solution a
c single step in the direction of tout. on return the parameters in
c the call list are set for continuing the integration. the user has
c only to call rkf45 again (and perhaps define a new value for tout).
c actually, rkf45 is an interfacing routine which calls subroutine
c rkfs for the solution. rkfs in turn calls subroutine fehl which
c computes an approximate solution over one step.
c
c rkf45 uses the runge-kutta-fehlberg (4,5) method described
c in the reference
c e.fehlberg , low-order classical runge-kutta formulas with stepsize
c control , nasa tr r-315
c
c the performance of rkf45 is illustrated in the reference
c l.f.shampine,h.a.watts,s.davenport, solving non-stiff ordinary
c differential equations-the state of the art ,
c sandia laboratories report sand75-0182 ,
c to appear in siam review.
c
c
c the parameters represent-
c f -- subroutine f(t,y,yp) to evaluate derivatives yp(i)=dy(i)/dt
c neqn -- number of equations to be integrated
c y(*) -- solution vector at t
c t -- independent variable
c tout -- output point at which solution is desired
c relerr,abserr -- relative and absolute error tolerances for local
c error test. at each step the code requires that
c abs(local error) .le. relerr*abs(y) + abserr
c for each component of the local error and solution vectors
c iflag -- indicator for status of integration
c work(*) -- array to hold information internal to rkf45 which is
c necessary for subsequent calls. must be dimensioned
c at least 3+6*neqn
c iwork(*) -- integer array used to hold information internal to
c rkf45 which is necessary for subsequent calls. must be
c dimensioned at least 5
c
c
c first call to rkf45
c
c the user must provide storage in his calling program for the arrays
c in the call list - y(neqn) , work(3+6*neqn) , iwork(5) ,
c declare f in an external statement, supply subroutine f(t,y,yp) and
c initialize the following parameters-
c
c neqn -- number of equations to be integrated. (neqn .ge. 1)
c y(*) -- vector of initial conditions
c t -- starting point of integration , must be a variable
c tout -- output point at which solution is desired.
c t=tout is allowed on the first call only, in which case
c rkf45 returns with iflag=2 if continuation is possible.
c relerr,abserr -- relative and absolute local error tolerances
c which must be non-negative. relerr must be a variable while
c abserr may be a constant. the code should normally not be
c used with relative error control smaller than about 1.e-8 .
c to avoid limiting precision difficulties the code requires
c relerr to be larger than an internally computed relative
c error parameter which is machine dependent. in particular,
c pure absolute error is not permitted. if a smaller than
c allowable value of relerr is attempted, rkf45 increases
c relerr appropriately and returns control to the user before
c continuing the integration.
c iflag -- +1,-1 indicator to initialize the code for each new
c problem. normal input is +1. the user should set iflag=-1
c only when one-step integrator control is essential. in this
c case, rkf45 attempts to advance the solution a single step
c in the direction of tout each time it is called. since this
c mode of operation results in extra computing overhead, it
c should be avoided unless needed.
c
c
c output from rkf45
c
c y(*) -- solution at t
c t -- last point reached in integration.
c iflag = 2 -- integration reached tout. indicates successful retur
c and is the normal mode for continuing integration.
c =-2 -- a single successful step in the direction of tout
c has been taken. normal mode for continuing
c integration one step at a time.
c = 3 -- integration was not completed because relative error
c tolerance was too small. relerr has been increased
c appropriately for continuing.
c = 4 -- integration was not completed because more than
c 3000 derivative evaluations were needed. this
c is approximately 500 steps.
c = 5 -- integration was not completed because solution
c vanished making a pure relative error test
c impossible. must use non-zero abserr to continue.
c using the one-step integration mode for one step
c is a good way to proceed.
c = 6 -- integration was not completed because requested
c accuracy could not be achieved using smallest
c allowable stepsize. user must increase the error
c tolerance before continued integration can be
c attempted.
c = 7 -- it is likely that rkf45 is inefficient for solving
c this problem. too much output is restricting the
c natural stepsize choice. use the one-step integrator
c mode.
c = 8 -- invalid input parameters
c this indicator occurs if any of the following is
c satisfied - neqn .le. 0
c t=tout and iflag .ne. +1 or -1
c relerr or abserr .lt. 0.
c iflag .eq. 0 or .lt. -2 or .gt. 8
c work(*),iwork(*) -- information which is usually of no interest
c to the user but necessary for subsequent calls.
c work(1),...,work(neqn) contain the first derivatives
c of the solution vector y at t. work(neqn+1) contains
c the stepsize h to be attempted on the next step.
c iwork(1) contains the derivative evaluation counter.
c
c
c subsequent calls to rkf45
c
c subroutine rkf45 returns with all information needed to continue
c the integration. if the integration reached tout, the user need onl
c define a new tout and call rkf45 again. in the one-step integrator
c mode (iflag=-2) the user must keep in mind that each step taken is
c in the direction of the current tout. upon reaching tout (indicated
c by changing iflag to 2),the user must then define a new tout and
c reset iflag to -2 to continue in the one-step integrator mode.
c
c if the integration was not completed but the user still wants to
c continue (iflag=3,4 cases), he just calls rkf45 again. with iflag=3
c the relerr parameter has been adjusted appropriately for continuing
c the integration. in the case of iflag=4 the function counter will
c be reset to 0 and another 3000 function evaluations are allowed.
c
c however,in the case iflag=5, the user must first alter the error
c criterion to use a positive value of abserr before integration can
c proceed. if he does not,execution is terminated.
c
c also,in the case iflag=6, it is necessary for the user to reset
c iflag to 2 (or -2 when the one-step integration mode is being used)
c as well as increasing either abserr,relerr or both before the
c integration can be continued. if this is not done, execution will
c be terminated. the occurrence of iflag=6 indicates a trouble spot
c (solution is changing rapidly,singularity may be present) and it
c often is inadvisable to continue.
c
c if iflag=7 is encountered, the user should use the one-step
c integration mode with the stepsize determined by the code or
c consider switching to the adams codes de/step,intrp. if the user
c insists upon continuing the integration with rkf45, he must reset
c iflag to 2 before calling rkf45 again. otherwise,execution will be
c terminated.
c
c if iflag=8 is obtained, integration can not be continued unless
c the invalid input parameters are corrected.
c
c it should be noted that the arrays work,iwork contain information
c required for subsequent integration. accordingly, work and iwork
c should not be altered.
c
c
integer neqn,iflag,iwork(5)
real y(neqn),t,tout,relerr,abserr,work(1)
c if compiler checks subscripts, change work(1) to work(3+6*neqn)
c
external f
c
integer k1,k2,k3,k4,k5,k6,k1m
c
c
c compute indices for the splitting of the work array
c
k1m=neqn+1
k1=k1m+1
k2=k1+neqn
k3=k2+neqn
k4=k3+neqn
k5=k4+neqn
k6=k5+neqn
c
c this interfacing routine merely relieves the user of a long
c calling list via the splitting apart of two working storage
c arrays. if this is not compatible with the users compiler,
c he must use rkfs directly.
c
call rkfs(f,neqn,y,t,tout,relerr,abserr,iflag,work(1),work(k1m),
1 work(k1),work(k2),work(k3),work(k4),work(k5),work(k6),
2 work(k6+1),iwork(1),iwork(2),iwork(3),iwork(4),iwork(5))
c
return
end