C ALGORITHM 771, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 23,NO. 3, September, 1997, P. 402--415. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc # Drivers # Src # This archive created: Sat May 23 19:49:22 1998 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'rksuite_90.doc' then echo shar: will not over-write existing file "'rksuite_90.doc'" else cat << \SHAR_EOF > 'rksuite_90.doc' ********************* Start of rksuite_90.doc ******************************** RKSUITE_90 Release 1.2 December 1995 written by R.W. Brankin(*) and I. Gladwell(**) (*) Numerical Algorithms Group Ltd. Wilkinson House Jordan Hill Road Oxford OX2 8DR U.K. email: richard@nag.co.uk na.brankin@na-net.ornl.gov International phone: + 44 1865 511245 International fax: + 44 1865 310139 (**) Department of Mathematics Southern Methodist University Dallas, Texas 75275-0156 U.S.A. email: gladwell@seas.smu.edu igladwel@sun.cis.smu.edu h5nr1001@vm.cis.smu.edu U.S. phone: (214) 768-2542 U.S. fax: (214) 768-4138 *** The authors would very much appreciate receiving notification *** *** of errors and constructive criticism of RKSUITE_90 *** *************************** RKSUITE_90 Overview ****************************** RKSUITE_90 is a module based on Runge-Kutta formulas that solves the initial value problem for ordinary differential equations. It integrates y' = f(t,y), y(t_start) = y_start. Here y is the solution (dependent variable) and t is the independent variable. The integration proceeds from t = t_start to t = t_end. The most commonly occurring case is where y and f are real vectors. At the end of this section, we discuss how to produce "automatically" a version of RKSUITE_90 for the cases where y and f are scalars, vectors or matrices and for the cases where y and f are real or complex. Algorithmically, RKSUITE_90 is taken almost entirely from the Fortran 77 package RKSUITE, designed and coded by L.F. Shampine and the authors of RKSUITE_90. RKSUITE_90 contains two recusive integration procedures, RANGE_INTEGRATE and STEP_INTEGRATE, and a number of associated auxiliary procedures. You define initial conditions and local error requirements, and set other options by a call to procedure SETUP. RANGE_INTEGRATE solves the commonly occurring problem of integrating the differential equations across a range obtaining answers at points you specify. STEP_INTEGRATE is designed for more complicated tasks and is a step oriented integrator. It is easy to change between RANGE_INTEGRATE and STEP_INTEGRATE: all the arguments accessible to you have the same meanings in both procedures. The distribution version of RKSUITE_90 sets the working precision, WP, using WP = SELECTED_REAL_KIND (10,50). When solving ordinary differential equations, it is advisable to use at least 10 decimal digits of working accuracy as selected in this setting (which also guarantees a decimal exponent range of 50 orders of magnitude). Edit the second line of the module RKSUITE_90_PREC for other choices of precision. The accuracy of the constants defining the Runge-Kutta formulas implies that you should not request more than 30 decimal digits of working accuracy. If you wish to work in standard DOUBLE PRECISION edit the precision line to WP = KIND(1.0D0) and for standard REAL (single precision) edit it to WP = KIND(1.0E0) Note that DOUBLE PRECISION on some machines (e.g. CRAYs) is about 30 decimal digits of working accuracy. To set various arguments in RKSUITE_90, you may need to know MACHEPS, the smallest number such that 1.0 + MACHEPS > 1.0 in the working precision, and DWARF, the smallest number available in the working precion. You may estimate these numbers by specifying a variable, DUMMY say, of kind WP and then making the calls MACHEPS = EPSILON (DUMMY) DWARF = TINY (DUMMY) Some errors that can arise in the use of RKSUITE_90 are catastrophic. Examples include input arguments that are meaningless in context or are inconsistent, and calls to procedures that have not been initialized or cannot be used in the context prevailing. A catastrophic error will automatically lead to an explanatory message on the standard output channel OUTPUT_CHANNEL and to the program STOPping. (It is possible to override this STOP so that the main program can go on to another task. This is achieved using procedures SET_STOP_ON_FATAL and GET_SAVED_STATE.) OUTPUT_CHANNEL is set (once) to 6 internally. To change this value you must edit the module RKSUITE_90. Module RKSUITE_90 needs to pass global information between its various procedures. This is achieved using a variable of a predetermined derived type, below called COMM of generic type "RK_COMM". The specification is (usually): COMM - type "RK_COMM", INTENT(INOUT) Used to communicate information between calls to the procedures in RKSUITE_90. The contents of COMM are private. ** In the base version of RKSUITE_90, "RK_COMM" is set to RK_COMM_REAL_1D. You may specify as many derived types of type "RK_COMM" as you need, hence permitting more than one simultaneous integration and/or recursive calls to the integrators RANGE_INTEGRATE and STEP_INTEGRATE in a single program. The use of this method of keeping global information also makes RKSUITE_90 "multiprocessor safe". If you wish to integrate a number of differential equations in succession in a single program, between calls you should deallocate the memory and nullify the pointers used by COMM. Procedure GARBAGE_COLLECT is provided for this purpose. In the specification below, a number of generic terms are used: 1. type("RK_COMM") 2. type(independent variable) 3. type(dependent variable) 4. shape(Y_START) In the base version these are to be interpreted as: 1. type(RK_COMM_REAL_1D) 2. REAL of KIND = WP 3. REAL of KIND = WP 4. One-dimensional array of length of the vector of initial conditions, Y_START = y_start With RKSUITE_90 we provide a UNIX script which will transform RKSUITE_90 to produce one or all of the combinations: A. type dependent variable = REAL or COMPLEX of KIND = WP B. shape(Y_START) can be zero (scalar), one (vector) or two (matrix) dimensional In these instances "RK_COMM" is REAL COMPLEX Zero dimensions RK_COMM_REAL_0D RK_COMM_COMPLEX_0D One dimension RK_COMM_REAL_1D RK_COMM_COMPLEX_1D Two dimensions RK_COMM_REAL_2D RK_COMM_COMPLEX_2D The other generic types are unchanged from the base version. There follows a template to which you should refer when reading the specification below. Note that it is for the base version of RKSUITE_90. module define_f use rksuite_90_prec, only:wp ! get precision of rksuite_90 ... ! define data for f here contains function f(t,y) real(kind=wp), intent(in) :: t real(kind=wp), dimension(:), intent(in) :: y real(kind=wp), dimension(size(y)) :: f f = (/ ... /) ! evaluate f end function f end module define_f program demo use rksuite_90 ! access rksuite_90 use define_f ! access f type(rk_comm_real_1d) :: comm ! the communication structure real(kind=wp) :: t_start=..., t_end=..., y_start(...)=(/.../), & tolerance=..., thresholds(size(y_start))=(/.../), & t_want, t_inc=..., y_maxvals(size(y_start)), & t_got, y_got(size(y_start)) integer :: flag ... ! f90 statements call setup(comm,t_start,y_start,t_end,tolerance,thresholds) ... ! f90 statements do t_want = t_want + t_inc if (t_want > t_end) exit call range_integrate(comm,f,t_want,t_got,y_got=y_got,flag=flag) if (flag /= 1) exit print*,' t = ',t_got,' y = ',y_got end do ... ! f90 statements call statistics(comm,y_maxvals=y_maxvals)! optional call to statistics print*,' y_maxvals = ',y_maxvals ... ! f90 statements end program demo ****************************************************************************** ---------------------Description of Subroutine SETUP------------------------ SUBROUTINE SETUP(COMM,T_START,Y_START,T_END,TOLERANCE,THRESHOLDS, METHOD,TASK,ERROR_ASSESS,H_START,MESSAGE) SETUP initializes the computation, so it is normally called only once. A call to SETUP must precede the first call to RANGE_INTEGRATE or STEP_INTEGRATE. Any subsequent call to SETUP reinitializes the computation. ARGUMENTS COMM - type("RK_COMM"), INTENT(OUT) Used to communicate information between calls to the procedures in RKSUITE_90. The contents of COMM are private. ** T_START - type(independent variable), INTENT(IN) The initial value of the independent variable. ** Y_START -type(dependent variable), INTENT(IN) The initial values of the solution Y. ** T_END - type(independent variable), INTENT(OUT) The integration proceeds from T_START in the direction of T_END. You can, and often will, terminate the integration before reaching T_END, but you cannot integrate past T_END. Constraint: T_END must be clearly distinguishable from T_START in the precision available. ** The integration proceeds by steps from T_START towards T_END. An approximate solution Y is computed at each step. The error in Y made in the step, i.e. the local error, is estimated. The step size is chosen automatically so that the integration will proceed efficiently while keeping this local error estimate smaller than a tolerance that you specify by means of the scalar argument TOLERANCE and THRESHOLDS which has the same shape as Y_START . You should consider carefully how you want the local error to be controlled. RKSUITE_90 basically uses relative local error control, with TOLERANCE being the desired relative accuracy. For reliable computation, integrators must work with approximate solutions that have some correct digits, so you are not allowed to specify TOLERANCE greater than 0.01. It is impossible to compute a numerical solution more accurate than the correctly rounded value of the true solution, so you are not allowed to specify a TOLERANCE that is too small for the precision you are using. Specifically, TOLERANCE must be greater than 10.0*MACHEPS. The magnitude of the local error in a component of Y on any step will not be greater than TOLERANCE * max("magnitude_y",THRESHOLDS) where "magnitude_y" is a average magnitude of the corresponding component of Y over the step. If THRESHOLDS is smaller than the current "magnitude_y", this is a relative error test and TOLERANCE indicates how many significant digits you want in Y. If THRESHOLDS is larger than the current "magnitude_y", this is an absolute error test with tolerance TOLERANCE*THRESHOLDS. Relative error control is recommended, but pure relative error control is not permitted. Specifically, no component of THRESHOLDS can be smaller than SQRT(DWARF). If certain solution components are of no interest when they are smaller in magnitude than a given threshold, you can inform RKSUITE_90 by setting corresponding components of THRESHOLDS to this threshold. Hence, you avoid the cost of computing unnecessary significant digits in a component of Y. This is important when a component of Y vanishes, and in particular, when a component of the initial value Y_START vanishes. Appropriate THRESHOLDS depend on the general size of Y in the course of the integration. Physical reasoning may help you select suitable threshold values. If you do not know what to expect of Y, you can find out by a preliminary integration using RANGE_INTEGRATE with nominal values of THRESHOLDS. As RANGE_INTEGRATE steps from T_START towards T_END, it forms Y_MAXVALS, containing the largest magnitude of each component of Y computed at any step in the integration so far. You can access Y_MAXVALS using procedure STATISTICS. Then, you can determine more appropriate values for THRESHOLDS for an accurate integration. You might, for example, take THRESHOLDS = 10.0*MACHEPS * Y_MAXVALS after completing the trial integration. TOLERANCE - REAL(KIND=WP), INTENT(IN) The relative error tolerance. Constraint: 0.01 >= TOLERANCE >= 10.0*MACHEPS ** THRESHOLDS - REAL(KIND=WP), shape(Y_START), INTENT(IN) Contains thresholds for the solution Y. Choose it so that a component of Y is not important when it is smaller in magnitude than the corresponding component of THRESHOLDS. Constraint: THRESHOLDS >= SQRT(DWARF) ** Like practically all IVP integrators, RANGE_INTEGRATE and STEP_INTEGRATE control the local error rather than the true (global) error, the difference between the numerical and true solution. Control of the local error controls the true error indirectly. Roughly speaking, RKSUITE_90 produces a solution that satisfies the differential equation with a discrepancy bounded in magnitude by the error tolerance. How close the numerical solution is to the true solution depends on the stability of the problem. Most practical problems are at least moderately stable, and the true error is then comparable to the error tolerance. To judge the accuracy of the numerical solution, you could reduce TOLERANCE substantially, e.g. by a factor of ten, and solve the problem again. This will usually result in a more accurate solution, and the true error of the first integration can be estimated by comparison. Alternatively, a global error assessment can be computed automatically by setting ERROR_ASSESS to .TRUE.. Because indirect control of the true error by controlling the local error is generally satisfactory and because both ways of assessing true errors cost twice, or more, the cost of the integration itself, such assessments should be used mainly for spot checks, for selecting appropriate tolerances for local error control, and for exploratory computations. RKSUITE_90 implements three Runge-Kutta formula pairs. You may want override the default selection depending on your requirements. METHOD - CHARACTER (LEN = *), OPTIONAL, INTENT(IN) Specifies which Runge-Kutta pair is to be used. Only the first character of METHOD is used. METHOD(1:1) = `L' or `l' - use the low order (2,3) pair = `M' or `m' - use the moderate order (4,5) pair (DEFAULT) = `H' or `h' - use the high order (7,8) pair (The order of accuracy is 3,5,8, respectively.) Constraint: METHOD(1:1) = `L', `l', `M', `m', `H' or `h'. ** The best choice for METHOD depends on the problem. As a rule, the smaller TOLERANCE, the higher order you should make METHOD. If the components of THRESHOLDS are small enough that you are effectively specifying relative error control, experience suggests TOLERANCE most efficient METHOD -4 0.01 -- 10 `L' or `l' -3 -6 10 -- 10 `M' or `m' -5 10 -- 10*MACHEPS `H' or `h' The overlap in the ranges of tolerances appropriate for a given METHOD reflects the dependence of efficiency on the problem. Making TOLERANCE smaller will normally make the integration more expensive. However, in the range of tolerances appropriate to a METHOD, the increase in cost is modest. Sometimes one METHOD, or even this kind of integrator, is a poor choice. You should not specify a very small THRESHOLDS component, like SQRT(DWARF), when the corresponding solution component might vanish. In particular, you should not do this when the component of Y_START is zero. If you do, the integrator will have to work hard with any METHOD to compute significant digits, but METHOD = `L' is particularly poor choice in this situation. All three METHODs are inefficient when the problem is "stiff". If it is only mildly stiff, you can solve it with acceptable efficiency with METHOD = `L', but if it is moderately or very stiff, an integrator designed specifically for such problems will be much more efficient. The higher the order selected by METHOD, the more smoothness is required of the solution for the METHOD to be efficient. You must decide which integrator, RANGE_INTEGRATE or STEP_INTEGRATE, to use. RANGE_INTEGRATE integrates the differential equation across a range to obtain answers at points you specify. STEP_INTEGRATE integrates across internally chosen steps and is used for all more complicated tasks. TASK - CHARACTER(LEN = *), OPTIONAL, INTENT(IN) Only the first character of TASK is significant. TASK(1:1) = `R' or `r' - RANGE_INTEGRATE is to be used (DEFAULT) = `S' or `s' - STEP_INTEGRATE is to be used Constraint: TASK(1:1) = `R' or `r' or `S' or `s' ** An assessment of the true (global) error is provided by setting ERROR_ASSESS = .TRUE.. The error assessment is updated at each step. Its value can be obtained at any time by a call to procedure GLOBAL_ERROR. Warning: RKSUITE_90 monitors the computation of the global error assessment and reports any doubts it has about the reliability of the results. The assessment scheme requires some smoothness of f(t,y), and it can be deceived if f is insufficiently smooth. At very crude tolerances the numerical solution can become so inaccurate that it is impossible to continue assessing the accuracy reliably. At very stringent tolerances the effects of finite precision arithmetic can make it impossible to assess the accuracy reliably. In either case the integrator returns with a message and a flag (FLAG = 6) reporting that global error assessment has broken down. ERROR_ASSESS - LOGICAL, OPTIONAL, INTENT(IN) = .FALSE. - do not attempt to assess the true error (DEFAULT) = .TRUE. - assess the true error, the difference between the numerical solution and the true solution. (The cost of this is roughly twice the cost of the integration itself with METHODs `M' and `H', and three times with METHOD = `L'.) ** The first step of the integration is critical because it sets the scale of the problem. The integrator will find a starting step size automatically if you set the variable H_START to 0.0. Automatic selection of the first step is so effective that you should normally use it. Nevertheless, you might want to specify a trial value for the first step to be certain that the integrator recognizes the scale on which phenomena occur near the initial point. Also, automatic computation of the first step size involves some cost, so supplying a good value for this step size will result in a less expensive start. If you are confident that you have a good value, provide it in the variable H_START. H_START - type(independent variable), OPTIONAL, INTENT(IN) = zero - The integrator will select automatically the first step size of the integration. (DEFAULT) = non-zero - The integrator will try ABS(H_START) for the first step size of the integration. ** RANGE_INTEGRATE and STEP_INTEGRATE communicate with your calling program by means of a argument FLAG. If you wish, you can also have messages written to the standard output channel OUTPUT_CHANNEL. The messages provide more detail, so it is advisable to permit them for all but production runs. MESSAGE - LOGICAL, OPTIONAL, INTENT(IN) Specifies whether you want informative messages written to OUTPUT_CHANNEL. = .TRUE. - provide messages (DEFAULT) = .FALSE. - do not provide messages ** The data input to SETUP is monitored. Any error detected is catastrophic. An error message is written to the output channel OUTPUT_CHANNEL (even if MESSAGE = .FALSE.), and the program STOPs. -------------------End of Description of Subroutine SETUP-------------------- ********************** Subroutine RANGE_INTEGRATE ************************** Subroutine RANGE_INTEGRATE is a member of the module, RKSUITE_90, for solving the initial value problem for ordinary differential equations by Runge-Kutta methods. RANGE_INTEGRATE is designed for computing an approximate solution at a sequence of points across a range of integration. First you call SETUP to specify the problem and how it is to be solved. Thereafter you (repeatedly) call RANGE_INTEGRATE to obtain answers at points you specify. Another integrator, STEP_INTEGRATE, in RKSUITE_90 is provided for more complicated tasks. RANGE_INTEGRATE requires you to specify a sequence of output points. They must be clearly distinguishable in the precision available, and the first one must be distinguishable from t_start as defined in the previous call to SETUP. You are extremely unlikely to specify points that are not clearly distinguishable except by mistake. Should this happen, the integrator will tell you how far apart the points must be. -------------------Description of Subroutine RANGE_INTEGRATE------------------ RECURSIVE SUBROUTINE RANGE_INTEGRATE(COMM,F,T_WANT,T_GOT,Y_GOT, YDERIV_GOT,FLAG) ARGUMENTS COMM - type("RK_COMM"), INTENT(IN)OUT Used to communicate information between calls to the procedures in RKSUITE_90. The contents of COMM are private. ** The differential equations are defined by the function F that you must provide. You may give this function any name you like but it must have the following interface: FUNCTION F(T,Y) T - type(independent variable), INTENT(IN) Y - type(dependent variable), shape(Y_START), INTENT(IN) F - type(dependent variable), shape(Y_START) Given input values of the independent variable T and the solution Y, evaluate the differential equations for the derivative and return the result in F. ** RANGE_INTEGRATE approximates the true solution of the initial value problem and its first derivative at the point T_WANT specified in the call to RANGE_INTEGRATE. On return, integration has been successful as far as T_GOT, and Y_GOT and YDERIV_GOT are accurate approximations to the solution and its first derivative at T_GOT. If all has gone well, T_GOT equals T_WANT and FLAG = 1. If the integrator did not reach T_WANT then T_GOT is not equal to T_WANT and an explanation is indicated by the value of FLAG. The integration proceeds by steps from T_START towards T_END (both specified in SETUP). For this reason, the specified T_WANT must be closer to T_END than the previous value of T_GOT (T_START on the first call to RANGE_INTEGRATE). T_WANT can equal T_END. T_WANT - type(independent variable), INTENT(IN) The next value of the independent variable where a solution is desired. Constraint: T_WANT must be closer to T_END than the previous value of T_GOT (T_START on the first call to RANGE_INTEGRATE). It can be equal to T_END. Unless exactly equal to T_END, it must be clearly distinguishable from T_END and from T_GOT, in the precision available. ** T_GOT - type(independent variable), INTENT(OUT) A solution has been computed at this value of the independent variable. If the task was completed successfully, it is the same as T_WANT. If the integrator did not reach T_WANT, an explanation is provided by the value of FLAG. ** Y_GOT - type(dependent variable), shape(Y_START), OPTIONAL, INTENT(OUT) Approximation to the true solution at T_GOT. At each step of the integration to T_GOT, the local error has been controlled as specified by TOLERANCE and THRESHOLDS. The local error in Y_GOT has been controlled even when T_GOT differs from T_WANT. ** YDERIV_GOT - type(dependent variable), shape(Y_START), OPTIONAL, INTENT(OUT) Approximation to the first derivative of the true solution at T_GOT. The local error has been controlled even when T_GOT differs from T_WANT. ** Normally, any difficulties encountered by RANGE_INTEGRATE will result with an error message written to the output channel OUTPUT_CHANNEL and the program STOPping unless the optional output argument FLAG is provided. If FLAG is present then RANGE_INTEGRATE reports success by setting FLAG = 1. Difficulties with the integration are reported by values FLAG > 1. If MESSAGE was set .TRUE. in the call to SETUP, details about the difficulties are written to the standard output channel OUTPUT_CHANNEL. If RANGE_INTEGRATE fails catastrophically, for example if values of its input variables are incompatible with those provided in SETUP, details are written to the output channel OUTPUT_CHANNEL (even if MESSAGE was set .FALSE.), and the program STOPs. FLAG - INTEGER, OPTIONAL, INTENT(OUT) "SUCCESS" T_GOT = T_WANT. = 1 - Successful call. To compute an approximation at a new T_WANT, just call RANGE_INTEGRATE again with the new value of T_WANT. "SOFT FAILURES" = 2 - This return is possible only when METHOD = `H'. The integrator is being used inefficiently because the step size has been reduced drastically many times to get answers at many points T_WANT. If you really need the solution at this many points, you should change to METHOD = `M' because it is (much) more efficient in this situation. To change METHOD, restart the integration from T_GOT, Y_GOT by a call to SETUP (note that you must b using the optional argument Y_GOT to be able to do this). Precede this call to SETUP by a call to GARBAGE_COLLECT to use memory efficiently. If you wish to continue on towards T_WANT with METHOD = `H', just call RANGE_INTEGRATE again. The monitor of this kind of inefficiency will be reset automatically so that the integration can proceed. = 3 - A considerable amount of work has been expended in the (primary) integration. This is measured by counting the number of calls to the function F. At least 5000 calls have been made since the last time this counter was reset. Calls to F in a secondary integration for global error assessment are not counted in this total. The integration was interrupted, so T_GOT is not equal to T_WANT. If you wish to continue on towards T_WANT, just call RANGE_INTEGRATE again. The counter measuring work will be reset to zero automatically. = 4 - It appears that this problem is "stiff". The methods implemented in RANGE_INTEGRATE can solve such problems, but they are inefficient. You should change to an integrator based on methods appropriate for stiff problems. If you want to continue on towards T_WANT, just call RANGE_INTEGRATE again. The stiffness monitor will be reset automatically. "HARD FAILURES" = 5 - It does not appear possible to achieve the accuracy specified by TOLERANCE and THRESHOLDS in the call to SETUP with the precision available on this computer and with this order of METHOD. You cannot continue integrating this problem. A higher order for METHOD, if possible, will permit greater accuracy with this precision. To increase the order of METHOD and/or continue with larger values of TOLERANCE and/or THRESHOLDS, restart the integration from T_GOT, Y_GOT by a call to SETUP (note that you must b using the optional argument Y_GOT to be able to do this). Precede this call to SETUP by a call to GARBAGE_COLLECT to use memory efficiently. = 6 - The global error assessment may not be reliable beyond the current integration point T_GOT. This may occur because either too little or too much accuracy has been requested or because f is not smooth enough near T_GOT and current values of the solution. The integration cannot be continued. This return does not mean that you cannot integrate past T_GOT, rather that you cannot do it with ERROR_ASSESS = .TRUE.. However, it may also indicate problems with the primary integration. You may continue the integration by restarting with a call to SETUP. Precede this call to SETUP by a call to GARBAGE_COLLECT to use memory efficiently. ** Performance statistics are available after any return from RANGE_INTEGRATE by a call to the procedure STATISTICS. IF ERROR_ASSESS was set .TRUE., global error assessment is available after any return from RANGE_INTEGRATE by a call to the procedure GLOBAL_ERROR. After a hard failure (FLAG = 5 or 6) the diagnostic procedures STATISTICS and GLOBAL_ERROR may be called only once. Other procedures from RKSUITE_90 may not be called at all, except SETUP to restart the integration. Precede this call to SETUP by a call to GARBAGE_COLLECT to use memory efficiently. -----------End of Description of Subroutine RANGE_INTEGRATE----------------- ------------Description of Subroutine GLOBAL_ERROR-------------------------- SUBROUTINE GLOBAL_ERROR(COMM,RMS_ERROR,MAX_ERROR,T_MAX_ERROR) To assess the true (global) error of the integration with STEP_INTEGRATE or RANGE_INTEGRATE, set ERROR_ASSESS = .TRUE. in the call to SETUP. After any call to STEP_INTEGRATE or RANGE_INTEGRATE, GLOBAL_ERROR may be called for information about the assessment. The solution Y is computed in the primary integration. (The values Y_GOT in RANGE_INTEGRATE and Y_NOW in STEP_INTEGRATE are instances of Y resulting from the primary integration.) A more accurate "true" solution Y_SECOND is computed in a secondary integration. The error is measured as specified in SETUP for local error control. At each step in the primary integration, an average magnitude "magnitude_y" of Y is computed, and the error is abs(Y - Y_SECOND)/max("magnitude_y", THRESHOLDS). It is difficult to estimate reliably the true error at a single point. For this reason the RMS (root-mean-square) average of the estimated global error is returned in RMS_ERROR. This average is taken over all steps from T_START through the current integration point. If all has gone well, the average errors reported in RMS_ERROR will be comparable to TOLERANCE. Other useful quantities are MAX_ERROR, the maximum error seen in any component in the integration so far, and T_MAX_ERROR, the point where the maximum error first occurred. You may call GLOBAL_ERROR only once after a hard failure in RANGE_INTEGRATE or STEP_INTEGRATE. ARGUMENTS COMM - type("RK_COMM"), INTENT(INOUT) Used to communicate information between calls to the procedures in RKSUITE_90. The contents of COMM are private. ** RMS_ERROR - REAL(KIND=WP), shape(Y_START), OPTIONAL, INTENT(OUT) RMS_ERROR approximates the RMS average of the true error of the numerical solution. The average is taken over all steps from T_START through the current integration point. If all has gone well, RMS_ERROR will be comparable to TOLERANCE. ** MAX_ERROR - REAL(KIND=WP), OPTIONAL, INTENT(OUT) The maximum weighted approximate true error taken over all solution components and all steps from T_START through the current integration point. If all has gone well, MAX_ERROR will be comparable to TOLERANCE. ** T_MAX_ERROR - type(independent variable), OPTIONAL, INTENT(OUT) First value of the independent variable where an approximate true error attains the maximum value MAX_ERROR. ** At least one optional argument must be present. The call to GLOBAL_ERROR is monitored. Any error is catastrophic. A message is written to OUTPUT_CHANNEL (even if MESSAGE = .FALSE.), and the program STOPs. -----------------End of Description of Subroutine GLOBAL_ERROR------------- -------------------Description of Subroutine STATISTICS------------------- SUBROUTINE STATISTICS(COMM,TOTAL_F_CALLS,STEP_COST,WASTE, NUM_SUCC_STEPS,H_NEXT,Y_MAXVALS) STATISTICS may be called after any call to RANGE_INTEGRATE or STEP_INTEGRATE to obtain details about the integration. You may call STATISTICS only once after a hard failure in RANGE_INTEGRATE or STEP_INTEGRATE. ARGUMENTS COMM - type("RK_COMM"), INTENT(INOUT) Used to communicate information between calls to the procedures in RKSUITE_90. The contents of COMM are private. ** TOTAL_F_CALLS - INTEGER, OPTIONAL, INTENT(OUT) Total number of calls to F in the integration so far -- a measure of the cost of the integration. ** STEP_COST - INTEGER, OPTIONAL, INTENT(OUT) Cost of a typical step with this METHOD measured in calls to F. ** WASTE - REAL(KIND=WP), OPTIONAL, INTENT(OUT) The number of attempted steps that failed to meet the local error requirement divided by the total number of steps attempted so far in the integration. A "large" fraction indicates that the integrator is having trouble with this problem. This can happen when the problem is "stiff" and also when the solution has discontinuities in a low order derivative. ** NUM_SUCC_STEPS - INTEGER, OPTIONAL, INTENT(OUT) The number of successful steps in the integratyion so far. ** H_NEXT - type(independent variable), OPTIONAL, INTENT(OUT) The step size the integrator plans to use for the next step. ** Y_MAXVALS - REAL(KIND=WP), shape(Y_START), OPTIONAL, INTENT(OUT) Y_MAXVALS is the componentwise largest value of ABS(Y) computed at any step in the integration so far. (With METHODs `L' and `M' in RANGE_INTEGRATE, Y_GOT is computed by interpolation, so Y_MAXVALS might be a little different than any value ABS(Y_GOT) reported so far.) ** At least one optional argument must be present. The call to STATISTICS is monitored. Any error is catastrophic. A message is written to OUTPUT_CHANNEL (even if MESSAGE = .FALSE.), and the program STOPs. --------------End of Description of Subroutine STATISTICS--------------------- ************ STEP_INTEGRATE: Subroutine for Complicated Tasks **************** STEP_INTEGRATE is a member of RKSUITE_90, for solving the initial value problem for ordinary differential equations by Runge-Kutta methods. SETUP is used to specify the problem and how it is to be solved. STEP_INTEGRATE is used to advance the integration one step. Another integrator, RANGE_INTEGRATE, in RKSUITE_90 is provided for the task of obtaining an approximate solution at a sequence of points. STEP_INTEGRATE is designed for more complicated tasks that require close monitoring of the integration and additional capabilities. To ease use, less common demands are handled by means of the auxiliary procedures INTERPOLATE and RESET_T_END. ----------------------------------------------------------------------------- ------------------Description of Subroutine STEP_INTEGRATE-------------------- RECURSIVE SUBROUTINE STEP_INTEGRATE(COMM,F,T_NOW,Y_NOW,YDERIV_NOW,FLAG) STEP_INTEGRATE is used to step from T_START in the direction of T_END as specified in SETUP. One way to use STEP_INTEGRATE involves repeatedly resetting T_END, so a procedure RESET_T_END is provided for this purpose. When STEP_INTEGRATE is called, Y_NOW and YDERIV_NOW approximate the solution its first derivative resectively at T_NOW. If the integrator encounters some difficulty in taking a step, it returns with these values unaltered and provides an explanation by means of the value of the optional argument FLAG, if supplied. If all goes well, STEP_INTEGRATE returns with FLAG = 1, and Y_NOW and YDERIV_NOW are the new values of the approximate solution and its first derivative at the end of a single step to the new T_NOW. STEP_INTEGRATE tries to advance the integration as far as possible subject to passing the test on the local error and not going past T_END. In the call to SETUP, you can specify that STEP_INTEGRATE try H_START as its first step size or that it compute automatically an appropriate value. Thereafter STEP_INTEGRATE estimates an appropriate step size for its next step. This value and other details of the integration can be obtained after any call to STEP_INTEGRATE by a call to procedure STATISTICS. The local error is controlled at every step as specified in SETUP. If you wish to assess the true error, you must set ERROR_ASSESS = .TRUE. in the call to SETUP. This assessment can be obtained after any call to STEP_INTEGRATE by a call to the procedure GLOBAL_ERROR. There are three ways to use STEP_INTEGRATE: 1. Step from T_START towards T_END, accepting answers at the points chosen by the integrator. This is often the best way to proceed when you want to see how the solution behaves throughout the interval of integration because the integrator tends to produce answers more frequently where the solution changes more rapidly (the step sizes are usually smaller there). ** If you want answers at specific points, two ways to proceed are: 2. The more efficient way is to step past the point where a solution is desired, and then call procedure INTERPOLATE to get an answer there. Within the span of the current step, you can get all the answers you want at very little cost by repeated calls to INTERPOLATE. This is very valuable when you want to find where something happens, e.g., where a particular solution component vanishes. You cannot proceed in this way with METHOD = `H'. 3. The other way to get an answer at a specific point is to set T_END to this value and integrate to T_END. STEP_INTEGRATE will not step past T_END. So, when a step would carry it past T_END, it will be reduced so as to produce an answer at T_END exactly. After getting an answer there (with T_NOW = T_END), you can reset T_END to the next point where you want an answer, and repeat. T_END could be reset by a call to SETUP, but you should not do this. You should use the procedure RESET_T_END because it is both easier to use and much more efficient. This way of getting answers at specific points can be used with any of the METHODs, but it is the only way with METHOD = `H'. It can be inefficient. Should this be the case, RKSUITE_90 will bring the matter to your attention. ARGUMENTS COMM - type("RK_COMM"), INTENT(INOUT) Used to communicate information between calls to the procedures in RKSUITE_90. The contents of COMM are private. ** The differential equations are defined by the function F that you must provide. You may give this function any name you like but it must have the following interface: FUNCTION F(T,Y) T - type(independent variable), INTENT(IN) Y - type(dependent variable), shape(Y_START), INTENT(IN) F - type(dependent variable), shape(Y_START) Given input values of the independent variable T and the solution Y, evaluate the differential equations for the derivative and return the result in F. ** On each call, STEP_INTEGRATE tries to take a step from the current value of T_NOW (initially T_START) in the direction of T_END. On the previous call, the integrator estimated an appropriate step size. (On the first call, either you provide the step size, or one is determined automatically.) If this step size is too big for the formula to achieve the specified accuracy, the integrator will adjust the step size and try again. It keeps trying until it produces a solution that is sufficiently accurate, or it decides to report that it has run into trouble via FLAG (and the standard output channel OUTPUT_CHANNEL if MESSAGE was set .TRUE. in SETUP). In any case, the values returned in Y_NOW and YDERIV_NOW satisfy the specified local accuracy requirement at the value T_NOW. T_NOW - type(independent variable), INTENT(OUT) Current value of the independent variable after a step towards T_END. ** Y_NOW - type(dependent variable), shape(Y_START), OPTIONAL, INTENT(OUT) Approximation to the solution at T_NOW. The local error of the step to T_NOW was no greater than permitted by the tolerances TOLERANCE and THRESHOLDS as specified in SETUP. ** YDERIV_NOW - type(dependent variable), shape(Y_START), OPTIONAL, INTENT(OUT) Approximation to the derivative of the solution at T_NOW. ** Normally, any difficulties encountered by STEP_INTEGRATE will result with an error message written to the output channel OUTPUT_CHANNEL and the program STOPping unless the optional output argument FLAG is provided. If FLAG is present then STEP_INTEGRATE reports success by setting FLAG = 1. Difficulties with the integration are reported by values of FLAG > 1. In such cases the integration has not advanced and T_NOW, Y_NOW, and YDERIV_NOW are unchanged. If MESSAGE was set .TRUE. in the call to SETUP, some details about the difficulty are written to the standard output channel OUTPUT_CHANNEL. The call to STEP_INTEGRATE is monitored. If a catastrophic error is detected, for example when STEP_INTEGRATE has been called out of context, then an error message is written to the output channel OUTPUT_CHANNEL (even if MESSAGE was set to .FALSE.), and the program STOPs. FLAG - INTEGER, INTENT(OUT) "SUCCESS" = 1 - Successful call. A step was taken to T_NOW. To continue the integration in the direction of T_END, just call STEP_INTEGRATE again. Do not alter any variables. "SOFT FAILURE" = 2 - The integrator is being used inefficiently because the step size has been reduced drastically many times to get answers at many values of T_END. If you really need the solution at this many specific points, you should use INTERPOLATE to get the answers inexpensively. If you need to change METHOD for this purpose, restart the integration from T_NOW, Y_NOW by a call to SETUP (note that you must b using the optional argument Y_NOW to be able to do this). Precede this call to SETUP by a call to GARBAGE_COLLECT to use memory efficiently. If you wish to continue toward T_END, just call STEP_INTEGRATE again. The monitor of this kind of inefficiency will be reset automatically so that the integration can proceed. = 3 - A considerable amount of work has been expended in the (primary) integration. This is measured by counting the number of calls to the procedure F. At least 5000 calls have been made since the last time this counter was reset. Calls to F in a secondary integration for global error assessment are not counted in this total. If you wish to continue towards T_END, just call STEP_INTEGRATE again. The counter measuring work will be reset to zero automatically. = 4 - It appears that this problem is "stiff". The methods implemented in STEP_INTEGRATE can solve such problems, but they are inefficient. You should change to an integrator based on methods appropriate for stiff problems. If you want to continue toward T_END, just call STEP_INTEGRATE again. The stiffness monitor will be reset automatically. "HARD FAILURE" = 5 - It does not appear possible to achieve the accuracy specified by TOLERANCE and THRESHOLDS in the call to SETUP with the precision available on this computer and with this order of METHOD. You cannot continue integrating this problem. A higher order for METHOD, if possible, will permit greater accuracy with this precision. To increase the order of METHOD and/or continue with larger values of TOLERANCE and/or THRESHOLDS, restart the integration from T_NOW, Y_NOW by a call to SETUP (note that you must b using the optional argument Y_NOW to be able to do this). Precede this call to SETUP by a call to GARBAGE_COLLECT to use memory efficiently. = 6 - The global error assessment may not be reliable beyond the current integration point T_NOW. This may occur because either too little or too much accuracy has been requested or because f is not smooth enough for values of t near T_NOW and current values of the solution. The integration cannot be continued. This return does not mean that you cannot integrate past T_NOW, rather that you cannot do it with ERROR_ASSESS = .TRUE.. However, it may also indicate problems with the primary integration. You may continue the integration by restarting with a call to SETUP. Precede this call to SETUP by a call to GARBAGE_COLLECT to use memory efficiently. ** Statistics of the performance of STEP_INTEGRATE are available after any return from STEP_INTEGRATE by a call to the procedure STATISTICS. If ERROR_ASSESS was set .TRUE., global error assessment is available after any return from STEP_INTEGRATE by a call to the procedure GLOBAL_ERROR. After a hard failure (FLAG = 5 or 6) the diagnostic procedures STATISTICS and GLOBAL_ERROR may be called only once. Other procedures from RKSUITE_90 may not be called at all, except SETUP to restart the integration. Precede this call to SETUP by a call to GARBAGE_COLLECT to use memory efficiently. ---------------Description of Subroutine RESET_T_END-------------------------- SUBROUTINE RESET_T_END(COMM,T_END_NEW) In the description of STEP_INTEGRATE it is explained how to get answers at specific values of the independent variable by resetting T_END. SETUP could be used to reset T_END, but there are good reasons for calling RESET_T_END for this specific task: * RESET_T_END is simpler to use. * RESET_T_END is much more efficient than SETUP because it only resets the value of a variable whereas SETUP completely restarts the integration. The integration proceeds from T_START in the direction of T_END, and at present is at T_NOW. To change T_END to a new value T_END_NEW, just call RESET_T_END with T_END_NEW as the argument. You must continue integrating in the same direction, so the sign of (T_END_NEW - T_NOW) must be the same as the sign of (T_END - T_START). To change direction of integration you must restart by a call to SETUP. Precede this call to SETUP by a call to GARBAGE_COLLECT to use memory efficiently. RESET_T_END cannot be called after a call to RANGE_INTEGRATE. ARGUMENTS COMM - type("RK_COMM"), INTENT(INOUT) Used to communicate information between calls to the procedures in RKSUITE_90. The contents of COMM are private. ** T_END_NEW - type(independent variable), INTENT(IN) The new value of T_END. Constraints: sign(T_END_NEW - T_NOW) = sign(T_END - T_START). T_END must be clearly distinguishable from T_NOW in the precision available. ** The call to RESET_T_END is monitored. Any error is catastrophic. An error message is written to the output channel OUTPUT_CHANNEL (even if MESSAGE = .FALSE.), and the program STOPs. --------------End of Description of Subroutine RESET_T_END-------------------- -------------------Description of Subroutine INTERPOLATE--------------------- RECURSIVE SUBROUTINE INTERPOLATE(COMM,F,T_WANT,Y_WANT,YDERIV_WANT) In the description of STEP_INTEGRATE it is explained that when integrating with METHOD = `L' or `M', answers can be obtained inexpensively by interpolation. Procedure INTERPOLATE is provided for this purpose. Within the span of each step the solution is approximated by a polynomial of degree 3 when METHOD = `L' and a polynomial of degree 6 when METHOD = `M'. The polynomials can be evaluated anywhere, but the theory assures accurate approximations for the solution and its first derivative only within the span of the current step, or very close to it. The interpolants for successive steps connect to form a piecewise polynomial approximation over the whole interval of integration that is continuous and has a continuous derivative (in the precision available). With METHOD = `L', the interpolant uses just solution and derivative information returned after calls to STEP_INTEGRATE. The matter is slightly more complicated and expensive with METHOD = `M'. In the latter case additional calls to function F are needed to initialize the computation. Although more expensive than for METHOD = `L', this extra cost is incurred only on those steps where you require interpolation, and then only once per step, no matter how many answers you require in the span of the step. In either case it is far more efficient to let the integrator work with the largest step size that will yield the desired accuracy and obtain answers by interpolation than to obtain answers by reducing the step size so as to step to the points where the answers are desired. INTERPOLATE is called after a successful step by STEP_INTEGRATE from a previous value of T_NOW, called T_OLD below, to the current value of T_NOW to get an answer at T_WANT. You can specify any value of T_WANT you wish, but specifying a value outside the interval [T_OLD,T_NOW], called "extrapolation", is likely to yield answers with unsatisfactory accuracy. Warning: You cannot call INTERPOLATE after a return from STEP_INTEGRATE with any kind of failure. You cannot call INTERPOLATE when you are using RANGE_INTEGRATE. ARGUMENTS COMM - type("RK_COMM"), INTENT(INOUT) Used to communicate information between calls to the procedures in RKSUITE_90. The contents of COMM are private. ** Because INTERPOLATE may need to evaluate the differential equations, it must be supplied with the name of the function you provided to STEP_INTEGRATE for evaluating the differential equations. FUNCTION F(T,Y) T - type(independent variable), INTENT(IN) Y - type(dependent variable), shape(Y_START), INTENT(IN) F - type(dependent variable), shape(Y_START) Given input values of the independent variable T and the solution Y, evaluate the differential equations for the derivative and return the result in F. ** T_WANT - type(independent variable), INTENT(IN) The value of the independent variable where a solution is desired. ** The interpolant is to be evaluated at T_WANT to approximate the solution and/or its first derivative there. Y_WANT - type(dependent variable), shape(Y_START), OPTIONAL, INTENT(OUT) Approximation the true solution at T_WANT. ** YDERIV_WANT - type(dependent variable), shape(Y_START), OPTIONAL, INTENT(OUT) Approximation to the first derivative of the true solution at T_WANT. ** If neither of Y_WANT and YDERIV_WANT is present, then the call to INTERPOLATE results in the coefficients of the interpolant being initialized. The call to INTERPOLATE and the data input are monitored. Any error is catastrophic. An error message is written to the output channel OUTPUT_CHANNEL (even if MESSAGE = .FALSE.), and the program STOPs. --------------End of Description of Subroutine INTERPOLATE------------------ ---------------Description of Subroutine GARBAGE_COLLECT-------------------- SUBROUTINE GARBAGE_COLLECT(COMM) If you wish to integrate a number of differential equations in succession in a single program, between calls you should deallocate the memory and nullify the pointers used by COMM. Procedure GARBAGE_COLLECT is provided for this purpose. COMM - type("RK_COMM"), INTENT(INOUT) Used to communicate information between calls to the procedures in RKSUITE_90. The contents of COMM are private. ** ----------------End of RKSUITE_90 Subroutine Documentation------------------- -----------------------------ACKNOWLEDGEMENTS------------------------------- P.J. Prince and J.R. Dormand of the University of Cleveland (UK) developed one of the pairs of Runge-Kutta formulas used in RKSUITE_90 and P. Bogacki of Old Dominion University and L.F. Shampine developed the others. We are grateful for the assistance these friends and colleagues gave us in implementing their formulas. We also wish to thank G. Kraut of the University of Texas at Tyler for her meticulous testing of the methods in the Fortran 77 version of RKSUITE. NATO Scientific Affairs Division (Grant 898/83) funded early joint work of I. Gladwell and L.F. Shampine that led to the development of RKSUITE. R.W. Brankin's involvement was entirely funded by the Numerical Algorithms Group Ltd. Some of L.F. Shampine's research on basic algorithms used later in RKSUITE was supported in part by the Applied Mathematical Sciences program of the Office of Energy Research of the U.S. Department of Energy. Development of RKSUITE_90 was partly funded by NATO grant CRG 920037. --------------------------------REFERENCES----------------------------------- Some references describe the formulas and algorithms used in RKSUITE. Others describe the design, implementation, and testing of codes based on explicit Runge-Kutta formulas and the development of additional capabilities that played a more-or-less direct role in the development of RKSUITE. Those publications directly pertaining to RKSUITE_90 are marked (**). Those that played a significant part in the design choices are marked (*). C.A. Addison, W.H. Enright, P.W. Gaffney, I. Gladwell and P.M. Hanson, "A Decision Tree for the Numerical Solution of Initial Value Ordinary Differential Equations", ACM Trans. on Math. Soft., 17 (1991), 1-10. (*) P. Bogacki and L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas", Appl. Math. Lett., 2 (1989) 321-325. Note: Derives the (2,3) pair implemented in the Fortran 77 RKSUITE and in RKSUITE_90 (*) P. Bogacki and L.F. Shampine, "An Efficient Runge-Kutta (4,5) Pair", Rept. 89-20, Math. Dept., Southern Methodist University, Dallas, Texas, USA, 1989. Note: Derives the (4,5) pair implemented in the Fortran 77 RKSUITE and in RKSUITE_90 R.W. Brankin, J.R. Dormand, I. Gladwell, P.J. Prince and W.L. Seward, "A Runge-Kutta-Nystrom Code", ACM Trans. on Math. Soft., 15 (1989) 31-40. (**) R.W. Brankin and I. Gladwell, "A Fortran 90 ODE Solver", pp. 363-376 of "Annals of Numer. Math., 1 (Proceedings of SCADE '93)" (ed. P. Sharp), J.C.Baltzer ag Science Publishers, 1994. Note: Describes Release 0 of RKSUITE_90 (**) R.W. Brankin, I. Gladwell and L.F. Shampine, "RKSUITE: a Suite of Runge-Kutta Codes for the Initial Value Problem for ODEs", Softreport 92-S1, Math. Dept., Southern Methodist University, Dallas, Texas, U.S.A, 1991, (also available by anonymous ftp and from netlib). Note: Contains the software and documentation for the Fortran 77 RKSUITE (**) R.W. Brankin, I. Gladwell and L.F. Shampine, "RKSUITE: A Suite of Explicit Runge-Kutta Codes", pp. 41-53 of "Contributions to Numerical Mathematics" (ed. R.P. Agarwal), WSSIAA 2, World Scientific Press, 1993. Note: Discusses the design philosophy and implementation of the Fortran 77 RKSUITE R.W. Brankin and I. Gladwell, "Using Shape Preserving Local Interpolants for Plotting Solutions of Ordinary Differential Equations", IMA J. of Numer. Anal., 9 (1989) 555-566. R.W. Brankin, I. Gladwell and L.F. Shampine, "Starting Adams and BDF Codes at Optimal Order", J. Comp. Appl. Math., 21 (1988) 357-368. I. Gladwell, "Initial Value Routines in the NAG Library", ACM Trans. on Math. Soft., 5 (1979) 386-400. I. Gladwell, J.A.I. Craigie and C.R. Crowther, "Testing Initial Value Routines as Black Boxes", Numer. Anal. Rept. 34, Math. Dept., Univ. of Manchester, U.K., 1979. I. Gladwell, M. Berzins and M. Brankin, "Design of Stiff Integrators in the NAG Library", SIGNUM Newsletter, 23 (1988) 16-23. (*) I. Gladwell, L.F. Shampine, L.S. Baca and R.W. Brankin, "Practical Aspects of Interpolation with Runge-Kutta Codes", SIAM J. Sci., Stat. Comp., 8 (1987) 322-341. Note: Discusses types of interpolation strategies employed in the Fortran 77 RKSUITE and in RKSUITE_90 (*) I. Gladwell, L.F. Shampine and R.W. Brankin, "Automatic Selection of the Initial Step Size for an ODE Solver", J. Comp. Appl. Math., 18 (1987) 175-192. Note: Presents early version of starting step size strategy implemented in the Fortran 77 RKSUITE and in RKSUITE I. Gladwell, L.F. Shampine and R.W. Brankin, "Locating Special Events When Solving ODEs", Appl. Math. Lett., 1 (1988) 153-156. (**) G. Kraut, "A Comparison of RKSUITE with Runge-Kutta Codes from the IMSL, NAG and SLATEC Libraries", Report 91-6, Math. Dept., Southern Methodist University, Dallas, Texas, U.S.A., 1991. Note: Presents efficiency tests of formulas implemented in the Fortran 77 RKSUITE and in RKSUITE_90 in comparison with formulas in the most widely used current ODE codes (*) P.J. Prince and J.R. Dormand, " High Order Embedded Runge-Kutta Formulae", J. Comput. Appl. Math., 7 (1981), 67-85. Note: Derives the (7,8) pair implemented in the Fortran 77 RKSUITE and in RKSUITE_90 L.F. Shampine, "Local Extrapolation in the Solution of Ordinary Differential Equations", Math. Comp., 27 (1973) 91-97. L.F. Shampine, "Limiting Precision in Differential Equation Solvers", Math. Comp., 28 (1974) 141-144. L.F. Shampine, "Storage Reduction for Runge-Kutta Codes", ACM Trans. on Math. Soft., 5 (1979) 245-250. L.F. Shampine, "Stiffness and Non-Stiff Differential Equation Solvers II: Detecting Stiffness with Runge-Kutta Methods", ACM Trans. on Math. Soft., 3 (1977) 44-53. L.F. Shampine, "Local Error Control in Codes for Ordinary Differential Equations", Appl. Math. Comp., 3 (1977) 189-210. L.F. Shampine, "Limiting Precision in Differential Equation Solvers II: Sources of Trouble and Starting a Code", Math. Comp., 32 (1978) 1115-1122. L.F. Shampine, "Lipschitz Constants and Robust ODE Codes", pp. 427-449 in J.T. Oden, ed., Computational Methods in Nonlinear Mechanics, North-Holland, Amsterdam, 1980. L.F. Shampine, "Estimating the Cost of Output in ODE Codes", Matematica Aplicada e Computacional, 2 (1983) 157-169. L.F. Shampine, "Stiffness and the Automatic Selection of ODE Code", J. Comp. Phys., 54 (1984) 74-86. L.F. Shampine, "Stability of Explicit Runge-Kutta Methods", Comp. & Maths. with Applics., 10 (1984) 419-432. L.F. Shampine, "Measuring Stiffness", Appl. Numer. Math., 1 (1985) 107-119. L.F. Shampine, "The Step Sizes Used by One-Step Codes for ODEs", Appl. Numer. Math., 1 (1985) 95-106. L.F. Shampine, "Interpolation for Runge-Kutta Methods", SIAM J. Numer. Anal., 22 (1985) 1014-1027. L.F. Shampine, "Global Error Estimation with One-Step Methods", Comp. & Maths. with Applics., 12A (1986) 885-894. L.F. Shampine, "Some Practical Runge-Kutta Formulas", Math. Comp., 46 (1986) 135-150. L.F. Shampine, "Interpolation for Variable Order Runge-Kutta Methods", Comp. & Maths. with Applics., 14 (1987) 255-260. L.F. Shampine, "Tolerance Proportionality in ODE Codes", pp. 118-136 in Numerical Methods in ODEs, A. Bellen et al., eds., Lecture Notes in Math., 1386, Springer, Berlin, 1989. (**) L.F. Shampine, "Diagnosing Stiffness for Runge-Kutta Methods", SIAM J. Sci., Stat. Comput., 12 (1991) 260-272. Note: Presents stiffness detection strategy used in the Fortran 77 RKSUITE and algorithmic base for the stiffness test in RKSUITE_90. L.F. Shampine and L.S. Baca, "Fixed vs. Variable Order Runge-Kutta", ACM Trans. on Math. Soft., 12 (1986) 1-23. L.F. Shampine, S.M. Davenport and H.A. Watts, " Comparison of Some Codes for the Initial Value Problem for Ordinary Differential Equations", pp. 349-353 in Numerical Solutions of Boundary Value Problems for Ordinary Differential Equations, A. K. Aziz, ed., Academic, New York,1975. L.F. Shampine, S.M. Davenport and H.A. Watts, "Solving Nonstiff Ordinary Differential Equations -- the State of the Art", SIAM Review, 18 (1976) 376-411. (*) L.F. Shampine and I. Gladwell, "The Next Generation of Runge-Kutta Codes", in Computational Ordinary Differential Equations, J.R.Cash and I.Gladwell, eds., IMA Conference Series, New Series 39, Clarendon Press, Oxford, 1992. Note: Outlines the requirements and design of the Fortran 77 RKSUITE L.F. Shampine, I. Gladwell and R.W. Brankin, "Reliable Solution of Special Root Finding Problems for ODE's", ACM Trans. on Math. Soft., 17 (1991) 11-25. L.F. Shampine and M.K. Gordon, "Numerical Solution of Ordinary Differential Equations: the Initial Value Problem", W. H. Freeman and Co., San Francisco, 1975. L.F. Shampine, M.K. Gordon and J.A. Wisniewski, "Variable order Runge-Kutta codes", pp. 83-101 in Computational Techniques for Ordinary Differential Equations, I. Gladwell and D.K. Sayers, eds., Academic, London, 1980. L.F. Shampine and H.A. Watts, "Comparing Error Estimators for Runge-Kutta methods", Math. Comp., 25 (1971) 445-455. L.F. Shampine and H.A. Watts, "Global Error Estimation for Ordinary Differential Equations", ACM Trans. on Math. Soft., 2 (1976) 172-186. L.F. Shampine and H.A. Watts, "The Art of Writing a Runge-Kutta Code, Part I", in Mathematical Software III, J. R. Rice, ed., Academic, New York, 1977. L.F. Shampine and H.A. Watts, "The Art of Writing a Runge-Kutta Code, II", Appl. Math. Comp., 5 (1979) 93-121. (*) H.A. Watts, "Step Size Control in Ordinary Differential Equation Solvers", Trans. Soc. for Computer Simulation, 1 (1984), 15-25. Note: Gives earlier version of the step size stragety implemented in Fortran 77 RKSUITE and in RKSUITE_90 ************************ End of rksuite_90.doc ******************************* SHAR_EOF fi # end of overwriting check if test -f 'README' then echo shar: will not over-write existing file "'README'" else cat << \SHAR_EOF > 'README' RKSUITE_90 Release 1.2 December 1995 The subdirectories are Src: containing 2 directories Base: containing 5 files rksuite_90.bas the base verion of the source make_rk a tool to create the desired version of rksuite_90 do_r_to_c a tool used by make_rk do_1_to_2 a tool used by make_rk do_1_to_0 a tool used by make_rk Derived: containing 7 files rksuite_90_r0.f90 Version with dependent variable real scalar rksuite_90_r1.f90 Version with dependent variable real 1d array rksuite_90_r2.f90 Version with dependent variable real 2d array rksuite_90_c0.f90 Version with dependent variable complex scalar rksuite_90_c1.f90 Version with dependent variable complex 1d array rksuite_90_c2.f90 Version with dependent variable complex 2d array rksuite_90_all.f90 Version treating all of the above Doc: containing 2 files README this file rksuite_90.doc documentation for rksuite_90 Drivers: containing 18 files eg_simple.f90 simple example for real vector dependent variable eg_simple_c0.f90 simple example for complex scalar :: :: eg_simple_c0_a.f90 simple example for complex scalar :: :: eg_stiff_r0.f90 simple stiff example for real scalar :: :: eg_two_body.f90 simple example for real 2-d array :: :: eg_invar_imbed.f90 invariant imbedding demonstrating recursion eg_quadrature.f90 multiple quadrature demonstrating recursion eg_orr_somm.f90 example for complex vector dependent variable eg_orr_somm_r.f90 real dependent variable equivalent of eg_orr_somm.f90 eg_simple.r eg_simple_c0.r eg_simple_c0_a.r eg_stiff_r0.r eg_two_body.r eg_invar_imbed.r eg_quadrature.r eg_orr_somm.r eg_orr_somm_r.r which are results files corresponding to the above examples. These have been generated by "NAGWare f90 compiler Version 2.1(557)" on an SGI Indigo under IRIX 4.0.5. Instructions for unix boxes... Change to directory Code Ensure that files make_rk, do_1_to_2, do_1_to_0 and do_r_to_c are executable Type make_rk opt to produce the desired version of rksuite_90.f90 opt = r0 treats case - dependent variable is a real scalar r1 treats case - .. .. .. .. 1d array r2 treats case - .. .. .. .. 2d array c0 treats case - .. .. .. complex scalar c1 treats case - .. .. .. .. 1d array c2 treats case - .. .. .. .. 2d array all - all the above are treated Use example(s) above to check the code is "ok" Instructions for non-unix boxes... Either: Use the codes available in the Src/Derived directory Or: Rename rksuite_90.bas to rksuite_90.f90 and you have a version which treats "dependent variable is a real 1d array". Try to emulate the script files make_rk, do_1_to_2, do_1_to_0 and do_r_to_c to create a version of rksuite_90 which treats required type for dependent variables. Use example(s) above to check the code is "ok" ------------end SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'eg_two_body.r' then echo shar: will not over-write existing file "'eg_two_body.r'" else cat << \SHAR_EOF > 'eg_two_body.r' ymax 1.90E+00 2.29E+00 4.36E-01 4.36E+00 The integration reached 2.00E+01 The cost of the integration in calls to F was 4578 The fraction of failed steps was 0.02 At t = 20, the true error is 2.25E-09 assessed error is 1.04E-07 maximum error was 1.16E-06 at 1.89E+01 SHAR_EOF fi # end of overwriting check if test -f 'eg_two_body.f90' then echo shar: will not over-write existing file "'eg_two_body.f90'" else cat << \SHAR_EOF > 'eg_two_body.f90' module two_body_f ! PROBLEM ! ! Integrate a two body problem. The equations for the coordinates ! (x(t),y(t)) of one body as functions of time t in a suitable frame of ! reference are ! ! x'' = - x/r**3, ! y'' = - y/r**3, where r = SQRT(x**2 + y**2). ! ! The initial conditions lead to elliptic motion with eccentricity ECC. ! This parameter will be taken to be 0.9. ! ! x(0) = 1-ECC, x'(0) = 0, ! y(0) = 0, y'(0) = SQRT((1+ECC)/(1-ECC)). ! ! An accurate solution that shows the general behavior of the orbit is ! desired. The coordinates will be returned at every time step in [0,20]. ! This is a standard test problem for which there is a semi-analytical ! solution. It will be compared to the computed solution at the end of ! the interval. ! ! SOLUTION ! ! For illustrative purposes we shall pose the system using a two dimensional ! array of unkowns ! ! ( x x' ) ! ( y y' ) ! ! to obtain ! ! ( x x' )' = ( x' -x/r**3 ) ! ( y y' ) ( y' -y/r**3 ) ! ! Since it is the general behavior of the solution that is desired, it is ! best to let the integrator choose where to provide answers. It will ! produce answers more frequently where the solution changes more rapidly. ! Because the solution is inspected at every step, the task is to use ! STEP_INTEGRATE. ! ! On physical grounds the solution is expected to be somewhat unstable ! when one body approaches the other. To obtain an accurate solution, a ! stringent relative error tolerance should be imposed -- TOL = 1.0D-10 ! will be used. At a tolerance this stringent the highest order pair, ! METHOD = 'H', is likely to be the most efficient choice. This method is ! inefficient when it is to produce answers at a great many specific ! points. It is most effective when used as in this template. The ! solution components are expected to be of order 1, so threshold values ! of THRESHOLDS = 1.0D-13 are reasonable. When a solution component is smaller ! in magnitude than this threshold, the code will control the local error ! to be no more than TOL*THRESHOLDS = 1.0D-23. The reasonableness of this ! choice will be monitored by printing out the maximum value seen for each ! solution component in the course of the integration. Error and warning ! messages, if any, will be printed out. ! ! This is the standard test problem D5 of T.E. Hull, W.H. Enright, ! B.M. Fellen, and A.E. Sedgwick, "Comparing Numerical Methods ! for Ordinary Differential Equations," SIAM Journal on Numerical ! Analysis, Vol. 9, pp. 603-637, 1972. The analytical solution ! in terms of the numerical solution of Kepler's equation can be ! found there as well as in most discussions of the two body ! problem. The results for the particular choice of eccentricity, ! initial conditions, and interval are providedin TRUE_Y_AT_TEND. ! ! Global error assessment has been selected as a check on the reliability ! of the results. From the results generated it will be seen that the ! global error at the end of the run is about 2.3E-9, rather bigger than ! the local error tolerance TOL. This illustrates the point that at best ! one can anticipate global errors comparable to the tolerance. In point ! of fact, this problem is unstable at some points of the integration and ! the global error assessment reveals that the worst global error is ! considerably worse than the error at the end -- an example of the value ! of the global error assessment capability. use rksuite_90_prec contains function f(t,y) real(kind=wp), intent(in) :: t real(kind=wp), dimension(:,:), intent(in) :: y real(kind=wp), dimension(size(y,1),size(y,2)) :: f real(kind=wp) :: r_cubed r_cubed = sqrt( sum ( y(:,1)**2 ) )**3 f(:,1) = y(:,2); f(:,2) = -y(:,1)/r_cubed end function f end module two_body_f program two_body use two_body_f use rksuite_90 implicit none integer :: flag, i, totf real(kind=wp) :: t_start=0.0_wp, t_end=20.0_wp, tol=1.0e-10_wp, t_now, & ecc=0.9_wp, true_error, waste, max_error, t_max_error real(kind=wp), dimension(2,2) :: y_start, thres = 1.0e-13_wp, y_now, & y_maxvals, true_y_at_t_end, assessed_error type(rk_comm_real_2d) :: comm true_y_at_t_end(:,1) = (/-1.29526625098758_wp, 0.400393896379232_wp /) true_y_at_t_end(:,2) = (/-0.67753909247075_wp, -0.127083815427869_wp /) y_start(:,1) = (/ 1.0_wp - ecc, 0.0_wp /) y_start(:,2) = (/ 0.0_wp, sqrt((1.0_wp+ecc)/(1.0_wp-ecc)) /) call setup(comm,t_start,y_start,t_end,tol,thres,method='high',task='step', & error_assess=.true.) do call step_integrate(comm,f,t_now,y_now=y_now,flag=flag) ! write (*,'(5e14.5') t_now, y_now if (t_now==t_end .or. flag > 3) exit end do call statistics(comm,y_maxvals=y_maxvals,total_f_calls=totf,waste=waste) write (*,'(/a,1p2e9.2/21x,1p2e9.2)') & ' ymax ', (y_maxvals(i,1:2),i=1,2) write (*,'(/a,1pe10.2,0p/a,i10/a,f10.2)') & ' The integration reached ', t_now, & ' The cost of the integration in calls to F was', totf, & ' The fraction of failed steps was ', waste if (t_now==t_end) then true_error = maxval(abs(true_y_at_t_end - y_now)/max(abs(y_now),thres)) call global_error(comm,rms_error=assessed_error,max_error=max_error, & t_max_error=t_max_error) write (*,'(/(a,1p,e9.2/a,e9.2/a,e9.2,a,e9.2))') & ' At t = 20, the true error is',true_error, & ' assessed error is',maxval(assessed_error), & ' maximum error was',max_error,' at ',t_max_error end if call collect_garbage(comm) end program two_body SHAR_EOF fi # end of overwriting check if test -f 'eg_stiff_r0.r' then echo shar: will not over-write existing file "'eg_stiff_r0.r'" else cat << \SHAR_EOF > 'eg_stiff_r0.r' Have integrated to t = 0.186 The problem has been diagnosed as stiff SHAR_EOF fi # end of overwriting check if test -f 'eg_stiff_r0.f90' then echo shar: will not over-write existing file "'eg_stiff_r0.f90'" else cat << \SHAR_EOF > 'eg_stiff_r0.f90' module stiff_r0_f ! PROBLEM ! ! Demonstrate that a scalar ODE can be stiff ! ! SOLUTION ! ! This problem is considered in detail in ! L.F. Shampine,"numerical solution of Ordinary Differential Equations", ! Chapman and Hall, 1994, p384 ! We use the example there ! ! y' = J ( y - p(x) ) + p'(x), y(0) = A ! ! which has solution y(x) = ( A - p(0) ) exp(Jx) + p(x). For p(x) ! slowly varying and J large and negative this problem is stiff. We choose ! p(x) = cos(x), A = 0 and J = -10**(n) with n=3. ! ! We use RANGE_INTEGRATE and switch of the printing of error messages by ! setting MESSAGE=.FALSE. in the call to SETUP. ! ! NOTE ! ! This problem is for illustrative purposes. We use the real scalar ! dependent variable facility of rksuite_90 use rksuite_90_prec real(kind=wp) :: lam=-1000.0_wp, a=0.0_wp contains function f(t,y) real(kind=wp), intent(in) :: t real(kind=wp), intent(in) :: y real(kind=wp) :: f f = lam*(y-cos(t)) - sin(t) end function f function exact(t) real(kind=wp), intent(in) :: t real(kind=wp) :: exact exact = (a-cos(0.0_wp))*exp(lam*t) + cos(t) end function exact end module stiff_r0_f program stiff_r0 use stiff_r0_f use rksuite_90 implicit none integer :: flag real(kind=wp) :: t_start=0.0_wp, t_end=10.0_wp, tol=1.0e-4_wp, t_got, & thres = 1.0_wp, y_start, y_got type(rk_comm_real_0d) :: comm y_start = a call setup(comm,t_start,y_start,t_end,tol,thres,message=.false.) call range_integrate(comm,f,t_end,t_got,y_got=y_got,flag=flag) write (*,'(a,f6.3)') ' Have integrated to t = ',t_got if (flag==4) then write (*,'(a)') ' The problem has been diagnosed as stiff' else if (flag/=1) then write (*,'(a/a)') ' The problem was not diagnosed as stiff !!! ', & ' did you change the program...????' end if call collect_garbage(comm) end program stiff_r0 SHAR_EOF fi # end of overwriting check if test -f 'eg_simple_c0_a.r' then echo shar: will not over-write existing file "'eg_simple_c0_a.r'" else cat << \SHAR_EOF > 'eg_simple_c0_a.r' 0.785 0.7071068 0.7071068 0.7071068 0.7071068 1.571 0.0000000 0.0000000 1.0000000 1.0000000 2.356 -0.7071068 -0.7071068 0.7071068 0.7071068 3.142 -1.0000000 -1.0000000 0.0000000 0.0000000 3.927 -0.7071068 -0.7071068 -0.7071068 -0.7071068 4.712 0.0000000 0.0000000 -1.0000000 -1.0000000 5.498 0.7071068 0.7071068 -0.7071068 -0.7071068 6.283 1.0000000 1.0000000 0.0000000 0.0000000 ymax 1.00E+00 The cost of the integration in calls to F was 117 SHAR_EOF fi # end of overwriting check if test -f 'eg_simple_c0_a.f90' then echo shar: will not over-write existing file "'eg_simple_c0_a.f90'" else cat << \SHAR_EOF > 'eg_simple_c0_a.f90' module simple_c0_a_f ! PROBLEM ! ! Compute about seven correct figures in the solution of ! ! y' = i y ! ! on the range [0,2*pi] at intervals of length pi/4, given that y(0)=(1,0). ! ! SOLUTION ! ! This is a variant of "eg_simple_c0.f90", where here somewhat more ! accuracy is required. As for "eg_simple_c0.f90", this is a "usual task" ! that is more appropriately solved with RANGE_INTEGRATE. However, for ! illustrative purposes, we shall use STEP_INTEGRATE in conjunction with ! RESET_T_END. ! Although the code controls the local error rather than ! the true error, it is "tuned" so that the true error will be comparable ! to the local error tolerance for typical problems. Thus a relative error ! tolerance TOL of 5.0D-8 is appropriate. In this range of tolerances, ! METHOD = 'H' is likely to be the most efficient choice. ! The solution components are expected to get as large as 1.0D0. With ! this in mind, solution components smaller than, say, 1.0D-15 are not ! very interesting, and requiring seven correct figures then is not worth ! the cost. For this reason, the threshold values are specified as ! THRESHOLDS = 1.0D-15. When a solution component is smaller than this ! threshold, the code will control the local error to be no more than ! TOL*THRESHOLDS = 5.0D-23. The true value of y will be printed out for ! comparison. ! ! NOTE ! ! This problem is similar to that in simple.f90. Here, for illustrative ! purposes, we use the complex scalar dependent variable facility of ! rksuite_90. use rksuite_90_prec contains function f(t,y) real(kind=wp), intent(in) :: t complex(kind=wp), intent(in) :: y complex(kind=wp) :: f f = cmplx(0.0_wp,1.0_wp) * y end function f end module simple_c0_a_f program simple_c0_a use simple_c0_a_f use rksuite_90 implicit none integer :: nout, totf, l real(kind=wp) :: t_start=0.0_wp, t_end_final, tol=5.0e-5_wp, t_now, & t_want, pi, t_inc, thres = 1.0e-10_wp, y_maxvals complex(kind=wp) :: y_start, y_now type(rk_comm_complex_0d) :: comm y_start = cmplx(1.0_wp,0.0_wp) pi = 4.0_wp*atan(1.0_wp) t_end_final = 2.0_wp*pi nout = 8; t_inc = (t_end_final-t_start)/nout l = -7; t_want = t_end_final + l*t_inc call setup(comm,t_start,y_start,t_want,tol,thres,task='step',method='high') do call step_integrate(comm,f,t_now,y_now=y_now) if (t_want == t_now) then write (*,'(1x,f6.3,3x,f12.7,3x,f12.7)') t_want, real(y_now), cos(t_want) write (*,'(1x,9x,f12.7,3x,f12.7/)') aimag(y_now), sin(t_want) l = l + 1; t_want = t_end_final + l*t_inc call reset_t_end(comm,t_want) end if if (t_now == t_end_final) exit end do call statistics(comm,y_maxvals=y_maxvals,total_f_calls=totf) write (*,'(/a,1pe9.2)') & ' ymax ', y_maxvals write (*,'(/a,i10)') & ' The cost of the integration in calls to F was', totf call collect_garbage(comm) end program simple_c0_a SHAR_EOF fi # end of overwriting check if test -f 'eg_simple_c0.r' then echo shar: will not over-write existing file "'eg_simple_c0.r'" else cat << \SHAR_EOF > 'eg_simple_c0.r' 0.785 0.7071 0.7071 0.7071 0.7071 1.571 0.0000 0.0000 1.0000 1.0000 2.356 -0.7071 -0.7071 0.7071 0.7071 3.142 -1.0000 -1.0000 0.0000 0.0000 3.927 -0.7071 -0.7071 -0.7071 -0.7071 4.712 0.0000 0.0000 -1.0000 -1.0000 5.498 0.7071 0.7071 -0.7071 -0.7071 6.283 1.0000 1.0000 0.0000 0.0000 ymax 1.00E+00 The cost of the integration in calls to F was 95 SHAR_EOF fi # end of overwriting check if test -f 'eg_simple_c0.f90' then echo shar: will not over-write existing file "'eg_simple_c0.f90'" else cat << \SHAR_EOF > 'eg_simple_c0.f90' module simple_c0_f ! PROBLEM ! ! Compute about four correct figures in the solution of ! ! y' = i y ! ! on the range [0,2*pi] at intervals of length pi/4, given that y(0)=(1,0). ! ! SOLUTION ! ! This is a "usual task" that is more appropriately solved with ! RANGE_INTEGRATE. However, for illustrative purposes, we shall ! STEP_INTEGRATE in conjunction with INTERPOLATE. ! Although the code controls the local error rather than ! the true error, it is "tuned" so that the true error will be comparable ! to the local error tolerance for typical problems. Thus a relative error ! tolerance TOL of 5.0D-5 is appropriate. In this range of tolerances, ! METHOD = 'M' (the default) is likely to be the most efficient choice. ! The solution components are expected to get as large as 1.0D0. With ! this in mind, solution components smaller than, say, 1.0D-10 are not ! very interesting, and requiring five correct figures then is not worth ! the cost. For this reason, the threshold values are specified as ! THRESHOLDS = 1.0D-10. When a solution component is smaller than this ! threshold, the code will control the local error to be no more than ! TOL*THRESHOLDS = 5.0D-15. Answers will be computed at equally spaced ! points, and the true value of y will be printed out for comparison. ! ! NOTE ! ! This problem is similar to that in simple.f90. Here, for illustrative ! purposes, we use the complex scalar dependent variable facility of ! rksuite_90. The number of function evaluations should differ to that ! required by simple.f90 since, there the error in two real components is ! controlled, and here the error in a single complex component is ! controlled. use rksuite_90_prec contains function f(t,y) real(kind=wp), intent(in) :: t complex(kind=wp), intent(in) :: y complex(kind=wp) :: f f = cmplx(0.0_wp,1.0_wp) * y end function f end module simple_c0_f program simple_c0 use simple_c0_f use rksuite_90 implicit none integer :: nout, totf, l real(kind=wp) :: t_start=0.0_wp, t_end, tol=5.0e-5_wp, t_now, & t_want, pi, t_inc, thres = 1.0e-10_wp, y_maxvals complex(kind=wp) :: y_start, y_now, y_want type(rk_comm_complex_0d) :: comm y_start = cmplx(1.0_wp,0.0_wp) pi = 4.0_wp*atan(1.0_wp) t_end = 2.0_wp*pi call setup(comm,t_start,y_start,t_end,tol,thres,task='step') nout = 8; t_inc = (t_end-t_start)/nout l = -7; t_want = t_end + l*t_inc do call step_integrate(comm,f,t_now,y_now=y_now) do if (t_want > t_now) exit call interpolate(comm,f,t_want,y_want=y_want) write (*,'(1x,f6.3,3x,f9.4,3x,f9.4)') t_want, real(y_want), cos(t_want) write (*,'(1x,9x,f9.4,3x,f9.4/)') aimag(y_want), sin(t_want) l = l + 1; t_want = t_end + l*t_inc end do if (t_now == t_end) exit end do call statistics(comm,y_maxvals=y_maxvals,total_f_calls=totf) write (*,'(/a,1pe9.2)') & ' ymax ', y_maxvals write (*,'(/a,i10)') & ' The cost of the integration in calls to F was', totf call collect_garbage(comm) end program simple_c0 SHAR_EOF fi # end of overwriting check if test -f 'eg_simple.r' then echo shar: will not over-write existing file "'eg_simple.r'" else cat << \SHAR_EOF > 'eg_simple.r' 0.785 0.7071 0.7071 0.7071 0.7071 1.571 1.0000 1.0000 0.0000 0.0000 2.356 0.7071 0.7071 -0.7071 -0.7071 3.142 0.0000 0.0000 -1.0000 -1.0000 3.927 -0.7071 -0.7071 -0.7071 -0.7071 4.712 -1.0000 -1.0000 0.0000 0.0000 5.498 -0.7071 -0.7071 0.7071 0.7071 6.283 0.0000 0.0000 1.0000 1.0000 ymax 9.90E-01 1.00E+00 The cost of the integration in calls to F was 102 SHAR_EOF fi # end of overwriting check if test -f 'eg_simple.f90' then echo shar: will not over-write existing file "'eg_simple.f90'" else cat << \SHAR_EOF > 'eg_simple.f90' module simple_f ! PROBLEM ! ! Compute about four correct figures in the solution of ! ! y'' = -y ! ! on the range [0,2*pi] at intervals of length pi/4, given that y(0)=0, ! y'(0)=1. ! ! SOLUTION ! ! Let y1 = y and y2 = y' to obtain the first order system ! ! y1' = y2 with initial values y1 = 0 ! y2' = - y1 y2 = 1 ! ! This is a "usual task" that is appropriately solved with ! RANGE_INTEGRATE. Although the code controls the local error rather than ! the true error, it is "tuned" so that the true error will be comparable ! to the local error tolerance for typical problems. Thus a relative error ! tolerance TOL of 5.0D-5 is appropriate. In this range of tolerances, ! METHOD = 'M' (the default) is likely to be the most efficient choice. ! The solution components are expected to get as large as 1.0D0. With ! this in mind, solution components smaller than, say, 1.0D-10 are not ! very interesting, and requiring five correct figures then is not worth ! the cost. For this reason, the threshold values are specified as ! THRESHOLDS = 1.0D-10. When a solution component is smaller than this ! threshold, the code will control the local error to be no more than ! TOL*THRESHOLDS = 5.0D-15. Answers will be computed at equally spaced ! points, and the true values of y and y' will be printed out for ! comparison. use rksuite_90_prec contains function f(t,y) real(kind=wp), intent(in) :: t real(kind=wp), dimension(:), intent(in) :: y real(kind=wp), dimension(size(y)) :: f f(:) = (/ y(2), -y(1) /) end function f end module simple_f program simple use simple_f use rksuite_90 implicit none integer :: nout, totf, l real(kind=wp) :: t_start=0.0_wp, t_end, tol=5.0e-5_wp, t_got, & t_want, pi, t_inc real(kind=wp), dimension(2) :: y_start = (/ 0.0_wp,1.0_wp /), & thres = 1.0e-10_wp, y_got, y_maxvals type(rk_comm_real_1d) :: comm pi = 4.0_wp*atan(1.0_wp) t_end = 2.0_wp*pi call setup(comm,t_start,y_start,t_end,tol,thres) nout = 8; t_inc = (t_end-t_start)/nout do l = 1, nout t_want = t_end + (l-nout)*t_inc call range_integrate(comm,f,t_want,t_got,y_got=y_got) write (*,'(1x,f6.3,3x,f9.4,3x,f9.4)') t_got, y_got(1), sin(t_got) write (*,'(1x,9x,f9.4,3x,f9.4/)') y_got(2), cos(t_got) end do call statistics(comm,y_maxvals=y_maxvals,total_f_calls=totf) write (*,'(/a,1p2e9.2)') & ' ymax ', y_maxvals write (*,'(/a,i10)') & ' The cost of the integration in calls to F was', totf call collect_garbage(comm) end program simple SHAR_EOF fi # end of overwriting check if test -f 'eg_quadrature.r' then echo shar: will not over-write existing file "'eg_quadrature.r'" else cat << \SHAR_EOF > 'eg_quadrature.r' ans = 3.9524 The cost of the outer integration in calls to gp was 119 SHAR_EOF fi # end of overwriting check if test -f 'eg_quadrature.f90' then echo shar: will not over-write existing file "'eg_quadrature.f90'" else cat << \SHAR_EOF > 'eg_quadrature.f90' module quadrature_inner_f use rksuite_90_prec implicit none contains function hp(x,h) real(kind=wp), intent(in) :: x real(kind=wp), intent(in) :: h real(kind=wp) :: hp hp = (1.0_wp + x)**3 end function hp end module quadrature_inner_f module quadrature_f ! PROBLEM ! ! Compute to five figures the integral ! ! 1 z ! I I (x+1)**3*(z+1)**2 dxdz ! 0 0 ! ! SOLUTION ! ! This simple nested quadrature problem can be solved by recursive use ! of the integrators provided in rksuite_90. First, we define the ! h'(x) ! ! h'(x) = (x+1)**3 with h(0) = 0.0 ! ! which is an initial value problem. The required integral can then be ! expressed as ! ! 1 ! I h(z)*(z+1)**2 dz ! 0 ! ! Similarly, we define g'(z) ! ! g'(z) = h(z)*(z+1)**2 with g(0) = 0.0, ! ! another initial value problem. Then, the answer to the problem is ! given by g(1) ! ! Hence we use RANGE_INTEGRATE to integrate g'(z) to compute the value ! g(1) and in the definition of g'(z) (the procedure GP) we use ! RANGE_INTEGRATE again (recursively) to integrate h'(x) (the procedure ! HP) to compute the required values of h(z). ! ! Given the nature of the integrand and the area of integration we expect ! the answer to be of order of magnitude 1. Hence, we use a relative ! TOLERANCE = 5.0e-5 and a value for THRESHOLD = 1.0e-10. use rksuite_90 use quadrature_inner_f implicit none real(kind=wp) :: x_lower=0.0_wp, x_inner_got, h_of_z_got, & tol=5.0e-6, thresh = 1.0e-10_wp type(rk_comm_real_0d) :: comm_inner contains function gp(z,g) real(kind=wp), intent(in) :: z real(kind=wp), intent(in) :: g real(kind=wp) :: gp if (z == 0.0_wp) then h_of_z_got = 0.0_wp else call setup(comm_inner,x_lower,0.0_wp,z,tol,thresh) call range_integrate(comm_inner,hp,z,x_inner_got,y_got=h_of_z_got) call collect_garbage(comm_inner) end if gp = h_of_z_got * (1.0_wp + z)**2 end function gp end module quadrature_f program quadrature use quadrature_f implicit none integer :: totf real(kind=wp) :: z_lower=0.0_wp, z_upper=1.0_wp, z_got, g_got type(rk_comm_real_0d) :: comm call setup(comm,z_lower,0.0_wp,z_upper,tol,thresh) call range_integrate(comm,gp,z_upper,z_got,y_got=g_got) write (*,'(/a,f9.4)') ' ans = ',g_got call statistics(comm,total_f_calls=totf) write (*,'(/a,i10)') & ' The cost of the outer integration in calls to gp was', totf call collect_garbage(comm) end program quadrature SHAR_EOF fi # end of overwriting check if test -f 'eg_orr_somm_r.r' then echo shar: will not over-write existing file "'eg_orr_somm_r.r'" else cat << \SHAR_EOF > 'eg_orr_somm_r.r' At t = 0.2367E-01 in 30 steps and imag(y'(t)) = 0.5398E+04 SHAR_EOF fi # end of overwriting check if test -f 'eg_orr_somm_r.f90' then echo shar: will not over-write existing file "'eg_orr_somm_r.f90'" else cat << \SHAR_EOF > 'eg_orr_somm_r.f90' module orr_somm_r_f ! PROBLEM ! ! Compute a trial solution for the Orr-Sommerfeld equation for plane ! Poiseuille flow ! ! y'''' - 2 k**2 y'' + k**4 y - i k R [(u(t) - c)(y'' - k**2 y) - u''(t)y] = 0 ! ! y(t) is the amplitude of the stream function, R is the Reynolds number, ! k is the wave number of the disturbance, i is sqrt(-1), c is the eigenvalue, ! and, in this instance, the laminar velocity profile is u(t) = 1 - t**2 ! ! The boundary conditions are y'(0) = y'''(0) = 0 = y(1) = y'(1). ! ! Use c = 0.06659252 - 0.01398327 i, k = 1, R = 1.0e+06, and the ! arbitrary initial conditions y(0) = y''(0) = 1. ! ! SOLUTION ! ! The real objective in solving this equation is to compute eigenvalues ! c such that the imaginary part of y'(1) = 0. See "Computational ! solution of two point BVPs via orthonormalization", MR Scott and HA Watts ! SIAM J. Num. Anal. 14, pp 40-40, 1977 and the references therein. ! ! We recast the problem in to a system of real first order ODEs using ! y1 = real(y), y3 = real(y'), y5 = real(y''), y7 = real(y'''), ! y2 = imag(y), y4 = imag(y'), y6 = imag(y''), y8 = imag(y'''), to obtain ! ! y1' = y3 ! y2' = y4 ! y3' = y5 ! y4' = y6 ! y5' = y7 ! y6' = y8 ! y7' = 2 k**2 y5 - k**4 y1 - k R * ! [(u(t) - cr)(y6 - k**2 y2) + ci (y5 - k**2 y1) - u''(t) y2 ] ! y8' = 2 k**2 y6 - k**4 y2 + k R * ! [(u(t) - cr)(y5 - k**2 y1) - ci (y6 - k**2 y2) - u''(t) y1 ] ! ! given y1(0) = 1, y3(0) = 0, y5(0) = 1, y7(0) = 0, ! y2(0) = 0, y4(0) = 0, y6(0) = 0, y8(0) = 0, ! ! where cr = real(c) and ci = imag(c). ! ! For such an arbitrary selection of initial conditions it is likely that ! the computed solution will "blow up". We use STEP_INTEGRATE to monitor ! the progress of the computation and if any solution component exceeds ! 1.0e+10 in magnitude the computation is abandoned. ! ! NOTE ! ! This example demonstrates the use of the real dependent variable ! facility offerred by rksuite_90 when in fact the complex dependent ! variable facility would have been more appropriate. The example ! eg_orr_somm.f90 illustrates how the same system can be handled more ! conveniently. use rksuite_90_prec real(kind=wp), parameter :: k=1.0_wp, k2=k*k, k4=k2*k2, r=1.0e+06_wp, kr=k*r, & cr=0.06659252_wp, ci=-0.01398327_wp contains function f(t,y) real(kind=wp), intent(in) :: t real(kind=wp), dimension(:), intent(in) :: y real(kind=wp), dimension(size(y)) :: f real(kind=wp) :: ucr, y2rk2yr, y2ik2yi f(1:6) = y(3:8) ucr = u(t) - cr y2rk2yr = y(5) - k2*y(1) y2ik2yi = y(6) - k2*y(2) f(7) = 2.0_wp*k2*y(5) - k4*y(1) - k*r* & ( ucr*y2ik2yi + ci*y2rk2yr - upp(t)*y(2) ) f(8) = 2.0_wp*k2*y(6) - k4*y(2) + k*r* & ( ucr*y2rk2yr - ci*y2ik2yi - upp(t)*y(1) ) end function f function u(t) real(kind=wp), intent(in) :: t real(kind=wp) :: u u = 1.0_wp - t**2 end function u function upp(t) real(kind=wp), intent(in) :: t real(kind=wp) :: upp upp = -2.0_wp end function upp end module orr_somm_r_f program orr_somm_r use orr_somm_r_f use rksuite_90 implicit none integer :: steps real(kind=wp) :: t_start=0.0_wp, t_end=1.0_wp, tol=5.0e-5_wp, t_now real(kind=wp), dimension(8) :: thres = 1.0e-10_wp real(kind=wp), dimension(8) :: y_start, y_now, yderiv_now type(rk_comm_real_1d) :: comm y_start(:) = 0.0_wp y_start(1) = 1.0_wp y_start(5) = 1.0_wp call setup(comm,t_start,y_start,t_end,tol,thres,task='step') do call step_integrate(comm,f,t_now,y_now=y_now,yderiv_now=yderiv_now) if (maxval(abs(y_now(:))) > 1.0e+10_wp .or. t_now == t_end) exit end do call statistics(comm,num_succ_steps=steps) write (*,'(a,e11.4,a,i4,a,e11.4)') & " At t = ",t_now," in ",steps," steps and imag(y'(t)) = ",y_now(4) call collect_garbage(comm) end program orr_somm_r SHAR_EOF fi # end of overwriting check if test -f 'eg_orr_somm.r' then echo shar: will not over-write existing file "'eg_orr_somm.r'" else cat << \SHAR_EOF > 'eg_orr_somm.r' At t = 0.2424E-01 in 27 steps with imag(y'(t)) = -0.2564E+04 SHAR_EOF fi # end of overwriting check if test -f 'eg_orr_somm.f90' then echo shar: will not over-write existing file "'eg_orr_somm.f90'" else cat << \SHAR_EOF > 'eg_orr_somm.f90' module orr_somm_f ! PROBLEM ! ! Compute a trial solution for the Orr-Sommerfeld equation for plane ! Poiseuille flow ! ! y'''' - 2 k**2 y'' + k**4 y - i k R [(u(t) - c)(y'' - k**2 y) - u''(t)y] = 0 ! ! y(t) is the amplitude of the stream function, R is the Reynolds number, ! k is the wave number of the disturbance, i is sqrt(-1), c is the eigenvalue, ! and, in this instance, the laminar velocity profile is u(t) = 1 - t**2 ! ! The boundary conditions are y'(0) = y'''(0) = 0 = y(1) = y'(1). ! ! Use c = 0.06659252 - 0.01398327 i, k = 1, R = 1.0e+06, and the ! arbitrary initial conditions y(0) = y''(0) = 1. ! ! SOLUTION ! ! The real objective in solving this equation is to compute eigenvalues ! c such that the imaginary part of y'(1) = 0. See "Computational ! solution of two point BVPs via orthonormalization", MR Scott and HA Watts ! SIAM J. Num. Anal. 14, pp 40-40, 1977 and the references therein. ! ! We recast the problem in to a system of complex first order ODEs using ! y1 = y, y2 = y', y3 = y'', y4 = y''', to obtain ! ! y1' = y2 ! y2' = y3 ! y3' = y4 ! y4' = 2 k**2 y3 - k**4 y1 + i k R [(u(t) - c)(y3 - k**2 y1) - u''(t)y1] ! ! given y1(0) = 1, y2(0) = 0, y3(0) = 1, y4(0) = 0, ! ! and use the complex dependent variable facility of rksuite_90. ! ! For such an arbitrary selection of initial conditions it is likely that ! the computed solution will "blow up". We use STEP_INTEGRATE to monitor ! the progress of the computation and if any solution component exceeds ! 1.0e+10 in magnitude the computation is abandoned. ! ! NOTE ! ! This example demonstrates the use of the complex dependent variable ! facility offerred by rksuite_90. The example eg_orr_somm_r.f90 ! illustrates how the same system can be handled in the more traditional ! (and more complicated) manner of treating separately the real imaginary ! parts of y as a larger system of real dependent variables. use rksuite_90_prec real(kind=wp), parameter :: k=1.0_wp, k2=k*k, k4=k2*k2, r=1.0e+06_wp, kr=k*r complex(kind=wp) :: c = (0.06659252_wp, -0.01398327_wp) contains function f(t,y) real(kind=wp), intent(in) :: t complex(kind=wp), dimension(:), intent(in) :: y complex(kind=wp), dimension(size(y)) :: f f(1:3) = y(2:4) f(4) = 2.0_wp*k2*y(3) - k4*y(1) + cmplx(0.0_wp,1.0_wp)*kr* & ( (u(t)-c)*(y(3)-k**2*y(1)) - upp(t)*y(1) ) end function f function u(t) real(kind=wp), intent(in) :: t real(kind=wp) :: u u = 1.0_wp - t**2 end function u function upp(t) real(kind=wp), intent(in) :: t real(kind=wp) :: upp upp = -2.0_wp end function upp end module orr_somm_f program orr_somm use orr_somm_f use rksuite_90 implicit none integer :: steps real(kind=wp) :: t_start=0.0_wp, t_end=1.0_wp, tol=5.0e-5_wp, t_now real(kind=wp), dimension(4) :: thres = 1.0e-10_wp complex(kind=wp), dimension(4) :: y_start, y_now, yderiv_now type(rk_comm_complex_1d) :: comm y_start(1) = 1.0_wp y_start(2) = 0.0_wp y_start(3) = 1.0_wp y_start(4) = 0.0_wp call setup(comm,t_start,y_start,t_end,tol,thres,task='step') do call step_integrate(comm,f,t_now,y_now=y_now,yderiv_now=yderiv_now) if (maxval(abs(y_now(:))) > 1.0e+10_wp .or. t_now == t_end) exit end do call statistics(comm,num_succ_steps=steps) write (*,'(a,e11.4,a,i4,a,e11.4)') & " At t = ",t_now," in ",steps," steps with imag(y'(t)) = ",aimag(y_now(2)) call collect_garbage(comm) end program orr_somm SHAR_EOF fi # end of overwriting check if test -f 'eg_invar_imbed.r' then echo shar: will not over-write existing file "'eg_invar_imbed.r'" else cat << \SHAR_EOF > 'eg_invar_imbed.r' t approx exact 4.362 0.12757E-01 0.12757E-01 3.807 0.22204E-01 0.22204E-01 3.257 0.38513E-01 0.38513E-01 2.706 0.66808E-01 0.66808E-01 2.154 0.11602E+00 0.11602E+00 1.599 0.20207E+00 0.20207E+00 1.035 0.35532E+00 0.35532E+00 0.517 0.59609E+00 0.59609E+00 0.000 0.10000E+01 0.10000E+01 SHAR_EOF fi # end of overwriting check if test -f 'eg_invar_imbed.f90' then echo shar: will not over-write existing file "'eg_invar_imbed.f90'" else cat << \SHAR_EOF > 'eg_invar_imbed.f90' module invar_imbed_fw use rksuite_90_prec contains function g(t,y) real(kind=wp), intent(in) :: t real(kind=wp), dimension(:), intent(in) :: y real(kind=wp), dimension(size(y)) :: g g(:) = (/ 1.0_wp - y(1)**2, -y(1)*y(2) /) end function g end module invar_imbed_fw module invar_imbed_bw use rksuite_90 use invar_imbed_fw real(kind=wp) :: u, v real(kind=wp) :: t_fw_start=0.0_wp, t_fw_end, t_fw_got, tol=1.0e-5_wp real(kind=wp), dimension(2) :: thres(2)=1.0e-20_wp, & y_fw_start, y_fw_got type(rk_comm_real_1d) :: comm_fw contains function f(t,y) real(kind=wp), intent(in) :: t real(kind=wp), intent(in) :: y real(kind=wp) :: f t_fw_end = t u = 0.0_wp; v = 1.0_wp; y_fw_start = (/u, v/) if (t_fw_end /= t_fw_start) then call setup(comm_fw,t_fw_start,y_fw_start,t_fw_end,tol,thres) call range_integrate(comm_fw,g,t_fw_end,t_fw_got,y_got=y_fw_got) call collect_garbage(comm_fw) u = y_fw_got(1); v = y_fw_got(2) end if f = u*y + v end function f end module invar_imbed_bw program invar_imbed use invar_imbed_bw implicit none real(kind=wp) :: t_bw_start=5.0_wp, t_bw_end, y_bw_start, & y_bw_now, t_bw_now type(rk_comm_real_0d) :: comm_bw t_fw_end = t_bw_start; y_fw_start = (/0.0_wp, 1.0_wp/) call setup(comm_fw,t_fw_start,y_fw_start,t_fw_end,tol,thres) call range_integrate(comm_fw,g,t_fw_end,t_fw_got,y_got=y_fw_got) call collect_garbage(comm_fw) u = y_fw_got(1); v = y_fw_got(2) y_bw_start = (exp(-t_fw_end) - v)/u t_bw_end = t_fw_start call setup(comm_bw,t_bw_start,y_bw_start,t_bw_end,tol,thres(1),task='s') write (*,'(a/)') ' t approx exact' do call step_integrate(comm_bw,f,t_bw_now,y_now=y_bw_now) write (*,'(1x,f6.3,3x,e12.5,3x,e12.5)') t_bw_now, u*y_bw_now+v, & exp(-t_bw_now) if (t_bw_now==t_bw_end) exit end do call collect_garbage(comm_bw) end program invar_imbed SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Derived' then mkdir 'Derived' fi cd 'Derived' if test -f 'rksuite_90_r2.f90' then echo shar: will not over-write existing file "'rksuite_90_r2.f90'" else cat << \SHAR_EOF > 'rksuite_90_r2.f90' module rksuite_90_prec ! ! Part of rksuite_90 v1.2 (December 1995) ! software for initial value problems in ODEs ! ! Authors: R.W. Brankin (NAG Ltd., Oxford, England) ! I. Gladwell (Math Dept., SMU, Dallas, TX, USA) ! see main doc for contact details ! integer, parameter :: wp = selected_real_kind(10,50) end module rksuite_90_prec module rksuite_90 ! ! Part of rksuite_90 v1.2 (December 1995) ! software for initial value problems in ODEs ! ! Authors: R.W. Brankin (NAG Ltd., Oxford, England) ! I. Gladwell (Math Dept., SMU, Dallas, TX, USA) ! see main doc for contact details ! use rksuite_90_prec, only:wp implicit none private public :: wp, setup, range_integrate, step_integrate, interpolate, & global_error, statistics, reset_t_end, collect_garbage, & set_stop_on_fatal, get_saved_fatal !starta! public :: rk_comm_real_2d type rk_comm_real_2d private real(kind=wp) :: t, t_old, t_start, t_end, dir !indep! real(kind=wp) :: h, h_old, h_start, h_average !indep! real(kind=wp) :: tol integer :: f_count, full_f_count, step_count, bad_step_count logical :: at_t_start, at_t_end real(kind=wp), dimension(:,:), pointer :: thresh, weights, ymax !shp-dep! real(kind=wp), dimension(:,:), pointer :: scratch, y, yp, y_new !dep! real(kind=wp), dimension(:,:), pointer :: y_old, yp_old, v0, v1 !dep! real(kind=wp), dimension(:,:), pointer :: err_estimates, v2, v3 !dep! real(kind=wp), dimension(:,:), pointer :: vtemp !dep! real(kind=wp), dimension(:,:,:), pointer :: stages !dep! real(kind=wp) :: a(13,13), b(13), c(13), bhat(13), r(11,6), e(7) integer :: ptr(13), no_of_stages, rk_method, intrp_degree logical :: intrp_able, intrp_needs_stages real(kind=wp) :: toosml, cost, safety, expon, stability_radius, tan_angle, & rs, rs1, rs2, rs3, rs4 integer :: order, last_stage, max_stiff_iters, no_of_ge_steps logical :: fsal real(kind=wp) :: ge_max_contrib real(kind=wp) :: t_ge_max_contrib !indep! integer :: ge_f_count real(kind=wp), dimension(:,:), pointer :: ge_assess !shp-dep! real(kind=wp), dimension(:,:), pointer :: ge_y, ge_yp, ge_y_new !dep! real(kind=wp), dimension(:,:), pointer :: ge_err_estimates !dep! real(kind=wp), dimension(:,:,:), pointer :: ge_stages !dep! logical :: erason, erasfl real(kind=wp) :: mcheps, dwarf, round_off, sqrrmc, cubrmc, sqtiny integer :: outch logical :: print_message, use_range character(len=80) :: rec(10) real(kind=wp) :: tlast, range_t_end !indep! real(kind=wp), dimension(:,:), pointer :: xstage, ytemp !dep! real(kind=wp), dimension(:,:,:), pointer :: p !dep! integer :: stiff_bad_step_count, hit_t_end_count real(kind=wp) :: errold logical :: chkeff, phase2 integer, dimension(7) :: save_states logical :: stop_on_fatal, saved_fatal_err end type rk_comm_real_2d !enda! interface setup module procedure setup_r2 end interface interface range_integrate module procedure range_integrate_r2 end interface interface step_integrate module procedure step_integrate_r2 end interface interface statistics module procedure statistics_r2 end interface interface global_error module procedure global_error_r2 end interface interface reset_t_end module procedure reset_t_end_r2 end interface interface interpolate module procedure interpolate_r2 end interface interface set_stop_on_fatal module procedure set_stop_on_fatal_r2 end interface interface get_saved_fatal module procedure get_saved_fatal_r2 end interface interface collect_garbage module procedure collect_garbage_r2 end interface contains !startb! subroutine machine_const(round_off,sqrrmc,cubrmc,sqtiny,outch) ! ! Part of rksuite_90 v1.2 (December 1995) ! software for initial value problems in ODEs ! ! Authors: R.W. Brankin (NAG Ltd., Oxford, England) ! I. Gladwell (Math Dept., SMU, Dallas, TX, USA) ! see main doc for contact details ! real(kind=wp), intent(out) :: round_off, sqrrmc, cubrmc, sqtiny integer, intent(out) :: outch ! real(kind=wp) :: dummy real(kind=wp), parameter :: third=1.0_wp/3.0_wp, ten=10.0_wp ! outch = 6 ! round_off = ten*epsilon(dummy) sqrrmc = sqrt(epsilon(dummy)) cubrmc = epsilon(dummy)**third sqtiny = sqrt(tiny(dummy)) ! end subroutine machine_const subroutine method_const(rk_method, a, b, c, bhat, r, e, ptr, no_of_stages, & intrp_degree, intrp_able, intrp_needs_stages, & cost, safety, expon, stability_radius, & tan_angle, rs, rs1, rs2, rs3, rs4, order, last_stage, & max_stiff_iters, no_of_ge_steps, fsal, cdiff) ! ! Part of rksuite_90 v1.2 (December 1995) ! software for initial value problems in ODEs ! ! Authors: R.W. Brankin (NAG Ltd., Oxford, England) ! I. Gladwell (Math Dept., SMU, Dallas, TX, USA) ! see main doc for contact details ! integer, intent(in) :: rk_method real(kind=wp), intent(out) :: a(13,13), b(13), c(13), bhat(13), r(11,6), e(7) integer, intent(out) :: ptr(13), no_of_stages, intrp_degree logical, intent(out) :: intrp_able, intrp_needs_stages real(kind=wp), intent(out) :: cost, safety, expon, stability_radius, & tan_angle, rs, rs1, rs2, rs3, rs4, cdiff integer, intent(out) :: order, last_stage, max_stiff_iters, no_of_ge_steps logical, intent(out) :: fsal ! integer :: i real(kind=wp), parameter :: fivepc=0.05_wp, one=1.0_wp, two=2.0_wp, & fifty=50.0_wp ! select case (rk_method) case(1) ! ! METHD = 1. ! This pair is from "A 3(2) Pair of Runge-Kutta Formulas" by P. Bogacki ! and L.F. Shampine, Appl. Math. Lett., 2, pp. 321-325, 1989. The authors ! are grateful to P. Bogacki for his assistance in implementing the pair. ! no_of_stages = 4; fsal = .true.; order = 2 tan_angle = 8.9_wp; stability_radius = 2.3_wp safety = 0.8_wp; intrp_able = .true.; intrp_degree = 3 intrp_needs_stages = .false.; no_of_ge_steps = 3 ! ptr(1:4) = (/ 0, 1, 2, 3 /) ! a(2,1) = 1.0_wp/2.0_wp a(3,1) = 0.0_wp a(3,2) = 3.0_wp/4.0_wp a(4,1) = 2.0_wp/9.0_wp a(4,2) = 1.0_wp/3.0_wp a(4,3) = 4.0_wp/9.0_wp ! ! The coefficients BHAT refer to the formula used to advance the ! integration, here the one of order 3. The coefficients B refer ! to the other formula, here the one of order 2. For this pair, BHAT ! is not needed since FSAL = .TRUE. ! b(1) = 7.0_wp/24.0_wp b(2) = 1.0_wp/4.0_wp b(3) = 1.0_wp/3.0_wp b(4) = 1.0_wp/8.0_wp ! c(1) = 0.0_wp c(2) = 1.0_wp/2.0_wp c(3) = 3.0_wp/4.0_wp c(4) = 1.0_wp ! case (2) ! ! METHD = 2 ! This pair is from "An Efficient Runge-Kutta (4,5) Pair" by P. Bogacki ! and L.F. Shampine, Rept. 89-20, Math. Dept., Southern Methodist ! University, Dallas, Texas, USA, 1989. The authors are grateful to ! P. Bogacki for his assistance in implementing the pair. Shampine and ! Bogacki subsequently modified the formula to enhance the reliability of ! the pair. The original fourth order formula is used in an estimate of ! the local error. If the step fails, the computation is broken off. If ! the step is acceptable, the first evaluation of the next step is done, ! i.e., the pair is implemented as FSAL and the local error of the step ! is again estimated with a fourth order formula using the additional data. ! The step must succeed with both estimators to be accepted. When the ! second estimate is formed, it is used for the subsequent adjustment of ! the step size because it is of higher quality. The two fourth order ! formulas are well matched to leading order, and only exceptionally do ! the estimators disagree -- problems with discontinuous coefficients are ! handled more reliably by using two estimators as is global error ! estimation. ! no_of_stages = 8; fsal = .true.; order = 4 tan_angle = 5.2_wp; stability_radius = 3.9_wp safety = 0.8_wp; intrp_able = .true. intrp_needs_stages = .true.; intrp_degree = 6 no_of_ge_steps = 2 ! ptr(1:8) = (/ 0, 1, 2, 3, 4, 5, 6, 7 /) ! a(2,1) = 1.0_wp/6.0_wp a(3,1) = 2.0_wp/27.0_wp a(3,2) = 4.0_wp/27.0_wp a(4,1) = 183.0_wp/1372.0_wp a(4,2) = -162.0_wp/343.0_wp a(4,3) = 1053.0_wp/1372.0_wp a(5,1) = 68.0_wp/297.0_wp a(5,2) = -4.0_wp/11.0_wp a(5,3) = 42.0_wp/143.0_wp a(5,4) = 1960.0_wp/3861.0_wp a(6,1) = 597.0_wp/22528.0_wp a(6,2) = 81.0_wp/352.0_wp a(6,3) = 63099.0_wp/585728.0_wp a(6,4) = 58653.0_wp/366080.0_wp a(6,5) = 4617.0_wp/20480.0_wp a(7,1) = 174197.0_wp/959244.0_wp a(7,2) = -30942.0_wp/79937.0_wp a(7,3) = 8152137.0_wp/19744439.0_wp a(7,4) = 666106.0_wp/1039181.0_wp a(7,5) = -29421.0_wp/29068.0_wp a(7,6) = 482048.0_wp/414219.0_wp a(8,1) = 587.0_wp/8064.0_wp a(8,2) = 0.0_wp a(8,3) = 4440339.0_wp/15491840.0_wp a(8,4) = 24353.0_wp/124800.0_wp a(8,5) = 387.0_wp/44800.0_wp a(8,6) = 2152.0_wp/5985.0_wp a(8,7) = 7267.0_wp/94080.0_wp ! ! The coefficients B refer to the formula of order 4. ! b(1) = 2479.0_wp/34992.0_wp b(2) = 0.0_wp b(3) = 123.0_wp/416.0_wp b(4) = 612941.0_wp/3411720.0_wp b(5) = 43.0_wp/1440.0_wp b(6) = 2272.0_wp/6561.0_wp b(7) = 79937.0_wp/1113912.0_wp b(8) = 3293.0_wp/556956.0_wp ! ! The coefficients E refer to an estimate of the local error based on ! the first formula of order 4. It is the difference of the fifth order ! result, here located in A(8,:), and the fourth order result. By ! construction both E(2) and E(7) are zero. ! e(1) = -3.0_wp/1280.0_wp e(2) = 0.0_wp e(3) = 6561.0_wp/632320.0_wp e(4) = -343.0_wp/20800.0_wp e(5) = 243.0_wp/12800.0_wp e(6) = -1.0_wp/95.0_wp e(7) = 0.0_wp ! c(1) = 0.0_wp c(2) = 1.0_wp/6.0_wp c(3) = 2.0_wp/9.0_wp c(4) = 3.0_wp/7.0_wp c(5) = 2.0_wp/3.0_wp c(6) = 3.0_wp/4.0_wp c(7) = 1.0_wp c(8) = 1.0_wp ! ! To do interpolation with this pair, some extra stages have to be computed. ! The following additional A and C coefficients are for this purpose. ! In addition there is an array R that plays a role for interpolation ! analogous to that of BHAT for the basic step. ! c(9) = 1.0_wp/2.0_wp c(10) = 5.0_wp/6.0_wp c(11) = 1.0_wp/9.0_wp ! a(9,1) = 455.0_wp/6144.0_wp a(10,1) = -837888343715.0_wp/13176988637184.0_wp a(11,1) = 98719073263.0_wp/1551965184000.0_wp a(9,2) = 0.0_wp a(10,2) = 30409415.0_wp/52955362.0_wp a(11,2) = 1307.0_wp/123552.0_wp a(9,3) = 10256301.0_wp/35409920.0_wp a(10,3) = -48321525963.0_wp/759168069632.0_wp a(11,3) = 4632066559387.0_wp/70181753241600.0_wp a(9,4) = 2307361.0_wp/17971200.0_wp a(10,4) = 8530738453321.0_wp/197654829557760.0_wp a(11,4) = 7828594302389.0_wp/382182512025600.0_wp a(9,5) = -387.0_wp/102400.0_wp a(10,5) = 1361640523001.0_wp/1626788720640.0_wp a(11,5) = 40763687.0_wp/11070259200.0_wp a(9,6) = 73.0_wp/5130.0_wp a(10,6) = -13143060689.0_wp/38604458898.0_wp a(11,6) = 34872732407.0_wp/224610586200.0_wp a(9,7) = -7267.0_wp/215040.0_wp a(10,7) = 18700221969.0_wp/379584034816.0_wp a(11,7) = -2561897.0_wp/30105600.0_wp a(9,8) = 1.0_wp/32.0_wp a(10,8) = -5831595.0_wp/847285792.0_wp a(11,8) = 1.0_wp/10.0_wp a(10,9) = -5183640.0_wp/26477681.0_wp a(11,9) = -1.0_wp/10.0_wp a(11,10) = -1403317093.0_wp/11371610250.0_wp ! r(1:11,1) = 0.0_wp; r(2,1:6) = 0.0_wp r(1,6) = -12134338393.0_wp/1050809760.0_wp r(1,5) = -1620741229.0_wp/50038560.0_wp r(1,4) = -2048058893.0_wp/59875200.0_wp r(1,3) = -87098480009.0_wp/5254048800.0_wp r(1,2) = -11513270273.0_wp/3502699200.0_wp ! r(3,6) = -33197340367.0_wp/1218433216.0_wp r(3,5) = -539868024987.0_wp/6092166080.0_wp r(3,4) = -39991188681.0_wp/374902528.0_wp r(3,3) = -69509738227.0_wp/1218433216.0_wp r(3,2) = -29327744613.0_wp/2436866432.0_wp ! r(4,6) = -284800997201.0_wp/19905339168.0_wp r(4,5) = -7896875450471.0_wp/165877826400.0_wp r(4,4) = -333945812879.0_wp/5671036800.0_wp r(4,3) = -16209923456237.0_wp/497633479200.0_wp r(4,2) = -2382590741699.0_wp/331755652800.0_wp ! r(5,6) = -540919.0_wp/741312.0_wp r(5,5) = -103626067.0_wp/43243200.0_wp r(5,4) = -633779.0_wp/211200.0_wp r(5,3) = -32406787.0_wp/18532800.0_wp r(5,2) = -36591193.0_wp/86486400.0_wp ! r(6,6) = 7157998304.0_wp/374350977.0_wp r(6,5) = 30405842464.0_wp/623918295.0_wp r(6,4) = 183022264.0_wp/5332635.0_wp r(6,3) = -3357024032.0_wp/1871754885.0_wp r(6,2) = -611586736.0_wp/89131185.0_wp ! r(7,6) = -138073.0_wp/9408.0_wp r(7,5) = -719433.0_wp/15680.0_wp r(7,4) = -1620541.0_wp/31360.0_wp r(7,3) = -385151.0_wp/15680.0_wp r(7,2) = -65403.0_wp/15680.0_wp ! r(8,6) = 1245.0_wp/64.0_wp r(8,5) = 3991.0_wp/64.0_wp r(8,4) = 4715.0_wp/64.0_wp r(8,3) = 2501.0_wp/64.0_wp r(8,2) = 149.0_wp/16.0_wp r(8,1) = 1.0_wp ! r(9,6) = 55.0_wp/3.0_wp r(9,5) = 71.0_wp r(9,4) = 103.0_wp r(9,3) = 199.0_wp/3.0_wp r(9,2) = 16.0d0 ! r(10,6) = -1774004627.0_wp/75810735.0_wp r(10,5) = -1774004627.0_wp/25270245.0_wp r(10,4) = -26477681.0_wp/359975.0_wp r(10,3) = -11411880511.0_wp/379053675.0_wp r(10,2) = -423642896.0_wp/126351225.0_wp ! r(11,6) = 35.0_wp r(11,5) = 105.0_wp r(11,4) = 117.0_wp r(11,3) = 59.0_wp r(11,2) = 12.0_wp ! case (3) ! ! METHD = 3 ! This pair is from "High Order Embedded Runge-Kutta Formulae" by P.J. ! Prince and J.R. Dormand, J. Comp. Appl. Math.,7, pp. 67-75, 1981. The ! authors are grateful to P. Prince and J. Dormand for their assistance in ! implementing the pair. ! no_of_stages = 13; fsal = .false.; order = 7 tan_angle = 11.0_wp; stability_radius = 5.2_wp safety = 0.8_wp; intrp_able = .false. intrp_needs_stages = .false.; intrp_degree = 0 no_of_ge_steps = 2 ! ptr(1:13) = (/ 0, 1, 2, 1, 3, 2, 4, 5, 6, 7, 8, 9, 1 /) ! a(2,1) = 5.55555555555555555555555555556e-2_wp a(3,1) = 2.08333333333333333333333333333e-2_wp a(3,2) = 6.25e-2_wp a(4,1) = 3.125e-2_wp a(4,2) = 0.0_wp a(4,3) = 9.375e-2_wp a(5,1) = 3.125e-1_wp a(5,2) = 0.0_wp a(5,3) = -1.171875_wp a(5,4) = 1.171875_wp a(6,1) = 3.75e-2_wp a(6,2) = 0.0_wp a(6,3) = 0.0_wp a(6,4) = 1.875e-1_wp a(6,5) = 1.5e-1_wp a(7,1) = 4.79101371111111111111111111111e-2_wp a(7,2) = 0.0_wp a(7,3) = 0.0_wp a(7,4) = 1.12248712777777777777777777778e-1_wp a(7,5) = -2.55056737777777777777777777778e-2_wp a(7,6) = 1.28468238888888888888888888889e-2_wp a(8,1) = 1.6917989787292281181431107136e-2_wp a(8,2) = 0.0_wp a(8,3) = 0.0_wp a(8,4) = 3.87848278486043169526545744159e-1_wp a(8,5) = 3.59773698515003278967008896348e-2_wp a(8,6) = 1.96970214215666060156715256072e-1_wp a(8,7) = -1.72713852340501838761392997002e-1_wp a(9,1) = 6.90957533591923006485645489846e-2_wp a(9,2) = 0.0_wp a(9,3) = 0.0_wp a(9,4) = -6.34247976728854151882807874972e-1_wp a(9,5) = -1.61197575224604080366876923982e-1_wp a(9,6) = 1.38650309458825255419866950133e-1_wp a(9,7) = 9.4092861403575626972423968413e-1_wp a(9,8) = 2.11636326481943981855372117132e-1_wp a(10,1) = 1.83556996839045385489806023537e-1_wp a(10,2) = 0.0_wp a(10,3) = 0.0_wp a(10,4) = -2.46876808431559245274431575997_wp a(10,5) = -2.91286887816300456388002572804e-1_wp a(10,6) = -2.6473020233117375688439799466e-2_wp a(10,7) = 2.84783876419280044916451825422_wp a(10,8) = 2.81387331469849792539403641827e-1_wp a(10,9) = 1.23744899863314657627030212664e-1_wp a(11,1) = -1.21542481739588805916051052503_wp a(11,2) = 0.0_wp a(11,3) = 0.0_wp a(11,4) = 1.66726086659457724322804132886e1_wp a(11,5) = 9.15741828416817960595718650451e-1_wp a(11,6) = -6.05660580435747094755450554309_wp a(11,7) = -1.60035735941561781118417064101e1_wp a(11,8) = 1.4849303086297662557545391898e1_wp a(11,9) = -1.33715757352898493182930413962e1_wp a(11,10) = 5.13418264817963793317325361166_wp a(12,1) = 2.58860916438264283815730932232e-1_wp a(12,2) = 0.0_wp a(12,3) = 0.0_wp a(12,4) = -4.77448578548920511231011750971_wp a(12,5) = -4.3509301377703250944070041181e-1_wp a(12,6) = -3.04948333207224150956051286631_wp a(12,7) = 5.57792003993609911742367663447_wp a(12,8) = 6.15583158986104009733868912669_wp a(12,9) = -5.06210458673693837007740643391_wp a(12,10) = 2.19392617318067906127491429047_wp a(12,11) = 1.34627998659334941535726237887e-1_wp a(13,1) = 8.22427599626507477963168204773e-1_wp a(13,2) = 0.0_wp a(13,3) = 0.0_wp a(13,4) = -1.16586732572776642839765530355e1_wp a(13,5) = -7.57622116690936195881116154088e-1_wp a(13,6) = 7.13973588159581527978269282765e-1_wp a(13,7) = 1.20757749868900567395661704486e1_wp a(13,8) = -2.12765911392040265639082085897_wp a(13,9) = 1.99016620704895541832807169835_wp a(13,10) = -2.34286471544040292660294691857e-1_wp a(13,11) = 1.7589857770794226507310510589e-1_wp a(13,12) = 0.0_wp ! ! The coefficients BHAT refer to the formula used to advance the ! integration, here the one of order 8. The coefficients B refer ! to the other formula, here the one of order 7. ! bhat(1) = 4.17474911415302462220859284685e-2_wp bhat(2) = 0.0_wp bhat(3) = 0.0_wp bhat(4) = 0.0_wp bhat(5) = 0.0_wp bhat(6) = -5.54523286112393089615218946547e-2_wp bhat(7) = 2.39312807201180097046747354249e-1_wp bhat(8) = 7.0351066940344302305804641089e-1_wp bhat(9) = -7.59759613814460929884487677085e-1_wp bhat(10) = 6.60563030922286341461378594838e-1_wp bhat(11) = 1.58187482510123335529614838601e-1_wp bhat(12) = -2.38109538752862804471863555306e-1_wp bhat(13) = 2.5e-1_wp ! b(1) = 2.9553213676353496981964883112e-2_wp b(2) = 0.0_wp b(3) = 0.0_wp b(4) = 0.0_wp b(5) = 0.0_wp b(6) = -8.28606276487797039766805612689e-1_wp b(7) = 3.11240900051118327929913751627e-1_wp b(8) = 2.46734519059988698196468570407_wp b(9) = -2.54694165184190873912738007542_wp b(10) = 1.44354858367677524030187495069_wp b(11) = 7.94155958811272872713019541622e-2_wp b(12) = 4.44444444444444444444444444445e-2_wp b(13) = 0.0_wp ! c(1) = 0.0_wp c(2) = 5.55555555555555555555555555556e-2_wp c(3) = 8.33333333333333333333333333334e-2_wp c(4) = 1.25e-1_wp c(5) = 3.125e-1_wp c(6) = 3.75e-1_wp c(7) = 1.475e-1_wp c(8) = 4.65e-1_wp c(9) = 5.64865451382259575398358501426e-1_wp c(10) = 6.5e-1_wp c(11) = 9.24656277640504446745013574318e-1_wp c(12) = 1.0_wp c(13) = c(12) ! end select ! ! The definitions of all pairs come here for the calculation of ! LAST_STAGE - the position of the last evaluated stage in a method ! RS1, RS2, RS3, RS4 - minimum and maximum rations used is step selection ! COST - the cost of a step ! MAX_STIFF_ITERS - the number of iterations permitted in stiffness detection ! There are at most Q = 3 function calls per iteration. MAX_STIFF_ITERS ! is determined so that Q*MAX_STIFF_ITERS <= 5% of the cost of ! 50 steps and 1 <= MAX_STIFF_ITERS <= 8. This limits the number of ! function calls in each diagnosis of stiffness to 24. ! EXPON - an exponent for use in step slection ! CDIFF - a coefficent used in determining the minimum permissible step ! last_stage = ptr(no_of_stages) if (fsal) then cost = real(no_of_stages-1,kind=wp) else cost = real(no_of_stages,kind=wp) end if ! max_stiff_iters = min(8,max(1,int(fivepc*cost*fifty))) ! expon = one/(order+one) ! ! In calculating CDIFF it is assumed that there will be a non-zero ! difference |C(I) - C(J)| less than one. If C(I) = C(J) for any I not ! equal to J, they should be made precisely equal by assignment. ! cdiff = one do i = 1, no_of_stages - 1 cdiff = min( cdiff, minval( & abs((c(i)-c(i+1:no_of_stages))), & mask=(c(i)-c(i+1:no_of_stages)/=0)) ) end do ! rs = two; rs1 = one/rs; rs2 = rs**2 rs3 = rs*rs2; rs4 = one/rs3 ! end subroutine method_const !endb! !startc! subroutine setup_r2(comm,t_start,y_start,t_end,tolerance,thresholds, & method,task,error_assess,h_start,message) ! ! Part of rksuite_90 v1.2 (December 1995) ! software for initial value problems in ODEs ! ! Authors: R.W. Brankin (NAG Ltd., Oxford, England) ! I. Gladwell (Math Dept., SMU, Dallas, TX, USA) ! see main doc for contact details ! real(kind=wp), intent(in) :: t_end, t_start !indep! real(kind=wp), intent(in) :: tolerance real(kind=wp), dimension(:,:), intent(in) :: y_start !dep! real(kind=wp), dimension(:,:), intent(in) :: thresholds !shp-dep! type(rk_comm_real_2d) :: comm real(kind=wp), intent(in), optional :: h_start !indep! logical, intent(in), optional :: error_assess, message character(len=*), intent(in), optional :: task, method ! character(len=*), parameter :: srname="SETUP" ! real(kind=wp) :: hmin !indep! real(kind=wp) :: cdiff integer :: ier, nrec, tr_dim_of_stages logical :: legalt character(len=1) :: task1, method1 ! integer, parameter :: not_ready=-1, fatal=911, just_fine=1 real(kind=wp), parameter :: zero=0.0_wp, pt01=0.01_wp, fivepc=0.05_wp, & third=1.0_wp/3.0_wp, one=1.0_wp, two=2.0_wp, ten=10.0_wp, fifty=50.0_wp ! ier = just_fine; nrec = 0 ! ! Clear previous state of the suite. ! call setup_global_stuff nullify(comm%thresh, comm%err_estimates, comm%weights, comm%y_old, & comm%scratch, & comm%y, comm%yp, comm%y_new, comm%yp_old, comm%stages, comm%ge_y, & comm%ge_yp, comm%ge_err_estimates, comm%ge_assess, comm%ge_y_new, & comm%ge_stages, comm%v0, comm%v1, comm%v2, comm%v3, comm%vtemp, & comm%xstage, comm%ytemp, comm%p) ! ! Fetch output channel and machine constants; ! call machine_const(comm%round_off,comm%sqrrmc,comm%cubrmc,comm%sqtiny, & comm%outch) ! body: do ! ! Check for valid shape if (size(shape(y_start))>0) then if (any(shape(y_start)==0)) then ier = fatal; nrec = 2; write (comm%rec,"(a)") & " ** An extent of Y_START has zero length. This is not permitted." exit body end if end if ! ! Check and process non-trivial optional arguments if (present(task)) then task1 = task(1:1); comm%use_range = task1 == "R" .or. task1 == "r" legalt = comm%use_range .or. task1 == "S" .or. task1 == "s" if (.not.legalt) then ier = fatal; nrec = 2; write (comm%rec,"(a,a,a/a)") & " ** You have set the first character of TASK to be '",TASK1,"'.", & " ** It must be one of 'R','r','S' or 's'." exit body end if end if if (present(method)) then method1 = method(1:1) select case (method1) case("L","l"); comm%rk_method = 1 case("M","m"); comm%rk_method = 2 case("H","h"); comm%rk_method = 3 case default ier = fatal; nrec = 2; write (comm%rec,"(a,a,a/a)") & " ** You have set the first character of METHOD to be '",METHOD1,"'.", & " ** It must be one of 'L','l','M','m','H' or 'h'." exit body end select end if if (present(message)) comm%print_message = message ! ! Check consistency of array arguments ! if (any(shape(y_start)/=shape(thresholds))) then ier = fatal; nrec = 1; write (comm%rec,"(a)") & " ** The shapes of Y_START and THRESHOLDS are not consistent." exit body end if ! ! Check and process compulsory arguments if (t_start == t_end) then ier = fatal; nrec = 1; write (comm%rec,"(a,e13.5,a)") & " ** You have set T_START = T_END = ",T_START,"." exit body else comm%t_end = t_end; comm%t_start = t_start comm%t_old = t_start; comm%t = t_start comm%dir = sign(one,t_end-t_start) end if if ((tolerance > pt01) .or. (tolerance < comm%round_off)) then ier = fatal; nrec = 2; write (comm%rec,"(a,e13.5,a/a,e13.5,a)") & " ** You have set TOLERANCE = ",tolerance," which is not permitted. The", & " ** range of permitted values is (",comm%round_off,",0.01)." exit body else comm%tol = tolerance end if if (minval(thresholds) < comm%sqtiny) then !spec-ar! ier = fatal; nrec = 2; write (comm%rec,"(a/a,e13.5,a)") & " ** You have set a component of THRESHOLDS to be less than the permitted", & " ** minimum,",comm%sqtiny,"." exit body end if ! ! Set formula definitions and characteristics call method_const(comm%rk_method, comm%a, comm%b, comm%c, comm%bhat, & comm%r, comm%e, comm%ptr, comm%no_of_stages, comm%intrp_degree, & comm%intrp_able, comm%intrp_needs_stages, comm%cost, & comm%safety, comm%expon, comm%stability_radius, comm%tan_angle, & comm%rs, comm%rs1, comm%rs2, comm%rs3, comm%rs4, comm%order, & comm%last_stage, comm%max_stiff_iters, comm%no_of_ge_steps, comm%fsal,& cdiff) ! tr_dim_of_stages = maxval(comm%ptr(2:comm%no_of_stages)) comm%toosml = comm%round_off/cdiff ! ! In STEP_INTEGRATE the first step taken will have magnitude H. If ! H_START = ABS(H_START) is not equal to zero, H = H_START. If H_START is ! equal to zero, the code is to find an on-scale initial step size H. To ! start this process, H is set here to an upper bound on the first step ! size that reflects the scale of the independent variable. ! RANGE_INTEGRATE has some additional information, namely the first output ! point, that is used to refine this bound in STEP_INTEGRATE when called ! from RANGE_INTEGRATE. If H_START is not zero, but it is either too big ! or too small, the input H_START is ignored and H_START is set to zero to ! activate the automatic determination of an on-scale initial step size. ! hmin = max(comm%sqtiny,comm%toosml*max(abs(t_start),abs(t_end))) if (abs(t_end-t_start) < hmin) then ier = fatal; nrec = 4; write (comm%rec,"(a/a/a,e13.5,a/a,e13.5,a)") & " ** You have set values for T_END and T_START that are not clearly", & " ** distinguishable for the method and the precision of the computer", & " ** being used. ABS(T_END-T_START) is ",ABS(T_END-T_START)," but should be", & " ** at least ",hmin,"." exit body end if if (present(h_start)) comm%h_start = abs(h_start) if (comm%h_start > abs(t_end-t_start) .or. comm%h_start < hmin) & comm%h_start = zero if (comm%h_start == zero) then comm%h = max(abs(t_end-t_start)/comm%rs3,hmin) else comm%h = comm%h_start end if ! ! Allocate a number of arrays using pointers. ! allocate(comm%thresh(size(y_start,1),size(y_start,2)), & !alloc! comm%err_estimates(size(y_start,1),size(y_start,2)), & !alloc! comm%weights(size(y_start,1),size(y_start,2)), & !alloc! comm%y_old(size(y_start,1),size(y_start,2)), & !alloc! comm%scratch(size(y_start,1),size(y_start,2)), & !alloc! comm%y(size(y_start,1),size(y_start,2)), & !alloc! comm%yp(size(y_start,1),size(y_start,2)), & !alloc! comm%stages(size(y_start,1),size(y_start,2),tr_dim_of_stages), & !alloc! comm%ymax(size(y_start,1),size(y_start,2)),stat=ier) !alloc! if (ier /= 0) then ier = fatal; nrec = 1 ; write (comm%rec,"(a)") & " ** Not enough storage available to create workspace required internally." exit body else comm%y = y_start; comm%ymax = abs(y_start) comm%thresh = thresholds; comm%y_new => comm%scratch; comm%yp_old => comm%scratch comm%v0 => comm%err_estimates; comm%vtemp => comm%scratch comm%v1 => comm%stages(:,:,1); comm%v2 => comm%stages(:,:,2) comm%v3 => comm%stages(:,:,3) end if ! ! Pre-allocate storage for interpolation if the TASK = `R' was specified. ! if (comm%use_range) then if (comm%intrp_able) then if (comm%rk_method==1) then comm%p => comm%stages(:,:,1:2) else if (comm%rk_method==2) then allocate(comm%p(size(y_start,1),size(y_start,2),5), & !alloc! comm%ytemp(size(y_start,1),size(y_start,2)), & !alloc! comm%xstage(size(y_start,1),size(y_start,2)),stat=ier) !alloc! if (ier /= 0) then ier = fatal; nrec = 1 ; write (comm%rec,"(a)") & " ** Not enough storage available to create workspace required internally." exit body end if end if end if end if ! ! Initialise state and allocate storage for global error assessment ! comm%t_ge_max_contrib = t_start if (present(error_assess)) comm%erason = error_assess if (comm%erason) then ! ! Storage is required for the stages of a secondary integration. The ! stages of the primary intergration can only be overwritten in the ! cases where there is no interpolant or the interpolant does not ! require information about the stages (e.g. METHOD 'H' and METHOD 'L', ! respectively). if (.not.comm%intrp_needs_stages) then comm%ge_stages => comm%stages else allocate(comm%ge_stages(size(y_start,1),size(y_start,2),tr_dim_of_stages),stat=ier) !alloc! if (ier /= 0) then ier = fatal; nrec = 1 ; write (comm%rec,"(a)") & " ** Not enough storage available to create workspace required internally." exit body end if end if allocate(comm%ge_y(size(y_start,1),size(y_start,2)), & !alloc! comm%ge_yp(size(y_start,1),size(y_start,2)), & !alloc! comm%ge_err_estimates(size(y_start,1),size(y_start,2)), & !alloc! comm%ge_assess(size(y_start,1),size(y_start,2)), & !alloc! comm%ge_y_new(size(y_start,1),size(y_start,2)),stat=ier) !alloc! if (ier /= 0) then ier = fatal; nrec = 1 ; write (comm%rec,"(a)") & " ** Not enough storage available to create workspace required internally." exit body else comm%ge_assess = 0.0_wp; comm%ge_y = y_start end if end if exit body end do body ! call rkmsg_r2(ier,srname,nrec,comm) ! contains subroutine setup_global_stuff ! comm%h_start=0.0_wp; comm%h_old=0.0_wp comm%f_count=0; comm%full_f_count=0; comm%step_count=0; comm%bad_step_count=0 comm%at_t_start=.true.; comm%at_t_end=.false. comm%rk_method=2; comm%ge_max_contrib=0.0_wp; comm%ge_f_count=0 comm%erason=.false.; comm%erasfl=.false. comm%print_message=.true.; comm%use_range=.true. comm%stiff_bad_step_count=0; comm%hit_t_end_count=0 comm%errold=0.0_wp; comm%h_average=0.0_wp comm%chkeff=.false.; comm%phase2=.true. comm%save_states(1:7)=not_ready comm%stop_on_fatal=.true.; comm%saved_fatal_err=.false. ! end subroutine setup_global_stuff end subroutine setup_r2 subroutine collect_garbage_r2(comm) ! ! Part of rksuite_90 v1.2 (December 1995) ! software for initial value problems in ODEs ! ! Authors: R.W. Brankin (NAG Ltd., Oxford, England) ! I. Gladwell (Math Dept., SMU, Dallas, TX, USA) ! see main doc for contact details ! type(rk_comm_real_2d) :: comm ! if (associated(comm%thresh)) then deallocate(comm%thresh); nullify(comm%thresh) end if if (associated(comm%y)) then deallocate(comm%y); nullify(comm%y) end if if (associated(comm%yp)) then deallocate(comm%yp); nullify(comm%yp) end if if (associated(comm%ymax)) then deallocate(comm%ymax); nullify(comm%ymax) end if if (associated(comm%scratch)) then deallocate(comm%scratch); nullify(comm%scratch) nullify(comm%y_new); nullify(comm%yp_old); nullify(comm%vtemp) end if if (associated(comm%weights)) then deallocate(comm%weights); nullify(comm%weights) end if if (associated(comm%ytemp)) then deallocate(comm%ytemp); nullify(comm%ytemp) end if if (associated(comm%y_old)) then deallocate(comm%y_old); nullify(comm%y_old) end if if (associated(comm%err_estimates)) then deallocate(comm%err_estimates); nullify(comm%err_estimates) nullify(comm%v0) end if if (associated(comm%p,comm%stages(:,:,1:2))) then nullify(comm%p) end if if (associated(comm%ge_stages,comm%stages)) then deallocate(comm%stages); nullify(comm%stages); nullify(comm%ge_stages); nullify(comm%v1,comm%v2,comm%v3) else if (associated(comm%ge_stages)) then deallocate(comm%ge_stages); nullify(comm%ge_stages) end if if (associated(comm%ge_y_new)) then deallocate(comm%ge_y_new); nullify(comm%ge_y_new) end if if (associated(comm%ge_assess)) then deallocate(comm%ge_assess); nullify(comm%ge_assess) end if if (associated(comm%ge_err_estimates)) then deallocate(comm%ge_err_estimates); nullify(comm%ge_err_e