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