Algorithm 771: rksuite_90: Fortran 90 software for ordinary differential equation initial-value problems.
Categories and Subject Descriptors: D.3.2 [Progamming Languages]: Language Classifications--Fortran 90; G.1.7 [Numerical Analysis]: Ordinary Differential Equations--initial value problems; GA [Mathematics of Computing]: Mathematical Software
General Terms: Algorithms
Additional Key Words and Phrases: Complex, recursion
MODULE rksuite_90 contains Fortran 90 software for the initial-value problem (IVP) in ordinary differential equations (ODEs)
y' = f(t, y), y(a) = [y.sub.a]
where a,t,b [element of] R, a [is less than] t [is less than] b (or b [is less than] t [is less than ] a), and y, f, [y.sub.a] [element of] R, [R.sup.n], [R.sup.m x n], G, [G.sup.n], or [G.sup.m x n]. In rksuite_90 we have deliberately adopted the integration algorithms of the Fortran 77 code RKSUITE though modified the code in many ways to exploit the efficiencies and facilities of Fortran 90; see Brankin et al.  for a fuller account of the design issues in Fortran 77 and Brankin et al.  for the Fortran 77 codes and documentation. In Brankin and Gladwell  we outlined our initial design of the Fortran 90 IVP software, concentrating on the language features exploited internally, e.g., long names, array syntax used internally, OPTIONAL arguments, the use of MODULES, INTERFACE blocks, and INTENT specifiers, etc. Here we discuss the software, its use, and its transformation in order to address different problem types and the algorithmic changes thus required. It is from features such as the OPTIONAL arguments and the transformations that the user will gain most benefits in comparison with the Fortran 77 software:
--an easier-to-use interface, once Fortran 90 is familiar;
--a somewhat more efficient code for large problems when Fortran 90 compilers produce code similarly efficient to that produced by Fortran 77 compilers;
--a more robust code, far more likely to trap subtle user-coding errors at compile time;
--a code able to deal directly with a wider range of problems types; and
--a code which the moderately sophisticated Unix user can transform further by modifying the scripts supplied.
As in RKSUITE, we provide two procedures to integrate across a range and compute approximations to y. One procedure is used to set up the integration with initial values, range, error requirements, and various options. The other integrates across the range and computes approximations to y at user-specified points t; for more complicated tasks we provide another procedure which integrates one step at a time. This is initialized using the same setup routine. Associated with the step integrator are procedures for interpolation and resetting the integration range. For use with either integrator are utility procedures which provide diagnostic information and memory deallocation.
In Section 2 we describe the interfaces provided in rksuite_90. Our emphasis is on the novel features in ODE software design permitted by the use of Fortran 90 (in contrast to Fortran 77). Section 3 discusses three design features. The first is a novel set of transformational tools to enable various formulations of the IVP to be solved directly without prior user transformation. The second concerns the method of managing global variables in a way which is "multiprocessor safe"; and this permits the third: the use of recursion. We provide an example using recursion for the invariant-embedding solution of boundary value problems. Finally, we conclude with a brief section discussing the efforts made to provide tags to identify independent-variable types for further transformation. An earlier, but almost identical, version of the code and examples has been available publicly for some time; we indicate where and how it was tested.
We use uppercase typewriter face letters for the names of fixed reserved Fortran 90 words, lowercase typewriter face letters for dummy and actual variables, and italicized lowercase letters for names of varying words a4d mathematical variables. We hope by these means to provide markers for those who are unfamiliar with Fortran 90 syntax.
2. THE FORTRAN 90 INTERFACES
We need a few preliminaries. Throughout, when computing with REAL machine numbers we work in a consistent precision, or KIND in Fortran 90 terminology. This KIND is set as an integer parameter wp used throughout MODULE rksuite_90 and set once in its own MODULE rksuite_90_prec. For procedures to communicate we define a special derived TYPE which we give the generic name rk_comm. Its contents are PRIVATE so as to prevent accidental modification. They include scalars, arrays whose extents depend on the structure of the dependent variables, character variables, and logicals.
We employ a number of generic terms:
(1) type(independent variable)
(2) type(dependent variable)
where type and dimension have different, but related, meanings to the Fortran 90 terms TYPE and DIMENSION. In the instances created by our Unix scripts, these terms can be interpreted as follows:
(1) type(independent variable) is equivalent to REAL of KIND wp.
(2) type(dependent variable) is equivalent to REAL or COMPLEX of KIND WP.
(3) dimension is absent for a scalar; DIMENSION(:) for a vector; and DIMENSION (: , :) for a matrix; this is used as an attribute for assumed-shape array arguments which are related to the dependent variable(s). When we specify the dimensionality of the result of a function (possibly array-valued) with respect to the dependent variable(s) y we write dimension (y), which is absent for a scalar; DIMENSION (SIZE (y)) for a vector; and DIMENSION (SIZE (y, 1), SIZE (y, 2)) for a matrix.
(4) for the types and dimensionalities of dependent variables currently available in rksuite_90, rk_comm is equivalent to those shown in Table I.
Table I. Actual Names for rk_comm 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
Otherwise the usual syntax of Fortran 90 applies.
Figure 1 is a simple template using the interfaces described below for the usual dependent-variable case, a real vector ODE system integrated over a real interval. The precision-defining MODULE rksuite_90_prec is USEd by MODULE rksuite_90. Hence, PROGRAM integrate_f derives its precision from USEing MODULE rksuite_90. When f is an array-valued FUNCTION, it must have an explicit INTERFACE, which enables checking at compile time. This INTERFACE can be provided by including the FUNCTION in a MODULE , as in Figure 1, or by including an INTERFACE block in the calling PROGRAM. To ensure use of the correct precision the MODULE define_f should also USE the MODULE rksuite_90_prec.
[Figure 1 ILLUSTRATION OMITTED]
2.1 The Procedures: setup and range_integrate
The procedure setup initializes the computation, so it is normally called only once; its specification is in Figure 2. A call to setup must precede the first call to range_integrate or step_integrate. Any subsequent call to setup reinitializes the computation. The argument comm is an instance of TYPE(rk_comm); t_start and t_end specify the range of integration [a, b]; and y_start specifies the initial conditions [y.sub.a]. tolerance AND thresholds define the error requirements. More precisely, the integration procedures ensure that the local error e in the computed solution y of size magnitude_y satisfies
ABS(e) / MAX(magnitude_y, thresholds) [is less than or equal to] tolerance
component wise. If the OPTIONAL arguments in setup are omitted, the default is to select a moderate-order integration method (equivalent to method= `m'), the integrator range_integrate (task= `r'), automatic computation of the initial step (h_start=0.0_wp), printing error messages (message=.true.), and no global error assessment (error_assess=.false.).
Fig. 2. Specification of setup. SUBROUTINE setup(comm,t_start,y_start,t_end,tolerance,thresholds, & method, task, error_assess, h_s tart,message TYPE(rk_comm), INTENT(OUT) :: comm type(independent variable), INTENT(IN) :: t_start, t_end type(dependent variable), dimension, INTENT(IN) :: y_start REAL(KIND=wp), INTENT(IN) :: tolerance REAL(KIND=wp), dimension(y), INTENT(IN) :: thresholds CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: method, task LOGICAL, OPTIONAL, INTENT(IN) :: error_assess, message type(independent variable), OPTIONAL, INTENT(IN) :: h_start
Procedure range_integrate computes approximations to the solution at user-specified points; its specification is in Figure 3. The array function argument f defines the ODEs as in f, t_want specifies where the solution is required, and y_got and yderiv_got return the solution and derivative at t_got. In a successful integration t_got = t_want. Otherwise t_got is the value of the independent variable where the integration halted; OPTIONALly the error message and the value of f lag indicate the difficulty. The program terminates on any error unless flag is provided when it returns with flag set to an error-reporting value. Problems that may be reported include the following: too small a step size was required to satisfy the error requirement; the global error assessment (where activated) ceased to be reliable; too many f evaluations have been used; or stiffness was detected.
Fig. 3. Specification of range_integrate. RECURSIVE SUBROUTINE range-integrate(comm,f,t_want,t_got,y_got, yderiv_got,flag) TYPE(rk_comm), INTENT(INOUT) :: comm INTERFACE FUNCTION f(t,y) type (independent variable), INTENT(IN) :: t type(dependent variable), dimension, INTENT(IN) :: y type(dependent variable), dimension(y) :: f END FUNCTION f END INTERFACE type(independent variable), INTENT(IN) :: t_want type(independent variable), INTENT(OUT) :: t_got type(dependent variable), dimension, OPTIONAL, INTENT(OUT) :: y_got, & yderiv_got INTEGER, OPTIONAL, INTENT(OUT):: flag
2.2 Utility Procedures
The utility procedures can be used in conjunction with either integrator. The procedure global_error reports the global error assessment if this option is selected in setup (via error_assess); its specification is in Figure 4. The argument rms_error reports, for each component of the solution, an assessment of the global error so far in the integration. For a successful integration, the components of rms_error will usually be a moderate factor of tolerance. The argument max_error reports the maximum contribution to this assessment, and the argument t_max_error reports the first point where max_error was observed.
Fig. 4. Specification of global error. SUBROUTINE global_error(comm,rms_error,max_error,t_max_error) TYPE(rk_comm), INTENT(INOUT) :: comm REAL(KIND=wp), dimension(y), OPTIONAL, INTENT(OUT) :: rms_error REAL(KIND=wp), OPTIONAL, INTENT(OUT) :: max_error type(independent variable), OPTIONAL, INTENT(OUT) :: t_max_error
The procedure statistics reports on the cost and performance of the integration; its specification is in Figure 5. total_f_calls, step_cost, and num_succ_steps report respectively the number of f evaluations used so far (excluding any used for global error assessment), the number of f evaluations required to take a single step with the selected method, and the number of successful steps taken so far. waste reports the fraction of attempted steps which failed to meet the local error requirement. A "large" fraction indicates the problem may be "stiff" or that the solution may have discontinuities in a low-order derivative. h_next reports the step size the integrator plans to use for the next step. y_maxvals reports the component wise largest value of ABS (y) computed at any step in the integration so far. It may be used in exploratory computations to determine a suitable setting for thresholds.
Fig. 5. Specification of statistics. SUBROUTINE statistics(comm,total_f_calls,step_cost,waste, & num_succ_steps, h_next,y_maxvals) TYPE(rk_comm), INTENT(INOUT) :: comm INTEGER, OPTIONAL, INTENT(OUT) :: total_f_Calls, step_cost, num_succ_steps REAL(KIND=wp), OPTIONAL, INTENT(OUT):: waste type(independent variable), OPTIONAL, INTENT(OUT) :: h_next REAL(KIND=wp), dimension(y), OPTIONAL, INTENT(OUT) :: y_maxvals
For each of global_error and statistics, at least one of the OPTIONAL arguments must be present. That is, we require that the call has an action, as a call without an action probably constitutes a user error. Either subroutine may be called at any time in the integration between calls to the integrators. Procedure garbage_collect has a single argument of TYPE rk_Comm and is provided to enable efficient use of memory. It should be called after completion of each completed or failed integration and prior to any reuse of the same instance of an argument in a new integration.
2.3 step-integrate and Associated Procedures
The procedure step_integrate is designed for tasks requiring close monitoring of the integration; its interface is in Figure 6. To ease use, less common demands are handled by the auxiliary procedures interpolate and reset_t_end. step_integrate integrates the solution one step at a time from t_start toward t_end. Arguments comm, f, and flag are as described for range_integrate. Arguments y_now and yderiv_now report the solution and derivative at t_now. If the integration has not been successful t_now, y_now, and yderiv_now retain their values at the end of the previous step.
Fig. 6. Specification of step_integrate. RECURSIVE SUBROUTINE step_integrate(comm,f,t_now,y_now,yderiv_now,flag) TYPE(rk_comm), INTENT(INOUT) :: comm INTERFACE FUNCTION f(t,y) type (independent variable), INTENT(IN) :: t type(dependent variable), dimension, INTENT(IN) :: y type(dependent variable), dimension(y) :: f END FUNCTION f END INTERFACE type (independent variable), INTENT(OUT) now type(dependent variable), dimension, OPTIONAL, INTENT(OUT) y-now, & yderiv-now INTEGER, OPTIONAL, INTENT(OUT)) flag
Procedure reset_t_end is used to reset the end of the integration range, for stepping to output points using step-integrate; its specification is in Figure 7. The argument t_end_new redefines the end of the range. Procedure reset_t_end should be used in preference to reinitializing using setup. It cannot be used to redefine the direction of integration.
Fig. 7. Specification of reset-t-end. SUBROUTINE reset_t_end(comm,t_end_new) TYPE(rk_comm), INTENT(INOUT) :: comm type(independent variable), INTENT(IN) :: t_end_new
In general, the most efficient way to compute the solution at specified points is to use interpolate; its specification is in Figure 8. It uses information calculated over the step just taken and may require additional evaluations of f in that step. (The number of additional evaluations required on any one step is fixed for any number of output points.) Argument t_want specifies the value of the independent variable where a solution value is required. Approximations to the solution and derivative are returned in y_want and yderiv_want. At least one of y_want or yderiv-want must be present so that the call has an action.
Fig. 8. Specification of interpolate. RECURSIVE SUBROUTINE interpolate(comm,f,t_want,y_want,yderiv_want) TYPE(rk_comm), INTENT(INOUT) :: comm INTERFACE FUNCTION f(t,y) type(independent variable), INTENT(IN) :: t type(dependent variable), dimension, INTENT(IN) :: y type(dependent variable), dimension(y) :: f END FUNCTION END INTERFACE type(independent variable), INTENT(IN) :: t_want type(dependent variable), dimension, OPTIONAL, INTENT(OUT) y_want, & yderiv_want
The availability in rksuite_90 of an interpolant with an internal power series storage mechanism permits a relatively straightforward attempt at event location. In Shampine et al. , we described a robust technique for a limited event location capability, for problems commonly solved by inverse interpolation. It describes a Fortran 77 event location code that can be used "side-by-side" with any Fortran 77 step-by-step integrator with the appropriate type of power series interpolant in the equivalent of a reverse communication scheme. Equally, it could be used side-by-side with step_integrate, since Fortran 77 is a subset of Fortran 90.
3. INTERNAL DESIGN OF RKSUITE_90
The major programming language and algorithmic design issues for the initial release of rksuite_90 were discussed in Brankin and Gladwell [19941. This final version contains three sets of changes: the treatment of a variety of dependent variables via transformations, the handling of global variables, and permitting recursion.
3.1 Dependent Variables
Almost all IVP solvers treat the case t [element of] R, y [element of] [R.sup.n]. As a result a user whose problem naturally has a solution y [element of] R or [R.sup.m x n], or a complex equivalent, is forced to recast the problem in [R.sup.n], a possibly error-prone task with a potentially inefficient result. We have treated these cases directly, aiming to maintain just one copy of the base software. We restrict the independent variable so that t [element of] R, but we have isolated the specifications of local variables of type (independent variable) to ease production of future codes with t in another field. We consider three aspects of this problem: the type of, the shape of, and the arithmetic operations involved with the dependent variable. We start from the software for y [element of] [R.sup.n]. In the code for y [element of] [R.sup.n] we have identified with a tag (a special in-line comment) all declarative sections of the code involving dependent variables. (We have used IMPLICIT NONE so that every entity has been declared.) This tag automatically identifies all dependent-variable quantities with an associated rank. However there are other quantities in the code which must have the same shape as the dependent variables. A further special tag has been inserted into the declarative sections to identify these. This tagging partially addresses the questions of type and shape of the variables, as mentioned above. For the arithmetic aspect, we have examined all assignments involving dependent variables. Throughout, we have ensured that array assignments
a = ... b + . . c
are used. That is, extents are not specified, since the arrays are always conformable. We find that all divisions involve real conformable divisors; multiplications of dependent-variable quantities involve a real scalar multiplier; the intrinsic ABS is used to return real quantities (magnitudes); the intrinsic mAx is used on real conformable quantities; and the intrinsics MINVAL and MAXVAL are used on real quantities. All operations involving y [element of] [R.sup.n] require no transformation to address y [element of] [C.sup.n]. Similarly, almost all assignments involving arrays require no transformation when considering y [element of] R and y [element of] [R.sup.m X n]. The exceptions are expressions using the intrinsics MINVAL, MAXVAL, and sum, which must take array arguments and are therefore inappropriate for y [element of] R. Instances of these intrinsics have been specially tagged so that they can be transformed.
We have created sed scripts exploiting the tags for use in a Unix environment to
(1) create a complex version from a real version,
(2) create a 2d version from a 1d version, and
(3) create a 0d version from a 1d version.
Two minor difficulties are addressed in the scripts. The first involves using the ALLOCATE statement to assign internal workspace via pointers. The exact extents of the arrays to be allocated must be supplied, and they depend on the shape of y. We have inserted special tags so that the "1d to 2d" script can add the second dimension of appropriate size and so that the "1d to 0d" script can remove the dimension information and allocate the equivalent scalar pointer. A residual problem in the scalar case is that some of the elemental array intrinsics used do not permit scalar arguments, requiring a workaround in the sed script. The second difficulty arises in the internal procedures implementing a stiffness check. A nonlinear power method is used to estimate eigenvalues of the Jacobian of f. It requires computation of inner products. For y,z [element of] [R.sup.n] this inner product is simply [y.sup.T]z. However for y,z [element of] [C.sup.n], the inner product is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The quantity
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]
is also required. Unfortunately the Fortran 90 elemental intrinsic CONJG requires an argument of type COMPLEX. In a complex arithmetic version of the code there are two places where CONJG is required. As before we have inserted tags so that the "real to complex" script can make the appropriate edit. Note that for y,z [element of] [R.sup.m x n] or [C.sup.m x n] the inner product remains well defined in Fortran 90, since the matrix-matrix multiplication is performed elementally. To derive a generic code, we restructured the stiffness check to use complex arithmetic even in the real arithmetic case (e.g., in a quadratic equation solver). After these changes, the results for y [element of] [R.sup.n] are identical to those computed using the codes described in Brankin et al. [1992; 1993]. When solving problems with y [element of] R or y [element of] C, the usual definition of stiffness (for systems) does not apply. However, it is still possible for a scalar ODE to be stiff. Consider
y'(t) = k (y(t) - p (t)) + p'(t), y (0) = A
with solution y(t) = (A - p(0))[e.sup.kt] + p(t); see Shampine [1994, p. 3841. If p (t) is smooth and slowly varying, and k is large and negative, this problem is stiff. We have solved this problem using rksuite_90 for p (t) = cos(t), A = 0 and k = - [10.sup.n] , n = 1,2,3, . . . . Stiffness is flagged for all n [is greater than] 2 on long-enough ranges of integration. So the default mechanisms in the stiffness check work here. Similarly, in RKSUITE and rksuite_90 these mechanisms work both for dependent-variable arrays of length 1 and for scalars.
The alternate versions of rksuite 90 are produced using a further Unix script provided with the Fortran 90 software. The user has the option of generating any individual version or all six. The versions are bundled in a generic interface which hides the number of interfaces. That is, the user sees just one each of setup, range_integrate, etc. The Fortran 90 compiler determines the specific version by the differences in types of the arguments in each generic interface. Particularly, the type of rk_comm is used to determine which instance to use if there are no other differences.
3.2 Global Variables
In RKSUITE global variables were passed in COMMON or, for those of length depending on n, in an unspecified work array passed as an argument between procedures. In the first version of rksuite_90 the same information was kept in PRIVATE (global variables) in a MODULE. This is natural in Fortran 90 but is not "multiprocessor safe" (nor is COMMON) and hence also unsuitable for multithreading environments, as separate copies of the global variables cannot be made. So, we have reverted to the traditional "packed work array" approach using the derived TYPE rk_comm for passing global information between procedures. A user (or system) may declare as many instances of TYPE (rk_comm) as required. A side effect of the change in handling of global variables is to permit simultaneous (side-by-side) integrations to be performed in the same program.
The approach adopted for global variables also enables recursive calls to the integrators. Hence a differential equation can be defined in terms of the solution of other differential equations. For example, recursion permits the direct evaluation of multidimensional quadrature using just the one-dimensional integration procedure defined by rksuite_90.
It also permits solving invariant-embedding problems in a natural (though inefficient) way. Consider the boundary value problem
y"(t) - y(t) = 0,y(0) = 1, y'(T) = [- e.sup.-T]
with general solution y(t) = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The boundary conditions force [C.sub.1] = 0, [C.sub.2] = 1. Forward or backward shooting will be unsuccessful for large T due to the exponential solution behavior, but invariant embedding may be used. As shown in Meyer [1973, p. 68], an invariant-embedding method consists of
(1) integrating together the ODEs
u'(t) = 1 - [u(t).sup.2], v'(t) = - u(t)v(t), t [element of] [0, T], u(O) = 0, v(0) = 1,
(2) finding x by solving
u(T)x = [e.sup.-T] - v(T),
(3) integrating ODE
x'(t) = u(t)x(t) + v(t), x(T) = x, t [element of] [0, T]
"backward" from t = T to t = 0.
The solution of the boundary value problem is y(t) = u(t)x(t) + v(t). To solve (3), we need the solutions of (1) everywhere. Probably the most efficient approach would be to save the interpolants for u(t) and v(t) everywhere on [0, T] and evaluate as necessary when solving (3). This may be achieved by saving copies of rk_comm on every step when using step_integrate. Then, the saved copies are used in reverse order in calls to interpolate when integrating (3). A more elegant, and much more expensive, approach is to use recursion calling step_integrate to solve (3) and calling range_integrate for (1) from the function defining the ODE for (3); a template is given in Figure 9. (Note that range_integrate calls step_integrate, so this is a recursive call.)
[Figure 9 ILLUSTRATION OMITTED]
4. CONCLUSIONS AND AVAILABILITY
We have developed software to solve initial-value problems when the dependent variables y [element of] R, [R.sup.n], [R.sup.m x n], C, [C.sup.n], or [C.sup.m x n] and the independent variable t [element of] R. It would be a simple matter to address y of dimension greater than two. We have identified and tagged those parts of the declarative sections which involve quantities of the type-dependent variable and type-independent variable. This should facilitate the creation of initial-value solvers with variables in fields different to those already considered. For example, it is relatively straightforward to address t [element of] C integrating along a complex straight line; indeed we have a version for this case produced by a further Unix script transformation, not included with rksuite_90. In contrast, modifying the software for integration along an are requires further transformational tools to incorporate the definition of the arc and to measure distance along it (corresponding to a Fortran 90 intrinsic in the simple straight line cases); hence it requires a variation of the interface. To address other fields, we must redefine various operations (such as addition, multiplication by a real scalar, a magnitude function, and an inner product). Implementation could be achieved by overloading the operators (using overloading procedures, with all the associated overhead). The interpretation of stiffness in these cases requires careful rethinking. Indeed, our complex straight line version has an automatically produced stiffness check, but it is not clear that it monitors what is usually considered to be stiffness.
The package comprises the base software (for y [element of] [R.sup.n]), the Unix scripts to generate alternate versions, the documentation, and several example programs, including some using the alternate versions. The first version of the software has been available on netlib (at email@example.com and its various mirror sites) in directory ode/rksuite and in the NAG Fortran 90 software repository at http://www.nag.co.uk. This software has been available and has been accessed frequently since mid-1994. The authors have received no error reports. The codes were developed using the NAGWare f90 compiler and tested on a variety of platforms using this compiler. They have also been tested using the latest versions available to us of the following: the DEC OSF/1 Fortran 90 compiler, the Cray CF90 compiler, and the IBM AIX XL Fortran Compiler/6000.
The authors thank Jeremy Du Croz and the anonymous referees for their suggestions on improving the presentation of this article.
BRANKIN, R. W. AND GLADWELL, I. 1994. A Fortran 90 version of RKSUITE: An ODE initial value solver. Ann. Num. Math. 1, 363-375.
BRANKIN, R. W., GLADWELL, I., AND SHAMPINE, L. F. 1992. RKSUITE: A suite of Runge-Kutta codes for the initial value problem for ODEs. Softreport 92-S1, Mathematics Dept., Southern Methodist Univ., Dallas, TX.
BRANKIN, R. W., GLADWELL, I., AND SHAMPINE, L. F. 1993. RKSUITE: A suite of explicit Runge-Kutta codes. In Contributions to Numerical Mathematics, R. P. Agarwal, Ed. WSSIAA, vol. 2. World Scientific Publishing Co., Inc., River Edge, NJ., 41-53.
MEYER, G. H. 1973. Initial Value Methods for Boundary Value Problems: Theory and Application of Invariant Imbedding. Academic Press Mathematics in Science and Engineering Series, vol. 100. Academic Press, Inc., Orlando, FL.
SHAMPINE, L. F. 1994. Numerical Solution of Ordinary Differential Equations. Chapman & Hall, Ltd., London, UK.
SHAMPINE, L. F., GLADWELL, I., AND BRANKIN, R. W. 1991. Reliable solution of special event location problems for ODEs. ACM Trans. Math. Softw. 17, 1 (Mar.), 11-25.
Received: July 1995; revised: July 1996; accepted: January 1997 STRIPACK is a
The work of R. W. Brankin was performed while at The Numerical Algorithms Group Ltd. Authors' addresses: R. W. Brankin, Systems and Software Engineering Centre, DERA Malvern, Saint Andrews Road, Malvern, WR14 3PS, England; e-mail: firstname.lastname@example.org; I. Gladwell, Department of Mathematics, Southern Methodist University, Dallas, TX 75275; email: email@example.com. Permission to make digital/hard copy of part or all of this work for personal or classroom use is granted without fee provided that the copies are not made or distributed for profit or commercial advantage, the copyright notice, the title of the publication, and its date appear, and notice is given that copying is by permission of the ACM, Inc. To copy otherwise, to republish, to post on servers, or to redistribute to lists, requires prior specific permission and/or a fee.
|Printer friendly Cite/link Email Feedback|
|Author:||Brankin, R.W.; Gladwell, I.|
|Publication:||ACM Transactions on Mathematical Software|
|Date:||Sep 1, 1997|
|Previous Article:||Level 3 Basic Linear Algebra Subprograms for sparse matrices: a user-level interface.|
|Next Article:||Algorithm 772: STRIPACK: Delaunay triangulation and Voronoi diagram on the surface of a sphere.|