c/EPA
United States
Environmental Protection
Agency
Environmental Research
Laboratory
Athens GA 30605
EPA-600/3-80-037
March 1980
Research and Development
Efficient
Algorithms for
Solving Systems of
Ordinary Differential
Equations for
Ecosystem Modeling
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development, U.S. Environmental
Protection Agency, have been grouped into nine series. These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
6. Scientific and Technical Assessment Reports (STAR)
7. Interagency Energy-Environment Research and Development
8. "Special" Reports
9. Miscellaneous Reports
This report has been assigned to the ECOLOGICAL RESEARCH series. This series
describes research on the effects of pollution on humans, plant and animal spe-
cies, and materials. Problems are assessed for their long- and short-term influ-
ences. Investigations include formation, transport, and pathway studies to deter-
mine the fate of pollutants and their effects. This work provides the technical basis
for setting standards to minimize undesirable changes in living organisms in the
aquatic, terrestrial, and atmospheric environments.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.
-------
EPA-600/3-80-037
March 1980
EFFICIENT ALGORITHMS FOR SOLVING SYSTEMS
OF ORDINARY DIFFERENTIAL EQUATIONS
FOR ECOSYSTEM MODELING
by
John Malanchuk, John Otis and Hubert Bouver*
Institute for Man and Environment
State University of New York
Plattsburgh, New York 12901
*Applied Physics Laboratory
Johns Hopkins University
Laurel, Maryland 20810
Grant No. R805452
Project Officer
Ray R. Lassiter
Environmental Systems Branch
Environmental Research Laboratory
Athens, Georgia 30605
ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ATHENS, GEORGIA 30605
-------
DISCLAIMER
This report has been reviewed by the Athens Environmental Research
Laboratory, U.S. Environmental Protection Agency, and approved for publi-
cation. Approval does not signify that the contents necessarily reflect the
views and policies of the U.S. Environmental Protection Agency, nor does
mention of trade names or commercial products constitute endorsement or
recommendation for use.
ii
-------
FOREWORD
Environmental protection efforts are increasingly directed towards
preventing adverse health and ecological effects associated with specific
compounds of natural or human origin. As part of this Laboratory's research
on the occurrence, movement, transformation, impact, and control of environ-
mental contaminants, the Environmental Systems Branch studies complexes of
environmental processes that control the transport, transformation, degrada-
tion, and impact of pollutants or other materials in soil and water, assesses
environmental factors that affect water quality, and develops mathematical
models for evaluating the environmental impact of pollutants.
Numerical integration packages for mathematical modeling, although
widely available, are usually documented in the language of the numerical
analyst. The unsophisticated user often has difficulty with the terminology
and frequently chooses a routine that is far from optimal. This report pro-
vides useful information and comprehensible guidance to one particular group
of users—mathematical modelers with backgrounds in subject areas such as
biology. The objective is to provide easily accessible routines capable of
efficiently addressing the class of problems often encountered by these
modelers. The documentation is extensive and rigorous yet couched in expli-
cit terms that are easily understood. Although nothing new appears here from
the point of view of a numerical analyst, it is hoped that modelers will find
the attempt at lucidity to have been worthwhile to them.
David W. Duttweiler
Director
Environmental Research Laboratory
Athens, Georgia
iii
-------
ABSTRACT
This report presents three packages of subroutines for
the numerical solution of systems of first order ordinary
differential equations. The three integration methods were
chosen to provide an efficient means of solving the full
range of initial value problems encountered in basin ecosys-
tem modeling. The subroutines have been designed to handle
all aspects of the numerical integration as automatically as
possible, so that their successful use does not require a
sophisticated knowledge of the techniques of numerical anal-
ysis.
After a general introduction to the principles and con-
cepts of numerical integration, the particular methods used
in the codes are discussed in some detail. Full instructions
for use are given, and some simple examples presented. The
complete code listings, in ANSI FORTRAN IV, are contained in
Appendices.
This report was submitted in fulfillment of Grant No.
R-805452 by Research Foundation of State University of New
York under the sponsorship of the U.S. Environmental
Protection Agency. This report covers the period July 1,
1977 to August 31, 1978, and work was completed as of
November 30, 1978.
-------
CONTENTS
Foreword ........................... ill
Abstract . iv
1. Introduction ...................... 1
2. Basic Ideas 4
3. Selection of Numerical Integration Programs 13
4. Runge-Kutta Methods ....... 18
5. Adams Methods . . 23
6. Stiff Methods 35
7. Using the Codes . . 44
8. Examples ....56
References .................62
Appendices
A. Runge-Kutta Programs ... ...... 64
B. Adams Method Programs ............ 84
C. Stiff Method Programs 112
D. Solving Partial Differential Equations . . 145
-------
SECTION 1
INTRODUCTION
The modeling of basin ecosystems is an immense undertaking which is most
efficiently conducted by a team of professionals who bring together a variety
of technical expertise. Mutual collaboration on such projects is aimed at
producing a unified representation of an entire system. Too often, however,
the greatest difficulties in such a "systems approach" are encountered not
within any single discipline, but at the interfaces between disciplines.
One such interface, where a large gap exists, is between ecologists and
mathematicians. The majority of ecologists have been unsuccessful in
articulating their problems in a form understandable to mathematicians.
Mathematicians, on the other hand, while equipped with the tools to solve
many ecological problems, have little or no familiarity with the systems
under study, and are equally unable to bridge the gap between problem
articulation and solution.
This report and the accompanying software are an attempt to ease this
problem for the ecologist by presenting methods for the numerical solution of
ordinary differential equations, methods which can be used with ease and
confidence without need for a sophisticated understanding of the techniques
of numerical analysis.
Other bridges need to be built. This work does not address the problem
of translating intimate knowledge 6f an ecosystem into a system of
differential equations. Fortunately, numerous examples of the necessary
techniques have appeared in the ecological literature as the importance of
these endeavors has become increasingly recognized. Patten(1971,1972,1975,
1976) has published several volumes which provide instruction in the methods
of ecological model construction, and which present examples from a variety
of ecosystems.
Figure 1-1. Simple Ecosystem Model.
The y. represent components of the
system and the F. . represent the
material fluxes between components.
PRODUCER
'21
I
INORGANIC
NUTRIENT
HERBIVORE
DECOMPOSER
-------
For illustrative purposes, a simple ecosystem model is presented.
Suppose that this ecosystem is adequately described by the existence of an
inorganic nutrient, a producer, an herbivore, and a decomposer. The
relationships among the various components are depicted in Figure 1-1.
Given the box and arrow diagram of Figure 1-1, the differential
equations describing the time dependence of the components may be written as
y\ = FHt ~ F21 US = F32 - F*3
(.1-1)
VZ = FZ1 - F32 ~ F*2. V* = FtZ * ^3 - F11
where the primes denote differentiation with respect to time.
The model is not yet complete, for the mechanism which controls the
material flows has not been specified. Suppose that each material flux is
regulated only by the donating component and, furthermore, is proportional to
that component. This assumption specifies the model described by
V2 = 211 32*22
(1-2)
9s =
where the a . . are positive constants
If this model is to describe a real system, the constants must be
determined. Typically, the ecologist will have observations of the values
of the components over some range of times, and will have some idea of the
approximate values of the constant parameters. His first task, then, is to
refine the values of these parameters in order to provide the best fit
between the model and the empirical data. The adequacy of the model will be
judged by the degree of agreement which can be obtained.
To carry out this fitting procedure, the ecologist must be able to solve
the equations (1—2) for a variety of choices of the parameters. The codes
presented in this report are intended to provide the tool needed for this
purpose. These routines have been designed to be easily implemented and
used, while automatically providing a sophisticated and efficient handling of
the numerical solution. They have been written in ANSI FORTRAN IV and,
except for the setting of a few machine dependent constants in DATA state-
ments, should require no modification to run on any system which provides
standard FORTRAN.
Since no single numerical integration method can solve all problems
efficiently, three packages of programs are presented: Runge-Kutta Programs
(listed in Appendix A), Adams Method Programs (Appendix B), and Stiff Method
Programs (Appendix C). Guidelines for choosing the method best suited to a
-------
particular problem are stated in Section 7, which constitutes a user's manual
for the codes. Furthermore, the three methods have been designed to be
easily interchangeable, so that a little experimentation to discover the best
choice should not be inconvenient.
For most users, it will suffice to read Section 2, which contains some
basic definitions, and Section 7, which explains the use of the codes.
Fairly complete instructions for use are also contained as comments within
the program listings in Appendices A, B, and C. Section 8 presents some
simple examples illustrating the performance of the codes.
Section 3 is preparatory for Sections 4, 5, and 6, which provide
detailed descriptions of, respectively, the Runge-Kutta Programs, the Adams
Method Programs, and the Stiff Method Programs. These sections should be
consulted when a more thorough understanding of the programs is necessary,
such as when modifications to the codes are being considered.
Appendix D contains a brief and elementary account of the manner in
which the programs provided here can be utilized for the solution of some
partial differential equation systems appearing in basin ecosystem models.
-------
SECTION 2
BASIC IDEAS
This section presents a number of definitions and concepts pertaining to
the numerical solution of ordinary differential equations. The aim is to
provide only that information which is truly necessary for understanding the
codes and for using them properly. Accordingly, a number of technical points
are ignored or glossed over. More rigorous and complete treatments of these
matters may be found in a number of textbooks, such as Gear (1971) and
Shampine( 1975) .
§2.1 The Initial Value Problem
The codes presented in this report are designed for the numerical
solution of the initial value problem for a system of first order ordinary
differential equations,
..... yn)
y2 =
(2-1)
yn= fn(t'yi'y2 yn>
W = ylO' y2(tO} = y20' ••" W = ynO
where the prime denotes differentiation with respect to the independent
variable t. The independent variable will often be referred to as "time",
which is its actual interpretation in many applications.
Such a problem is compactly described using a vector notation
y = (yj/y2/.--/yn)
in terms of which (2-1) may be written as
l(t,yj y(tQ) = }Q (2-2)
-------
In discussing numerical integration procedures, it occasionally will be
necessary to refer to the Jacobian matrix of f. This is the n x n matrix
jt(t,y) =
3y
3y
n
*2
3y
n
3y
n
n
3y.
(2-3)
where all the partial derivatives are evaluated at (t,y).
The numerical methods presented here apply only to the solution of first
order equations. To solve the initial value problem for a higher order
differential equation,
y(n)(t) = f(t,y,y'/y",...,y
-------
It should be understood that this procedure for reducing (2-4) to a
system of first order equations is not unique. For a particular problem,
physical or mathematical considerations may suggest a different choice of
auxiliary variables which is more meaningful and appropriate.
§2.2 Numerical Solution of the Initial Value Problem
For purposes of analysis and explication, it is generally sufficient to
discuss the single-equation problem
y' = f(t,y) y(t0) = yQ (2-5)
Because of the formal resemblance between (2-5) and (2-2), insights gained
from discussion of (2-5) can usually be taken over in obvious fashion to the
n-equation problem (2-1).
A numerical solution of (2-5) consists of finding approximations
yo» yi» ya» ••• to the solution y(t) on a sequence of discrete mesh points
t0, t}, t2, ... in an interval [t0/^]' -"-he spacing between adjacent mesh
points is called the stepsize and is denoted by h. The mesh points are
usually not equally spaced; when the stepsize referred to is not clear from
the context, it is distinguished by use of a subscript:
h. = t. - t.
i i i-l
Starting with the given initial point (t0,y0), the numerical solution proceeds
step by step, with the solution at each mesh point being calculated by some
numerical integration methodt which specifies the solution in terms of f(t,y)
and solution values already found at previous mesh points. Integration
methods may be classified by the number of previous values which they
require: an #-step method utilizes the solution values at the fc preceding
mesh points. A single-step method (&=1) is also called self -starting , since
it can make the first step of the solution using only the information given
in the problem statement (2-5). A multistep method (&>1) cannot be used to
begin the solution, but it may be used to continue after a sufficient number
of points have been calculated by a different method. The simplest single-
-step methods are
Euler method: y.^ = y^j + hif(ti_}'yi_l'> (2-6a)
backward Euler method: .y. = y. + fi..f(t.,y.) (2-6b)
Integration methods may be further classified as explicit or implicit. An
explicit method, such as (2-6a), expresses the solution directly in terms of
previously calculated values, while an implicit method, such as (2-6b),
specifies the solution implicitly by stating an algebraic equation which it
must satisfy.
-------
The problem (2-5) is well-posed over the interval Q:0,zQ if small
changes in the equation or in the initial condition produce only small
changes in the solution. This is a more stringent condition than the mere
existence of a unique solution.
Since any numerical computation must inevitably introduce small errors,
well-posedness is necessary for (2-5) to be amenable to numerical solution.
Otherwise, these small errors would completely change the character of the
computed solution. Although it is possible to artificially construct
problems which are not well-posed, such problems are unlikely to arise as
models of real systems, since such models must display continuous behavior
within the uncertainites of the empirical data upon which they are con-
structed. In any case, well-posedness need not be a serious a. priori
concern of the modeler, since its lack will become apparent from the anom-
alous behavior of the solution as the parameters of the model are varied.
§2.3 Errors in Numerical Solution
The error (or global error) in a point of the numerical solution is the
amount by which it differs from the correct value,
ej = y(ti) - y. (2-7)
where y(t) is the exact solution of (2-5).
Errors in the numerical solution arise from two sources: truncation
errors and roundoff errors . Roundoff errors are those which occur because
arithmetic is not performed exactly. Even if the arithmetic were exact,
truncation errors would occur because the integration method itself is only
an approximation. All the methods discussed here are convergent, meaning
that (assuming exact arithmetic) the error can be kept arbitrarily small by
choosing the stepsizes sufficiently small.
The local truncation error for a step of the integration is the error
which would be introduced at that step if the solution to the beginning of
the step were exact:
T. = y(t.) - y.
given: y^_l = y(ti_])/
(2-8)
The form of the local truncation error is characteristic of the integration
method being employed. A method is said to be of order k if the local
truncation error is of order h^+'. (This definition of the order of the
method assumes that y(t) possesses at least k continuous derivatives; if this
is not the case for a particular problem, then the effective order of the
method is reduces accordingly when applied to that problem.)
If the numerical solution is to be a useful approximation of y(t), the
global error must be kept within reasonable bounds throughout the interval of
-------
integration. In the codes, the error is controlled indirectly by choosing
each stepsize small enough to keep a quantity called the local error smaller
than limits supplied by the user. To understand the meaning of local error,
consider all the solution curves in the (t,y) plane, of the differential
equation in (2-5). For a well-posed problem, these curves do not intersect
in any region of interest. The initial condition in (2-5) selects as the
solution the unique curve y(t) which contains the point (t0,y0). Suppose
that the numerical solution has proceeded to t^.p Because y£_j contains
some error, it dees not lie on y(t), but rather on the nearby solution curve
u(t) specified by u(t£_j)=yi_j. In making the step to tj_, the numerical
method will then be following u(t) instead of y(t). The local error in the
step from t^_j to t^ is defined to be
n. = u(t.) - y^ (2-9)
and is a measure of the amount by which the step fails to follow the solution
curve on which it starts: ni may be thought of as the truncation error of the
method as applied to u(t) instead of to y(t).
Formally, the local errors may be made arbitrarily small by using
sufficiently small stepsizes. However, this cannot be done in practice,
since'too small a stepsize will produce roundoff errors which exceed the
local error. Thus, the accuracy achievable in the computed solution is
ultimately limited by the precision of the arithmetic.
Although limiting the local error as just described does place a bound
on the global error, the size of this bound depends on the characteristics of
the problem being solved. If 3r~/3y>0, then the curves u(t) and y(t) diverge
with increasing t, so that the global errors e^, ££+1, ... will increase even
if the local errors n^, rii+l» ••• vanish. A problem with this behavior is
termed unstable. If the nearby solution curves diverge rapidly from the
desired solution, the global error at the end of the integration interval
may be expected to be substantially larger than the limit placed on the local
error. Confining the global error within a given bound will then require
that the local errors be limited to much smaller values, in order to keep the
computed solution away from those curves which move too far away from y(t)
during the integration interval. The need to solve a severely unstable
problem, 3f/3y»0, is unlikely to arise in practice, since any uncertainty in
empirical initial conditions necessarily encompasses widely differing
solutions: in a practical sense, such a problem is not well-posed.
If, on the other hand, 3-f/3y<0, then the nearby solution curves converge
toward y(t). Under these circumstances, the effect of previous errors is
damped out as integration proceeds, and the global errors may be expected to
be comparable in magnitude to the limit placed on local errors. Such
problems are termed stable.
For a given method, the expense of integrating (2-5) is proportional to
the number of steps required to span the integration interval Ct0,r]. For
most problems, the solution curves neither diverge rapidly nor converge
rapidly. For these problems, ri^-T^ and the stepsize necessary for an
accurate integration may be plausibly estimated from the form of y(t), i.e.,
on the need for keeping the truncation errors small. As discussed above, a
8
-------
severely unstable problem requires a much smaller stepsize, and so is much
more expensive to integrate accurately. It might be supposed that the
opposite case of a highly stable problem, 3f/3y«0, would present no partic-
ular difficulty. Unfortunately, this is not true, for reasons which are
discussed next.
§2.4 Stability of an Integration Method
If 3.f/3y<0, the problem (2-5) is stable, so that errors in the numerical
solution are normally expected to be damped out as the integration proceeds.
The point to be made here is that this is untrue unless the stepsize is kept
smaller than a maximum value which depends both upon the integration method
and upon the problem being solved. l
For purposes of illustration, advance the solution of (2-5) from t£ to
method (2-6a),
y. = y(t.) - e.
To determine the fate of an error e£ in y^, substitute
(2_n)
. . u n)
into (2-10), obtaining
ei+l = ei * ^'i+l* " y
-------
exceeded, the method is said to be unstable: errors occurring in one step
tend to grow in subsequent steps, even though the problem itself is stable.
For the solution of a system of equations (2-?), the role of 3-f/3y is
played by the eigenvalues of the Jacobian matrix Jr. Since these eigenvalues
are, in general, complex, it is customary to describe the stability of an
integration method by reference to the standard problem
y'(t) = \y
<+ ^ (2
y(to) = y<>
where y is complex and X is a complex constant. When ReX<0, this problem is
stable. For the Euler method (2-10), the analysis proceeds precisely as
before, with X in place of 3f/3y, except that the approximation in (2-12) is
now exact. The stability condition for the Euler method is thus
| 1 + hX | < 1 (2-14)
The region of stability of an integration method is that portion of the
left half of the complex plane in which hX must lie if the method is to be
stable when applied to (2-13). Each integration method has its own charac-
teristic stability region. The stability region of the Euler method, as
given by (2-14), is the interior of the unit circle centered at -1+Oi.
The stepsize is also bounded by the need to keep the local truncation
error smaller than an acceptable limit. If, as is normally the case, this
bound is more stringent than that imposed by the stability condition, the
stepsize is limited for reasons of accuracy. If, on the other hand, the
stability condition imposes the more stringent bound, the stepsize is
limited for reasons oiT stability.
Stiff problems are those for which most integration methods limit h for
reasons of stability. Typically, the solution is smooth and slowly varying,
and so would appear easy to follow; but 3-f/3y«0, requiring a much smaller
stepsize than is needed to keep the local truncation error within bounds.
Thus, a curve which should be easy to compute requires much more work than
expected. In extreme cases, the stepsize limitation may be so severe that
roundoff errors make it impossible to proceed at all.
The difficulty presented by stiff problems may be pictured as follows:
The computed points do not lie exactly on y(t), but rather lie on nearby
solution curves which, because of their rapid convergence toward y(t), vary
much more rapidly than y(t) itself. Consequently, for most integration
methods,
»
|t.| (2-15)
and the stepsize must be small encough to follow these rapid variations, and
cannot have the much larger size which would be sufficient to follow the
slowly varying y(t). Although controlling the local error closely enough
10
-------
will prevent the method from becoming unstable, the integration will be very
expensive unless it is possible to employ a method for which the inequality
(2-15) is not so severe.
It is rare that an initial value problem is stiff over the entire inter-
val of integration. Usually, the solution has an initial transient which
decays rapidly to approach a slowly varying curve. Until this transient has
decayed to an insignificant value, the steps ize limitation is one of accuracy
rather than one of stability, and the problem is not stiff over this initial
interval .
In order to solve stiff problems efficiently, i£ is necessary to seek
methods which do not limit h for reasons of stability. An example of such a
method is the backward Euler method (2-6b),
- y. + Wi+1,y) (2-16)
For the standard problem (2-13), this implicit method yields a linear
equation which is easily solved to give
- h\
The substitution (2-11) then gives
Thus, the stability condition for the backward Euler method is
| 1 - h\ \ > 1
and the stability region occupies the entire left half of the complex plane.
Such a method is suitable for solving stiff problems, since it imposes
no stability limitation on the stepsize, which is then limited only for
reasons of accuracy. No explicit method has such favorable stability proper-
ties, but applying an implicit method entails the additional labor of solving
a (generally nonlinear) algebraic equation at each step. For a sufficiently
stiff problem, the ability to use a larger stepsize will more than compensate
for this additional work.
One approach to solving (2-16) is to refine an initial approximation
y. . . by successive iterations
Unfortunately, this simple iteration will not converge for stiff problems.
Define
11
-------
the amount by which the mth iterate differs from the solution. Then it is
easily found that, to a first order approximation in 6m,
Therefore, convergence is not obtained unless |M3^/3y)| < !• For stiff
problems, 3f/3y«0, this implies a limitation on h of the sort which must be
avoided.
To solve (2-16) without this limitation, it is necessary to employ a
Newton iteration,
(2-17)
For nearly all practical problems 3f/3y varies only slowly, and the Newton
iteration will converge for fairly large values of h.
If an ^-equation system (2-2) is being solved, then (2-17) is a vector
equation in which ( l-MS-f/Sy))"1 is replaced by the n x n matrix
r1
where I is the identity matrix and Jf is the Jacobian matrix of f. The need
to evaluate the n2 components of the Jacobian and to invert an n x n matrix,
which is common to all methods suitable for stiff equations, makes the
numerical solution of very stiff problems generally intractable for large n •
Nonetheless, these methods can be very useful for systems of moderate size.
12
-------
SECTION 3
SELECTION OF NUMERICAL INTEGRATION PROGRAMS
The foremost consideration in the construction of a numerical integra-
tion routine is that of efficiency: of two programs, the more efficient
requires less computing to solve the same problem. The efficiency of a
program depends not only on the basic integration method used, but also on
the particular coding which is used to implement that method. Thus, compar-
isons of efficiency can generally be made only between the complete codes,
and not between the integration methods themselves.
Efficiency, however, cannot be the only criterion in constructing a
practical code. A user will certainly desire that the code be reliable:
it should not fail on routine problems, and it should not return erroneous
results while appearing to operate normally. Furthermore, the code should
be easy to use, requiring little more than a statement of the problem to be
solved. While the user must be concerned with the nature of his problem, a
good code should automatically handle all concerns which arise in the
solution process. In selecting the codes to be incorporated in this report,
attention was restricted to those existing codes which meet these criteria
of reliability and ease of use. While no code is equally efficient for all
problems, each of the three codes presented here may be regarded as repre-
senting state-of-the-art performance when applied to those problems for which
it is intended.
§3.1 General Considerations and Terminology
The cost of solving (2-1) numerically may be viewed as
cost = (number of steps required) x (amount of calculation/step)
For a given integration method, the amount of calculation per step is fixed,
and the cost is minimized by choosing each stepsize no smaller than is
necessary to achieve the desired accuracy in the solution. If the program is
to be reliable and easy to use, it must contain provision for choosing this
optimum stepsize automatically. As mentioned in §2.3, the programs given
here select the stepsize to keep the local error at each step smaller than
some limit requested by the user. As the integration proceeds, this stepsize
is increased or decreased to reflect the local requirements of accuracy.
These programs are thus called variable-stepsize codes.
Of two methods, that having the smaller truncation error will, ceteris
paribus, be able to solve the problem at lower cost by using larger steps.
However, the construction of an efficient variable-stepsize code requires the
ability to form a reliable estimate of the local error without greatly
13
-------
increasing the amount of calculation per step. Consequently, a number of
methods which are excellent in the sense of having small truncation errors,
are unsuitable because they do not allow a cheap estimation of the local
error.
If the solution is sufficiently smooth, the use of a higher order inte-
gration method will, because of the smaller truncation error, allow use of a
larger stepsize. Although higher order methods require more computation for
each step, the stepsize can usually be sufficiently larger to produce an
overall increase in efficiency. However, the use of a high order method for
a problem with discontinuities in a low order derivative will be less
efficient, since no increase of stepsize will then be possible. For this
reason, high order methods are best incorporated into programs which can
automatically choose an order which is suitable for the particular problem
being solved. Programs which can do this are termed variable order codes.
Part of the calculation performed in each step of the integration is in
the form of function evaluations, the evaluation of the derivatives in (2-1)
at a given (t,y). The remaining calculation performed in the step is termed
overhead. The division of labor between function evaluations and overhead
varies widely among integration methods. Some methods require many function
evaluations, but relatively small overhead; while others employ fewer function
evaluations at a cost of increased overhead. Methods of the former type will
be more efficient for problems whose function evaluations are cheap, while
those of the latter type are preferable for problems whose function eval-
uations are expensive.
Increasing the stringency of the local error tolerances will always
increase the expense of the integration, since smaller stepsizes must be
used. However, the increase in cost will occur at different rates for dif-
ferent methods, so that the method of choice for use at moderate accuracy
may be different from that which is best for extreme accuracy.
To monitor the accuracy of the solution and to choose an appropriate
stepsize, the code must estimate the local error at each step. In a typical
case, the local error may be estimated as follows: using the same stepsize h,
compute y. by a method of order k (less accurate solution)
compute 0. by a method of order k+] (more accurate solution)
estimate local error in y. as T). = y. - y.
k+1
Estimated in this way, r)i=0(h ) and so is the local error in the less
accurate solution. The stepsize control is thus based on the less accurate
solution. However, since the more accurate solution y^ has already been
found, it is reasonable to accept'it as the computed value and to use this
value when continuing the integration. In effect, one improves the accuracy
of y^by adding the estimated local error to it. This procedure, which is
called local extrapolation, may be used whenever there.is good reason to
believe that the local error estimate is quite accurate. It may be noted
that the above method for estimating n£ is practical only when the calcu-
lation of y^ can be done using very little computation beyond that already
needed for y^: integration methods which do not allow this cannot economically
estimate the local error in this way.
14
-------
§3.2 Control of Stepsize
As discussed above, achieving maximum efficiency for a given integration
method requires that each stepsize be .chosen as large as is permitted by
the constraint
n < E
where n is the absolute value of the local error in the step, and E is an
error tolerance specified by the user of the code.
k+1
Suppose that a method of order k is being used, so that n=0(fc ). Then
the basic procedure for error control and stepsize adjustment is as follows:
1.) Perform the step with the (previously chosen) stepsize h and
estimate f|.
2.) If r\£, the stepsize was too large. Since n=0(^k+l),
the stepsize h1 which will exactly meet the error criterion may be
estimated as
1
rk+1
• h (3-1)
3.) If r]>E, then the step did not meet the required error tolerance,
so repeat it with the new (smaller) stepsize h'. If r\E.
Since each failed step must be repeated, efficiency will be lost unless the
number of such failures is kept small by a conservative choice of stepsize.
This problem is handled somewhat differently in the various codes, so the
details are reserved for the detailed code descriptions in Sections 4, 5,
and 6. The choice of an initial stepsize for beginning the integration will
also be described there.
§3.3 Choice of the Codes
The codes chosen for inclusion in this report are incorporated into
three packages of subroutines:
Rjwge-Kutta programs: A fixed-order, variable-stepsize code employing
the single-step Runge-Kutta-Fehlberg fourth-fifth order method.
Local extrapolation is used, the computed solution coming from the
fifth order formula. This program is most efficient for nonstiff
problems whose function evaluations are cheap, and for which
extreme accuracy is not requested.
Adams method programs: A variable-order, variable-stepsize code
employing the multistep Adams PECE methods of orders 1 through 12.
15
-------
Local extrapolation is used, so that the effective orders are 2
through 13. This program is most efficient for nonstiff problems
whose function evaluations are fairly expensive or for which high
accuracy is requested.
Stiff method programs: A variable-order, variable-stepsize code
employing the multistep backward differentiation methods of orders
1 through 6. This program is intended only for solving stiff
problems. When applied to nonstiff problems, it will be less
efficient, and often less accurate, than the other two methods.
The integration routine for the Runge-Kutta programs, RKFS, is a
modified version of that given in Shampine(1976). On the basis of tests
reported there, and after further tests comparing RKFS with some Runge-Kutta
methods recently proposed by Bettis(1978), RKFS was selected as the most
efficient of those presently existing programs which utilize Runge-Kutta
methods of fifth order or less. The Runge-Kutta-Fehlberg method is described
in Section 4, which also describes the coding of the Runge-Kutta programs.
The power and efficiency of variable-order, variable-stepsize codes
employing Adams methods was first recognized in codes written by Krogh(1969)
and Gear(1971,197 la). The code STEP given here is taken from Shampine(1975),
and represents the most sophisticated level of development of such codes.
Section 5 explains the Adams PECE methods, and also describes the coding
of the Adams method programs.
As discussed in §2.4, stiff problems require special methods for their
efficient solution. The backward differentiation methods which are used for
the stiff method programs are described in Section 6. The program GEAR
which incorporates these methods is a modified version of an algorithm
originally developed by Gear(1971,197 la), the chief modification being a
more conservative stepsize selection procedure. While this is presented as
a general-purpose code for stiff problems, it cannot, for reasons mentioned
in §2.4, be expected to solve very large systems of equations. The handling
of very large stiff problems generally requires much special-purpose coding
tailored to the particular problem, and so falls outside the purview of this
report.
One '•las.: of methods not represented here is that of the extrapolation
methods, whic1. are single-step methods suitable for incorporation into a
variable-order, variable-stepsize code. The highest state of their develop-
ment in a practical code is represented by the program DODES described in
Schryer( 1975). Such codesi.are best suited for problems whose function
evaluations are cheap, but for which high accuracy is required. For these
problems, extrapolation codes can exceed the efficiency of Adams codes,
which are, in turn, more efficient than Runge-Kutta codes for work at high
accuracies. However, the overall reliability of extrapolation codes is
somewhat uncertain because of unresolved theoretical and practical questions
about their behavior upon encountering stiffness or lack of smoothness in the
solution. The chief reason that they have been excluded from further con-
sideration here is that the efficiency of all extrapolation methods tends to
be very sensitive to the choice of output points: since the average user
desires to specify the particular points at which solution values are to be
returned, he is likely to find this feature a serious inconvenience.
16
-------
§3.4 la/natations of Precision
As was seen in §2.3, the accuracy which can be achieved is limited by
the precision of the arithmetic used to perform the computations. The
precision of the arithmetic is most conveniently described by the computer
unit roundoff error u, which is the smallest positive value for which
l.+u > 1. in floating point arithmetic. Typical values of u for single
precision arithmetic range from 9.54xlO~7 for the IBM 360 to 7.11xlO~15 for
the CDC 6600.
During the numerical integration, it is clear that one cannot make
progress with a stepsize so small that t+h=t in floating point arithmetic.
Thus, each stepsize must be chosen to be at least as large as tu, where t is
the current value of the independent variable. In fact, h must be chosen
somewhat larger than this in order to allow the stepsize adjustment mechanism
some range in which to operate. All of the codes in this report require h to
be at least as large as 4ut. In the case of the Runge-Kutta code, the
properties of the integration formula require a much larger limit, and h is
limited to a minimum value of 26ut.
Similarly, it is not possible to require the solution to have an error
smaller than yu, which is the error inherent in its correctly rounded value.
The error tolerance E must be somewhat larger than this ultimate limit if
stepsize adjustment is to work properly. The codes all require that E be at
least as large as 4uy, where y is the current value of the solution.
If, at any point during the integration, the codes find that the
user's requested error tolerance is too small, or that meeting that tolerance
requires too small a stepsize, the integration is stopped and the user is
informed of what (larger) tolerance will allow a successful continuation.
This provision is necessary for the reliability of the codes, so that they
will not return meaningless or inaccurate values heavily contaminated by
roundoff errors.
17
-------
SECTION 4
RUNGE-KUTTA METHODS
An r-stage explicit Runge-Kutta method advances the solution of (2-5)
according to
where
(4-1)
**-
i+C\ ' Vj=1 Aj j
£-1
I
J-l
r
I Y,
A=l *
2j 3) •
The number of stages, r, is the number of evaluations of f needed to perform
one step. If the g and y coefficients are chosen so that T^+i is of order
h8**, then the method is of order s.
A Runge-Kutta is conventionally specified by presenting the coefficients
in a tableau:
a
32
18
-------
The program RKFS employs the following six-stage, fifth order method derived
by Fehlberg(l969):
0
1
4
3
8
12
13
1
1
2
1
4
3
32
1932
2197
439
216
-8
27
16
135
9
32
-7200
2197
-8
2
0
7296
2197
3680
513
-3544
2565
6656
12825
-845
4104
1859 -11
4 104 40
2856 1 -9
56430 50
2
55
(4-2)
The utility of this method derives from the fact that it has embedded within
it the five-stage, fourth order method given by the first five rows of the
above table together with the y coefficients
(4-3)
25
216
0
1408
2565
2197
4104
-1
5
Thus, at a cost of six function evaluations, both a fifth order solution and
a fourth order solution are obtained. The difference of these two results
provides an estimate of the local error in the fourth order solution, and
this estimate is used for error control and stepsize selection. However, the
fifth order solution is the value actually accepted; i.e., local extrapo-
lation is used.
The remainder of this section describes the implementation of this
method in the code RKFS, and should be read in conjunction with the program
listing in Appendix A and the instructions for use in Section,7. RKFS is
very similar to the program in Shampine(1976), which reference should be
consulted if more detailed information is needed.
Applying (4-1) to a system of n equations (2-2) is accomplished simply
by interpreting the y's, fs, and &'s as n-dimensional vectors. At each step
of the integration, there is then a vector of local error estimates n(j) to
be compared with a vector E(j) of local error tolerances, j = 1, 2, ..., n.
The £(j) are computed from user-specified relative and absolute error
tolerances, as described in Section 7. For a successful step, it is required
that n(j) < E(j) for all components of the solution: j = 1, 2, ..., n. For
stepsize control, it is necessary to form a single number which characterizes
the maximum error occurring in any of the components. This is done by
19
-------
calculating
max
n --- (4~4)
• n
so that p<1 for a successful step and p>1 for a failed step. The stepsize to
be used for the next step, or to repeat the last step if it failed, is
then computed as
h = h x 0.9 x pI/5 (4-5)
new old
where the factor 0.9 causes the new stepsize to be chosen conservatively in
order to avoid unnecessary step failures. Since the assumptions underlying
(4-5) are not valid if the solution is not sufficiently smooth, possible
difficulties are avoided through limiting its operation by three constraints:
1.) Stepsize increase is limited to a factor of 5.
2.) Stepsize decrease is limited to a factor of 0.1.
3.) The stepsize is not allowed to increase for the step following
one which failed on the first attempt.
While (4-5) specifies how the stepsize is to be adjusted during the inte-
gration, a separate procedure is needed to select a stepsize for beginning
the integration. The starting stepsize need not be highly accurate, since
it will be adjusted to an appropriate value by (4-5) after the first step
is attempted. If the first step were taken with a constant approximation,
yi=yo» then the error would be estimated as #tXt0,yo). Estimating the error
in a fourth order method as h* times the error in this "zeroth order" method,
an appropriate starting stepsize may be estimated as
mm
n
1/5
(4-6)
where Eo(j)=#r|yo(j)\+Ea» ET and Sa being the user-specified relative and
absolute error tolerances RELERR and ABSERR, and y0(j) the initial value of
the j-th component of the solution. (4-6) is supplemented by two rules:
1.) h is not allowed to be smaller than 26u*max(|t0|,At), where At
is the distance from t0 to the first output point. This is
the value used if all the E0(j) vanish.
2.) h cannot be larger than At. This deals with the case in/
which all the initial derivatives -fi(to,yo) vanish.
RKFS performs a number of checks designed to protect the user from
improper use of the code. Attempts at impossible precision are prevented
by requiring that the relative error tolerance RELERR be at least as large as
the machine dependent constant REMIN, whose value is set in a DATA statement.
20
-------
A request for greater accuracy results in an error return signaled by
IFLAG=3, and with RELERR set to the value of REMIN. The value of REMIN
should not be set smaller than 4", where u is the computer unit roundoff
error. It is further recommended that REMIN not be smaller than 10~12, since
it is more efficient to use Adams methods if greater accuracy than this is
required.
To distinguish the values of the arguments t+a^h = t+(i2/13)£ and
t+ash m t+h, which appear in the evaluations of k^ and k5 in (4-0, the
stepsize is never allowed to be smaller than 26ut. If the code finds it
impossible to satisfy the error test using the smallest permissible stepsize,
it returns with IFLAG=4 and with RELERR and ABSERR increased to values which
should be achievable.
Unlike Adams codes, which can interpolate to obtain the solution at an
arbitrary point, Runge-Kutta methods must adjust the stepsize so that the
output point will be a mesh point of the integration. Unless some foresight
is applied, the last step to the output point may be exceedingly small. This
would impair the efficiency of the stepsize selection procedure (4-5), and
might even require a step smaller than 26"t. RKFS avoids this difficulty by
looking ahead two steps of the current stepsize. If the output point lies
within this interval, then the stepsize is adjusted as follows: If the
output point is within A/0.9 if the current value of t, the stepsize is set
to complete the integration in one step. Otherwise, the stepsize is set to
half the distance to the output point, so that the integration will be com-
pleted in two steps.
While this procedure greatly reduces the adverse impact of output points
on stepsize selection, inefficiency may still occur if the user requests a
sequence of output points which are spaced much more closely than the natural
stepsize which would be selected by (4-5). To inform the user of this, RKFS
counts, in the variable KOP, the number of output intervals which are less
than half the natural stepsize, and returns with IFLAG=7 if KOP exceeds 100.
Another possible difficulty, which should happen very seldom, occurs
when the user requests an output point at distance less than 26ut from the
previous output point. In this case, RKFS simply extrapolates the solution
by the Euler method. No error check is made, and there is no special
indication that this has been done.
In order to inform the user when a problem is unusually expensive to
integrate, RKFS counts, in the variable NFE, the number of function evalu-
ations which have been performed since the start of the integration. If NFE
exceeds the value of MAXNFE, which is set in a DATA statement, RKFS returns
with IFLAG=5 or IFLAG=6.
Since a common cause of an expensive integration is stiffness of the
problem, RKFS incorporates a procedure to detect stiffness, as described in
Shampine(1977). This test for stiffness is based on a first-second order
pair of Runge-Kutta formulae which use the same function evaluations as (4-2)
with the Y coefficients given by:
1094951
13000000
-2120820
13000000
9893169
13000000
5275998
13000000
-1715610
13000000
572312
13000000
(4-7)
WWW UVSVUUUU I J\J\J\J\J\J
-------
sei;uuu ULUCI. uiciuuu.
1815846
13000000
-2582209
13000000
9417746
13000000
5576389
13000000
-1839305
13000000
611533
13000000
(4-8)
This pair may be used to estimate the local error in the first order method,
in precisely the manner described above for the fourth-fifth order pair
(4-3) and (4-2). Although (4-7) and (4-8) are less accurate than (4-3) and
(4-2), their stability regions are substantially larger than those of the
fourth-fifth order methods. Thus, (4-7) and (4-8) are unlikely to meet the
error test unless (4-2) and (4-3) have severely restricted the stepsize for
reasons of stability. During the integration, RKFS examines successive
blocks of 50 successful steps. In each block, the number of steps which
pass the first-second order error test is counted in KST. If KST reaches 25,
the problem is labeled as stiff, and it is unnecessary to test subsequent
blocks. Since stiffness is not a problem unless it makes the integration too
expensive, a finding of stiffness is not reported to the user unless there is
an error return for too much work. Then, IFLAG=6 indicates that stiffness
was found. While tests indicate that this method is quite effective in
detecting severe stiffness, moderate stiffness may not be reported, partic-
ularly at stringent error tolerances.
22
-------
SECTION 5
ADAMS METHODS
The solution of (2-5) at t may be written as
y(V + J n+1J? 1
and ^ is in the interval
23
-------
The Adams-Moulton method of order k differs from the Adams-Bashforth
method in replacing f(t,y(t)) with P. -H^' so that yn+l is 8iven
/•n+l
y , = y + I P, ,(t) dt (5-5)
n+l yn J k,n+l
n
When the stepsize is constant, the truncation error in (5-5) is
(5-6)
where y* = 1 and y* = Y- ~ Y-_i
Since P^ n+l^) depends on yn+l. (5-5) is an implicit method. For k=l, it is
the backward Euler method (2-6b). Although applying the implicit Adams-
Moulton method requires additional labor in solving (5-5) for yn+i» the
truncation error is smaller, and the stability region larger, than for the
same order Adams-Bashforth method. Thus, the Adams-Moulton method is able to
take larger steps while achieving the same accuracy, so that it is usually
more efficient in spite of requiring more work for each step. Rather than
attempting to solve (5-5) exactly, the code STEP uses the Adams-Bashforth and
Adams-Moulton methods in combination as follows:
1. The Adams-Bashforth method of order k is used as a predictor to
obtain a predicted solution
/•
J
Pkn(t)
n
2. The derivative r~p =r~(t ,,p ,) is evaluated using this predicted
, . . n+l n+1 n+1 e
solution.
3. The Adams-Moulton method of order k+1 is used as a corrector, with
the k-degree interpolating polynomial ^k+l,n+l being approximated by
***• given by
j j = 1, 2, ..... k
P* (t ) = fP (5-8)
k+l,n n+l n+l
so that the corrected solution is given by
//"
n
4. To complete the step, the derivative fn+\= f(tn+i,yn+i) is evaluated
using the corrected solution. This derivative is needed to begin the
next step.
24
-------
Because of the successive operations of prediction, derivative
evaluation, correction, and final derivative evaluation, this method is
conventionally referred to as a PECE method. Since it begins with the Adams-
-Bashforth method of order k, it will be called a PECE method of order k.
However, the use of a corrector whose order is one higher than that of the
predictor, which may be viewed as local extrapolation, causes this PECE
method of order k to have an effective order k+1. Subroutine STEP employs
the Adams PECE methods of orders 1 through 12, with local extrapolation
giving effective orders of 2 through 13. During the integration, the
stepsize and the order of the method are adjusted in order to proceed
efficiently while keeping the local errors within the bounds specified by the
user of the code.
The Adams methods may be expressed in many different equivalent forms,
corresponding to different algebraic representations of the interpolating
polynomials in (5-7) and (5-9). For practical computing, preference must be
given to those forms which facilitate changes in stepsize and order, and
which are economical of storage. The derivations of the formulae given below
may be found in Shampine( 1975) . These derivations are quite lengthy and will
not be repeated here .
The divided difference (or Newtonian) form of the polynomial
which interpolates to f , f , . . . , f _, is
Pkn(t) =
where
a
,t ..,..., t . .3 = i f. . n
-
. . . t t
* i=0 J s=0 j-i j-s
The quantities ^[•••3 are called divided differences because they obey the
relation
fCt./t. ,...,t._j D - r"[t. j,... ft. £+ j,t. ^3
r~Ct./t. ....ft. j —
J J-l J"* t. - t.
(5-11)
Starting with fn»^n_! fn-k+l' (5~H) may be used to 8enerate a11 the
divided differences~needed for the interpolating polynomial (5-10).
In the code STEP, the interpolating polynomials are represented in terms
of modified divided differences defined by
25
-------
4>. 3
Then the predicted solution (5-7) may be written as
l
k
I
26
-------
where (j>*(n) = BjdH-1)^(11) and the coefficients
c. (s)
0 ln
may be calculated recursively from
where
The corrector polynomial P* interpolates to the same data as P, plus the
value k+1'n kn
and so is given by
where the superscript p indicates that the divided difference is generated
using fP . rather than f ."^(t^,, »_._,). The corrected solution (5-9) is
0 n+ 1 n+ 1 n+ 1 n+ 1
then
where < .(n+1) may be calculated from
To complete the step, it is necessary to compute the divided differences
<|>.(n+l) to be used in beginning the next step. They are given by
27
-------
A careful arrangement of the computation avoids the need for separate
arrays to store the 4>.(n) and the *(n). To achieve this economy of storage,
STEP introduces intermediate quantities (|>f(nH-l), and successively overwrites
the divided differences,
4>?(n-H)
The calculations performed by STEP in advancing the solution from t to t
are as follows:
computing 3-(n+l) for i=l,2,...,k and g. . for i=l ,2, . . . ,k+ 1
predicting <(>*(n) = B-(n+l)4>.(n) i=l,2 ..... k
k
Z, g
-i- (|)*(n) i=k,k-l
evaluating
correcting y = P + ( - *>+!»
evaluating
i=k,k-l, . . ., 1
This formulation is particularly convenient for changing the order of
the method. Since the step taken with order k evaluates (^k+iCn+l), this
difference may be used to raise the order to k+1 for the next step. Lowering
the order from k to k-1 is accomplished simply by not using ^^jCn+l) and
in starting the next step.
When the stepsize is constant, it is not necessary to recompute all the
coefficients. Let ng be the number of successive steps taken with the
28
-------
constant stepsize h, including the current step. Then
= ih i = 1, 2 ..... ng-l
a.(n+l) = a.(n) =4- 1*1,2 ..... n -1
-L 1 1 &
B.(n+l) = 3.(n) = 1 i = 1, 2 ..... n -1
11 S
and the remaining values are computed from
\|j (n+1) = n h
Yns s
4). (n+1) = \)j. (n) + h i = n +1, ..., k
1 1~* J S
a (n+1) = —
ns ns
£ = V1 ..... k
8 (n+1) = 1
ns
*.
g.(n+l) = 3. ,(n+l) -=— = - i = n +1, ..., k
*•'
In STEP, the coefficients g^j are stored in the array G(i), i=l ..... k+1.
For i=l ..... ns, the G(i) are unchanged from the previous step. The infor-
mation necessary for generating the remaining G(i) is stored by the previous
step in the array V(q) according to
v(2) = gn 2
•
•
•
V(k+2-ns) = gn
V(k+3-ng) = gn
V(k+4-ns) = srng.2fk+4^l8
•
•
*
V(k) = 92k
29
-------
For this step, V(q) is updated according to
V(q) «-V(q) - a (n+l)V(q+l) q=l ,2, . . . ,k+l-n (5-19)
ng s
The G(i) are generated using a working vector W(q) which is initialized by
W(q) «-V(q) q = 1, 2, ..., k+l-n (5-20)
S
Then G(ns+l)=W( 1), and the remaining G(i) are obtained by calculating in the
order i = n +2, . .., k+1:
S
W(q) «-W(q) - ou (n+l)W(q+l) q = 1, 2, ..., k+2-i
(5-21)
If the order is lowered, no special action is required in computing the
coefficients, since the only effect is that some of the stored values are not
used. If the order is raised, k=k -,.+ !» V must be updated according to
V(k) f
V(k-j) + V(k-j) - a. ,(n+l)V(k-j+l) j-1,2 n -2
j + j s
after which the computation of G(i) proceeds as before starting with (5-19).
If the stepsize is changed, then ns=l and all the coefficients must be
recomputed. In this case, V and W are initialized as
W(q) «-V(q) * q(q+1) q = 1, 2, ..., k
and the G(i) are then computed according to (5-21).
The PECE method of order k produces the k+1 order solution yn+i given
by (5-17). If the corrector formula were of order k rather than order k+1,
then the k order solution
T\
y , (k) = p , + h g. (b, , (n+1)
n+1 n+1 n+1 klYk+l
would be obtained. STEP uses the difference between these two solutions as
an estimate of the local error in the step:
ERR" ^^^^0! - <5-22>
The step is accepted or rejected depending on whether or not ERR
-------
test does not require the final function evaluation, so this evaluation is
not done if the step fails.
Regardless of whether the step is successful, the following quantities
are used to test whether the order should be lowered:
ERKM2 = |H-^_20K.,(n+l)<|{_I(n+l)|
ERKM1 = |fcY_a(n-H)<|>(n+l) (5-23)
ERK =
where the a's are computed along with the other coefficients by using
•a-(n-H) = 1 i=l,2,...,n +1
8 (5-24)
= (i-l)a. .(n+l)a. ,(n+l) i=n +2,..
- 1~ 1 1" 1 S
ERKM2, ERKM1, and ERK are estimates of the errors which would occur at,
respectively, orders k-2, k-1, and k, if a constant stepsize were being used.
The order is lowered to k-1 if
k > 2 and max(ERKMl ,ERKM2) & ERK
(5-25)
k = 2 and ERKM1 £ 0.5 ERK
A test for raising the order is made only if (1) the step is successful,
(2) the stepsize is actually constant (ns > k+1), and (3) it has not already
been decided to lower the order. If these three conditions are fulfilled,
the quantity
ERKP1 = |^Yk+,k+2l (5"26)
where
Wn+1) = Wn+1) - Wn)
is calculated as an estimate of the error at order k+1. The order is raised
if
1 < k < 12 and ERKP1 < ERK
(5-27)
k = 1 and ERKP1 < 0.5 ERK
After the order is selected, it is renamed k, and the corresponding
error estimate for constant stepsize is renamed ERK. To economize in the
evaluation of the coefficients, the stepsize selection procedure is biased in
31
-------
favor of using a constant stepsize. Stepsize increase is made only by
doubling, which is performed only if
k+l
0.5 EPS > 2 ERK
If the stepsize cannot be doubled, the current stepsize is retained provided
0.5 EPS > ERK
If this condition is not satisfied, the stepsize is reduced by the factor
0.5 EPS^*1
ERK /
except that r is not permitted to be greater than 0.9 nor smaller than 0.5.
By aiming for an error of 0.5 EPS rather than EPS, the stepsize is
chosen conservatively so as to prevent unnecessary step failures.
To start the integration, the method of order 1, which requires no
memorized values, is used. The user provides an upper bound -&max for the
size of the first step, and STEP computes a starting stepsize from
in
\
0.5 EPS
min I n , 7-
1 max 4
During the first phase of the integration, signaled in the code by the
logical variable PHASE 1 s,et .TRUE., the stepsize is doubled and the order is
raised after every step. The purpose of this device is to reach the stepsize
and order appropriate for the problem as rapidly as possible. This first
phase is terminated when there is a step failure, when (5-25) says to lower
the order, or when the maximum order of 12 is reached.
After leaving the starting phase, STEP follows the previously stated
rules for order and stepsize selection with one exception: if three succes-
sive failures occur in the same step, all memorized values are discarded and
the order is set to 1. This is done because such a sequence of failures
indicates a drastic change in the character of the solution, so that the use
of higher order methods cannot be justified.
The interpolating polynomial £|!+j n in (5-9) can be used not only to
approximate the solution at tn+j, but also at any t near tn+j according to
y(t) e u + / P* (T) dT
n j. K+l,n
n
This allows the solution to be obtained at points which are not mesh points
if the integration, so that it is not necessary of adjust the stepsize to
32
-------
make mesh points and output points coincide, as must be done for Runge-Kutta
methods. The subroutine SINTRP, when used in conjunction with STEP, performs
this interpolation of the solution.
Let t be an output point such that
out
t < t < t ,
n - out - n+1
Then the solution and its derivative at t are calculated by SINTRP as
k+1
yout ' yn+l + hl £ *i,4
k+1
yout ' .
where the coefficients are generated from
giq = q
and
Pi' '
where
ril
r
n
Tl
t
out
By interpreting the foregoing formulae in terms of vectors, their
application to a system of n equations (2-2) is straightforward. The only
modification required is in the error estimates (5-22), (5-23) and (5-26),
which are computed in STEP as
33
-------
error
max
*n
error(L)
WT(L)
where error(L) is the error estimated for the L-th component of the solution,
and WT(L) is a vector of nonzero weights specified by the user of the code,
as described below in §7.7.
STEP is functionally organized into five blocks:
Block 0: Checks that the stepsize and error tolerance are not too
small for the machine precision, and performs the
initialization necessary for the first step of the
integration.
Block 1: Computes the a, 3, a, 4», and g coefficients needed for the
step.
Block 2: Computes the predicted solution, evaluates the derivatives
using the predicted solution, tests whether the order
should be lowered, and tests whether the step is
successful. If the step succeeds", control is passed to
Block 4. If the step fails, control passes to Block 3.
Block 3: Restores variables for repeating the failed step, sets
the order selected in Block 2, halves the stepsize, and
transfers control to Block 1 to repeat the step.
Block 4: Corrects the predicted solution, evaluates the derivatives,
updates the divided differences, sets the order and
stepsize appropriate for the next step, and returns to the
calling program.
As a means of controlling propagated roundoff error, STEP incorporates
an alternate route of computation in the DO loops terminating at statements
25, 235, and 405. The operation of this device is described in Chapter 9 of
Shampine(1975). Since roundoff control is not important except near the
limits of the machine precision, these loops are not used if the error
tolerance is greater than 100 times limiting precision on the first call to
STEP.
The code ADAM is an interfacing routine designed to provide the user
with an access to the Adams methods which is similar to that provided by the
Runge-Kutta code RKFS. The test for stiffness which is used in ADAM is that
described in Shampine(1975): the problem is labeled as stiff if a sequence of
50 consecutive steps are taken with order less than or equal to 4. Since
the lower order PECE methods have larger stability regions than higher order
methods, they will tend to be selected when the problem is stiff, as their
larger stability regions will allow larger steps to be taken. While this
procedure will also label as stiff some nonstiff problems which*happen to be
most efficiently integrated by low order methods, most false indications of
stiffness are avoided because a finding of stiffness is not reported to the
user unless the code also finds that too much work is required to perform the
integration.
34
-------
SECTION 6
STIFF METHODS
Among methods suitable for solving stiff problems, the backward
differentiation methods of orders 1 through 6 are the best understood and
most widely used. Their use was popularized through a code written by
Gear(1971,197 la), for which reason they are sometimes referred to as Gear's
method.
The solution of (2-5) at t. satisfies
L = f(t y(t )) (6-1)
'i x *
The backward differentiation method of order k approximates (6-1) by
replacing y(t^) with y^, and by replacing y(t) with the (unique) k degree
polynomial ^^(t) which interpolates to y^ and to the computed solution
values at the k preceding points. This yields an implicit integration method
in which y^ is obtained by solving the (generally nonlinear) algebraic
equation
0/.(t.) - ^(t.,y.) (6-2)
For k=l, this is the backward Euler method (2-6b).
Suppose that the solution has been obtained up to tj_j, so that
is known. When the solution is advanced to t£, the new polynomial
is obtained, which satisfies (6-2),
0ki (6-3)
and which interpolates to y. ,, y. 0, ...» y. t»
/ 1— I l—^r !"•»,
^ki'W ' yi-* = Ck,i-.(ti-^ (6-4)
A/ = 1 y Z y •••} K
Define the k degree polynomial pki(fc) by
35
-------
From (6-4), P, . must satisfy
A/ = 1 y
y K
(6-6)
The condition (6-6) completely determines P^i up to a constant factor. To
simplify the discussion, assume that the stepsize has been constant for at
least k steps. Then P . (t) may be written as
where
(6-7)
The coefficients of fi. (x) = a + ax + ax + •
f-ii'j-i-i K I t, j
following table,
+ a x are given in the
JC+ 1
k =
-a,
"a2
-a3
-a4
-a5
-6
-37
. 2 6 12
3 11 25
1111
1 6 7
3 11 10
1 1
11 5
1
50
60
137
1
225
274
85
274
15
274
274
20
49
1
58
63
5
12
25
252
1
84
1
1764
Rewriting (6-5) as
t-t.
(6-8)
36
-------
and substituting this into (6-3) gives, after a little rearrangement,
yi(0)
The algebraic equation (6-9) may be solved for e^ by a Newton iteration
starting with ei(0)=0 and proceeding according to
' O^aM) .
= yi(m) + Vi
y
= e
i(m) + Ci
where 3f/3y is evaluated at (^.y., J. If h is sufficiently small, this
iteration will always converge,
ciU)
(6-11)
and the polynomial G, .(t) is obtained from (6-8).
tCl
For k = 1, 2, ..., 6, the stability region of the backward differen-
tiation method of of order k includes that portion of the left half of the
complex plane which lies between the two lines passing through the origin
and making angles ±6^ with the negative real axis, where
61 = 90° 6* = 73.5°
62 = 90° 65 = 51.8°
03 = 86.0° 66 = 17.2°
Plots of the full stability regions are shown in Gear(197l). If the problem
37
-------
being integrated is such that all the eigenvalues X of the Jacobian matrix
fall within the wedge-shaped region specified by'6^, then h\ will fall within
this region for all ht and the backward differentiation method of order k
will not limit the stepsize for reasons of stability. It is this feature
which makes the backward differentiation methods suitable for solving stiff
problems, although the occurrence of an eigenvalue close to the imaginary
axis may preclude use of the higher order methods.
The truncation error of the backward differentiation method of order k
is given by
(6-12)
where £ is some value in the interval between *£_j and t£. Thus, for k>l,
the backward differentiation method of order k is less accurate than the
same order Adams-Moulton method. For this reason, the backward differen-
tiation methods are more efficient only for stiff problems, for which Adams
methods must limit the stepsize for reasons of stability rather than for
reasons of accuracy.
Let the interpolating polynomial ^^(t) be represented by the vector
of its coefficients when it is expressed in terms of powers of (t
•s2 +
k
t-t.
where s =
These coefficients may be written as
b. ,a> =
(6-13)
where y^J approximates the j-th derivative of the solution at £«
polynomial S^ is similarly represented by the vector a = (a ,a ,...,a )
of its coefficients. With this notation, (6-8) may be written as
The
a
where I is the Pascal triangle matrix
1
1
(6-14)
0
1
2
1
1
3
3
1
1 ••• l"
4 *•• k
6 • •'• •
4 ... .
1 ... .
* •
38
-------
The formulation (6-13) and (6-14) suggest a convenient means of altering the
stepsize: to make the next step with stepsize h , it is only necessary to
scale the components of 6(i) according to
'.(i)
j = 1,
.,
k+1
(6-15)
Since (6-14) was derived under the assumption of a constant stepsize, this
procedure can be justified only if the stepsize changes are neither too large
nor too frequent. Accordingly, in the code GEAR,
1.) Stepsize increases are limited to a factor no more than 2,
and decreases to a factor not smaller than 1/10.
2.) No stepsize or order change is made until there has been a
sequence of k+1 successful steps since the last change, unless
a step failure necessitates a reduction in stepsize and/or
order.
The truncation error (6-12) may be written as
T.
X
y(k+1)(t.) + higher order terms
1
(6-16)
By ignoring the higher order terms in (6-16), the magnitude of the local
error is estimated as
with
(k+1)
(6-17)
,OO ,,(k)
k!
T+T
k!
"k+T ei ak+i
(6-18)
where the last equality follows from the observation that multiplication by
T in (6-14) does not alter the value of ^k+l^-1). In GEAR, the constants
(k+l)/(k!ak+1) are stored in the array ERTST, and the quantity
(6-19)
Tl
EPS
=
e.
i
ERTST(k)*EPS
39
-------
is calculated, where EPS is the local error tolerance specified by the user
of the code. Thus, the step succeeds if D<1 and fails if D>1.
If the rule (2) above allows a change in stepsize or order, the optimum
stepsize for continuing integration with the present order k is estimated as
h/R, where
R = 1.2(2D)
1/U+l)
(6-20)
The factor 2 is equivalent to aiming for an error of 0.5 EPS, in order to
prevent unnecessary step failures, while the overall factor 1.2 provides
an allowance for the approximations made in estimating the local error by
(6-17).
The stepsize appropriate for order k-1 is estimated as h/RDWN, where
RDWN
1.3
1.3
= 1.3
EPS
1/k
2 (k-0! ak+1(i)
EPS
1/k
ERDWN(k)*EPS
1/k
(6-21)
where the array ERDWN contains the quantities l/(k-l)!.
Raising the order is considered only if the step was successful. Then,
the stepsize appropriate for order k+1 is estimated as 'h/RUP, where
with
RUP - 1.4
,(k+2)
_L
k+2
(k+2)
EPS
I/(k+2)
k! a.
k+1
uk+2
.-e, ,)
(6-22)
where the last equality follows from use of (6-18). The constants
) are stored in the array ERUP, and RUP is calculated by
40
-------
RUP = 1.4
ERUP(k)*EPS
l/U+2)
(6-23)
The order for the next step is set to that which corresponds to the smallest
of the values R, RDWN, and RUP. The factors 1.3 and 1.4 in (6-21) and (6-23)
bias the choice, first, in favor of retaining the same order, and, second, in
favor of lowering the order. Thus, the order is raised only if there appears
to be a substantial advantage in doing so.
If the order is unchanged, the stepsize is not changed unless
JR-0.95| < 0.06. This avoids the need for reevaluating the Jacobian, as
described below, in response to an insignificant stepsize change. Otherwise,
the stepsize is changed to A/R, subject to the constraint of rule (1) above.
If the order is changed to k-1 or k+1 , the Jacobian must be reevaluated
anyway, and the new stepsize is set to fc/RDWN or A/RUP, or as limited by rule
(1). The only exception to this is when the order is lowered after a step
failure, in which case the new stepsize is not permitted to be larger than
.1 .
Lowering the order is accomplished simply by ignoring £k+i(i) in
starting the next step. When the order is raised, the quantity
(k+1)! *T* k+1
as given by (6-18), is appended to b(i) for use in starting the next step.
The integration is started by using the method of order 1, which
requires no memorized values. The initial stepsize is chosen in the same
manner as for the Adams code STEP, as described in Section 5.
Applying the formulae given above to a system of n equations (2-2) is
a fairly straightforward matter of interpreting these formulae as vector
equations, with a separate component for each of the n components of the
solution. Since (6-19), (6-21), and (6-23) must produce single numbers, the
absolute value signs in these equations are interpreted as vector norms,
according to
Z,
1*1 +11*11- ,
i-\,<
WT(fc)
where WT is a vector of nonzero weights specified by the user of the code,
as described in §7.7.
The most significant complication occurs in (6-10), where the quantity
(1+ha^f/ay) must be replaced by (I+hajjt), where Jf is the Jacobian matrix
(2-3) and I is the nxn identity matrix. In GEAR, this matrix is not actually
41
-------
inverted. Instead, the equation
(6-24)
is solved numerically for C£/m+j). The most expensive part of this solution
is performed by the subroutine DECOMP, which is taken from Forsy the ( 1967) .
DECOMP factors the matrix on the left hand side of (6-24) into a product of
upper and lower triangular matrices. After this is done, it is easy and
inexpensive to solve the equation. DECOMP performs a number of operations
which increases as n3, so that this may be the most expensive part of the
entire integration when the number of equations is large. To economize on
this expense, this matrix is reevaluated only when the iteration (6-10) will
not converge, or when the stepsize or order are changed. Up to three
iterations (6-10) are taken; the iteration is stopped when practical
convergence is indicated by
I" i(m)'1- 2«n-(k+2)
If convergence does not occur in three iterations, the matrix in (6-24) is
first reevaluated and the calculation is repeated. If this does not suffice
to produce convergence, the stepsize is reduced to A/4 for another attempt.
In GEAR, the polynomial coefficients (6-13) are in the array Y(*,*),
where the first index denotes the component of the solution and the second
denotes the order of the coefficient. GEAR is functionally organized into
five blocks:
Block 0: Checks whether the stepsize and error tolerance are too
small for the machine precision. Performs necessary
initializations the first time GEAR is called.
Block 1: Computes coefficients to correspond to the order of the
method, and scales Y(*,*) in accordance with any stepsize
change.
Block 2: Advances the solution one step, and tests whether the
error criterion is met. If the step succeeds, control
is passed to Block 4. If the step fails, control passes
to Block 3.
Block 3: Restores values to be used for repeating the failed step.
Tests whether to lower the order, reduces the stepsize,
and transfers control to Block 1 to repeat the step.
Block 4: Completes the successful step by updating those components
of Y(*,*) which were not computed in Block 2. Sets
order and stepsize appropriate for the next step, and
returns control to the calling program.
Subroutine GINTRP, when used in conjunction with GEAR, interpolates the
solution in order to provide output at points which are not mesh points of
the integration. It is used in the same way as the subroutine SINTRP for
42
-------
the Adams method code STEP, as described in Section 5.
The code STIFF is an interfacing routine designed to provide the user
with an access to the methods of GEAR which is similar in format to the
Runge-Kutta code RKFS and the Adams method code ADAM.
43
-------
SECTION 7
USING THE CODES
§7.1 Implementing the Codes
The three packages of subroutines, Runge-Kutta Programs, Adams Method
Programs, and Stiff Method Programs, are independent of one another and may
be implemented separately. The actual integration is performed by one of the
subroutines RKFS, STEP, and GEAR. The remaining programs are interfacing
routines designed to give the user a more convenient, though less flexible,
mode of access. The relationships among the various subroutines are shown
in Figure 7-1. The boxes connected by dotted lines represent named COMMON
blocks used to pass parameter values between the subroutines. The subrou-
tines denoted by OUTP, FCT, and FDER are not included in the packages, but
are subroutines supplied by the user. Complete listings of the programs are
contained in Appendices A, B, and C.
All the subroutines are written to conform to the specifications of ANSI
FORTRAN IV. The only known departure from this standard is in the form of
those DATA statements in STEP and GEAR which involve arrays. However, the
form of these DATA statements is believed to be acceptable to all current
implementations of FORTRAN IV.
The subroutines RKFS, ADAM, STEP, STFINT, GEAR, and QJACOB contain
machine dependent constants which are used to prevent attempts to compute
beyond the limits set by the machine precision. Before the programs are
compiled, the DATA statements which set these constants must be modified to
provide values appropriate for the computer being used. To avoid logical
conflicts when integrating at high precision, it is essential that these
constants be set identically for all subroutines in the same package. The
recommended values are
FOURU = 4u
U26 = 26u
REMIN = larger of 4u and 10~12
\ ,-
:V-
where u is the computer unit roundoff error, the smallest positive value such
that I.H-U is greater than 1 in floating point arithmetic. It is given by
1— s
u = 3 for truncated arithmetic
or (7-1)
" = k& for rounded arithmetic
44
-------
RKFINT
OUTP
RKDATA
RKFS
FCT
a) Runge-Kutta programs
ADMINT |
\
»
*
•
1
»
*
1 •
OUTP
ADAM
SINTRP
STDATA
STEP
FCT
b) Adams method programs
c) Stiff method programs
Figure 7-1. Organization of the Subroutine Packages.
Arrow direction is from calling program to called program. Those
blocks connected only by dotted lines represent COMMON blocks.
OUTP, FCT, and FDER represent subroutines supplied by the user.
45
-------
where 3 is the base of the floating point representation, and s is the number
of base 3 digits in the mantissa. (If 3=2, the most significant digit of the
mantissa may be not actually stored, since it is always 1 in normalized form.
It should still be counted in determining s.) If u is unknown, a sufficient-
ly accurate value may be obtained by running the following simple program.
DOUBLE PRECISION DU
U=1.
10 U=U/2.
IF (C1.+U).GT.1.) GO TO 10
U=2.*U
DU=1.DO
20 DU=DU/2.DO
IF (d.DO+DU).GT.1.DO) GO TO 20
DU=2.DO*DU
WRITE(6*100) U/DU
STOP
100 FORMAT (9H S.P. U =,E12.4/9H D.P. U'=/El2.4)
END
Variables in common blocks must be ordered so as to force proper align-
ment between the variables and word boundaries. The order used here is real
variables first, integer variables second, and logical variables last. With
most installations, this should produce correct alignment for both single
and double precision versions of the programs.
The programs as given are written in single precision. This should be
adequate for most purposes, except possibly on the IBM 360/370 machines,
whose short precision is quite poor. In the programs, all constants whose
values are critical are representable exactly in the single precision of the
IBM 360/370, and hence are exact for all machines on which these programs are
likely to be implemented. Since all arithmetic operations have been coded so
as to force double precision arithmetic whenever all the real variables are
double precision, the programs may be converted to double precision simply by
declaring all real variables double precision, and by replacing all library
functions with their double precision counterparts. In each of the programs
there are several contiguous lines of comments which will accomplish this
conversion if they are turned into program statements by removal of the C
in the first column. The machine dependent constants must also be set to the
values correct for double precision. Any conversion to double precision must
be done for all the programs in the same package, for else they cannot
interact properly.
§7.2 Specifying the Equation
The user specifies the differential equations to be solved by supplying
f subroutine, denoted by FCT in Figure 7-1, which evaluates the derivatives
(t,y) in (2-2). A call to FCT will be referred to as a function evaluation.
FCT must have the argument list (T,Y,YP), where T is a real variable and Y
and YP are real arrays each of dimension n^ whe^re n is the number of
equations. When called, FCr must compute f(t,y) for (t,y)=(T,Y), and place
these derivative values in YP. FCT must not alter the values of T and Y.
46
-------
Frequently, f will depend upon a number of parameters which the user
wishes to vary in order to observe the effect upon the solution. This is
most easily accomplished by coding these parameters as variables in FCT, and
then passing their values from the main program through a COMMON block.
As an example, the FCT subroutine for the simple ecosystem model (1-2)
might be coded as follows:
SUBROUTINE SIMPLE(T,Y,YP)
DIMENSION Y(4)*YP(4)
COMMON A14,A21,A32,A42,A43
F14=A14*Y(4)
F21=A21*Y(1)
F32=A32*Y(2)
F42=A42*Y(2)
F43=A43*Y(3)
YP(1)=F14-F21
YP(2)=F21-F32-F42
YP(3)=F32-F43
YP(4)=F42+F43-F14
RETURN
END
During the integration, many function evaluations are performed. To
make the integration as inexpensive as possible, FCT must be carefully coded
to eliminate redundant and unnecessary operations. This is particularly
important if f(t,y) is expensive to evaluate, since function evaluations will
then contribute a large portion of the integration cost.
The stiff method programs additionally require a subroutine, denoted by
FDER in Figure 7-l(c), to evaluate the Jacobian matrix Jf(t,y) defined by
(2-3). For FDER, the user may specify the subroutine QJACOB, which employs
FCT to compute the Jacobian by numerical differencing. However, unless the
form of the Jacobian is very complex,. it is preferable to supply an FDER
which evaluates the Jacobian directly. FDER must have the argument list
(T,Y,W,YP,DER), where Y and YP are arrays of dimension n, and W is a two-
dimensional array of dimensions (n,n), DER is a subroutine name which is
included in the argument list only to accommodate the requirements of QJACOB:
when FDER is called by GEAR, DER is the subroutine which the user has sup-
plied for FCT.
FDER must compute all the partial derivatives of f(t,y) at (t,y)=(T,Y)
and place them in W according to
, 2, ..., n
3yJ J = 1, 2, ..., n
When FDER is called by GEAR, YP contains the derivatives f(t,y) at (T,Y).
For some problems, these values may be useful in simplifying the evaluation
of the Jacobian. FDER must not alter the values of T, Y, and YP.
For the simple ecosystem model (1-2). the Jacobian is constant, and FDER
might be coded as follows:
47
-------
SUBROUTINE SIMPJ
COMMON A14*A21,A32*A42,A43
W(1/1)=-A21
W(2/1)=A21
W(3/-1)=0
W(4x1)=0
WC1,2)=0
W(2,2)=-A32-A42
W(3/2)=A32
W(4/-2)=A42
W(1,3>=0
W(2/3)=0
W(3/3)=-A43
W(4/3)=A43
W(1,4)=A14
W<2,4)=0
W(3,4)=0
W(4,4)=-A14
RETURN
END
§ 7 . 3 The Local Error Tolerance
Except for direct use of STEP or GEAR, whose requirements are described
in §7.7 below, the user specifies the local error test by supplying nonneg-
ative values for two parameters: a relative error tolerance Br and an
absolute error tolerance Ea. In the programs, these parameters are denoted
by RELERR and ABSERR, respectively. If y(fc) is the current value of the fcth
component of the solution, the local error tolerance for that component is
computed as
+ E& (7-2)
and the codes require that
& = 1, 2, ..., n (7-3)
for a successful step, where T|(£) is the magnitude of the estimated local
error in the &th component of the solution.
At least one of ET and E& must be nonzero, for else it will not be pos-
sible to satisfy (7-3). For a pure relative error test (#a=0 and £r>0) , £(&)
will vanish if y(&)=0. The Runge-Kutta programs compute |y(£)| as the
average magnitude of the solution at the beginning and end of the step, so
that this difficulty is unlikely to arise. However, the Adams and stiff
method programs use for y(&) the solution value at the beginning of the step,
so that a pure relative error test -cannot be used for starting a solution
which has any component with vanishing initial value.
48
-------
To avoid attempts at impossible accuracy, the programs also require that
> 4u|yU)| A = i, 2, ..., n (7-4)
where u is the computer unit roundoff error (7-1). For the Runge-Kutta
programs, (7-4) is automatically satisfied by requiring that RELERR be
greater than or equal to REMIN. The Adams and stiff method programs do allow
a pure absolute error test (RELERR=0, ABSERtyO), but a check is made at each
step to ensure that (7-4) is satisfied.
In choosing the error tolerances, the user should consider the accuracy
of the model which he is investigating. It is pointless to request that
errors in the integration be very much smaller than the empirical uncertain-
ties in the model itself. Doing so will simply increase the cost of solution
without producing any additional usable information.
The user must be fully aware that his specified error tolerances are
used only to control the local error in each step of the integration. While
it is true for most problems that the resulting global error is the same
order of magnitude as the local error tolerances, this cannot be guaranteed,
for the reasons discussed in §2.3 above. A reliable estimate of the global
error may be obtained by repeating the integration with the error tolerances
reduced by, say, a factor of 1/10. The difference of the two solutions thus
obtained may then be taken as an estimate of the global error in the less
accurate solution.
§7.4 Choice of Method and Choice of Program
The best method for a particular problem is that which computes the
solution most efficiently, i.e., with the least possible expense. Factors
which affect the efficiency of the codes are the degree of stiffness of the
problem, the expense of the function evaluations, and the stringency of the
error tolerances.
The stiff method programs are designed specifically for the solution of
stiff problems. Although they will integrate nonstiff problems, they will do
so less efficiently than the Runge-Kutta or Adams methods. Normally they are
used only when the problem is already known to be stiff, or after stiffness
is suspected as the cause of an expensive integration with the Runge-Kutta or
Adams programs.
Unless stiffness is fairly severe, the choice is between the Runge-Kutta
and Adams methods. The Runge-Kutta methods are best when the function eval-
uations are cheap and only moderate accuracy is required. The Adams methods
are preferred when function evaluations are expensive or when high accuracy
is required.
The expense of a function evaluation is judged by the number of oper-
ations per component required to evaluate f(t,y). If, as in the subroutine
SIMPLE in §7.2 above, each component is computed by only a few arithmetic
operations, the function evaluation is classified as cheap. On the other
hand, if the evaluation of each component requires, on the average, many
arithmetic operations, or the evaluation of a transcendental function, the
49
-------
function evaluation must be regarded as expensive. The effect of error
tolerance on the relative efficiencies of the Ruhge-Kutta and Adams methods
is illustrated in Tables 8-2 and 8-3 in Section 8.
If, as is often the case, the same problem is to be integrated repeated-
ly with only minor adjustment of parameters in the definition of ?(t,y), it
may be worthwhile to perform a trial solution with the different methods in
order to determine which is best for the production run.
If the number n of equations is very large, storage requirements may
constrain choice of a method. To hold the solution and intermediate results,
the programs require the following:
Runge-Kutta programs: 8n floating point locations
Adams method programs: 2Qn floating point locations
Stiff method programs: \8n+n2 floating point locations
and n fixed point (integer) locations
Also, while the amount of computation increases linearly with n for
Runge-Kutta and Adams methods, the amount of computation required by the
stiff methods has a term which increases as n3. Thus, when n is large, a
moderately stiff problem may be handled more efficiently by the Runge-Kutta
or Adams methods than by the stiff method programs. Problems which are both
very large and very stiff cannot be integrated by any of these codes. Such
problems are tractable only when the Jacobian matrix is sparse (i.e., contains
few nonzero elements), in which case a specially coded integration routine
may be used to economize on storage requirements and on the cost of matrix
inversion.
Once the basic method has been selected, the choice of a particular
program depends on the amount of control which the user must have over the
integration process. The easiest programs to use are RKFINT, ADMINT, and
STFINT. A single call to one of';these subroutines calculates the solution
over the entire interval of integration, using the mixed relative-absolute
error test described in §7.3. Output of the solution at equally spaced points
in the interval is accomplished by calls to a user-supplied output subroutine.
When an error is encountered, integration is continued if it is possible to
take corrective action. If an error cannot be corrected, execution of the
program is terminated. Informative and diagnostic messages are written into
a data set specified by the user.
Greater flexibility in the choice of output points and in the handling
of error conditions may be achieved through direct use of RKFS, ADAM, or
STIFF. A call to one of these programs advances the solution to the output
point specified by the user, using a mixed relative-absolute error tolerance.
Return to the calling program occurs when the output point is reached, or
sooner, if an error condition occurs. The complete solution requires a
sequence of calls to the integration routine, each call carrying the inte-
gration from one output point to the next. Using these programs is consider-
ably more demanding of the user, since all provisions for output and error
handling must be written into the calling program.
For Adams and stiff methods, the greatest control over the integration
can be achieved by direct use of STEP or GEAR. Each of these subroutines
50
-------
returns control to the calling program after advancing the solution over one
step. They also allow specification of a different error test for each com-
ponent of the solution. For the Runge-Kutta method, similar step-by-step
monitoring of the solution can be achieved by using RKFS in single-step mode.
In this case, however, the error test is limited to the simple relative-
-absolute specification.
§7.5 Using RKFINT, ADMINT, and STFINT
A single call to RKFINT, ADMINT, or STFINT will produce a complete
integration of the user's problem. The user specifies the beginning and end
of^the integration interval by TINIT and TFINAL, and the spacing of output
points by TINCR. To provide the output, the program calls a user-supplied
output subroutine, denoted by OUTP in Figure 7-1, at intervals of TINCR in the
independent variable. The OUTP subroutine must have an argument list of the
for (T,Y,YP), where T is a real variable, and Y and YP are real arrays each
of dimension NEQN, where NEQN is the number of equations being integrated.
When OUTP is called, T contains the current value of the independent variable,
Y contains the solution at T, and YP contains the derivatives of the solution
at T. (YP is actually the working storage array W described below, but only
the first NEQN locations of this array, which contain the derivative values,
are likely to be of interest to the user.)
As an example, the following output subroutine may be used to print out
the solution and derivative values for the problem specified by the subrou-
tine SIMPLE in §7.2. For this problem, the sum of the components, Y(1)+Y(2)
+Y(3)+Y(4), should be constant. As a check against possible coding errors,
and as a rough indication of the solution's accuracy, this output routine
also computes and prints out this sum.
SUBROUTINE SIMOUT(T,Y,YP)
DIMENSION Y(4),YP(4)
SUM=Y(1)+Y(2)+Y(3)+Y(4>
WRITE(6x10) T/Y/SUM/YP
RETURN
10 FORMATCIH ,F7.2/5G12.4/5X,4G12.4>
END
All information necessary for integration by RKFINT, ADMINT, and STFINT
is passed to these programs through their call lists, which are
RKFINT(NEQN,Y,W,TINIT,TFINAL,TINCR,RELERR,ABSERR,IUNIT,FCr,Ot/rP)
ADMINT (NEQN, Y ,W, TINIT, TFINAL, TINCR, RELERR, ABSERR, IUNIT, FCT, OUTP )
STFINT(NEQN,Y,W,IW, TINIT, TFINAL, TINCR,RELERR,ABSERR,IUNIT,FCr,FDH?,Oy:rP)
The quantities in the call list have the following meanings:
NEQN the number of equations
Y a real array of dimension NEQN, into which the user places the
initial values of the solution
W a real array used for working storage. Its minimum dimensions
must be 7*NEQN for RKFINT, 19*NEQN for ADMINT, and 17*NEQN+NEQN**2
51
-------
for STFINT
IW (used by STFINT only) an integer array of dimension NEQN, used
for working storage
TINIT the starting value of the independent variable
TFINAL the value of the independent variable at the end of the
integration interval
TINCR the spacing between the output points
RELERR the relative error tolerance for the local error test
ABSERR the absolute error tolerance for the local error test
IUNIT the logical unit number of the data set to be used for writing
informative and diagnostic messages. If IUNIT is less than or
equal to zero, no messages will be written.
FCT the name of the user-supplied subroutine to be used for calcu-
lating the derivatives, as described in §7.2
FDER (used by STFINT only) the name of the user-supplied subroutine to
be used for calculating the Jacobian matrix, as described in §7.2
OUTB the name of the user-supplied subroutine used to handle output
To use one of these programs, the user must
1. Allocate storage for the arrays Y and W (and also for IW, if STFINT
is being used) in a DIMENSION statement.
/
2. Declare the actual names of FCT and OUTP (and. also of FDER, if
STFINT is being used) in an EXTERNAL statement, and place those names
in the call list.
3. Place the initial values of the solution in Y.
4. Provide appropriate values for NEQN, TINIT, TFINAL, TINCR, RELERR,
ABSERR, and IUNIT. Since these quantities are not altered by the
programs, they may be specified by writing constants in the call list,
As an example, the following main program uses RKFINT to integrate the
simple ecosystem model (1-2) over the interval CO,1003 with output at steps
of 1 and error tolerances RELERR=ABSERR=10~I*.
DIMENSION Y(4)/W(28)
EXTERNAL SIMPLE,SIMOUT
COMMON A14/A21/A32/A42/A43
C GET VALUES OF PARAMETERS
READ(5,10) A14'A21,A32,A42/A43
C READ INITIAL VALUES
READ(5/20) Y
C WRITE TABLE HEADINGS FOR OUTPUT
WRITE(6/30)
C INTEGRATE. MESSAGES TO BE WRITTEN TO UNIT 4.
CALL RKFINT(4/-Y/W/0.,100.,1 ./1 .E-4/1 .E-4/4/SIMPLE/SIMOUT)
STOP
52
-------
10 FORMAT(5F12.2)
20 FORMATC4F15.3)
30 FORMAT(3X,4HTIME,5X,4HYC1),12X,4HY(2),12X*4HY(3),12Xx3HSUM/
1 27X/11HDERIVATIVES)
END
To integrate with ADMINT, it is only necessary to change the DIMENSION
statement to
DIMENSION Y(4)/W(76)
and to replace "RKFINT" with "ADMINT" in the CALL statement.
To integrate with STFINT, the DIMENSION, EXTERNAL, and CALL statements
must be changed to
DIMENSION Y(4)*W(84),IW(4)
EXTERNAL SIMPLE/-SIMOUT/SIMPJ
CALL STFINT (4*Y,W/IW,0./100./'1 .^1 .E-4x1 .E-4/4/SIMPLE/>SIMPJ,SIMOUT)
Integration may be done in the reverse direction; i.e., TFINAL < TINIT
is permitted. In this connection, it should be noted that the reverse inte-
gration of a stiff problem is a highly unstable problem, and vice versa,
The error messages written into the data set specified by IUNIT are self
explanatory. Each message contains the value of the independent variable at
which the error occurred, and states what action was taken. As an example,
a request for too much accuracy is handled by increasing RELERR and ABSERR to
achievable values, printing these new values in the error message, and then
continuing integration with the increased tolerances.
The OUTP subroutine need not be limited to merely writing out the
solution. Depending on the user's needs, it may be written to perform calcu-
lations requiring the solution values as input, to prepare the output for
graphical display, to compare the solution with expected values, or, in an
interactive environment, to force a pause in integration to give the user
time to assess the quality of the results already obtained. A limited amount
of control over the integration can be achieved by using OUTP to manipulate
some of the values in the COMMON blocks shown in Figure 7-1. For example, if
the problem is known to be very expensive to integrate, error messages for
"too much work" may be suppressed by setting the function evaluation counter
NFE to a large negative value the first time OUTP is called, which is at
T = TINIT. NFE is in the COMMON block RKDATA, STDATA, or GEDATA, for Runge-
-Kutta, Adams, or stiff methods, respectively. When using RKFINT, the test
for stiffness may be suppressed by referring to the COMMON block RKDATA and
setting KST > 25 the first time OUTP is called. This will reduce the over-
head of integration slightly, but any report of "too much work" will then
give a false indication of stiffness.
§7.6 Using RKFS, ADAM, and STIFF
If unequally spaced output points are needed, or a different handling of
errors is required, either the user may write a modified driver program on
the model of RKFINT, ADMINT, or STFINT, or he may use RKFS, ADAM, or STIFF
53
-------
directly. Each of these three programs returns control to the calling
program after advancing the solution from one output point to the next.
To start the integration, the user supplies the initial conditions (T,Y),
sets TOUT to the value of the first output point, and sets IFLAG=1 to
indicate that this is the beginning of the problem. Assuming that the inte-
gration is successful, the subroutine returns with T = TOUT, with Y the
solution at TOUT, and with IFLAG=2. To continue to the next output point,
it is only necessary to set the new value of TOUT and call the subroutine
again. This process is continued until the entire solution is generated.
Although RKFS, ADAM, and STIFF are similar in their requirements for
use, their call lists are quite different. Full instructions for using each
of these programs are contained as comments in the program listings in the
Appendices. Values of IFLAG other than 2 are used to inform the user of
errors which occur during the integration: the significance of each value of
IFLAG is explained in the comments to the appropriate program.
RKFS may be used in single-step mode by starting the integration with
IFLAG=-1 rather than IFLAG=1. In single-step mode, RKFS returns to the
calling program after each step of the integration with IFLAG=-2. Calling
RKFS again then advances the integration another step, and so forth. Arrival
at TOUT is signaled by IFLAG=2. To continue in single-step mode to a new
TOUT then requires the user to change IFLAG from 2 to -2.
For reasons of efficiency, ADAM and STIFF integrate beyond TOUT inter-
nally and interpolate the solution to TOUT by using the subroutines SINTRP
and GINTRP. For some problems, it may not be possible to integrate beyond
TOUT because the derivatives are not defined there. Negative values of
IFLAG are used to signal ADAM and STIFF that the integration cannot be
carried past TOUT. This option should never be used unless there is a
genuine difficulty at TOUT, since any subsequent continuation to a new TOUT
then requires restarting the integration with a first order method, so that
much efficiency will be lost.
§7.7 Using STEP and GEAR
STEP and GEAR are the step-by-step integration routines for the Adams
and STIFF methods, respectively. Full instructions for their use are given
as comments in the program listings in Appendices B and C.
These two codes provide for an error tolerance specification which is
much more flexible than the relative-absolute test described in §7.3. In
STEP and GEAR, the user supplies a basic error tolerance EPS, and specifies
how this tolerance is to be applied to each component of the solution by
giving a vector of nonzero weights WT(L), L=l,2,...,NEQN. For example,
WT(I) = A
specifies an absolute error test of A*EPS for the I-th component of the
solution,
WT(I) = ABS(Y(I))
54
-------
specifies an error test of EPS relative to the most recent value of Y(I), and
WT(I) «- AMAX1(WT(I),ABS(Y(I)))
specifies an error test of EPS relative to the largest magnitude so far
achieved by the I-th component. This last example is effectively a relative
error test if the I-th component is increasing, and an absolute error test
if the I-th component is decreasing.
Rather than manipulating the stepsize to cause exact arrival at an
output point, STEP and GEAR are best used in conjunction with the interpo-
lation routines SINTRP and GINTRP. In general, any alteration of the step-
size automatically chosen by STEP or GEAR will cause inefficient operation of
the program, and will, consequently, cause the integration to be more expen-
sive than necessary.
STEP and GEAR do not check for errors in input other than requests for
a stepsize or an error tolerance which is too small for the precision of the
machine. Thus, the user must be quite cautious when using these programs
directly. An examination of the listings of ADAM and STIFF may be a useful
aid to understanding the necessary precautions.
§7.8 Miscellanea
The COMMON blocks shown in Figure 7-1 are used for two purposes: (1) to
transmit values from and to return values to the calling program, and (2) to
hold values needed for use on subsequent calls to the subroutine. The
use of a COMMON block to transmit arguments to a subroutine is usually more
efficient than placing the values in the call list, particularly when the
number of variables is large. If, however, the use of these COMMON block is
inconvenient, it is a simple matter to modify the programs to place these
variables in the call lists.
Since no values are retained internally between calls to the subroutines,
these programs are easily used in an overlay structure if care is taken to
preserve the COMMON blocks from being overwritten. If such a use is not
contemplated, then the programs may be modified to allow internal retention
of values. Although this is not standard ANSI FORTRAN, it is a commonly used
device which allows call lists and interfacing requirements to be greatly
simplified.
It is sometimes necessary to integrate problems (2-2) for which f(t,y)
has discontinuities in t. While the integration routines given here can
usually pass such discontinuities, they can do so only by locating the point
of discontinuity through a succession of step failures. It is more efficient
to treat such a problem as two distinct problems: the first problem is
integrated up to the point of discontinuity, and the solution there is used
as initial condition for the second problem, which provides the solution
beyond the discontinuity.
55
-------
SECTION 8
EXAMPLES
This section presents a few simple examples to illustrate how the
relative efficiencies of the codes vary with the characteristics of the
problem and with the accuracy requested for the integration. These examples
were run on the IBM 370 at the Washington Computer Center, using double
precision versions of the codes. The integration routines were compiled with
the IBM FORTRAN H compiler, OFT=2. Main programs, function subroutines, and
output subroutines were compiled with the FORTRAN G compiler, which is faster
but produces less efficient object code. The timings quoted are for inte-
gration and output only, and do not include the time needed for compilation
and linkage of the programs. These timings should not be regarded as highly
accurate, since they are somewhat dependent on the number of tasks being
handled by the operating system: differences in the least significant digit
quoted are not necessarily meaningful.
The first example illustrates the response of the codes to increasing
stiffness by integrating the problem
y' = -X(y-t2) + 2t
(8-1)
y(0) = 0
This problem has the exact solution y(t) = t2, regardless of the value of X.
However, the problem becomes increasingly stiff for increasingly large
positive values of X. Table 8-1 shows the performance of the codes RKFINT,
ADMINT, and STFINT in integrating this problem on the interval CO,503 with
output at steps of 1, and error tolerances of RELERR=ABSERR=10~5. The error
quoted is the maximum relative error which was found at the output points
t = 1, 2 50.
In this example, the behavior of the errors is unusual because of the
simple form of the solution. For X = 0, the solution curves of the differ-
ential equation in (8-1) are of the form u(t) = t2 +constant. An inte-
gration method of order greater than or equal to 2 can follow these curves
exactly. Since RKFINT uses a method of order 5, and ADMINT uses a method of
order at least 2 (because of local extrapolation), the errors in these cases
come only from conputer roundoff errors. STFINT, however, starts with a
method of order 1, which is not exact for this problem, so that the error
tolerance permits an error of about 10~5. For X > 0, the stability of the
problem (8-1) causes the starting error in STFINT to be damped out after the
order of the method is raised. When X is large, the error which remains at
56
-------
X = 0
X = 1
X = 10
X = 102
X = 103
X = 10*
RKFINT
301
ixlO-15
.020
461
.4xlO~5
.038
1625
.5xlO-5
.107
9573
.5xlO"5
.58
82745
.5xlO-5
4.86
not
attempted
ADMIKT
42
7xlQ-16
.015
55
HxlO~5
.018
687
9xlO-5
.135
6402
lOxlO-5
1.03
64849
lOxlO-5
10.89
not
attempted
STPINT
47(15)
IxlO-5
.019
49(15)
4xlO-6
.018
48(15)
8xlQ-8
.021
52(15)
4xlO~13
.02
54(16)
3xlO-16
.02
51(16)
4xlO~16
.02
NF(NJ)
error
time
NF(NJ)
error
time
NF(NJ)
error
time
NF(NJ)
error
time
NF(NJ)
error
time
NF(NJ)
error
time
Table 8-1. Performance of the codes in integrating problem (8-1),
illustrating the effect of increasing stiffness. The problem was
integrated on the interval E0,503 with output at steps of 1, and
error tolerances RELERR=ABSERR=10~5. NF is the number of function
evaluations, and NJ the number of Jacobian evaluations, needed for
the integration. Error is the maximum relative error actually
occurring at any of the output points, and time is the integration
time in seconds.
57
-------
t = 1 will be much smaller than 10~5. For RKFINT and ADMINT, errors are not
damped out in this fashion because the methods are operating near the
boundaries of their stability regions.
Since RKFINT requires that output points be mesh points of the inte-
gration, the minimum cost of solving (8-1) with 50 output points is
50 steps x 6 faction evaluations + , f£nal evaluation
step
= 301 function evaluations
Thus, the cost of using RKFINT for X = 0 and X = 1 is largely determined by
the number of output points: the integration of these cases would be
significantly cheaper if fewer output points were requested.
The second example is the system of three differential equations
The constants were set to u = 0.04, k = 200, Z = 0.5, and the problem was
integrated on the interval CO,503, with output at steps of 5, starting with
RKFINT
ADMINT
STFINT
e = 10-"
73
.009
38
.016
60(12)
.032
e = 10~6
91
.012
52
.020
89(15)
.044
e = 10~8
199
.022
78
.027
140(20)
.065
e = 10-10
433
.047
104
.036
209(28)
.096
e = ID"12
1051
.110
128
.046
338(26)
.142
NF
time
NF
time
NF(NJ)
time
Table 8-2. Illustrating the effect of error tolerance on the cost
of integrating the problem (8-2), with constants and initial
conditions as given in the text. Integration was on the interval
110,503 with output at steps of 5, and error tolerances of
RELERR=ABSERR=e. NF and NJ are the number of function evaluations
and the number of Jacobian evaluations, respectively. Integration
times are in seconds.
58
-------
the initial conditions yi(0) = 10, y2(0) = 1000, and y3(0) = 0. This problem
is not stiff, and the function evaluations are cheap. Thus, the most
efficient integrator is expected to be RKFINT for moderate tolerances and
ADMINT for more stringent tolerances, while STFINT should not be competitive.
These expectations are verified by the timings shown in Table 8-2.
The final example is the problem
- by z + (-
(-a-i+Oe (8-3)
yi(0) = 2 y2(0) = 1
for which the exact solution is
yi(t) = e cos bt + e
at t
y2(t) = e sin-bt + e
For this problem, the function evaluations must be regarded as moderately
expensive because of the need to evaluate the exponential function e . The
problem was integrated over the interval CO, 103 with output at steps of 0.5 .
Table 8-3 presents the results of the integration for a =» -1, b = 2
using pure relative error tolerances RELERR = e. This problem is not stiff.
Because of the cost of the function evaluations, ADMINT performs more
efficiently than RKFINT, except possibly for e = 10""*. On this machine, the
smallest tolerances permitted are e = 10~12 for RKFINT and 8.89xlO~16 for
ADMINT and STFINT. The errors quoted are the maximum relative errors which
occur at any of the output points. Observe that ADMINT is somewhat more
successful than STFINT in achieving accuracy near limiting precision.
Table 8-4 shows the results of integrating (8-3) for some different
values of a and Jb. These problems are all stiff to some extent, and STFINT
generally performs best. The case a = -20 and b = 70 provides an excellent
example of the fact that the degree of stiffness depends not only on the
problem, but also on the amount of accuracy which is sought: for c - 10~"*,
STFINT is clearly superior; but at e = 10~8, ADMINT appears to have a slight
advantage .
59
-------
e
10-*
1
-------
e
10-*
10"6
io~8
10"*
io-6
ID'8
a = -20 b - 70
RKFINT
1686
.59x10"*
.22
4212
.29xlO~6
.53
10439
.50xlO~8
1.31
a =
2093
.24x10"*
.30
4901
.20xlO~6
.64
11615
.43xlO~8
1.48
ADMINT
1298
2.7x10-*
.33
1482
5.6xlO"6
.36
1937
1.3xlO"8
.47
= -100 b --
1384
4.4x10"*
.37
1540
6.9xlO"6
.40
1787
4.2xlO~8
.48
STFINT
344(28)
.26x10-"
.14
766(37)
.59xlO"6
.30
1571(63)
1.2xlO"8
.56
= 0
206(19)
.08x10-*
.10
319(25)
.OSxlO"6
.16
599(33)
.OSxlO"8
.26
a = -50 b - 50
RKFINT
1734
. SSxlO"1*
.24
4188
.18xlO"6
.59
10086
.25xlQ-8
1.25
a =
6540
.17x10"*
.77
10565
.34xlQ-6
1.20
24863
.38xlO"8
2.94
ADMINT
1208
8.3x10"*
.30
1368
5.8xlO"6
.36
1642
l.SxlO"8
.41
= -200 b --
3649
4.0x10-*
.97
3705
2.2xlO"6
.96
4069
4.2x10-"
1.14
STFINT
223(17)
.07x10-*
.13
420(27)
.OSxlO"6
.17
802(50)
.24xlQ-8
.30
= 100
236(24)
.02x10"*
.10
439(30)
.03xlO"6
.17
665(36)
.OSxlO"8
.25
NF(NJ)
error
time
NF(NJ)
error
time
NF(NJ)
error
time
NF(NJ)
error
time
NF(NJ)
error
time
NF(NJ)
error
time
Table 8-4. Performance of the codes in integrating problem (8-3) with the parameters a,Jb as given.
-------
REFERENCES
Ames, W. F. (1969) Numerical Methods for Partial Differential Equations,
Barnes and Noble, New York
Bettis, D. G. (1978) Embedded Runge-Kutta Methods of Order Four and Five.
TICOM Report 78-2, The Texas Institute for Computational Mechanics,
University of Texas, Austin
Fehlberg, E. (1969) Low-Order Classical Runge-Kutta Formulas With Stepsize
Control and their Application to some Heat Transfer Problems.
NASA Technical Report NASA TR R-315
Forsythe, G.,and C. Moler (1967) Computer Solution of Linear Algebraic
Systems, Prentice-Hall, Englewood Cliffs, New Jersey
Gear, C. W. (1971) Numerical Initial Value Problems in Ordinary Differential
Equations, Prentice-Hall, Englewood Cliffs, New Jersey
Gear, C. W. (1971a) The Automatic Integration of Ordinary Differential
Equations. Communications of the ACM 14:176-179
Algorithm 407: DIFSUB for Solution of Ordinary Differential
Equations. Communications of the ACM 14:185-190
Krogh, F. T. (1969) VODQ/SVDQ/DVDQ—Variable Order Integrators for the
Numerical Solution of Ordinary Differential Equations.
TU Doc. CP-2308, NPO-11643, Jet Propulsion Laboratory, Pasadena
Lassiter, R. R., J. L. Malanchuk, and G. L. Baughman (1976) Comparison of
Processes Determining the Fate of Mercury in Aquatic Ecosystems.
Environmental Modeling and Simulation, U. S. Environmental
Protection Agency, Washington, pp. 619-623
Patten, B. C., ed. Systems Analysis and Simulation in Ecology, Academic
Press, New York
(1971) vol. 1
(1972) vol. 2
(1975) vol. 3
(1976) vol. 4
Schryer, N. L. (1975) A User's Guide to DODES, a Double Precision Ordinary
Differential Equation Solver. Bell Laboratories, Murray Hill,
New Jersey
62
-------
Shampine, L. F., and M. K. Gordon (1975) Computer Solution of Ordinary
Differential Equations: The Initial Value Probleme W. H. Freeman,
San Francisco
Shampine, L. F., and H. A. Watts (1976) Practical Solution of Ordinary
Differential Equations by Runge-Kutta Methods. Sandia Laboratories
Report SAND76-0585, Albuquerque, New Mexico
Shampine, L. F., and K. L. Hiebert (1977) Detecting Stiffness with the
Fehlberg (4,5) Formulas. Comp. and Maths, with Appls. 3:41-46,
Pergamon Press
63
-------
APPENDIX A
RUNGE-KUTTA PROGRAMS
This appendix contains the complete listings of the subroutines RKFINT
and RKFS.
SUBROUTINE RKFINTCNEQN/Y/W/TINIT/TFINAL/TINCR/RELERR/
1 ABSERR/IUNIT/FCT/OUTP)
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
SUBROUTINE RKFINT USES RKFS TO INTEGRATE A SYSTEM OF NEQN
FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS
DY/DT = FCT/Y)
Y =
-------
OUTP(T/Y/W> USED TO HANDLE OUTPUT.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c****************************************************************c
C
C
C
C
C
C
C
C
C
C
C
TO USE RKFINT/ THE CALLING PROGRAM MUST
1. ALLOCATE SPACE FOR THE ARRAYS Y(*> AND W(*) IN A
DIMENSION STATEMENT. (NOTE THAT Y MUST BE AN ARRAY/
EVEN IF NEQN=1.)
2. SET Y(1)/Y(2)/.../Y(NEQN) TO THE INITIAL VALUES TO BE
USED FOR THE SOLUTION.
3. PROVIDE APPROPRIATE VALUES FOR NEQN/ TINIT/ TFINAL/
TINCR/ RELERR/ ABSERR/ AND IUNIT. (SINCE THESE
QUANTITIES ARE NOT ALTERED BY RKFINT/ THEY MAY/ IF
DESIRED/ BE WRITTEN AS CONSTANTS IN THE CALL LIST.)
4. PLACE THE ACTUAL NAMES OF FCT AND OUTP IN THE CALL
LIST/ AND DECLARE THOSE NAMES IN AN EXTERNAL STATEMENT.
AT EACH CALL TO OUTP(T/Y/W)/ THE CURRENT VALUES OF THE
DERIVATIVES DY(I)/DT ARE IN W(1)/W(2)/.../W(NEQN). THE
FIRST CALL TO OUTP IS MADE AT TINIT.
IT IS NOT NECESSARY THAT THE INTEGRATION INTERVAL BE EVENLY
DIVISIBLE BY THE OUTPUT MESH SIZE TINCR. HOWEVER/ IF IT IS
NOT/ INTEGRATION WILL PROCEED TO THE FIRST MESH POINT
BEYOND TFINAL/ RATHER THAN STOPPING AT TFINAL.
UNDER CERTAIN CIRCUMSTANCES/ THE OUTPUT ROUTINE OUTP MAY
BE USED TO ALTER SOME OF THE PARAMETERS IN THE COMMON
BLOCK RKDATA. SEE THE PROGRAM DOCUMENTATION FOR A
DISCUSSION OF THIS.
THIS SUBROUTINE IS WRITTEN IN SINGLE PRECISION.
A DOUBLE PRECISION VERSION FOR USE WITH A DOUBLE PRECISION
VERSION OF RKFS MAY BE OBTAINED SIMPLY BY REMOVING THE C
FROM COLUMN 1 OF THE NEXT FIVE CARDS.
DOUBLE PRECISION ABS/ABSER/ABSERR/DABS/DSIGN/DT/D1/D2/
1 EPSDT/H/RELER/RELERR/SIGN/T/TFINAL/TINCR/
2 TINIT/TMAXIN/TOUT/W/Y
ABS(D1) = DABS(D1)
SIGN(D1/D2) = DSIGN(D1/D2)
COMMON /RKDATA/ T/TOUT/RELER/ABSER/H/NEQ/IFLAG/NFE/KOP/
1 INIT/NSS/KST
DIMENSION Y(NEQN)/W(1)
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
13700
13800
13900
14000
14100
14200
14300
14400
14500
14600
14700
14800
14900
15000
15100
15200
15300
15400
15500
15600
15700
15800
15900
16000
16100
16200
16300
16400
16500
16600
16700
16800
16900
17000
17100
17200
17300
17400
17500
17600
17700
17800
17900
18000
18100
18200
18300
18400
18500
18600
18700
18800
65
-------
C THE ACTUAL DIMENSION OF W, AS DEFINED BY THE CALLING PROGRAM/ C
C MUST BE AT LEAST 7*NEQN.
C
C
EXTERNAL FCT
C
C WRITE STARTING DATA INTO MESSAGE FILE
C
IF (IUNIT.GT.O) WRITE(IUNIT,100) NEQN/TINIT/TFINAL/TINCR/
, 1 RELERRsABSERR
C
C COMPUTE INDICES FOR SPLITTING THE WORK ARRAY W(*>.
C
K1=NEQN+1
K2=K1+NEQN
K3=K2+NEQN
K4=K3+NEQN
K5=K4+NEQN
K6=K5+NEQN
C
C INITIALIZE VARIABLES TO START INTEGRATION WITH RKFS.
C
NEQ=NEQN
RELER=RELERR
ABSER=ABSERR
IFLAG=1
T=TINIT
TOUT=TINIT
C
C SAVE INFORMATION NEEDED TO LOCATE END OF INTEGRATION.
C
TMAXIN=ABS(TFINAL-TINIT)
EPSDT=0.001*ABS(TINCR)
C
C SET SIGN OF OUTPUT INTERVAL TO PROCEED FROM TINIT TO TFINAL.
C
DT=SIGN(TINCR/TFINAL-TINIT)
C
C INITIALIZE COUNTER OF ERROR RETURNS.
C
IER=0
C
C CALL RKFS TO INTEGRATE TO THE NEXT OUTPUT POINT.
C
5 CALL RKFS(FCT,Y,W(1)/W(K1),WCK2),WCK3),W(K4),W
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
CHECK WHETHER TFINAL HAS BEEN REACHED OR PASSED.
IF (ABS(T-TINIT).GE.TMAXIN) GO TO 95
INTEGRATION IS NOT YET COMPLETE/ SO THE NEXT OUTPUT POINT IS
SET. IF THE NEW POINT DIFFERS FROM TFINAL BY AN AMOUNT
ATTRIBUTABLE TO ROUNDOFF ERROR/ IT IS RESET TO TFINAL EXACTLY,
INTEGRATION IS THEN CONTINUED.
TOUT=T+DT
IF (ABS(TFINAL-TOUT).LT.EPSDT) TOUT=TFINAL
GO TO 5
INTEGRATION WAS INTERRUPTED BEFORE REACHING THE OUTPUT POINT.
INCREMENT THE ERROR COUNTER AND BRANCH FOR WRITING THE
APPROPRIATE DIAGNOSTIC MESSAGE. INTEGRATION WILL THEN BE
CONTINUED IF THE ERROR IS NOT TOO SEVERE.
8 IER=IER+1
GO TO (20/20/30/40/50/60/70/80/90)/IFLAG
IFLAG HAS AN ILLEGAL VALUE/ SO THAT A PRECISE DEFINITION OF
THE ERROR IS NOT POSSIBLE/ AND EXECUTION MUST BE TERMINATED.
20 IF (IUNIT.GT.O) WRITECIUNIT/200) T/IFLAG
STOP
IFLAG = 3 -- THE REQUESTED RELATIVE ERROR TOLERANCE WAS TOO
SMALL FOR THE MACHINE PRECISION AND HAS BEEN INCREASED TO
THE MINIMUM PERMISSIBLE VALUE.
30 IF (IUNIT.GT.O) WRITEUUNIT/300) T/RELER
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 4 — THE REQUESTED ERROR TOLERANCES WERE TOO
STRINGENT FOR THIS PROBLEM AND HAVE BEEN INCREASED TO
VALUES WHICH SHOULD BE ACCEPTABLE.
40 IF (IUNIT.GT.O) WRITE(IUNIT/400) T/RELER/ABSER
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 5 — THE COUNTER OF CALLS TO FCT WAS RESET AFTER
EXCEEDING THE LIMIT IMPOSED IN RKFS.
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
24100
24200
24300
24400
24500
24600
24700
24800
24900
25000
25100
25200
25300
25400
25500
25600
25700
25800
25900
26000
26100
26200
26300
26400
26500
26600
26700
26800
26900
27000
27100
27200
27300
27400
27500
27600
27700
27800
27900
28000
28100
28200
28300
28400
28500
28600
28700
28800
28900
29000
29100
29200
67
-------
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
50 IF (IUNIT.GT.O) WRITE(IUNIT/500) T
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 6 — THE COUNTER OF CALLS TO FCT WAS RESET AFTER
STIFFNESS CAUSED IT TO EXCEED THE LIMIT IMPOSED IN RKFS.
60 IF (IUNIT.GT.O) WRITE(IUNIT/600) T
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 7 — THE MONITOR OF OUTPUT FREQUENCY WAS RESET
AFTER INDICATING THAT THE OUTPUT INTERVAL IS TOO SMALL
TO ALLOW EFFICIENT OPERATION OF RKFS.
70 IF (IUNIT.GT.O) WRITE(IUNIT/700) T
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 8 — FATAL ERROR — A COMPONENT OF THE SOLUTION
VANISHED/ MAKING A PURE RELATIVE ERROR TEST IMPOSSIBLE.
THE USER MUST SUPPLY AN APPROPRIATE NONZERO VALUE OF
ABSERR FOR SUCCESSFUL INTEGRATION.
80 IF (IUNIT.GT.O) WRITE(IUNIT/800) T
STOP
IFLAG = 9 — FATAL ERROR — AN ILLEGAL PARAMETER VALUE
WAS PASSED TO RKFS. THE USER MUST LOCATE AND CORRECT THE
ERROR.
90 IF (IUNIT.GT.O) WRITE(IUNIT/900) T/NEQ/RELER/ABSER/TINCR
STOP
THE PROGRAM IS STOPPED BECAUSE THERE HAVE BEEN FIVE ERROR
RETURNS FROM RKFS/ INDICATING SERIOUS PROBLEMS WHICH
REQUIRE INTERVENTION BY THE USER.
92 IF (IUNIT.GT.O) WRITEUUNIT/920) T
STOP
THE INTEGRATION HAS BEEN COMPLETED TO OR BEYOND TFINAL.
95 IF (IUNIT.GT.O) WRITE(IUNIT/950) T
RETURN
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
29300
29400
29500
29600
29700
29800
29900
30000
30100
30200
30300
30400
30500
30600
30700
30800
30900
31000
31100
31200
31300
31400
31500
31600
31700
31800
31900
32000
32100
32200
32300
32400
32500
32600
32700
32800
32900
33000
33100
33200
33300
33400
33500
33600
33700
33800
33900
34000
34100
34200
34300
34400
68
-------
C****************************************************************c 34500
C C 34600
C THE REMAINDER OF THE PROGRAM CONTAINS THE FORMAT STATEMENTS C 34700
C FOR THE DIAGNOSTIC MESSAGES. C 34800
C C 34900
100 FORMAT(36H RKFINT HAS BEEN CALLED TO INTEGRATE/15/ 35000
1 10H EQUATIONS/18H FROM INITIAL TIME/612.5/ 35100
2 14H TO FINAL TIME/612.5/ 35200
3 24H WITH OUTPUT AT STEPS OF/G12.5/ 35300
4 32H AND INITIAL ERROR TOLERANCES OF/ 35400
5 9H RELERR =/G12.5/9H ABSERR =/612.5/) 35500
C C 35600
200 FORMAT(4H ***/12(4H****)/ 35700
1 37H *** THE PROGRAM IS TERMINATED AT T =/G12.5/3H***/ 35800
2 52H *** BECAUSE AN UNIDENTIFIABLE ERROR HAS CAUSED ***/ 35900
3 52H *** RKFS TO RETURN WITH THE ILLEGAL FLAG VALUE ***/ 36000
4 12H *** IFLAG =/I11/26X/3H***/4H ***/12(4H****)) 36100
C C 36200
300 FORMAT(17H *WARNING—AT T =/G12.5/21H/ THE VALUE OF RELERR/ 36300
1 49H PASSED TO RKFS WAS FOUND TO BE TOO SMALL FOR THE/ 36400
2 45H MACHINE PRECISION. IT HAS BEEN INCREASED TO/ 36500
3 9H RELERR =/G12.5/23H TO CONTINUE EXECUTION./) 36600
C C 36700
400 FORMAT(17H *WARNING—AT T =/612.5/20H/ RKFS WAS UNABLE TO/ 36800
1 52H MEET THE REQUESTED ERROR TOLERANCES AT THE SMALLEST/ 36900
2 52H ALLOWABLE STEPSIZE. INTEGRATION IS CONTINUING WITH/ 37000
3 37H THE TOLERANCES INCREASED TO RELERR =/G12.5/ 37100
4 13H AND ABSERR =/G12.5/) 37200
C C 37300
500 FORMAT(17H *WARNING—AT T =/G12.5/24H/ THE NUMBER OF FUNCTION/ 37400
1 51H CALLS EXCEEDED THE MAXIMUM NUMBER/ MAXNFE/ ALLOWED/ 37500
2 51H BY RKFS. THE FUNCTION EVALUATION COUNTER HAS BEEN/ 37600
3 51H REDUCED BY MAXNFE TO ALLOW INTEGRATION TO PROCEED./) 37700
C C 37800
600 FORMAT(17H *WARNING—AT T =/G12.5/22H/ STIFFNESS HAS CAUSED/ 37900
1 51H THE NUMBER OF FUNCTION CALLS TO EXCEED THE MAXIMUM/ ' 38000
2 47H NUMBER/ MAXNFE/ ALLOWED BY RKFS. THE FUNCTION/ 38100
3 49H EVALUATION COUNTER HAS BEEN REDUCED BY MAXNFE TO/ 38200
4 30H ALLOW INTEGRATION TO PROCEED./) 38300
C C 38400
700 FORMATC17H *WARNING—AT T =/G12.5/23H/ RKFS RETURNED WITH AN/ 38500
1 50H INDICATION THAT THIS INTEGRATION IS UNNECESSARILY/ 38600
2 48H COSTLY BECAUSE THE REQUESTED OUTPUT INTERVAL IS/ 38700
3 49H SUBSTANTIALLY SMALLER THAN THE NATURAL STEPSIZE./ 38800
4 48H (REFER TO PROGRAM DOCUMENTATION FOR ALTERNATIVE/ 38900
5 52H APPROACHES.) THE OUTPUT FREQUENCY MONITOR HAS BEEN/ 39000
6 39H RESET TO ALLOW INTEGRATION TO PROCEED./) 39100
C C 39200
800 FORMAT(4H ***/12<4H****)/ 39300
1 37H *** THE PROGRAM IS TERMINATED AT T =/G12.5/3H***/ 39400
2 40H *** BECAUSE A COMPONENT OF THE SOLUTION/9X/3H***/ 39500
3 52H *** VANISHED WHILE ABSERR = 0. IT WILL BE ***/ 39600
69
-------
4 52H *** NECESSARY TO PROVIDE A NONZERO ABSOLUTE ***/
5 52H *** ERROR TOLERANCE TO PERFORM THIS INTEGRATION.***/
6 4H ***/12(4H****))
900 FORMAT(4H ***/12(4H****)/
1 37H *** THE PROGRAM IS TERMINATED AT T =/Gl2.5/3H***/
52H *** BECAUSE INVALID INPUT WAS RECEIVED BY RKFS. ***/
37H *** CHECK THE FOLLOWING PARAMETERS—/12X/3H***/
12H *** NEQN =/I11/26X/3H***/
14H *** RELERR =/G12.5/23X/3H***/
14H *** ABSERR =/G12.5/23X/3H***/
13H *** TINCR =/G12.5/24X/3H***/
52H *** ALSO, CHECK THAT THE OUTPUT ROUTINE DOES ***/
2
3
4
5
6
7
8
9
A
41H *** NOT RESET IFLAG TO AN ILLEGAL VALUE./8X/3H***/
4H ***/12<4H****))
920 FORMAT(4H ***,12(4H****)/
1 37H *** THE PROGRAM IS TERMINATED AT T =/G12.5/3H***/
2 52H *** BECAUSE RKFS HAS RETURNED IRREGULARLY FIVE ***/
3 52H *** TIMES/ WHICH INDICATES SOME SORT OF EXTREME ***/
4 52H *** DIFFICULTY IN THE INTEGRATION. THE USER ***/
5 52H *** SHOULD CAREFULLY CHECK THE ACCURACY OF HIS ***/
6 52H *** PROBLEM STATEMENT AND CODING. IT IS ALSO ***/
7 52H *** POSSIBLE THAT THIS PROBLEM IS TOO STIFF TO ***/
8 27H *** BE INTEGRATED BY RKFS./22X/3H***/
9 4H ***/12(4H****)>
C C
950 FORMATC4H ***,12<4H****>/
1 32H *** INTEGRATION COMPLETE AT T =/G12.5/5X/3H***/
2 4H ***/12(4H****)///>
C C
C****************************************************************c
C C
C END OF RKFINT C
C C
END
39700
39800
39900
40000
40100
40200
40300
40400
40500
40600
40700
40800
40900
41000
41100
41200
41300
41400
41500
41600
41700
41800
41900
42000
42100
42200
42300
42400
42500
42600
42700
42800
42900
43000
43100
43200
SUBROUTINE RKFSCF/Y/YP/F1/F2/F3/F4/F5/W)
C C
c****************************************************************c
C
C
C
C
C
C
C
C
C
C
C
SUBROUTINE RKFS USES THE FEHLBERG FOURTH-FIFTH ORDER RUNGE-
-KUTTA METHOD TO INTEGRATE A SYSTEM OF NEQN FIRST ORDER
ORDINARY DIFFERENTIAL EQUATIONS
DY/DT = FCT/Y) / Y = /.../Y
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
THE PARAMETERS IN THE ARGUMENT LIST ARE
F — THE NAME OF THE USER-SUPPLIED SUBROUTINE F(T/Y/YP)/
WHICH EVALUATES THE DERIVATIVES YPU)=DY(I)/DT/
I=1/2/.../NEQN.
Y — AN ARRAY OF DIMENSION NEQN/ WHICH CONTAINS THE VALUE
OF THE SOLUTION.
YP — AN ARRAY OF DIMENSION NEQN/ WHICH CONTAINS THE VALUE
OF THE DERIVATIVE DY/DT.
F1/F2/F3/F4/F5/W — ARRAYS EACH OF DIMENSION NEQN/ WHICH
ARE USED FOR WORKING STORAGE BY RKFS.
>
THE PARAMETERS IN THE COMMON BLOCK RKDATA ARE
T — THE INDEPENDENT VARIABLE
TOUT — THE VALUE OF T AT WHICH OUTPUT IS DESIRED
RELERR/ABSERR — RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR
THE LOCAL ERROR TEST. EACH STEPSIZE IS CHOSEN SO THAT
ABSCLOCAL ERROR(D) .LE. RELERR*ABS(Y(I))-l-ABSERR FOR
EACH COMPONENT OF THE SOLUTION.
NEQN — THE NUMBER OF EQUATIONS
IFLAG — INDICATOR OF THE STATUS OF INTEGRATION
H/NFE/KOP/INIT/NSS/KST — THOSE QUANTITIES USED FOR
INTERNAL BOOKKEEPING BY RKFS WHICH ARE NEEDED FOR
CONTINUING INTEGRATION/ BUT WHICH ARE NORMALLY OF NO
INTEREST TO THE USER. THEY ARE PLACED IN THE COMMON
BLOCK TO AVOID LOCAL RETENTION BY RKFS BETWEEN CALLS.
SEE THE PROGRAM DOCUMENTATION FOR A DESCRIPTION OF THEIR
FUNCTIONING.
TO USE RKFS/ THE CALLING PROGRAM MUST
1. ALLOCATE SPACE FOR THE ARRAYS Y(*)/ YPC*)/ F1(*)/ F2(*)/
F3(*)/ F4C*)/ F5(*)/ W<*) IN A DIMENSION STATEMENT.
(THESE MUST BE ARRAYS EVEN IF NEQN=1.)
2. DECLARE THE NAME OF F IN AN EXTERNAL STATEMENT/ AND
PLACE THAT NAME IN THE CALL LIST.
FIRST CALL TO RKFS
TO BEGIN INTEGRATION/ THE CALLING PROGRAM MUST
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
V
c
c
c
11400
11500
11600
11700
11800
11900
12000
12100
12200
12300
12400
12500
12600
12700
12800
12900
13000
13100
13200
13300
13400
13500
13600
13700
13800
13900
14000
14100
14200
14300
14400
14500
14600
14700
14800
15000
15100
15200
15300
15400
15500
15600
15700
15800
15900
16000
16100
16200
1 WtWw
16300
16400
16500
71
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
1. INITIALIZE THE FOLLOWING VARIABLES IN THE COMMON BLOCK
RKDATA
T — THE STARTING VALUE OF THE INDEPENDENT VARIABLE
TOUT — THE VALUE OF THE INDEPENDENT VARIABLE AT WHICH
OUTPUT IS DESIRED. (TOUT=T IS PERMITTED ON THE
FIRST CALL ONLY/ IN WHICH CASE RKFS EVALUATES THE
DERIVATIVES YP AND RETURNS WITH IFLAG=2.)
RELERR — THE DESIRED RELATIVE ERROR TOLERANCE
ABSERR -- THE DESIRED ABSOLUTE ERROR TOLERANCE
IFLAG — -1-1 FOR INTERVAL MODE OR -1 FOR SINGLE-STEP
MODE.
2. INITIALIZE YC*> TO CONTAIN THE STARTING VALUE OF THE
SOLUTION.
UPON RETURN/ T IS THE LAST VALUE OF THE INDEPENDENT VARIABLE
SUCCESSFULLY REACHED DURING THE INTEGRATION/ Y CONTAINS THE
SOLUTION AT T/ YP CONTAINS THE DERIVATIVES AT T/ AND IFLAG
INDICATES THE STATUS OF THE INTEGRATION AS DESCRIBED BELOW.
SUBSEQUENT CALLS TO RKFS
RKFS NORMALLY RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE
INTEGRATION.
INTERVAL MODE —
IF THE INTEGRATION REACHED TOUT (IFLAG=2)/ THE USER NEED
ONLY DEFINE A NEW TOUT AND CALL RKFS AGAIN.
IF THE INTEGRATION WAS INTERRUPTED WITH IFLAG = 3/ 4/ 5/ 6/
OR 7/ INTEGRATION MAY BE CONTINUED SIMPLY BY CALLING RKFS
AGAIN. HOWEVER/ SUCH A RETURN INDICATES THAT THE
INTEGRATION IS NOT PROCEEDING SO SMOOTHLY AS MIGHT BE
EXPECTED/ SO THE USER SHOULD EXERCISE SOME CAUTION.
IF THE INTEGRATION WAS INTERRUPTED WITH IFLAG = 8/ THE USER
MUST GIVE ABSERR A NONZERO POSITIVE VALUE AND SET IFLAG = 2
IN ORDER TO CONTINUE.
IF RKFS RETURNED WITH IFLAG = 9/ THE USER MUST CORRECT THE
INVALID INPUT AND SET IFLAG = 2 (OR IFLAG = 1/ IF
INTEGRATION HAS NOT YET STARTED) IN ORDER TO CONTINUE.
THE USER MAY AT ANY TIME SWITCH TO SINGLE-STEP MODE BY
REVERSING THE SIGN OF IFLAG.
SINGLE-STEP MODE —
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
16600
16700
16800
16900
17000
17100
17200
17300
17400
17500
17600
17700
17800
17900
18000
18100
18200
18300
18400
18500
18600
18700
18800
A conn
loVUU
19000
19100
19200
19300
19400
19500
19600
19700
19800
19900
20000
20100
20200
20300
20400
20500
20600
20700
20800
20900
21000
21100
21200
21300
21400
21500
21600
21700
72
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
IF THE INTEGRATION HAS PROCEEDED ONE STEP (IFLAG=-2)/. THE
USER NEED ONLY CALL RKFS AGAIN TO PROCEED ANOTHER STEP.
IF THE INTEGRATION HAS REACHED TOUT UFLAG=2>/- THE USER
MUST DEFINE A NEW TOUT AND SET IFLAG = -2 IN ORDER TO
CONTINUE IN SINGLE-STEP MODE. IF IFLAG IS NOT RESET/.
INTEGRATION WILL PROCEED IN INTERVAL MODE.
IF THE INTEGRATION WAS INTERRUPTED WITH IFLAG = -3/ -4/ -5,
-6, OR -7, INTEGRATION MAY BE CONTINUED SIMPLY BY CALLING
RKFS AGAIN.
IF THE INTEGRATION WAS INTERRUPTED WITH IFLAG = 8, THE USER
MUST GIVE ABSERR A NONZERO POSITIVE VALUE AND SET IFLAG =-2
IN ORDER TO CONTINUE.
IF RKFS RETURNED WITH IFLAG = 9, THE USER MUST CORRECT THE
INVALID INPUT AND SET IFLAG =-2 (OR IFLAG =-1* IF
INTEGRATION HAS NOT YET STARTED) IN ORDER TO CONTINUE.
THE USER MAY AT ANY TIME SWITCH TO INTERVAL MODE BY
REVERSING THE SIGN OF A NEGATIVE IFLAG.
UPON RETURN/. THE VALUE OF IFLAG INDICATES THE STATUS OF THE
INTEGRATION AS FOLLOWS/.
IFLAG = 2 ~ T=TOUT, SO THAT INTEGRATION WAS SUCCESSFUL TO
THE OUTPUT POINT.
IFLAG = 3 — INTEGRATION WAS NOT ATTEMPTED BECAUSE THE
REQUESTED ERROR TOLERANCE RELERR WAS SMALLER
THAN THE MINIMUM PERMISSIBLE VALUE REMIN.
RELERR HAS BEEN RESET TO REMIN.
IFLAG = 4 — INTEGRATION WAS INTERRUPTED BECAUSE THE
REQUESTED ERROR TOLERANCES COULD NOT BE
ATTAINED AT THE SMALLEST ALLOWABLE STEPSIZE.
THE TOLERANCES HAVE BEEN INCREASED TO VALUES
WHICH SHOULD ALLOW SUCCESSFUL CONTINUATION.
IFLAG = 5 — INTEGRATION WAS INTERRUPTED AFTER THE COUNTER
OF CALLS TO F EXCEEDED MAXNFE. THE COUNTER HAS
BEEN REDUCED BY MAXNFE.
IFLAG = 6 — SAME AS IFLAG=5/ EXCEPT THAT THE PROBLEM HAS
BEEN ADDITIONALLY DIAGNOSED AS STIFF.
IFLAG = 7 — INTEGRATION WAS NOT CONTINUED BECAUSE THE
DISTANCE FROM T TO TOUT HAS BEEN LESS THAN HALF
THE NATURAL STEPSIZE FOR 100 CALLS TO RKFS/
INDICATING THAT THE CODE IS BEING USED IN A WAY
WHICH GREATLY IMPAIRS ITS EFFICIENCY. THE
MONITOR OF THIS EFFECT WAS RESET TO ALLOW
CONTINUATION IF DESIRED.
IFLAG = 8 — INTEGRATION WAS INTERRUPTED WHEN A COMPONENT OF
THE SOLUTION VANISHED WHILE A PURE RELATIVE
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
21800
21900
22000
22100
22200
22300
22400
22500
22600
22700
22800
22900
23000
23100
23200
23300
23400
23500
23600
23700
23800
23900
24000
o/ 1 nn
CHlUU
24200
24300
24400
24500
24600
24700
24800
24900
25000
25100
25200
25300
25400
25500
25600
25700
25800
25900
26000
26100
26200
26300
26400
26500
26600
26700
26800
26900
73
-------
ERROR TEST WAS SPECIFIED. THE USER MUST SUPPLY
A NONZERO ABSERR TO CONTINUE INTEGRATION.
IFLAG = 9 — INVALID INPUT WAS PASSED TO RKFS. THE USER MUST
LOCATE AND CORRECT THE ERROR IN ORDER TO
CONTINUE INTEGRATION.
IFLAG =-2 — A SINGLE SUCCESSFUL STEP HAS BEEN TAKEN IN THE
DIRECTION OF TOUT/. BUT TOUT HAS NOT YET BEEN
REACHED. THIS OCCURS ONLY IF RKFS WAS CALLED
WITH IFLAG NEGATIVE.
IFLAG =-3/ -4/ -5/ -6/ OR -7 — SAME AS IFLAG = 3/ 4/ 5/ 6/
OR 7/ RESPECTIVELY/ EXCEPT THAT RKFS WAS CALLED
WITH IFLAG NEGATIVE AND IS OPERATING IN SINGLE-
-STEP MODE.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c****************************************************************c
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
THIS SUBROUTINE IS WRITTEN IN SINGLE PRECISION.
WITH MOST COMPILERS/ A DOUBLE PRECISION VERSION MAY BE
OBTAINED SIMPLY BY SETTING APPROPRIATE VALUES FOR U26
AND REMIN IN THE DATA STATEMENT BELOW AND BY REMOVING
THE C FROM COLUMN 1 OF THE NEXT EIGHT CARDS.
DOUBLE PRECISION ABS/ABSERR/AE/AMAX1/CH/DABS/DMAX1/DSIGN/
1 DT/D1/D2/EE/EEC/EECOET/EEOET/ESTTOL/ET/
2 F1/F2/F3/F4/F5/GFACT/H/HMIN/RELERR/REMIN/
3 S/SCALE/SIGN/T/TOL/TOLN/TOUT/U26/
4 Y/YP/YPK/W
ABS(D1) = DABSCD1)
SIGNCD1/D2) = DSIGN(D1/D2)
AMAXKD1/D2) = DMAXKD1/D2)
LOGICAL HFAILD/OUTPUT/STIFF/NOTEST
1
DIMENSION Y(NEQN)/YP(NEQN)/F1(NEQN)/F2(NEQN)/F3(NEQN)/
F4(NEQN)/F5(NEQN)/W(NEQN)
COMMON /RKDATA/ T/TOUT/RELERR/ABSERR/H/NEQN/IFLAG/NFE/KOP/
1 INIT/NSS/KST
IN THE FOLLOWING DATA STATEMENT/ U26 AND REMIN MUST BE SET
TO VALUES APPROPRIATE TO THE MACHINE BEING USED.
U26 IS 26.*U/ WHERE THE COMPUTER UNIT ROUNDOFF ERROR U IS
THE SMALLEST POSITIVE VALUE REPRESENTABLE IN THE MACHINE
SUCH THAT (1.+U).GT.1.
REMIN IS A TOLERANCE THRESHOLD WHICH SHOULD NOT BE SMALLER
THAN 4.*U. SINCE THIS INTEGRATION METHOD CANNOT OPERATE
EFFICIENTLY AT EXTREME TOLERANCES/ REMIN SHOULD BE CHOSEN
AS THE LARGER OF 4.*U AND 1.E-12
SOME REPRESENTATIVE VALUES ARE
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
27000
27100
27200
27300
27400
27500
27600
27700
27800
27900
28000
28100
28200
28300
28400
28500
28600
28700
28800
28900
29000
29100
29200
29300
29400
29500
29600
29700
29800
29900
30000
30100
30200
30300
30400
30500
30600
30700
30800
30900
31000
31100
31200
31300
31400
31500
31600
31700
31800
31900
32000
32100
74
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
U=9.54E-7
U=5.96E-8
U=1.49E-8
U=7.45E-9
U=7.11E-15
U=2.22D-16
U=1.390-17
U26=2.48E-5
U26=1.55E-6
U26=3.88E-7
U26=1.94E-7
U26=1.85E-13
U26=5.780-15
U26=3.61D-16
REMIN=3.82E-6
REMIN=2.39E-7
REMIN=5.97E-8
REMIN=2.99E-8
REMIN=1.E-12
REMIN=1.D-12
REMIN=1.D-12
IBM 360/370 S.P.
PDP-11 S.P.
UNIVAC 1108
PDP-10
CDC 6000 SERIES
IBM 360/370 D.P.
PDP-11 D.P.
DATA U26/1.55E-6//REMIN/2.39E-7/
THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER
OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE.
AS SET/ THIS CORRESPONDS TO ABOUT 1000 STEPS.
DATA MAXNFE/6000/
C CHECK INPUT PARAMETERS
C
MFLAG=IABS(IFLAG)
DT=TOUT-T
C
C THE INPUT IS INVALID IF NEQN IS LESS THAN 1/ IF IFLAG IS ZERO C
C OR HAS ABSOLUTE VALUE GREATER THAN 7/ OR IF EITHER OF RELERR
C AND ABSERR IS NEGATIVE.
C
IF ((MFLAG.GT.7).OR.(MFLAG.EQ.O).OR.(NEQN.LT.D) GO TO 430
IF ((RELERR.LT.0.0).OR.(ABSERR.LT.0.0)) GO TO 430
C
C IF THIS IS THE FIRST CALL/ GO TO INITIALIZATION SECTION.
C
IF (MFLAG.EQ.1) GO TO 400
C
C SINCE THIS IS NOT THE FIRST CALL/ T = TOUT IS NOT VALID INPUT.C
C
IF (T.EQ.TOUT) GO TO 430
C
C IF NO STEP HAS YET BEEN ATTEMPTED/ GO TO INITIALIZATION
C SECTION TO COMPUTE THE STARTING STEPSIZE.
IF (INIT.EQ.O) GO TO 410
C BEFORE CONTINUING INTEGRATION/ CHECK WHETHER REQUESTED ERROR C
C TOLERANCE IS TOO SMALL FOR THE MACHINE PRECISION.
C
IF (RELERR.LT.REMIN) GO TO 440
C
C SET STEPSIZE FOR INTEGRATION IN THE DIRECTION FROM T TO TOUT C
C
80 H=SIGN(H/DT)
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c •
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
32200
32300
32400
32500
32600
32700
32800
32900
33000
33100
33200
33300
33400
33500
33600
33700
33800
33900
34000
34100
34200
34300
34400
34500
34600
34700
34800
34900
35000
35100
35200
35300
35400
35500
35600
35700
35800
35900
360QO
36100
36200
36300
36400
36500
36600
36700
36800
36900
37000
37100
37200
37300
75
-------
C
C TEST WHETHER THE EFFICIENCY OF INTEGRATION IS BEING SEVERELY
C IMPAIRED BY TOO FREQUENT OUTPUT POINTS.
C
IF (ABS(H).GE.2.0*ABS(DT)) KOP=KOP+1
IF (KOP.GT.100) GO TO 450
C
C TEST WHETHER THE INTERVAL FROM T TO TOUT IS TOO SMALL FOR
C MINIMUM ALLOWABLE STEPSIZE.
C
IF (ABS(DT).LT.U26*ABS(T)) GO TO 460
C
C SET INDICATORS FOR STIFFNESS TEST
C
STIFF=. FALSE.
IF (KST.GE.25) STIFF=.TRUE.
NOTEST=STIFF
IF
-------
H=DT
GO TO 200
150 H=0.5*DT
C CORE INTEGRATOR FOR TAKING A SINGCE STEP
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
THE TOLERANCES HAVE BEEN SCALED TO AVOID PREMATURE UNDERFLOW
IN COMPUTING THE ERROR TOLERANCE FUNCTION ET.
TO AVOID PROBLEMS WITH ZERO CROSSINGS/ RELATIVE ERROR IS
MEASURED USING THE AVERAGE OF THE MAGNITUDES OF THE SOLUTION
AT THE BEGINNING AND END OF A STEP.
TO DISTINGUISH THE VARIOUS ARGUMENTS/ H IS NOT PERMITTED TO
BECOME SMALLER THAN 26 UNITS OF ROUNDOFF IN T.
PRACTICAL LIMITS ON THE CHANGE IN THE STEPSIZE ARE ENFORCED
TO SMOOTH THE STEPSIZE SELECTION PROCESS AND TO AVOID
EXCESSIVE CHATTERING ON PROBLEMS HAVING DISCONTINUITIES.
TO PREVENT UNNECESSARY FAILURES/ THE CODE USES 9/10 THE
STEPSIZE IT ESTIMATES WILL SUCCEED.
AFTER A STEP FAILURE/ THE STEPSIZE IS NOT ALLOWED TO INCREASE
FOR THE NEXT ATTEMPTED STEP. THIS MAKES THE CODE MORE
EFFICIENT ON PROBLEMS HAVING DISCONTINUITIES AND MORE
EFFECTIVE IN GENERAL SINCE LOCAL EXTRAPOLATION IS BEING USED
AND EXTRA CAUTION SEEMS WARRANTED.
TEST WHETHER THE NUMBER OF FUNCTION EVALUATIONS HAS EXCEEDED
THE ALLOWED LIMIT.
200 IF (NFE.GT.MAXNFE) GO TO 500
ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H
THE FORMULAS HAVE BEEN GROUPED TO CONTROL LOSS OF
SIGNIFICANCE.
CH=0.25*H
DO 221 I=1/NEQN
221 W(I)=Y(I)+CH*YP(I)
CALL F(T+0.25*H/W/F1)
CH=0.09375*H
DO 222 I=1/NEQN
222 W(I)=Y(I)+CH*(YP(I)+3.*F1(I))
CALL F(T+0.375*H/W/F2)
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
42600
42700
42800
42900
43000
43100
43200
43300
43400
43500
43600
43700
43800
43900
44000
44100
44200
44300
44400
44500
44600
44700
44800
44900
45000
45100
45200
45300
45400
45500
45600
45700
45800
45900
46000
46100
46200
46300
46400
46500
46600
46700
46800
46900
47000
47100
47200
47300
47400
47500
47600
47700
77
-------
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
CH=H/2197.
DO 223 I=1/NEQN
223 W(I)=Y(I)+CH*(1932.*YP(I)+(7296.*F2(I)-7200.*F1(I)))
CALL F(T+12.*H/13./W/F3>
CH=H/4104.
DO 224 I=1/NEQN
224 WU>=Y(mCH*((8341.*YP(I)-845.*F3(I» +
1 (29440.*F2(I)-32832.*F1(I)»
CALL F(T+H/W/F4)
CH=H/20520.
DO 225 I=1,NEQN
225 W(I>=Y(I)KH*((-6080.*YP(I) + (9295.*F3U)-
1 5643.*F4CI)» + (41040.*F1(I)-28352.*F2(I)»
CALL F(T+0.5*H,W,F5)
NFE=NFE+5
COMPUTE SOLUTION W AT T+H
COMPUTE SCALED RATIO EEOET OF LOCAL ERROR ESTIMATE TO
ALLOWABLE TOLERANCE.
CH=H/7618050.
EEOET=0.0
DO 230 I=1,NEQN
W(I)=Y(I)+CH*((902880.*YPU) + (3855735.*F3(I>-
1 1371249. *F4(I»> + (3953664.*F2(I)+277020.*F5(I»)
ET=ABS(Y(I))+ABS(W(I))+AE
IF ET VANISHES/1 IT IS IMPOSSIBLE TO PROCEED. THE USER
MUST PROVIDE A NONZERO VALUE OF ABSERR.
IF (ET.EQ.0.0) 60 TO 480
EE=ABS«-2090.*YPU)-K21970.*F3(I)-15048.*F4(I)» +
1 (22528. *F2(I>-27360.*F5(I>»
230 EEOET=AMAX1(EEOET/EE/ET)
REMOVE SCALING OF THE ERROR ESTIMATE RATIO. THE STEP IS
ACCEPTED IF THE RESULT DOES NOT EXCEED 1.
ESTTOL = ABS(H)*EEOET*SCALE/752400.
IF (ESTTOL.LE.1.) GO TO 260
THE STEP TO T+H DID NOT SUCCEED.
REDUCE THE STEPSIZE AND MAKE ANOTHER ATTEMPT.
THE DECREASE IS LIMITED TO A FACTOR OF 1/10
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
47800
47900
48000
48100
48200
48300
48400
48500
48600
48700
48800
48900
49000
49100
49200
49300
49400
49500
49600
49700
49800
49900
50000
50100
50200
50300
50400
50500
50600
50700
50800
50900
51000
51100
51200
51300
51400
51500
51600
51700
51800
51900
52000
52100
52200
52300
52400
52500
52600
52700
52800
52900
78
-------
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
HFAILD=.TRUE.
OUTPUT=.FALSE.
S=0.1
IF (ESTTOL.LT.59049.) S=0.9/ESTTOL**0.2
H=S*H
IF (ABS(H).GE.HMIN) GO TO 200
THE REQUESTED ERROR TOLERANCES RELERR AND ABSERR ARE TOO
SMALL TO BE MET AT THE SMALLEST ALLOWABLE STEPSIZE. THEY
ARE INCREASED TO VALUES WHICH SHOULD ALLOW SUCCESS.
GFACT=(ABS(HMIN/H))**5
IF (GFACT.LT.2.) GFACT=2.
RELERR=GFACT*RELERR
ABSERR=GFACT*ABSERR
H=SIGN(HMIN/H)
IFLAG=ISIGN(4/IFLAG)
RETURN
THE STEP TO T+H WAS SUCCESSFUL.
IF THE PROBLEM HAS ALREADY BEEN FOUND STIFF/ BYPASS FURTHER
TESTING.
260 IF (STIFF) GO TO 290
INCREMENT COUNTER OF SUCCESSFUL STEPS. IF THE CURRENT BLOCK
OF 50 STEPS HAS ALREADY BEEN FOUND NOT STIFF/ DO NOT TEST
THIS STEP.
NSS=NSS+1
IF (NOTEST) GO TO 285
TEST WHETHER THIS STEP SHOWS STIFFNESS.
EECOETO.O
DO 275 I=1/NEQN
ET=ABS(Y(imABS(W(I))+AE
EEC=ABS((-720895.*YP(I)+475423.*F2(I))+(461389.*F1(I)
1 -300391 .*F3(I))+123695.*F4(I)-39221.*F5(D)
275 EECOET=AMAX1(EECOET/EEC/ET)
IF (ABS(H)*EECOET*SCALE.GT.13.E6) GO TO 280
C
C
C
C
C
C
THIS STEP SHOWS STIFFNESS. IF THE NUMBER OF SUCH STEPS IN
THE CURRENT BLOCK HAS REACHED 25/ SET STIFF AND NOTEST TO
INDICATE THAT THE PROBLEM IS STIFF/ AND TO BYPASS ALL
TESTING OF SUBSEQUENT STEPS.
KST=KST+1
IF (KST.LT.25) GO TO 285
STIFF=.TRUE.
NOTEST=.TRUE.
GO TO 290
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
53000
53100
53200
53300
53400
53500
53600
53700
53800
53900
54000
54100
54200
54300
54400
54500
54600
54700
54800
54900
55000
55100
55200
55300
55400
55500
55600
55700
55800
55900
56000
56100
56200
56300
56400
56500
56600
56700
56800
56900
57000
57100
57200
57300
57400
57500
57600
57700
57800
57900
58000
58100
79
-------
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
C
C
C
C
C
C
THIS STEP DOES NOT SHOW STIFFNESS. IF THE CURRENT BLOCK IS
THEREBY ESTABLISHED AS NOT STIFF, SET NOTEST=.TRUE. TO
SUSPEND TESTING FOR THE REMAINDER OF THIS BLOCK.
280 IF UNSS-KST).GT.25> NOTEST=.TRUE.
THE PROBLEM HAS NOT BEEN ESTABLISHED AS STIFF. IF THE
CURRENT BLOCK IS COMPLETE, REINITIALIZE INDICATORS TO
TEST NEXT BLOCK.
285 IF (NSS.LT.50) GO TO 290
NSS=0
KST=0
NOTEST=.FALSE.
REPLACE T WITH T+H
STORE SOLUTION IN Y
EVALUATE DERIVATIVES YP
290 T=T+H
DO 295 K=1,NEQN
295 Y(K)=W(K)
CALL F(T,Y,YP)
NFE=NFE+1
ESTIMATE OPTIMUM STEPSIZE FOR NEXT STEP.
THE INCREASE IS LIMITED TO A FACTOR OF 5
IF STEP FAILURE HAS JUST OCCURRED, THE STEPSIZE IS NOT
ALLOWED TO INCREASE.
S=5.0
IF (ESTTOL.GT.1.889568E-4) S=0.9/ESTTOL**0.2
IF (HFAILD.AND.(S.GT.I.O)) S=1.0
H=SIGN(AMAX1
C END OF CORE INTEGRATOR
c****************************************************************c
c
C DETERMINE WHETHER TO TAKE ANOTHER STEP BEFORE RETURN.
C
IF (OUTPUT) GO TO 300
IF (IFLAG.GT.O) GO TO 100
C
c****************************************************************c
C****************************************************************c
c
C INTEGRATION SUCCESSFULLY COMPLETED
C
C SINGLE-STEP MODE
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
rC
c
fC
c
c
c
c
trC
*c
c
c
c
c
c
58200
58300
58400
58500
58600
58700
58800
58900
59000
59100
59200
59300
59400
59500
59600
59700
59800
59900
60000
60100
60200
60300
60400
60500
60600
60700
60800
60900
61000
61100
61200
61300
61400
61500
61600
61700
61800
61900
62000
62100
62200
62300
62400
62500
62600
62700
62800
62900
63000
63100
63200
63300
80
-------
C
C
C
IFLAG=-2
RETURN
INTERVAL MODE
300 T=TOUT
IFLAG=2
RETURN
c****************************************************************c
C
C INITIALIZATION SECTION —
C
C
C
C
C
C
C
C
C
C
ON FIRST CALL (IFLA6 = +1 OR -1),
SET INITIALIZATION COMPLETION INDICATOR INIT
SET EXCESSIVE OUTPUT FREQUENCY INDICATOR KOP
SET FUNCTION EVALUATION COUNTER NFE
SET STIFFNESS TEST COUNTERS NSS/KST
EVALUATE INITIAL DERIVATIVES
IF TOUT = T, RETURN TO CALLING PROGRAM
400 INIT=0
KOP=0
NFE=1
NSS=0
KST=0
CALL F(T,Y,YP)
IF (T.NE.TOUT) GO TO 410
IFLAG=2
RETURN
C
C
C
C
C
C
C
C
TEST WHETHER REQUESTED ERROR TOLERANCE IS TOO SMALL FOR THE
MACHINE PRECISION.
410 IF (RELERR.LT.REMIN) GO TO 440
ESTIMATE STARTING STEPSIZE, SET INIT = 1 TO INDICATE THAT
THIS HAS BEEN DONE/ AND TRANSFER TO BEGIN INTEGRATION.
INIT=1
H=ABS(DT)
TOLN=0.0
DO 420 K=1/NEQN
TOL=RELERR*ABS(YCK))+ABSERR
IF (TOL.LE.0.0) GO TO 420
TOLN=TOL
YPK=ABS
IF (YPK*H**5.GT.TOL> H=(TOL/YPK)**0.2
420 CONTINUE
IF (TOLN.LE.0.0) H=0.0
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
63400
63500
63600
63700
63800
63900
64000
64100
64200
64300
64400
64500
64600
64700
64800
64900
65000
65100
65200
65300
65400
65500
65600
65700
65800
65900
66000
66100
66200
66300
66400
66500
66600
66700
66800
66900
67000
67100
67200
67300
67400
67500
67600
67700
67800
67900
68000
68100
68200
68300
68400
68500
81
-------
H=AMAX1 (H,U26*AMAX1 (ABS (T) , ABS
-------
KST=0 73800
IF (STIFF) IFLAG=ISIGN(6/IFLAG) 73900
RETURN 74000
C C 74100
C****************************************************************C 74200
C****************************************************************C 74300
C : C 74400
C END OF RKFS C 74500
C C 74600
END 74700
83
-------
APPENDIX B
ADAMS METHOD PROGRAMS
This appendix contains the complete listings of the subroutines ADMINT,
ADAM, STEP, and SINTRP.
SUBROUTINE ADMINT USES ADAM TO INTEGRATE A SYSTEM OF NEQN
FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS
DY/DT = FCT/Y) / Y =
-------
OUTP(T/Y/W) USED TO HANDLE OUTPUT.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c****************************************************************c
C
C
C
C
C
C
C
C
C
C
C
TO USE ADMINT/ THE CALLING PROGRAM MUST
1. ALLOCATE SPACE FOR THE ARRAYS Y(*) AND W(*) IN A
DIMENSION STATEMENT. (NOTE THAT Y MUST BE AN ARRAY/
EVEN IF NEQN=1.)
2. SET Y(1)/Y(2)/.../Y(NEQN> TO THE INITIAL VALUES TO BE
USED FOR THE SOLUTION.
3. PROVIDE APPROPRIATE VALUES FOR NEQN/ TINIT/ TFINAL/
TINCR/ RELERR/ ABSERR/ AND IUNIT. (SINCE THESE
QUANTITIES ARE NOT ALTERED BY ADMINT/ THEY MAY/ IF
DESIRED/ BE WRITTEN AS CONSTANTS IN THE CALL LIST.)
4. PLACE THE ACTUAL NAMES OF FCT AND OUTP IN THE CALL
LIST/ AND DECLARE THOSE NAMES IN AN EXTERNAL STATEMENT.
AT EACH CALL TO OUTP(T/Y/W)/ THE CURRENT VALUES OF THE
DERIVATIVES DY(I)/DT ARE IN W(1)/W(2)/.../W(NEQN). THE
FIRST CALL TO OUTP IS MADE AT TINIT.
IT IS NOT NECESSARY THAT THE INTEGRATION INTERVAL BE EVENLY
DIVISIBLE BY THE OUTPUT MESH SIZE TINCR. HOWEVER/ IF IT IS
NOT/ INTEGRATION WILL PROCEED TO THE FIRST MESH POINT
BEYOND TFINAL/ RATHER THAN STOPPING AT TFINAL.
C
C
C
C
C
THIS SUBROUTINE IS WRITTEN IN SINGLE PRECISION.
A DOUBLE PRECISION VERSION FOR USE WITH A DOUBLE PRECISION
VERSION OF ADAM MAY BE OBTAINED SIMPLY BY REMOVING THE C
FROM COLUMN 1 OF THE NEXT FIVE CARDS.
DOUBLE PRECISION ABS/ABSER/ABSERR/DABS/DSIGN/DT/DTSGN/
1 D1/D2/EPSDT/RELER/RELERR/SIGN/T/TFINAL/
2 TINCR/TINIT/TMAXIN/TOUT/W/Y
ABS(D1) = DABSCD1)
SIGN(D1/D2) = DSIGNCD1/D2)
COMMON /ADDATA/ T/TOUT/RELER/ABSER/DTSGN/NEQ/IFLAG/
1 IOLD/KLE4/INIT/KOP
DIMENSION Y(NEQN)/W(1)
THE ACTUAL DIMENSION OF W/ AS DEFINED BY THE CALLING PROGRAM/
MUST BE AT LEAST 19*NEQN.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
13700
13800
13900
14000
14100
14200
14300
14400
14500
14600
14700
14800
14900
15000
15100
15200
15300
15400
15500
15600
15700
15800
15900
16000
16100
16200
16300
16400
16500
16600
16700
16800
16900
17000
17100
17200
17300
17400
17500
17600
17700
17800
17900
18000
18100
18200
18300
18400
18500
18600
18700
18800
85
-------
EXTERNAL FCT
C
C WRITE STARTING DATA INTO MESSAGE FILE
C
IF (IUNIT.GT.O) WRITEUUNIT/100) NEQNrTINIT/TFINALxTINCR/'
1 RELERR/ABSERR
C
C COMPUTE INDICES FOR SPLITTING WORK ARRAY W<*).
C
K1=NEQN+1
K2=K1+NEQN
K3=K2+NEQN
C
C INITIALIZE VARIABLES TO START INTEGRATION WITH ADAM.
C
NEQ=NEQN
RELER=RELERR
ABSER=ABSERR
IFLAG=1
T=TINIT
TOUT=TINIT
C
C SAVE INFORMATION NEEDED TO LOCATE END OF INTEGRATION.
C
TMAXIN=ABS(TFINAL-TINIT)
EPSDT=0.001*ABSCTINCR)
C SET SIGN OF OUTPUT INTERVAL TO PROCEED FROM TINIT TO TFINAL.
DT=SIGN (TINCR^TFINAL-TINIT)
C
C
C
C
C INITIALIZE COUNTER OF ERROR RETURNS.
C
IER=0
C
C CALL ADAM TO INTEGRATE TO THE NEXT OUTPUT POINT.
C
5 CALL ADAM(FCT,Y,W(1)/W(K1)/W(K2),W(K3)>
IF (IFLAG.NE.2) GO TO 8
C
C SUCCESSFUL INTEGRATION TO OUTPUT POINT.
C
CALL OUTP(T/Y,W)
C
C CHECK WHETHER TFINAL HAS BEEN REACHED OR PASSED.
C
IF (ABS(T-TINIT).GE.TMAXIN) GO TO 95
C
C INTEGRATION IS NOT YET COMPLETE/ SO THE NEXT OUTPUT POINT IS C
C SET. IF THE NEW POINT DIFFERS FROM TFINAL BY AN AMOUNT
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
•c
C
C
C
C
C
C
C
C
C
C
C
C
18900
19000
19100
19200
19300
19400
19500
19600
19700
19800
19900
20000
20100
20200
20300
20400
20500
20600
20700
20800
20900
21000
21100
21200
21300
21400
21500
21600
21700
21800
21900
22000
22100
22200
22300
22400
22500
22600
22700
22800
22900
23000
23100
23200
23300
23400
23500
23600
23700
23800
23900
24000
86
-------
C ATTRIBUTABLE TO ROUNDOFF ERROR/ IT IS RESET TO TFINAL EXACTLY.C
C INTEGRATION IS THEN CONTINUED.
C
TOUT=T+DT
IF (ABS(TFINAL-TOUT).LT.EPSDT) TOUT=TFINAL
GO TO 5
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
INTEGRATION WAS INTERRUPTED BEFORE REACHING THE OUTPUT POINT.
INCREMENT THE ERROR COUNTER AND BRANCH FOR WRITING THE
APPROPRIATE DIAGNOSTIC MESSAGE. INTEGRATION WILL THEN BE
CONTINUED IF THE ERROR IS NOT TOO SEVERE.
8 IER=IER+1
GO TO (20/20/30/40/50/60/70/80/90)/IFLAG
IFLAG HAS AN ILLEGAL VALUE/ SO THAT A PRECISE DEFINITION OF
THE ERROR IS NOT POSSIBLE/ AND EXECUTION MUST BE TERMINATED.
20 IF (IUNIT.GT.O) WRITECIUNIT/200) T/IFLAG
STOP
IFLAG = 3 — THE REQUESTED ERROR TOLERANCES WERE TOO
SMALL FOR THE MACHINE PRECISION AND HAVE BEEN INCREASED TO
MINIMUM PERMISSIBLE VALUES.
30 IF (IUNIT.GT.O) WRITE(IUNIT/300) T/RELER/ABSER
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 4 — THE REQUESTED ERROR TOLERANCES WERE TOO
STRINGENT FOR THIS PROBLEM AND HAVE BEEN INCREASED TO
VALUES WHICH SHOULD BE ACCEPTABLE.
40 IF (IUNIT.GT.O) WRITEUUNIT/400) T/RELER/ABSER
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 5 — THE COUNTER OF CALLS TO FCT WAS RESET AFTER
EXCEEDING THE LIMIT IMPOSED IN ADAM.
50 IF UUNIT.GT.O) WRITEUUNIT/500) T
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG .= 6 — THE COUNTER OF CALLS TO FCT WAS RESET AFTER
STIFFNESS CAUSED IT TO EXCEED THE LIMIT IMPOSED IN ADAM.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
24100
24200
24300
24400
24500
24600
24700
24800
24900
25000
25100
25200
25300
25400
25500
25600
25700
25800
25900
26000
26100
26200
26300
26400
26500
26600
26700
26800
26900
27000
27100
27200
27300
27400
27500
27600
27700
27800
27900
28000
28100
28200
28300
28400
28500
28600
28700
28800
28900
29000
29100
29200
87
-------
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
60 IF (IUNIT.GT.O) WRITEUUNIT/600) T
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 7 — THE MONITOR OF OUTPUT FREQUENCY WAS RESET
AFTER INDICATING THAT THE OUTPUT INTERVAL IS TOO SMALL
TO ALLOW EFFICIENT OPERATION OF ADAM.
70 IF (IUNIT.GT.O) WRITEUUNIT/700) T
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 8 — FATAL ERROR — A COMPONENT OF THE SOLUTION
VANISHED/ MAKING A PURE RELATIVE ERROR TEST IMPOSSIBLE.
THE USER MUST SUPPLY AN APPROPRIATE NONZERO VALUE OF
ABSERR FOR SUCCESSFUL INTEGRATION.
80 IF (IUNIT.GT.O) WRITE(IUNIT/800) T
STOP
IFLAG = 9 — FATAL ERROR — AN ILLEGAL PARAMETER VALUE
WAS PASSED TO ADAM. THE USER MUST LOCATE AND CORRECT THE
ERROR.
90 IF (IUNIT.GT.O) WRITE(IUNIT/900) T/NEQ/RELER/ABSER/TINCR
STOP
THE PROGRAM IS STOPPED BECAUSE THERE HAVE BEEN FIVE ERROR
RETURNS FROM ADAM, INDICATING SERIOUS PROBLEMS WHICH
REQUIRE INTERVENTION-BY THE USER.
92 IF (IUNIT.GT.O) WRITE(IUNIT/920) T
STOP
THE INTEGRATION HAS BEEN CONPLETED TO OR BEYOND TFINAL.
95 IF (IUNIT.GT.O) WRITE(IUNIT/950) T
RETURN
C
C
C
C
THE REMAINDER OF THE PROGRAM CONTAINS THE FORMAT STATEMENTS
FOR THE IDAGNOSTIC MESSAGES.
100 FORMAT(36H ADMINT HAS BEEN CALLED TO INTEGRATE/15/
1 10H EQUATIONS/18H FROM INITIAL TIME/612. 5/
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
rC
c
c
c
c
29300
29400
29500
29600
29700
29800
29900
30000
30100
30200
30300
30400
30500
30600
30700
30800
30900
31000
31100
31200
31300
31400
31500
31600
31700
31800
31900
32000
32100
32200
32300
32400
32500
32600
32700
32800
32900
33000
33100
33200
33300
33400
33500
33600
33700
33800
33900
34000
34100
34200
34300
34400
88
-------
2 14H TO FINAL TIME/G12.5/ 34500
3 24H WITH OUTPUT AT STEPS OF/G12.5/ 34600
4 32H AND INITIAL ERROR TOLERANCES OF/ 34700
5 9H RELERR =/G12.5/9H ABSERR,=/G12.5/) 34800
C C 34900
200 FORMAT(4H ***/12(4H****)/ 35000
1 37H *** THE PROGRAM IS TERMINATED AT T =/G12.5/3H***/ 35100
2 52H *** BECAUSE AN UNIDENTIFIABLE ERROR HAS CAUSED ***/ 35200
3 52H *** ADAM TO RETURN WITH THE ILLEGAL FLAG VALUE ***/ 35300
4 12H *** IFLAG =/I11/26X/3H***/4H ***,12(4H****)) 35400
C C 35500
300 FORMATC17H *WARNING—AT T =/G12.5/22H/ THE ERROR TOLERANCES/ 35600
1 50H PASSED TO ADAM WERE FOUND TO BE TOO SMALL FOR THE/ 35700
2 48H MACHINE PRECISION. THEY HAVE BEEN INCREASED TO/ 35800
3 9H RELERR =/G12.5/13H AND ABSERR =/G12.5/ 35900
4 25H TO CONTINUE INTEGRATION./) 36000
C C 36100
400 FORMATC17H *WARNING—AT T =/G12.5/20H/ ADAM WAS UNABLE TO/ 36200
1 52H MEET THE REQUESTED ERROR TOLERANCES AT THE SMALLEST/ 36300
2 52H ALLOWABLE STEPSIZE. INTEGRATION IS CONTINUING WITH/ 36400
3 37H THE TOLERANCES INCREASED TO RELERR =/G12.5/ 36500
4 13H AND ABSERR =/G12.5/> 36600
C C 36700
500 FORMAT(17H *WARNING—AT T =/G12.5/24H/ THE NUMBER OF FUNCTION/ 36800
1 51H CALLS EXCEEDED THE MAXIMUM NUMBER/ MAXNFE/ ALLOWED/ 36900
2 51H BY ADAM. THE FUNCTION EVALUATION COUNTER HAS BEEN/ 37000
3 51H REDUCED BY MAXNFE TO ALLOW INTEGRATION TO PROCEED./) 37100
C C 37200
600 FORMATC17H *WARNING—AT T =/G12.5/22H/ STIFFNESS HAS CAUSED/ 37300
1 51H THE NUMBER OF FUNCTION CALLS TO EXCEED THE MAXIMUM/ 37400
2 47H NUMBER/ MAXNFE/ ALLOWED BY ADAM. THE FUNCTION/ 37500
3 49H EVALUATION COUNTER HAS BEEN REDUCED BY MAXNFE TO/ 37600
4 30H ALLOW INTEGRATION TO PROCEED./) 37700
C C 37800
700 FORMAT(17H *WARNING—AT T =/G12.5/23H/ ADAM RETURNED WITH AN/ 37900
1 50H INDICATION THAT THIS INTEGRATION IS UNNECESSARILY/ 38000
2 48H COSTLY BECAUSE THE REQUESTED OUTPUT INTERVAL IS/ 38100
3 49H SUBSTANTIALLY SMALLER THAN THE NATURAL STEPSIZE./ 38200
4 48H (REFER TO PROGRAM DOCUMENTATION FOR ALTERNATIVE/ 38300
5 52H APPROACHES.) THE OUTPUT FREQUENCY MONITOR HAS BEEN/ 38400
6 39H RESET TO ALLOW INTEGRATION TO PROCEED./) 38500
C C 38600
800 FORMAT(4H ***/12<4H****)/ 38700
1 37H *** THE PROGRAM IS TERMINATED AT T =/G12.5/3H***/ 38800
2 40H *** BECAUSE A COMPONENT OF THE SOLUTION/9X/3H***/ 38900
3 52H *** VANISHED WHILE ABSERR = 0. IT WILL BE ***/ 39000
4 52H *** NECESSARY TO PROVIDE A NONZERO ABSOLUTE ***/ 39100
5 52H *** ERROR TOLERANCE TO PERFORM THIS INTEGRATION.***/ 39200
6 4H ***/12<4H****)) 39300
C C 39400
900 FORMAT<4H ***/12(4H****)/ 39500
1 37H *** THE PROGRAM IS TERMINATED AT T =/G12.5/3H***/ 39600
89
-------
2
3
4
5
6
7
8
9
A
52H *** BECAUSE INVALID INPUT WAS RECEIVED BY ADAM. ***/
37H *** CHECK THE FOLLOWING PARAMETERS— ,12X,3H***/
NEQN =,I11/26X/3H***/
RELERR =,G12.5/23X/3H***/
ABSERR =/G12.5/23X/3H***/
TINCR =,G12.5/24X/3H***/
52H *** ALSO, CHECK THAT THE OUTPUT ROUTINE DOES ***/
41 H *** NOT RESET IFLAG TO AN ILLEGAL VALUE. /8X/3H***/
4H ***,12(4H****)>
12H ***
HH ***
14H ***
13H ***
920 FORMAT (4H ***/12<4H****)/
1 37H *** THE PROGRAM IS TERMINATED AT T =,G12.5,3H***/
2 52H *** BECAUSE ADAM HAS RETURNED IRREGULARLY FIVE ***/
3 52H *** TIMES/ WHICH INDICATES SOME SORT OF EXTREME ***/
4 52H *** DIFFICULTY IN THE INTEGRATION. THE USER ***/
5 52H *** SHOULD CAREFULLY CHECK THE ACCURACY OF HIS ***/
6 52H *** PROBLEM STATEMENT AND CODING. IT IS ALSO ***/
7 52H *** POSSIBLE THAT THIS PROBLEM IS TOO STIFF TO ***/
8 27H *** BE INTEGRATED BY ADAM./22X/3H***/
9 4H ***,12(4H****)>
950 FORMAT <4H ***,12<4H****)/
1 32H *** INTEGRATION COMPLETE AT T =/G12.5,5X/3H***/
2 4H ***/12<4H****>///>
C
C
C
END OF ADMINT
END
C
C
C
C
C
C
C
39700
39800
39900
40000
40100
40200
40300
40400
40500
40600
40700
40800
40900
41000
41100
41200
41300
41400
41500
41600
41700
41800
41900
42000
42100
42200
42300
42400
42500
42600
SUBROUTINE ADAM USES SUBROUTINE STEP TO INTEGRATE A YSYTEM OF
NEQN FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS
DY/DT = F(T/Y) / Y =
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
WHICH EVALUATES THE DERIVATIVES YP(I)=DY(I)/DT/
I=1/2/.../NEQN.
Y — AN ARRAY OF DIMENSION NEQN/ WHICH CONTAINS THE VALUE
OF THE SOLUTION.
YP — AN ARRAY OF DIMENSION NEQN/ WHICH CONTAINS THE VALUE
OF THE DERIVATIVE DY/DT.
WT/YY/PHI — ARRAYS OF DIMENSIONS NEQN/ NEQN/ AND 16*NEQN/
RESPECTIVELY/ WHICH ARE USED FOR WORKING
STORAGE BY ADAM AND STEP. THE VALUES STORED
IN YY AND PHI ARE NEEDED FOR CONTINUING
INTEGRATION ON SUBSEQUENT CALLS.
THE PARAMETERS IN THE COMMON BLOCK ADDATA ARE
T — THE INDEPENDENT VARIABLE
TOUT — THE VALUE OF T AT WHICH OUTPUT IS DESIRED
RELERR/ABSERR — RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR
THE LOCAL ERROR TEST. EACH STEPSIZE IS CHOSEN SO THAT
ABS(LOCAL ERROR(D) .LE. RELERR*ABS(YU))+ABSERR FOR
EACH COMPONENT OF THE SOLUTION.
NEQN — THE NUMBER OF EQUATIONS
IFLAG — INDICATOR OF THE STATUS OF INTEGRATION
DTSGN/IOLD/KLE4/INIT/KOP — THOSE QUANTITIES USED FOR
INTERNAL BOOKKEEPING BY ADAM WHICH ARE NEEDED FOR
CONTINUING INTEGRATION/ BUT WHICH ARE NORMALLY OF NO
INTEREST TO THE USER. THEY ARE PLACED IN THE COMMON
BLOCK TO AVOID LOCAL RETENTION BY ADAM BETWEEN CALLS.
SEE THE PROGRAM DOCUMENTATION FOR A DESCRIPTION OF THEIR
FUNCTIONING.
TO USE ADAM/ THE CALLING PROGRAM MUST
1. ALLOCATE SPACE FOR THE ARRAYS YCNEQN)/ YP(NEQN)/
WT(NEQN)/ YY(NEQN)/ AND PHKNEQN/16) .
CY/ YP/ WT/ AND YY MUST BE ARRAYS EVEN IF NEQN=1.)
2. DECLARE THE NAME OF F IN AN EXTERNAL STATEMENT/ AND
PLACE THAT NAME IN THE CALL LIST.
FIRST CALL TO ADAM
TO BEGIN INTEGRATION/ THE CALLING PROGRAM MUST
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
V
c
c
c
12000
12100
12200
12300
12400
12500
12600
12700
12800
12900
13000
13100
13200
13300
13400
13500
13600
13700
13800
13900
14000
14100
14200
14300
14400
14500
14600
14700
14800
14900
15000
15100
15200
15300
15400
•jccnn
1 JJUU
15600
15700
15800
15900
16000
16100
16200
16300
16400
16500
16600
16700
16800
1 WWWW
16900
17000
17100
91
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
1. INITIALIZE THE FOLLOWING VARIABLES IN THE COMMON BLOCK
ADDATA
T — THE STARTING VALUE OF THE INDEPENDENT VARIABLE
TOUT — THE VALUE OF THE INDEPENDENT VARIABLE AT WHICH
OUTPUT IS DESIRED. (TOUT=T IS PERMITTED ON THE
FIRST CALL ONLY/ IN WHICH CASE ADAM EVALUATES THE
DERIVATIVES YP AND RETURNS WITH IFLAG=2.)
RELERR — THE DESIRED RELATIVE ERROR TOLERANCE
ABSERR — THE DESIRED ABSOLUTE ERROR TOLERANCE
I FLAG — NORMALLY +1. IF THE INTEGRATION CANNOT
PROCEED BEYOND TOUT/ SET IFLAG = -1.
2. INITIALIZE Y(*) TO CONTAIN THE STARTING VALUE OF THE
SOLUTION.
UPON SUCCESSFUL RETURN (IFLAG=2)/ T=TOUT/ Y CONTAINS THE
SOLUTION AT TOUT/ AND YP CONTAINS THE DERIVATIVES AT TOUT.
OTHERWISE/ T IS THE LAST VALUE OF THE INDEPENDENT VARIABLE
SUCCESSFULLY REACHED DURING THE INTEGRATION/ Y CONTAINS THE
SOLUTION AT T/ YP CONTAINS THE DERIVATIVES AT T/ AND IFLAG
INDICATES THE STATUS OF THE INTEGRATION AS DESCRIBED BELOW.
SUBSEQUENT CALLS TO ADAM
ADAM NORMALLY RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE
INTEGRATION.
IF THE INTEGRATION REACHED TOUT (IFLAG=2)/ THE USER NEED
ONLY DEFINE A NEW TOUT AND CALL ADAM AGAIN. IF THE INTE-
GRATION CANNOT CONTINUE INTERNALLY BEYOND THE NEW TOUT/
THE USER NUST ALSO SET IFLAG=-2. IF THE NEW TOUT REPRESENTS
A REVERSAL OF THE DIRECTION OF INTEGRATION/ IT IS NECESSARY
THAT Y NOT HAVE BEEN ALTERED.
IF THE INTEGRATION WAS INTERRUPTED WITH ABS(IFLAG) = 3/ 4/ 5/
6/ OR 7/ INTEGRATION MAY BE CONTINUED SIMPLY BY CALLING ADAM
AGAIN. HOWEVER/ SUCH A RETURN INDICATES THAT THE
INTEGRATION IS NOT PROCEEDING SO SMOOTHLY AS MIGHT BE
EXPECTED/ SO THE USER SHOULD EXERCISE SOME CAUTION.
IF THE INTEGRATION WAS INTERRUPTED WITH IFLAG = 8/ THE USER
MUST GIVE ABSERR A NONZERO POSITIVE VALUE AND SET IFLAG = 2
(OR IFLAG=-2 IF IT IS NECESSARY TO STOP AT TOUT) IN ORDER TO
CONTINUE.
IF ADAM RETURNED WITH IFLAG = 9/ THE USER MUST CORRECT THE
INVALID INPUT AND SET IFLAG = 2 OR -2 (OR IFLAG = 1 OR -1/
IF INTEGRATION HAS NOT YET STARTED) IN ORDER TO CONTINUE.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
17200
17300
17400
17500
17600
17700
17800
17900
18000
18100
18200
18300
18400
18500
18600
18700
18800
18900
19000
19100
19200
19300
19400
19500
19600
19700
1 oonn
IVoUU
19900
20000
20100
20200
20300
20400
20500
20600
20700
20800
20900
21000
21100
21200
21300
21400
21500
21600
21700
21800
21900
22000
22100
22200
22300
92
-------
c c
c c
C UPON RETURN/- THE VALUE OF IFLA6 INDICATES THE STATUS OF THE C
C INTEGRATION AS FOLLOWS, C
c c
C IFLAG = 2 — T=TOUT/- SO THAT INTEGRATION WAS SUCCESSFUL TO C
C THE OUTPUT POINT. C
C IFLAG = 3 — INTEGRATION WAS INTERRUPTED BECAUSE THE C
C REQUESTED ERROR TOLERANCES WERE TOO SMALL FOR C
C THE MACHINE PRECISION. THEY HAVE BEEN INCREASED C
C TO MINIMUM ACCEPTABLE VALUES. C
C IFLAG = 4 — INTEGRATION WAS INTERRUPTED BECAUSE THE C
C REQUESTED ERROR TOLERANCES COULD NOT BE C
C ATTAINED AT THE SMALLEST ALLOWABLE STEPSIZE. C
C THE TOLERANCES HAVE BEEN INCREASED TO VALUES C
C WHICH SHOULD ALLOW SUCCESSFUL CONTINUATION. C
C IFLAG = 5 — INTEGRATION WAS INTERRUPTED AFTER THE COUNTER C
C OF CALLS TO F EXCEEDED MAXNFE. THE COUNTER HAS C
C BEEN REDUCED BY MAXNFE. C
C IFLAG = 6 — SAME AS IFLAG=5/ EXCEPT THAT THE PROBLEM HAS C
C BEEN ADDITIONALLY DIAGNOSED AS STIFF. C
C IFLAG = 7 — INTEGRATION WAS NOT CONTINUED BECAUSE THE C
C SPACING OF OUTPUT POINTS HAS BEEN SO SMALL AS TO C
C RESTRICT THE STEPSIZE FOR 50 CALLS TO STEP/ C
C INDICATING THAT THE CODE IS BEING USED IN A WAY C
C WHICH GREATLY IMPAIRS ITS EFFICIENCY. THE C
C MONITOR OF THIS EFFECT WAS RESET TO ALLOW C
C CONTINUATION IF DESIRED. C
C IFLAG = 8 — INTEGRATION WAS INTERRUPTED WHEN A COMPONENT OF C
C THE SOLUTION VANISHED WHILE A PURE RELATIVE C
C ERROR TEST WAS SPECIFIED. THE USER MUST SUPPLY C
C A NONZERO ABSERR TO CONTINUE INTEGRATION. C
C IFLAG = 9 — INVALID INPUT WAS PASSED TO ADAM. THE USER MUST C
C LOCATE AND CORRECT THE ERROR IN ORDER TO C
C CONTINUE INTEGRATION. C
C IFLAG =-3/ -4/ -5, -6/ OR -7 — SAME AS IFLAG = 3/ 4, 5, 6/ C
C OR 7/ RESPECTIVELY, EXCEPT THAT ADAM WAS CALLED C
C WITH IFLAG NEGATIVE. C
C C
c****************************************************************c
c c
C THIS SUBROUTINE IS WRITTEN IN SINGLE PRECISION. C
C WITH MOST COMPILERS/ A DOUBLE PRECISION VERSION/ FOR USE WITH C
C DOUBLE PRECISION VERSIONS OF SUBROUTINES STEP AND SINTRP/ MAY C
C BE OBTAINED SIMPLY BY SETTING THE APPROPRIATE VALUE FOR FOURU C
C IN THE DATA STATEMENT BELOW AND BY REMOVING THE C FROM C
C COLUMN 1 OF THE NEXT SEVEN CARDS. C
C DOUBLE PRECISION ABS/ABSDT/ABSEPS/ABSERR/AMAX1/DABS/DMAX1/
C 1 DSIGN/DT/DTSGN/D1/D2/EPS/FOURU/H/HOLD/
C 2 PHI/RELEPS/RELERR/RVAR/SIGN/T/TEND/TOUT/
C 3 WT/X/Y/YP/YY
22400
O^ c nn
22500
22600
22700
22800
22900
23000
23100
23200
23300
23400
23500
23600
23700
23800
23900
24000
24100
24200
24300
24400
24500
24600
24700
24800
24900
25000
25100
25200
25300
25400
25500
25600
25700
25800
25900
26000
26100
26200
26300
26400
26500
26600
26700
26800
26900
27000
27100
27200
27300
27400
27500
93
-------
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c*
C
c
c
c
ABS(D1) = DABS(D1)
SIGN(D1/D2) = DSIGN(D1/D2)
AMAX1CD1/D2) = DMAXKD1/D2)
LOGICAL NORND/PHASE1 /START
DIMENSION Y(NEQN)/YP(NEQN)/YY(NEQN)/WT(NEQN)/PHI(NEQN/16>
COMMON /ADDATA/ T/TOUT/RELERR/ABSERR/DTSGN/NEQN/IFLAG/
1 IOLD/KLE4/INIT/KOP
COMMON /STDATA/ X/H/HOLD/EPS/RVAR<86)/NEQ,KOLD/ICRASH/K/
1 NFE/NS/NORND/PHASE1/START
EXTERNAL F
IN THE FOLLOWING DATA STATEMENT/ FOURU MUST BE SET TO A
VALUE APPROPRIATE TO THE MACHINE BEING USED.
FOURU IS 4.*U/ WHERE THE COMPUTER UNIT ROUNDOFF ERROR U IS
THE SMALLEST POSITIVE VALUE REPRESENTABLE IN THE MACHINE
SUCH THAT (1.+U).GT.1.
SOME REPRESENTATIVE VALUES ARE
U = 9.54E-7 FOURU = 3.82E-6 IBM 360/370 S.P.
U = 5.96E-8 FOURU = 2.39E-7 PDP-11 S.P.
U = 1.49E-8 FOURU = 5.97E-8 UNIVAC 1108
U = 7.45E-9 FOURU = 2.99E-8 PDP-10
U = 7.11E-15 FOURU = 2.85E-14 CDC 6000 SERIES
U = 2.22D-16 FOURU = 8.89D-16 IBM 360/370 D.P.
U = 1.39D-17 FOURU = 5.56D-17 PDP-11 D.P.
DATA FOURU/2.39E-7/
THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER
OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE.
AS SET/ THIS CORRESPONDS TO ABOUT 1500 STEPS.
DATA MAXNFE/3000/
**************************************************************i
TEST FOR IMPROPER PARAMETERS
MFLAG=IABSUFLAG)
IF
-------
NEQ=NEQN
DT=TOUT-T
ABSDT=ABS(DT)
TENO=T + 10.*DT
IF(IFLAG.LT.O) TEND=TOUT
RELEPS=RELERR/EPS
ABSEPS=ABSERR/EPS
C
C IF THIS IS THE FIRST CALL, GO TO INITIALIZATION SECTION.
C
IF (MFLAG.EQ.1) GO TO 150
C
C SINCE THIS IS NOT THE FIRST CALL/ T=TOUT IS NOT VALID
C INPUT.
C
IF (T.EQ.TOUT) GO TO 200
C
C IF INITIALIZATION WAS NOT COMPLETED ON THE FIRST CALL/ DO
C IT NOW.
C
IF (INIT.EQ.O) GO TO 180
C
C IF THE LAST STEP WAS MADE WITH IFLAG NEGATIVE/ OR IF THE
C DIRECTION OF INTEGRATION IS CHANGED/ RESTART THE INTEGRATION. C
C
IF (IOLD .LT. 0) GO TO 160
IF (DTSGN*DT .LT. 0.0) GO TO 160
C
C IF ALREADY PAST OUTPUT POINT/ INTERPOLATE AND RETURN.
C
50 IF(ABS.LT.FOURU*ABS(X)» GO TO 210
C
C TEST FOR TOO MUCH WORK.
C
IFCNFE .GT. MAXNFE) GO TO 230
C
C LIMIT STEPSIZE/ SET WEIGHT VECTOR/ AND TAKE A STEP.
C
IF (ABS(H).LT.ABS(TEND-X)) GO TO 100
KOP=KOP+1
IF (KOP.GT.50) GO TO 290
H=TEND-X
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
110
C
C
C
C
C
C
32800
32900
33000
33100
33200
33300
33400
33500
33600
33700
33800
33900
34000
34100
34200
34300
34400
34500
34600
34700
34800
34900
35000
35100
35200
35300
35400
35500
35600
35700
35800
35900
36000
36100
36200
36300
36400
36500
36600
36700
36800
36900
37000
37100
37200
37300
37400
37500
37600
37700
37800
37900
95
-------
100 DO 110 L=1/NEQN
WT(L)=RELEPS*ABS(YY(L))+ABSEPS
IF (WT(L).EQ.O.O) GO TO 270
100 CONTINUE
C
CALL STEPCF/YY/YP/WT/PHI/Y)
C
C TEST WHETHER STEP WAS SUCCESSFUL.
C
IF (ICRASH.NE.O) GO TO 250
C
C STEP WAS SUCCESSFUL. TEST FOR STIFFNESS.
C
IF (KLE4.GE.50) GO TO 50
KLE4=KLE4+1
IF(KOLD .GT. 4) KLE4=0
GO TO 50
C
C**************************************************************
C
C INITIALIZATION SECTION —
C
C
C ON START AND RESTART/ SET WORK VARIABLES X AND YY(*>/ STORE
C THE DIRECTION OF INTEGRATION/ AND INITIALIZE THE STEPSIZE.
C
150 NFE=0
160 START=.TRUE.
KLE4=0
INIT=0
KOP=0
X=T
DO 170 L=1/NEQN
170 YY(L)=Y(L)
C
C IF T=TOUT ON FIRST CALL/ EVALUATE DERIVATIVES AND RETURN.
C
IF (T.NE.TOUT) GO TO 180
CALL F(X/YY/YP)
NFE=NFE+1
IFLAG=2
RETURN
C
C SET STEPSIZE AND TRANSFER TO BEGIN INTEGRATION.
C
180 INIT=1
DTSGN=DT
H=SIGN CAMAX1(ABS(TOUT-X)/FOURU*ABS(X))/TOUT-X)
GO TO 50
C
C****************************************************************c
96
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
.
c
c
c
c
*c
38000
38100
38200
38300
38400
38500
38600
38700
38800
38900
39000
39100
39200
39300
39400
39500
39600
39700
39800
39900
40000
40100
40200
40300
40400
40500
40600
40700
40800
40900
41000
41100
41200
41300
41400
41500
41600
41700
41800
41900
42000
42100
42200
42300
42400
42500
42600
42700
42800
42900
43000
43100
-------
c
C THE REMAINDER OF THE PROGRAM CONTAINS STATEMENTS TO HANDLE
C VARIOUS EXCEPTIONAL OCCURRENCES.
C
C
C INVALID INPUT WAS PASSED TO ADAM.
C
200 IFLAG=9
RETURN
C
C CANNOT PASS TOUT, BUT TOO CLOSE TO CALL STEP. EXTRAPOLATE
C BY EULER METHOD AND RETURN.
C
210 H=TOUT-X
CALL FCX,YY,YP)
NFE=NFE+1
DO 220 L=1,NEQN
220 Y(L)=YY(L)+H*YP(L)
IFLAG=2
T=TOUT
CALL F(T,Y/YP)
NFE=NFE+1
IOLD=-1
RETURN
C
C NUMBER OF FUNCTION CALLS EXCEEDED MAXNFE. SET IFLAG TO
C INDICATE WHETHER THIS RESULTS FROM STIFFNESS, AND RESET
C COUNTER NFE.
C
230 IFLAG=ISIGN(5,IFLAG)
IFCKLE4.GE.50) IFLAG=ISIGN(6/IFLAG)
NFE=NFE-MAXNFE
KLE4=0
GO TO 300
C
C ERROR TOLERANCES WERE TOO SMALL TO PERMIT A SUCCESSFUL STEP.
C INCREASE RELERR AND ABSERR.
C
250 IFLAG=ISIGN(4,IFLAG>
IF (ICRASH.EQ.2) IFLAG=ISIGN(3/IFLAG)
RELERR=EPS*RELEPS
ABSERR=EPS*ABSEPS
GO TO 300
C
C A COMPONENT OF THE SOLUTION VANISHED WHEN A PURE RELATIVE
C ERROR TEST WAS SPECIFIED. USER MUST SET ABSERR POSITIVE
C IN ORDER TO PROCEED.
C
270 IFLAG=8
GO TO 300
C
C THE OUTPUT INTERVAL HAS BEEN TOO SMALL TO ALLOW EFFICIENT
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
43200
43300
43400
43500
43600
43700
43800
43900
44000
44100
44200
44300
44400
44500
44600
44700
44800
44900
45000
45100
45200
45300
45400
45500
45600
45700
45800
45900
46000
46100
46200
46300
46400
46500
46600
46700
46800
46900
47000
47100
47200
47300
47400
47500
47600
47700
47800
47900
48000
48100
48200
48300
97
-------
c
c
c
c
c
c
c
c
c
c
c
c
OPERATION OF ADAM.
290 KOP=0
IFLAG=ISIGN(7/IFLAG)
FOR ERROR RETURN/ SET T/ Y<*)/ AND YP(*) TO THE MOST RECENT
VALUES REACHED IN THE INTEGRATION.
300 T=X
IOLD=1
DO 310 L=1/NEQN
Y(L)=YY(L)
310 YP(L)=PHI(L/1)
RETURN
C C
c****************************************************************c
c c
C END OF ADAM C
C C
END
48400
48500
48600
48700
48800
48900
49000
49100
49200
49300
49400
49500
49600
49700
49800
49900
50000
50100
50200
50300
SUBROUTINE STEP USES A MODIFIED DIVIDED DIFFERENCE FORM OF
AN ADAMS PECE METHOD (OF ORDER.LE.12) TO PERFORM ONE STEP
IN THE INTEGRATION OF A SYSTEM OF NEQN FIRST ORDER ORDINARY
DIFFERENTIAL EQUATIONS
DY/DX = F(X/Y) / Y = (Y(1)/Y(2)/.../Y(NEQN))
LOCAL EXTRAPOLATION IS USED TO IMPROVE STABILITY AND ACCURACY,
TO SPECIFY THE DIFFERENTIAL EQUATIONS TO BE SOLVED/ THE USER
MUST SUPPLY A SUBROUTINE/ DENOTED HERE BY F/ TO
CALCULATE THE NEQN DERIVATIVES DY(I)/DX.
F MUST BE OF THE FORM FCX/Y/YP)/ WHERE Y AND YP ARE
ARRAYS EACH OF DIMENSION NEQN. F MUST CALCULATE THE
VALUES OF THE DERIVATIVES DYCD/DX AT (X/Y) AND PLACE
THEM IN YP. X AND Y MUST NOT BE ALTERED BY F.
INFORMATION IS EXCHANGED WITH STEP BOTH THROUGH THE ARGUMENT
LIST AND THROUGH THE COMMON BLOCK STDATA.
THE PARAMETERS IN THE ARGUMENT LIST ARE
F -- THE NAME OF THE USER-SUPPLIED SUBROUTINE TO
CALCULATE THE DERIVATIVE
Y — AN ARRAY OF DIMENSION NEQN WHICH CONTAINS THE
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
10000
10100
10200
10300
10400
10500
10600
10700
10800
10900
11000
11100
11200
11300
11400
11500
11600
11700
11800
11900
12000
12100
12200
12300
12400
12500
12600
12700
12800
12900
98
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
SOLUTION. ITS VALUE MUST NOT BE ALTERED BETWEEN
CALLS TO STEP.
YP — AN ARRAY OF DIMENSION NEQN WHICH CONTAINS THE
DERIVATIVE OF THE SOLUTION AFTER A SUCCESSFUL STEP.
ITS VALUE NEED NOT BE PRESERVED BETWEEN STEPS.
WT — AN ARRAY OF DIMENSION NEQN USED IN SPECIFYING THE
ERROR TEST.
PHI — AN ARRAY OF DIMENSIONS (NEQN/16) WHICH CONTAINS
THE MODIFIED DIVIDED DIFFERENCES. PHI(*/15) AND
PHI<*,16) ARE USED FOR PROPAGATED ROUNDOFF CONTROL.
PHI MUST NOT BE CHANGED BETWEEN CALLS.
P — AN ARRAY OF DIMENSION NEQN USED FOR WORKING
STORAGE. ITS VALUE IS NOT NEEDED FOR SUBSEQUENT
CALLS.
THE PARAMETERS IN THE COMMON BLOCK STDATA ARE
X — THE INDEPENDENT VARIABLE
H — THE STEPSIZE TO BE USED IN ATTEMPTING THE NEXT STEP
HOLD — THE STEPSIZE USED FOR THE LAST SUCCESSFUL STEP
EPS — THE DESIRED TOLERANCE FOR THE ERROR TEST
NEQN — THE NUMBER OF EQUATIONS BEING INTEGRATED
KOLD — THE ORDER OF THE METHOD USED FOR THE LAST
SUCCESSFUL STEP
ICRASH — AN INDICATOR OF ERROR CONDITIONS UPON RETURN
K __ THE ORDER OF THE METHOD TO BE USED FOR THE NEXT STEP
NFE — A COUNTER WHICH IS INCREMENTED BY 1 AFTER EACH
CALL TO F
START — A LOGICAL VARIABLE WHICH MUST BE SET .TRUE. ON
THE FIRST CALL TO STEP
ALPHA(12),BETA(12),G<13),SIG(13),V(12),W(12)/
PSIC12),NS/NORND,PHASE1 — QUANTITIES USED FOR INTERNAL
BOOKKEEPING BY STEP WHOSE VALUES ARE NEEDED FOR
CONTINUING' INTEGRATION. THEIR VALUES ARE NORMALLY OF
NO INTEREST TO THE USER. THEY ARE PLACED IN THE
COMMON BLOCK TO AVOID LOCAL RETENTION BETWEEN CALLS.
NORND AND PHASE1 ARE LOGICAL VARIABLES.
EXCEPT FOR H/ EPS, AND NFE* THE VALUES IN THE COMMON BLOCK
MUST NOT BE ALTERED BETWEEN CALLS TO STEP.
THE COUNTER NFE IS PROVIDED SO THAT THE USER CAN MONITOR THE
AMOUNT OF WORK BEING DONE. ITS VALUE HAS NO MEANING TO STEP.
STEP AUTOMATICALLY ADJUSTS THE STEPSIZE AND THE ORDER OF THE
METHOD SO THAT THE MAX NORM OF THE VECTOR WITH COMPONENTS
LOCAL ERROR(I)/WT(I) IS LESS THAN EPS FOR EACH STEP. THE
ARRAY WT ALLOWS THE USER TO SPECIFY AN ERROR TEST WHICH IS
APPROPRIATE FOR HIS PROBLEM.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
13000
13100
13200
13300
13400
13500
13600
13700
13800
13900
14000
14100
14200
14300
14400
14500
14600
14700
14800
14900
15000
15100
15200
15300
15400
15500
15600
15700
15800
15900
16000
16100
16200
16300
16400
16500
16600
16700
16800
16900
17000
17100
17200
17300
i?Ann
1 f *rUw
17500
17600
17700
17800
17900
18000
18100
99
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
FOR EXAMPLE,
WTU) = 1.0 SPECIFIES ABSOLUTE ERROR
WTU) = ABS(YU)) SPECIFIES ERROR RELATIVE TO THE MOST
RECENT VALUE OF THE I-TH COMPONENT OF THE
SOLUTION
WTU) = AMAX1(WTU),ABS(YU))) SPECIFIES ERROR RELATIVE
TO THE LARGEST MAGNITUDE YET ACHIEVED BY THE
I-TH COMPONENT OF THE SOLUTION
WTU) = ABS(YU))*RELERR/EPS+ABSERR/EPS, WHERE EPS =?
AMAX1(RELERR,ABSERR), SPECIFIES A MIXED ERROR
TEST WHERE RELERR IS RELATIVE ERROR AND
ABSERR IS ABSOLUTE ERROR
ALL COMPONENTS OF WT MUST BE NONZERO.
UPON RETURN, THE VALUE OF ICRASH INDICATES THE ACTION WHICH
WAS TAKEN.
ICRASH = 0 — A STEP OF LENGTH HOLD WAS SUCCESSFULLY
TAKEN TO THE VALUE NOW IN X. THE SOLUTION AT X
IS IN Y, AND YP CONTAINS THE DERIVATIVES. THE
DERIVATIVES ARE ALSO IN PHI(*,1).
ICRASH = 1 — NO STEP WAS TAKEN BECAUSE H WAS SMALLER THAN
PERMITTED BY THE MACHINE PRECISION. H WAS INCREASED
TO AN ACCEPTABLE VALUE.
ICRASH = 2 — NO STEP WAS TAKEN BECAUSE EPS WAS TOO SMALL
FOR THE MACHINE PRECISION. EPS WAS INCREASED TO AN
ACCEPTABLE VALUE.
ICRASH = 3 -- NO STEP WAS TAKEN BECAUSE THE ERROR TEST
COULD NOT BE SATISFIED AT THE SMALLEST ALLOWED STEPSIZE.
EPS IS INCREASED TO A VALUE WHICH IS ACHIEVABLE.
FIRST CALL TO STEP
TO BEGIN THE SOLUTION, THE CALLING PROGRAM MUST
1. PROVIDE STORAGE FOR THE ARRAYS Y(NEQN), YP(NEQN),
WT(NEQN), PHI(NEQN,16), AND P(NEQN).
2. PLACE THE INITIAL VALUES OF THE SOLUTION IN Y.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
18200
18300
18400
18500
18600
18700
18800
18900
19000
19100
19200
19300
19400
19500
19600
19700
19800
19900
20000
20200
20300
20400
20500
20600
20700
20800
20900
21000
21100
21200
21300
21400
21500
21600
21700
21800
21900
22000
22100
22200
oo^nn
bbOUU
22400
22500
22700
22800
22900
23000
23100
23200
23300
100
-------
c
C 3. SET WTC1)/ .../ WT(NEQN) TO NONZERO VALUES TO PROVIDE
C THE ERROR TEST DESIRED.
C
C 4. INITIALIZE THE FOLLOWING VARIABLES IN THE COMMON BLOCK
C STDATA
C X — THE STARTING VALUE OF THE INDEPENDENT VARIABLE
C H — A NOMINAL STEPSIZE FOR THE FIRST STEP/
C INDICATING THE DIRECTION OF INTEGRATION AND THE
C MAXIMUM STEPSIZE. H WILL BE REDUCED AS NECESSARY
C TO MEET THE ERROR CRITERION.
C EPS — THE DESIRED ERROR TOLERANCE
C NEQN — THE NUMBER OF EQUATIONS '
C NFE — SET TO 0 TO BEGIN A COUNT OF FUNCTION CALLS
C START — SET .TRUE.
C
C 5. PLACE THE ACTUAL NAME OF THE USER SUPPLIED SUBROUTINE
C F(X/Y/YP) IN THE CALL LIST/ AND DECLARE THAT NAME
C IN AN EXTERNAL STATEMENT.
C
C
C SUBSEQUENT CALLS TO STEP
C
C STEP RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE THE
C INTEGRATION. TO MAINTAIN RELATIVE ERROR TESTS AS DESCRIBED
C ABOVE/ WT(*) MUST BE UPDATED AFTER EACH STEP. CALLING STEP
C AGAIN THEN ADVANCES THE SOLUTION ANOTHER STEP. TO CONTINUE
C INTEGRATION USING THE REVISED STEPSIZE OR TOLERANCE
C PROVIDED WITH AN ERROR RETURN (ICRASH.NE.O)/ IT IS ONLY
C NECESSARY TO CALL STEP AGAIN.
C
C NORMALLY/ THE INTEGRATION IS CONTINUED TO THE FIRST STEP
C BEYOND THE DESIRED OUTPUT POINT/ AND THE SOLUTION IS INTER-
C POLATED TO THE OUTPUT POINT USING THE SUBROUTINE SINTRP.
C IF IT IS IMPOSSIBLE TO INTEGRATE BEYOND THE OUTPUT POINT/ H
C MAY BE REDUCED TO HIT THE OUTPUT POINT. OTHERWISE/ H SHOULD
C NOT BE ALTERED/ SINCE CHANGING H IMPAIRS THE EFFICIENCY OF
C THE CODE. H MUST NEVER BE INCREASED. EPS MAY BE CHANGED
C BETWEEN CALLS/ BUT LARGE DECREASES SHOULD BE AVOIDED/ SINCE
C THIS WILL INVALIDATE THE STEP SELECTION MECHANISM AND LEAD
C TO STEP FAILURES.
C
c***************************************************************
C
C THIS SUBROUTINE IS WRITTEN IN SINGLE PRECISION.
C WITH MOST COMPILERS/ A DOUBLE PRECISION VERSION MAY BE
C OBTAINED BY SETTING THE APPROPRIATE VALUE FOR FOURU IN THE
C DATA STATEMENT BELOW/ AND BY REMOVING THE C FROM COLUMN 1 OF
C THE NEXT ELEVEN CARDS.
C DOUBLE PRECISION ABS/ABSH/ALPHA/AMAX1/AMIN1 /BETA/DABS/
C 1 DMAX1/DMIN1/DSIGN/DSQRT/DUM1/DUM2/EPS/
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
•c
c
c
c
c
c
*c
c
c
c
c
c
c
23400
23500
23600
23700
23800
23900
24000
24100
24200
24300
24400
24500
24600
24700
24800
24900
25000
25100
25200
25300
25400
25500
"%c^nn
cbOUU
25700
25800
25900
26000
26100
26200
26300
26400
26500
26600
26700
26800
26900
27000
27100
27200
27300
27400
27500
27600
27700
27800
27900
28000
28100
28200
28300
28400
28500
101
-------
C 2 ERK,ERKM1,ERKM2,ERKP1, ERR, FOURU, G,GSTR,
C 3 H,HNEW,HOLD,P,PHI,PSI,R,REALNS,RHO,ROUND,
C 4 SIG,SIGN,SQRT,SUM,TAU,TEMP1,TEMP2,TEMP3,
C 5 TEMP4,TEMP5,TEMP6,TWO,V,W,WT,X,XOLD,Y,YP
C ABS(DUM1)=DABS(DUM1)
C AMAX1(DUM1,DUM2)=DMAX1(DUM1,DUM2)
C AMIN1 (DUM1,DUM2)=DMIN1 (DUM1,DUM2)
C SIGNCDUM1,DUM2)=DSIGN(DUM1,DUM2)
C SQRT(DUM1)=DSQRT(DUM1)
C
LOGICAL START, PHASE1,NORND
C
EXTERNAL F
C
COMMON /STDATA/ X,H,HOLD,EPS,ALPHA(12),BETA(12),GC13>,
1 SIG(13),V(12),W(12),PSI(12),NEQN,KOLD,
2 ICRASH,K,NFE,NS,NORND,PHASE1 , START
C
DIMENSION Y(NEQN),WTCNEQN),PHI(NEQN,16),P(NEQN),
1 YP(NEQN),GSTR(13),TWO(12)
C
DATA TWO/8.,16.,32.,64.,128.,256.,512.,1024.,2048.,4096.,
1 8192., 16384. /
C
DATA GSTR/. 5,. 0833,. 0417,. 0264,. 0188,. 0143,. 01 14,
1 .00936, .00789, .00697, .00592, .00524, .004687
C
C
C IN THE FOLLOWING DATA STATEMENT, FOURU MUST BE SET TO A
C VALUE APPROPRIATE TO THE MACHINE BEING USED.
C
C FOURU IS 4.*U, WHERE THE COMPUTER UNIT ROUNDOFF ERROR U IS
C THE SMALLEST POSITIVE VALUE REPRESENTABLE IN THE MACHINE
C SUCH THAT (1,+U).GT.1.
C
C SOME REPRESENTATIVE VALUES ARE
C U = 9.54E-7 FOURU = 3.82E-6 IBM 360/370 S.P.
C U t 5.96E-8 FOURU = 2.39E-7 PDP-11 S.P.
C U = 1.49E-8 FOURU = 5.97E-8 UNIVAC 1108
C i U = 7.45E-9 FOURU = 2.99E-8 PDP-10
C U = 7.11E-15 FOURU = 2.85E-14 CDC 6000 SERIES
C U = 2.22D-16 FOURU = 8.89D-16 IBM 360/370 D.P.
C U = 1.39D-17 FOURU = 5.56D-17 PDP-11 D.P.
C
DATA FOURU/2.39E-7/
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
C****************************************************************C
c ******
C BEGIN BLOCK 0 ******
c ******
c***********************
C
c
c
c
c
c
28600
28700
28800
28900
29000
29100
29200
29300
29400
29500
29600
29700
29800
29900
30000
30100
30200
30300
30400
30500
30600
30700
30800
30900
31000
31100
31200
31300
31400
31500
31600
31700
31800
31900
32000
32100
32200
32300
32400
32500
32600
32700
32800
32900
33000
33100
33200
33300
33400
33500
33600
33700
102
-------
C CHECK WHETHER STEPSIZE OR ERROR TOLERANCE IS TOO SMALL FOR C 33800
C MACHINE PRECISION. IF FIRST STEP/ INITIALIZE PHI ARRAY AND C 33900
C ESTIMATE A STARTING STEPSIZE. C 34000
C C 34100
C C 34200
C IF STEPSIZE IS TOO SMALL/ DETERMINE AN ACCEPTABLE ONE. C 34300
C C 34400
IF (ABS(H) .GE. FOURU*ABS(X)) GO TO 5 34500
H=SIGN(FOURU*ABS(X),H) 34600
ICRASH=1 34700
RETURN 34800
C C 34900
C IF ERROR TOLERANCE IS TOO SMALL/ INCREASE IT TO AN ACCEPTABLE C 35000
C VALUE. C 35100
C C 35200
5 ROUND=0. 35300
DO 10 L=1/NEQN 35400
10 ROUND=AMAX1(ROUND,ABS(Y(L)/WT(L))) 35500
ROUND=FOURU*ROUND 35600
IF (EPS .GE. ROUND) GO TO 15 35700
EPS=ROUND*(1.+FOURU) 35800
ICRASH=2 35900
RETURN 36000
C C 36100
15 ICRASH=0 36200
G(1)=1.0 36300
G(2)=0.5 36400
SIG(1)=1.0 36500
IFAIL=0 36600
IF (.NOT. START) GO TO 100 36700
C C 36800
C INITIALIZE. COMPUTE APPROPRIATE STEPSIZE FOR FIRST STEP. C 36900
C C 37000
CALL F(X,Y,YP) 37100
NFE=NFE+1 37200
SUM=0. 37300
DO 20 L=1/NEQN 37400
PHI(L/1)=YP(L) 37500
PHI(L,2)=0. 37600
20 SUM=AMAX1(SUM/ABS(YP(L)/WT(L))) 37700
C C 37800
ABSH=ABS(H) 37900
IF (EPS .LT. 16.*SUM*H*H) ABSH=0.25*SQRT(EPS/SUM) 38000
H=SIGN(AMAX1(ABSH,FOURU*ABS(X)),H) 38100
HOLD=0. 38200
K=1 38300
KOLD=0 38400
START=.FALSE. 38500
PHASE1=.TRUE. 38600
NORND=.TRUE. 38700
IF (EPS .GT. 200.*ROUND) GO TO 100 38800
NORND=.FALSE. 38900
103
-------
DO 25 L=1/NEQN
25 PHI(L,15)=0.
C
c***********************
C ******
C END BLOCK 0 ******
C ******
C
C
C
C
C
c****************************************************************c
c ******
C BEGIN BLOCK 1 ******
C ******
c***********************
c
C COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP. AVOID
C COMPUTING THOSE QUANTITIES NOT CHANGED WHEN STEPSIZE IS NOT
C CHANGED.
C
100 KP1=K+1
KP2=K+2
KM1=K-1
KM2=K-2
C .*
C NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, INCLUDING THE
C CURRENT ONE. WHEN J<.LT.NS/ NO COEFFICIENTS CHANGE.
C
IF (H .NE. HOLD) NS=0
NS=MINO(NS+1xKOLD+1)
NSP1=NS+1
IF (K .LT. NS) GO TO 200
C
C COMPUTE THOSE COMPONENTS OF ALPHA(*)/ BETA(*), PSK*), AND
C SIGO) WHICH ARE CHANGED.
C
BETA(NS)=1.
REALNS=NS
ALPHA(NS)=1./REALNS
TEMPI =H*REALNS
SIG(NSP1)=1. v
IF (K .LT. NSP1) GO TO 110
DO 105 I=NSP1/K
IM1=I-1
TEMP2=PSI(IM1)
PSI(IM1)=TEMP1
BETA(I)=BETA(IM1)*PSI(IM1)/TEMP2
TEMPI =TEMP2+H
ALPHA(I)=H/TEMP1
105 516(1+1 )=FLOAT(I)*ALPHA(I)*SIG(I>
110 PSI(K)=TEMP1
C
C COMPUTE COEFFICIENTS G(*)
C
C INITIALIZE V(*) AND SET W<*).
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
39000
39100
39200
39300
39400
39500
39600
39700
39800
39900
40000
40100
40200
40300
40400
40500
40600
40700
40800
40900
41000
41100
41200
41300
41400
41500
41600
41700
41800
41900
42000
42100
42200
42300
42400
42500
42600
42700
42800
42900
43000
43100
43200
43300
43400
43500
43600
43700
43800
43900
44000
44100
104
-------
115
IF (NS .GT. 1) GO TO 120
DO 115 IQ=1,K
TEMP3=IQ*(IQ+1)
V(IQ)=1./TEMP3
WUQ)=V(IQ)
GO TO 140
C
C
C
C
C
C
IF ORDER WAS RAISED/ UPDATE DIAGONAL PART OF B(*),
120 IF (K .LE. KOLD) GO TO 130
TEMP4=K*KP1
V(K)=1./TEMP4
NSM2=NS-2
IF (NSM2 .LT. 1) GO TO 130
DO 125 J=1,NSM2
I=K-J
125 V(I)=V
G(NSP1)=W(1)
C
C
C
COMPUTE THE G<*) IN THE WORK VECTOR WC*),
140 NSP2=NS+2
IF (KP1 .LT. NSP2) GO TO 200
DO 150 I=NSP2/KP1
LIMIT2=KP2-I
TEMP6=ALPHA(I-1)
DO 145 IQ=1/LIMIT2
145 W(IQ)=W(IQ)-TEMP6*W(IQ+1)
150 G(I)=W(1)
C***********************
c ******
C END BLOCK 1 ******
C ******
C****************************************************************C
C ******
C BEGIN BLOCK 2 ******
C ******
c***********************
c
C PREDICT A SOLUTION P(*>/ EVALUATE DERIVATIVES USING PREDICTED
C SOLUTION, AND ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT
C ORDERS K/ K-1, K-2 AS IF CONSTANT STEPSIZE WERE USED.
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
44200
44300
44400
44500
44600
44700
44800
44900
45000
45100
45200
45300
45400
45500
45600
45700
45800
45900
46000
46100
46200
46300
46400
46500
46600
46700
46800
46900
47000
47100
47200
47300
47400
47500
47600
47700
47800
47900
48000
48100
48200
48300
48400
48500
48600
48700
48800
48900
49000
49100
49200
49300
105
-------
C C 49400
C C 49500
C CHANGE PHI TO PHI STAR. C 49600
C C 49700
200 IF (K .LT. NSP1) GO TO 215 49800
DO 210 I=NSP1/K 49900
TEMP1=BETA(I) 50000
DO 205 L=1/NEQN 50100
205 PHI(L/I)=TEMP1*PHI(L/I) 50200
210 CONTINUE 50300
C C 50400
C PREDICT SOLUTION AND DIFFERENCES. C 50500
C C 50600
215 DO 220 L=1/NEQN 50700
PHI(L/KP2)=PHI(L/KP1> 50800
PHI(L/KP1>=0. 50900
220 P(L)=0. 51000
DO 230 J=1/K 51100
I=KP1-J 51200
IP1=I+1 51300
TEMP2=G(I) 51400
DO 225 L=1/NEQN 51500
P(L)=P(L)+TEMP2*PHI(L/I) 51600
225 PHI(L/I)=PHI(L/I)+PHI(L/IP1> 51700
230 CONTINUE 51800
IF (NORND) GO TO 240 51900
DO 235 L=1/NEQN 52000
TAU=H*P(L)-PHKL/15) 52100
P(L)=YCL)+TAU 52200
235 PHI(L/16) = (P 54000
IF (KM2) 265/260/255 54100
255 ERKM2=AMAX1(ERKM2/ABS() 54200
260 ERKM1=AMAX1(ERKM1/ABS«PHI(L/K)+TEMP4)*TEMP3» 54300
265 ERK=AMAX1(ERK/ABSCTEMP4*TEMP3)> 54400
IF (KM2) 280/275/270 54500
106
-------
C
c
c
c
c
c
270 ERKM2=ABSH*SIG(KM1)*GSTR(KM2)*ERKM2
275 ERKM1=ABSH*SIG(K)*GSTR(KM1)*ERKM1
280 TEMP5=ABSH*ERK
ERR=TEMP5*(G(K)-G(KP1))
ERK=TEMP5*SIG(KP1)*GSTR(K>
KNEW=K
TEST WHETHER ORDER SHOULD BE LOWERED.
IF (KM2) 299x290/285
285 IF (AMAXKERKM1/ERKM2) .LE. ERK) KNEW=KM1
GO TO 299
290 IF (ERKM1 .LE. 0.5*ERK) KNEW=KM1
TEST WHETHER STEP WAS SUCCESSFUL.
299 IF (ERR .LE. EPS) GO TO 400
c ******
C END BLOCK 2 ******
C ******
c****************************************************************c
c ******
C BEGIN BLOCK 3 ******
C ******
c***********************
c
THE STEP IS UNSUCCESSFUL. RESTORE X/ PHK*/*)/ PSK*).
C
C
C
C
C
C
C
C
C
IF THIRD CONSECUTIVE FAILURE SET ORDER TO ONE. IF STEP FAILS
MORE THAN THREE TIMES/ CONSIDER AN OPTIMAL STEPSIZE. DOUBLE
ERROR TOLERANCE AND RETURN IF ESTIMATED STEPSIZE IS TOO SMALL
FOR MACHINE PRECISION.
RESTORE X/ PHK*/*)/ AND PSI(*).
PHASE1=. FALSE.
, X=XOLD
DO 310 1=1 /K
TEMPI =1./BETAU>
305
310
315
DO 305 L=1/NEQN
PHI(L/I)=TEMP1*(PHI(L/I)-PHI(L/IP1)>
CONTINUE
IF (K .LT. 2) GO TO 320
DO 315 I=2/K
PSIU-1)=PSI(I)-H
C
C
C
C
ON THIRD FAILURE/ SET ORDER TO ONE.
STEPSIZE.
THEREAFTER/ USE OPTIMAL
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
54600
54700
54800
54900
55000
55100
55200
55300
55400
55500
55600
55700
55800
55900
56000
56100
56200
56300
56400
56500
56600
56700
56800
56900
57000
57100
57200
57300
57400
57500
57600
57700
57800
57900
58000
58100
58200
58300
58400
58500
58600
58700
58800
58900
59000
59100
59200
59300
59400
59500
59600
59700
107
-------
320 IFAIL=IFAIL+1
TEMP2=0.5
IF (IFAIL-3) 335/330/325
325 IF (EPS .LT. 0.5*ERK) TEMP2=SQRT(0.5*EPS/ERK)
330 KNEW=1
335 H=TEMP2*H
K=KNEW
IF (ABS(H) .GE. FOURU*ABSCX» GO TO 100
C
ICRASH=3
H=SIGN ( FOURU*ABS (X) /H>
EPS=EPS+EPS
RETURN
C
c***********************
c ******
C END BLOCK 3 ******
C ******
C
C
C
C
C
C
c****************************************************************c
c ******
C BEGIN BLOCK 4 ******
C ******
c***********************
c
C THE STEP IS SUCCESSFUL. CORRECT THE PREDICTED SOLUTION/
C EVALUATE THE DERIVATIVES USING THE CORRECTED SOLUTION/ AND
C UPDATE THE DIFFERENCES/ DETERMINE BEST ORDER AND STEPSIZE
C FOR THE NEXT STEP.
C
400 KOLD=K
HOLD=H
C
C CORRECT AND EVALUATE.
C
TEMPI =H*G(KP1)
IF (NORND) GO TO 410
DO 405 L=1/NEQN
RHO=TEMP1*(YP(L)-PHI(L/1))-PHI(L/16)
Y(L)=P(L)+RHO
405 PHI(L/15)=(Y(L)-P(L))-RHO
GO TO 420
410 DO 415 L=1/NEQN
415 Y(L)=P
-------
C
c
c
c
c
c
DO 430 L=1/NEQN
430 PHI(L/I)=PHI(L/I)+PHI(L/KPD
ESTIMATE ERROR AT ORDER K+1 UNLESS/
IN STARTING PHASE WHEN ORDER IS ALWAYS RAISED/
ALREADY DECIDED IN BLOCK 2 TO LOWER ORDER/
STEPSIZE IS NOT CONSTANT SO ESTIMATE WOULD BE UNRELIABLE.
ERKP1=0.
IF «KNEW.EQ.KM1).OR.(K.EQ.12)> PHASE1=.FALSE.
IF (PHASED GO TO 450
IF (KNEW .EQ. KM1) GO TO 455
IF (KP1 .GT. NS) GO TO 460
DO 440 L=1/NEQN
ERKP1=AMAX1(ERKP1/ABS(PHI(L/KP2)/WT(L)))
440
ERKP1=ABSH*GSTR(KP1)*ERKP1
C
C USING ESTIMATED ERROR AT ORDER K+1/ DETERMINE
C ORDER FOR NEXT STEP.
C
IF (K .GT. 1) GO TO 445
IF (ERKP1 .GE. 0.5*ERK) GO TO 460
GO TO 450
445 IF (ERKM1 .LE. AMIN1 (ERK/ERKPD) GO TO 455
IF ((ERKP1.GE.ERK).OR.(K.EQ.12» GO TO 460
C
C HERE ERKP1.LT.ERK.LT.AMAXKERKM1/ERKM2)/ ELSE
C HAVE BEEN LOWERED IN BLOCK 2. THUS/ ORDER IS
C
450 K=KP1
ERK=ERKP1
GO TO 460
C
C LOWER ORDER.
C
455 K=KM1
ERK=ERKM1
C
C WITH NEW ORDER DETERMINE APPROPRIATE STEPSIZE
C
APPROPRIATE
ORDER WOULD
RAISED.
FOR NEXT ST
460 HNEW=H+H
IF (PHASED GO TO 465
IF (EPS .GE. ERK*TWO(K» GO TO 465
HNEW=H
IF (EPS .GE. 2.*ERK) GO TO 465
R=(0.5*EPS/ERK)**(1./FLOAT(K+1»
IF (R.GT.0.9) R=0.9
IF (R.LT.0.5) R=0.5
HNEW=SIGN(AMAX1(ABSH*R/FOURU*ABS(X))/H)
465 H=HNEW
RETURN
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
65000
65100
65200
65300
65400
65500
65600
65700
65800
65900
66000
66100
66200
66300
66400
66500
66600
66700
66800
66900
67000
67100
67200
67300
67400
67500
67600
67700
67800
67900
68000
68100
68200
68300
68400
68500
68600
68700
68800
68900
69000
69100
69200
69300
69400
69500
69600
69700
69800
- 69900
70000
70100
109
-------
C*********************** C 70200
C ****** C 70300
C END BLOCK 4 ****** C 70400
C ****** C 70500
C****************************************************************C 70600
C C 70700
C END OF STEP C 70800
C C 70900
END 71000
SUBROUTINE SINTRP(XOUT,YOUT/YPOUT,Y/PHI> 10000
C C 10100
C****************************************************************C 10200
C C 10300
C SUBROUTINE SINTRP MAY BE USED IN CONJUNCTION WITH SUBROUTINE C 10400
C STEP TO INTERPOLATE THE SOLUTION TO THE NONMESH POINT XOUT. C 10500
C SINTRP SHOULD BE CALLED WHEN THE INTEGRATION HAS PROCEEDED TO C 10600
C THE FIRST MESH POINT BEYOND XOUT. C 10700
C C 10800
C THE METHODS IN SUBROUTINE STEP APPROXIMATE THE SOLUTION NEAR C 10900
C X BY A POLYNOMIAL. SINTRP APPROXIMATES THE SOLUTION AT C 11000
C XOUT BY EVALUATING THE POLYNOMIAL THERE. C 11100
C C 11200
C THE PARAMETERS IN THE ARGUMENT LIST ARE C 11300
C C 11400
C XOUT — THE VALUE OF THE INDEPENDENT VARIABLE FOR WHICH C 11500
C THE SOLUTION IS TO BE FOUND C 11600
C YOUT — AN ARRAY OF DIMENSION NEQN TO HOLD THE VALUE OF C 11700
C THE SOLUTION AT XOUT. C 11800
C YPOUT — AN ARRAY OF DIMENSION NEQN TO HOLD THE VALUE OF C 11900
C DERIVATIVE AT XOUT. C 12000
C Y — THE ARRAY OF DIMENSION NEQN WHICH CONTAINS THE C 12100
C SOLUTION RETURNED AT X BY STEP C 12200
C PHI — THE ARRAY OF DIMENSIONS (NEQN/16) WHICH HOLDS THE C 12300
C MODIFIED DIVIDED DIFFERENCES RETURNED AT X BY STEP C 12400
C C 12500
C THE REMAINING PARAMETERS X, PSK12)/ NEQN/ AND KOLD WHICH C 12600
C ARE NEEDED FOR DEFINING THE POLYNOMIAL/ ARE PASSED FROM C 12700
C SUBROUTINE STEP THROUGH THE COMMON BLOCK STDATA. C 12800
C C 12900
C****************************************************************C 13000
C C 13100
C THIS PROGRAM IS WRITTEN IN SINGLE PRECISION. TO OBTAIN A C 13200
C DOUBLE PRECISION VERSION/ REMOVE THE C FROM COLUMN 1 OF C 13300
C THE NEXT THREE CARDS. C 13400
C DOUBLE PRECISION XOUT/YOUT/YPOUT/Y/PHI/G/W/RHO/X/RVAR/ 13500
C 1 PSI/HI/TEMP1/TERM/PSIJM1/GAMMA/ETA/ 13600
C 2 TEMP2/TEMP3 13700
C C 13800
LOGICAL LVAR 13900
C C 14000
110
-------
DIMENSION Y(NEQN)/YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16)/ 14100
1 G(13)/WC13),RHO(13) 14200
C C 14300
COMMON /STDATA/ X,RVARC77),PSI(12),NEQN,KOLD,IVAR(4),LVAR(3) 14400
C C 14500
DATA G(1)/1.0//RHO(1>/1.0/ 14600
C C 14700
HI=XOUT-X 14800
KI=KOLD+1 14900
KIP1=KI+1 150QO
C C 15100
C INITIALIZE W(*) FOR COMPUTING G(*). C 15200
C C 15300
DO 5 1=1,KI 15400
TEMPI=1 15500
5 W(I)=1.0/TEMP1 15600
TERM=0.0 15700
C C 15800
C COMPUTE G(*). C 15900
C C 16000
DO 15 J=2,KI 16100
JM1=J-1 16200
PSIJM1=PSI(JM1) 16300
GAMMA=(HI+TERM)/PSIJM1 16400
ETA=HI/PSIJM1 16500
LIMIT1=KIP1-J 16600
DO 10 I=1,LIMIT1 16700
10 WU)=GAMMA*W(I)-ETA*W(I-H) 16800
G(J)=W(1) 16900
RHO(J)=GAMMA*RHO(JM1) 17000
15 TERM=PSIJM1 17100
C C 17200
C INTERPOLATE. C 17300
C C 17400
DO 20 L=1/NEQN 17500
YPOUT(L)=0.0 17600
20 YOUT(L)=0.0 17700
DO 25 J=1/KI 17800
I=KIP1-J 17900
TEMP2=G(I) 18000
TEMP3=RHO(I) 18100
DO 25 L=1/NEQN 18200
YOUT(L)=YOUT(L)+TEMP2*PHKL/I) 18300
25 YPOUTCL)=YPOUT(L)+TEMP3*PHI(L,I) 18400
DO 35 L=1,NEQN 18500
35 YOUT(L)=Y(L)+HI*YOUT(L) 18600
C C 18700
RETURN 18800
C C 18900
C END OF SINTRP C 19000
C C 19100
END 19200
111
-------
APPENDIX C
STIFF METHOD PROGRAMS
This appendix contains the complete listings of the subroutines STFINT,
STIFF, GEAR, DECOMP, GINTRP, and QJACOB.
SUBROUTINE STFINT(NEQN/Y/W/IW/TINIT/TFINAL/TINCR/RELERR/
1 ABSERR/IUNIT/FCT/-FDER/OUTP)
C
c****************************************************************c
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
SUBROUTINE STFINT USES STIFF TO INTEGRATE A SYSTEM OF NEQN
FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS
DY/DT = F(T/Y) , Y =
-------
FCT — THE NAME OF THE USER-SUPPLIED SUBROUTINE
FCTCT/Y/YP) WHICH EVALUATES THE DERIVATIVES
YP(I)=DY(I)/DT/ I=1/2/.../NEQN
FDER — THE NAME OF THE USER SUPPLIED SUBROUTINE WHICH
EVALUATES THE JACOBIAN MATRIX AS DESCRIBED IN THE
COMMENTS TO GEAR.
OUTP ~ THE NAME OF THE USER-SUPPLIED SUBROUTINE
OUTPCT/Y/W) USED TO HANDLE OUTPUT.
TO USE STFINT/ THE CALLING PROGRAM MUST
1. ALLOCATE SPACE FOR THE ARRAYS Y(*)/ IW(*)/ AND W(*) IN
DIMENSION STATEMENT. (NOTE THAT Y MUST BE AN ARRAY/
EVEN IF NEQN=1.)
2. SET Y(1)/Y(2)/.../Y(NEQN) TO THE INITIAL VALUES TO BE
USED FOR THE SOLUTION.
3. PROVIDE APPROPRIATE VALUES FOR NEQN/ TINIT/ TFINAL/
TINCR/ RELERR/ ABSERR/ AND IUNIT. (SINCE THESE
QUANTITIES ARE NOT ALTERED BY STFINT/ THEY MAY/ IF
DESIRED/ BE WRITTEN AS CONSTANTS IN THE CALL LIST.)
4. PLACE THE ACTUAL NAMES OF FCT/ FDER/ AND OUTP IN THE
CALL LIST/ AND DECLARE THOSE NAMES IN AN EXTERNAL
STATEMENT.
AT EACH CALL TO OUTP(T/Y/W)/ THE CURRENT VALUES OF THE
DERIVATIVES DY(I)/DT ARE IN W(1)/W(2)/.../W(NEQN). THE
FIRST CALL TO OUTP IS MADE AT TINIT.
IT IS NOT NECESSARY THAT THE INTEGRATION INTERVAL BE EVENLY
DIVISIBLE BY THE OUTPUT MESH SIZE TINCR. HOWEVER/ IF IT IS
NOT/ INTEGRATION WILL PROCEED TO THE FIRST MESH POINT
BEYOND TFINAL/ RATHER THAN STOPPING AT TFINAL.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c****************************************************************c
C
C
C
C
C
C
C
C
C
C
C
THIS SUBROUTINE IS WRITTEN IN SINGLE PRECISION.
A DOUBLE PRECISION VERSION FOR USE WITH A DOUBLE PRECISION
VERSION OF STIFF MAY BE OBTAINED SIMPLY BY REMOVING THE C
FROM COLUMN 1 OF THE NEXT FIVE CARDS.
DOUBLE PRECISION ABS/ABSER/ABSERR/DABS/DSIGN/DT/DTSGN/
1 D1/D2/EPSDT/RELER/RELERR/SIGN/T/TFINAL/
2 TINCR/TINIT/TMAXIN/TOUT/W/Y
ABSCD1) = DABSCD1)
SIGNCD1/D2) = DSIGN(D1/D2)
COMMON /SFDATA/ T/TOUT/RELER/ABSER/DTSGN/NEQ/IFLAG/IOLD/
C
C
C
C
C
C
C
C
C
C
C
C
A C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
**c
C
C
C
C
C
C
13700
13800
13900
14000
14100
14200
14300
14400
14500
14600
14700
14800
14900
15000
15100
15200
15300
15400
15500
15600
15700
15800
15900
16000
16100
16200
16300
16400
16500
16600
16700
16800
16900
17000
17100
17200
17300
17400
17500
17600
17700
17800
17900
18000
18100
18200
18300
18400
18500
18600
18700
18800
113
-------
1 INIT/KOP
C
DIMENSION Y(NEQN)/IW(NEQN)/W(1)
C
C THE ACTUAL DIMENSION OF W/ AS DEFINED BY THE CALLING PROGRAM/ C
C MUST BE AT LEAST NEQN**2+17*NEQN.
C
C
EXTERNAL FCT/FDER
C
C WRITE STARTING DATA INTO MESSAGE FILE
C
IF (IUNIT.GT.O) WRITEUUNIT/100) NEQN/TINIT/TFINAL/TINCR/
1 RELERR/ABSERR
C
C COMPUTE INDICES FOR SPLITTING THE WORK ARRAY W<*>.
C
K1=NEQN+1
K2=K1+9*NEQN
K3=K2+7*NEQN
C
C
INITIALIZE VARIABLES TO START INTEGRATION WITH STIFF.
NEQ=NEQN
RELER=RELERR
ABSER=ABSERR
IFLAG=1
T=TINIT
TOUT=TINIT
C SAVE INFORMATION NEEDED TO LOCATE END OF INTEGRATION.
C
TMAXIN=ABS(TFINAL-TINIT)
EPSDT=0.001*ABS(TINCR>
C
C SET SIGN OF OUTPUT INTERVAL TO PROCEED FROM TINIT TO TFINAL.
C
DT=SIGN(TINCR/TFINAL-TINIT)
C
C
C
INITIALIZE COUNTER OF ERROR RETURNS.
IER=0
C**************************************************************
C
C CALL STIFF TO INTEGRATE TO THE NEXT OUTPUT POINT.
C
5 CALL STIFF(FCT/FDER,Y,W(1)/IW/.W(K2),W
-------
C
c
c
c
c
c
c
c
c
CALL OUTP(T/Y/W)
CHECK WHETHER TFINAL HAS BEEN REACHED OR PASSED.
IF (ABS(T-TINIT).GE.TMAXIN) GO TO 95
INTEGRATION IS NOT YET COMPLETE/ SO THE NEXT OUTPUT POINT IS
SET. IF THE NEW POINT DIFFERS FROM TFINAL BY AN AMOUNT
ATTRIBUTABLE TO ROUNDOFF ERROR/ IT IS RESET TO TFINAL EXACTLY,
INTEGRATION IS THEN CONTINUED.
TOUT=T+DT
IF (ABS(TFINAL-TOUT).LT.EPSDT) TOUT=TFINAL
GO TO 5
c****************************************************************c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
INTEGRATION WAS INTERRUPTED BEFORE REACHING THE OUTPUT POINT.
INCREMENT THE ERROR COUNTER AND BRANCH FOR WRITING THE
APPROPRIATE DIAGNOSTIC MESSAGE. INTEGRATION WILL THEN BE
CONTINUED IF THE ERROR IS NOT TOO SEVERE.
8 IER=IER+1
GO TO (20/20/30/40/50/60/70/80/90)/IFLAG
IFLAG HAS AN ILLEGAL VALUE/ SO THAT A PRECISE DEFINITION OF
THE ERROR IS NOT POSSIBLE/ AND EXECUTION MUST BE TERMINATED.
20 IF (IUNIT.GT.O) WRITEUUNIT/200) T/IFLAG
STOP
IFLAG = 3 — THE REQUESTED ERROR TOLERANCES WERE TOO
SMALL FOR THE MACHINE PRECISION AND HAVE BEEN INCREASED TO
MINIMUM PERMISSIBLE VALUES.
30 IF (IUNIT.GT.O) WRITEUUNIT/300) T/RELER/ABSER
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 4 — THE REQUESTED ERROR TOLERANCES WERE TOO
STRINGENT FOR THIS PROBLEM AND HAVE BEEN INCREASED TO
VALUES WHICH SHOULD BE ACCEPTABLE.
40 IF (IUNIT.GT.O) WRITE(IUNIT/400) T/RELER/ABSER
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG f 5 — THE COUNTER OF CALLS TO FCT WAS RESET AFTER
EXCEEDING THE LIMIT IMPOSED IN STIFF.
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
24100
24200
24300
24400
24500
24600
24700
24800
24900
25000
25100
25200
25300
25400
25500
25600
25700
25800
25900
26000
26100
26200
26300
26400
26500
26600*
26700
26800
26900
27000
27100
27200
27300
27400
27500
27600
27700
27800
27900
28000
28100
28200
28300
28400
28500
28600
28700
28800
28900
29000
29100
29200
115
-------
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
50 IF UUNIT.GT.O) WRITE(IUNIT/500) T
IF (IER.6T.4) GO TO 92
GO TO 5
IFLAG = 6 — THE CONVERGENCE TEST IN GEAR COULD NOT BE MET
USING THE SMALLEST ALLOWABLE STEPSIZE. THE ERROR TOLERANCES
HAVE BEEN INCREASED TO VALUES WHICH WILL ALLOW THE
CONVERGENCE TEST TO BE PASSED.
60 IF UUNIT.GT.O) WRITEUUNIT/600) T/RELER/ABSER
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 7 — THE MONITOR OF OUTPUT FREQUENCY WAS RESET
AFTER INDICATING THAT THE OUTPUT INTERVAL IS TOO SMALL
TO ALLOW EFFICIENT OPERATION OF STIFF.
70 IF (IUNIT.GT.O) WRITE(IUNIT/700) T
IF (IER.GT.4) GO TO 92
GO TO 5
IFLAG = 8 — FATAL ERROR — A COMPONENT OF THE SOLUTION
VANISHED/ MAKING A PURE RELATIVE ERROR TEST IMPOSSIBLE.
THE USER MUST SUPPLY AN APPROPRIATE NONZERO VALUE OF
ABSERR FOR SUCCESSFUL INTEGRATION.
80 IF (IUNIT.GT.O) WRITE(IUNIT,800) T
STOP
IFLAG = 9 — FATAL ERROR — AN ILLEGAL PARAMETER VALUE
WAS PASSED TO STIFF. THE USER MUST LOCATE AND CORRECT THE
ERROR.
90 IF (IUNIT.GT.O) WRITEUUNIT/900) T/NEQ/RELER/ABSER/TINCR
STOP
THE PROGRAM IS STOPPED BECAUSE THERE HAVE BEEN FIVE ERROR
RETURNS FROM STIFF/ INDICATING SERIOUS PROBLEMS WHICH
REQUIRE INTERVENTION BY THE USER.
92 IF (IUNIT.GT.O) WRITE(IUNIT,920) T
STOP
THE INTEGRATION HAS BEEN COMPLETED TO OR BEYOND TFINAL.
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
29300
29400
29500
29600
29700
29800
29900
30000
30100
30200
30300
30400
30500
30600
30700
30800
30900
31000
31100
31200
31300
31400
31500
31600
31700
31800
31900
32000
32100
32200
32300
32400
32500
32600
32700
32800
32900
33000
33100
33200
33300
33400
33500
33600
33700
33800
33900
34000
34100
34200
34300
34400
116
-------
95 IF (IUNIT.6T.O) WRITE(IUNIT,950) T 34500
RETURN 34600
c C 34700
C*******************************************A 35800
C C 35900
200 FORMAT(4H ***,12(4H****)/ 36000
1 37H *** THE PROGRAM IS TERMINATED AT T =,G12.5,3H***/ 36100
2 52H *** BECAUSE AN UNIDENTIFIABLE ERROR HAS CAUSED ***/ 36200
3 52H *** STIFF TO RETURN WITH THE ILLEGAL FLAG VALUE ***/ 36300
4 12H *** IFLAG =,I11,26X,3H***/4H ***,12<4H****)) 36400
C C 36500
300 FORMAT(17H *WARNING—AT T =,G12.5,22H, THE ERROR TOLERANCES/ 36600
1 51H PASSED TO STIFF WERE FOUND TO BE TOO SMALL FOR THE/ 36700
2 48H MACHINE PRECISION. THEY HAVE BEEN INCREASED TO/ 36800
3 9H RELERR =,G12.5,13H AND ABSERR =,G12.5/ 36900
4 25H TO CONTINUE INTEGRATION./) 37000
C C 37100
400 FORMATC17H *WARNING—AT T =,G12.5,21H, STIFF WAS UNABLE TO/ 37200
1 52H MEET THE REQUESTED ERROR TOLERANCES AT THE SMALLEST/ 37300
2 52H ALLOWABLE STEPSIZE. INTEGRATION IS CONTINUING WITH/ 37400
3 37H THE TOLERANCES INCREASED TO RELERR =,G12.5/ 37500
4 13H AND ABSERR =,G12.5/> 37600
C C 37700
500 FORMATC17H *WARNING—AT T =,G12.5,24tt, THE NUMBER OF FUNCTION/ 37800
1 51H CALLS EXCEEDED THE MAXIMUM NUMBER, MAXNFE, ALLOWED/ 37900
2 52H BY STIFF. THE FUNCTION EVALUATION COUNTER HAS BEEN/ 38000
3 51H REDUCED BY MAXNFE TO ALLOW INTEGRATION TO PROCEED./) 38100
C C 38200
600 FORMAT(17H *WARNING—AT T =,G12.5,22H, THE CONVERGENCE TEST/ 38300
1 41H IN GEAR COULD NOT BE MET AT THE SMALLEST/ 38400
2 47H ALLOWABLE STEPSIZE. INTEGRATION IS CONTINUING/ 38500
3 48H WITH THE ERROR TOLERANCES INCREASED TO RELERR =/ 38600
4 1H ,G12.5,13H AND ABSERR =,G12.5/) 38700
C C 38800
700 FORMAT(17H *WARNING—AT T =,G12.5,24H, STIFF RETURNED WITH AN/ 38900
1 50H INDICATION THAT THIS INTEGRATION IS UNNECESSARILY/ 39000
2 48H COSTLY BECAUSE THE REQUESTED OUTPUT INTERVAL IS/ 39100
3 49H SUBSTANTIALLY SMALLER THAN THE NATURAL STEPSIZE./ 39200
4 48H (REFER TO PROGRAM DOCUMENTATION FOR ALTERNATIVE/ 39300
5 52H APPROACHES.) THE OUTPUT FREQUENCY MONITOR HAS BEEN/ 39400
6 39H RESET TO ALLOW INTEGRATION TO PROCEED./) 39500
C C 39600
117
-------
800 FORMAT(4H ***/12(4H****)/
1 37H *** THE PROGRAM IS TERMINATED AT T =,G12.5/3H***/
2 40H *** BECAUSE A COMPONENT OF THE SOLUTION/9X/3H***/
3 52H *** VANISHED WHILE ABSERR = 0. IT WILL BE ***/
4 52H *** NECESSARY TO PROVIDE A NONZERO ABSOLUTE ***/
5 52H *** ERROR TOLERANCE TO PERFORM THIS INTEGRATION.***/
6 4H ***,12<4H****)>
900 FORMAT(4H ***,12(4H****)/
1 37H *** THE,PROGRAM IS TERMINATED AT T =/G12.5/3H***/
52H *** BECAUSE INVALID INPUT WAS RECEIVED BY STIFF.***/
37H *** CHECK THE FOLLOWING PARAMETERS—/12X/3H***/
12H *** NEON =>I11,26X,3H***/
14H *** RELERR =,G12.5/23X,3H***/
14H *** ABSERR =/.G12.5/'23X/3H***/
13H *** TINCR =/G12.5/24X/3H***/
52H *** ALSO, CHECK THAT THE OUTPUT ROUTINE DOES ***/
2
3
4
5
6
7
8
9
A
41H *** NOT RESET IFLAG TO AN ILLEGAL VALUE./8X/3H***/
4H ***,12(4H****»
920 FORMAT(4H ***/12<4H****)/
1 37H *** THE PROGRAM IS TERMINATED AT T =/G12.5/3H***/
2 52H *** BECAUSE STIFF HAS RETURNED IRREGULARLY FIVE ***/
3 52H *** TIMES/ WHICH INDICATES SOME SORT OF EXTREME ***/
4 52H *** DIFFICULTY IN THE INTEGRATION. THE USER ***/
5 52H *** SHOULD CAREFULLY CHECK THE ACCURACY OF HIS ***/
6 52H *** PROBLEM STATEMENT AND CODING. ***/
7 4H ***,12(4H****»
C C
950 FORMAT(4H ***/12(4H****)/
1 32H *** INTEGRATION COMPLETE AT T =/G12.5/5X/3H***/
2 4H ***/12(4H****>///)
C C
c****************************************************************c
C C
C END OF STFINT C
C C
END
39700
39800
39900
40000
40100
40200
40300
40400
40500
40600
40700
40800
40900
41000
41100
41200
41300
41400
41500
41600
41700
41800
41900
42000
42100
42200
42300
42400
42500
42600
42700
42800
42900
43000
43100
43200
43300
43400
SUBROUTINE STIFF(F/FDER/Y/YP/IPS/YY/W/WORK)
C C
c****************************************************************c
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
BEYOND TOUT INTERNALLY/ THOUGH NEVER BEYOND T + 10.*(TOUT-T)/ C
AND CALLS SUBROUTINE GINTRP TO INTERPOLATE THE SOLUTION AT C
SUBROUTINE STIFF USES SUBROUTINE GEAR TO INTEGRATE A SYSTEM
OF NEQN FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS
DY/DT = F(T/Y) / Y = (Y(1)/Y(2)/.../Y(NEQN»
FROM T TO TOUT. FOR REASONS OF EFFICIENCY/ STIFF INTEGRATES
10000
10100
10200
10300
10400
10500
10600
10700
10800
10900
11000
11100
118
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
TOUT. AN OPTION IS PROVIDED TO STOP THE INTEGRATION AT TOUT/
BUT IT SHOULD BE USED ONLY IF THE CHARACTERISTICS OF THE
PROBLEM MAKE IT IMPOSSIBLE TO CONTINUE THE INTEGRATION BEYOND
TOUT.
STIFF IS INTENDED FOR SOLVING STIFF PROBLEMS. IT IS NOT
AN EFFICIENT MEANS FOR SOLVING NONSTIFF PROBLEMS.
THE PARAMETERS IN THE ARGUMENT LIST ARE
F — THE NAME OF THE USER-SUPPLIED SUBROUTINE FCT/Y/YP)/
WHICH EVALUATES THE DERIVATIVES YPU)=DY(I)/DT/
1=1/2, .../NEQN.
FDER — THE NAME OF THE USER SUPPLIED SUBROUTINE WHICH
EVALUATES THE JACOBIAN MATRIX AS DESCRIBED IN THE
COMMENTS TO GEAR.
Y — AN ARRAY OF DIMENSION NEQN/ WHICH CONTAINS THE VALUE
OF THE SOLUTION.
YP — AN ARRAY OF DIMENSION NEQN/ WHICH CONTAINS THE VALUE
OF THE DERIVATIVE DY/DT.
IPS/YY/W — AN INTEGER ARRAY OF DIMENSION NEQN/ AND
ARRAYS OF DIMENSIONS (NEQN/7) AND (NEQN/NEQN)/
RESPECTIVELY. THEIR VALUES ARE NEEDED FOR
CONTINUING INTEGRATION ON SUBSEQUENT CALLS.
WORK — AN ARRAY OF DIMENSIONS (NEQN/9) WHICH IS USED
FOR WORKING STORAGE BY GEAR.
THE PARAMETERS IN THE COMMON BLOCK SFDATA ARE
T -- THE INDEPENDENT VARIABLE
TOUT — THE VALUE OF T AT WHICH OUTPUT IS DESIRED
RELERR/ABSERR — RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR
THE LOCAL ERROR TEST. EACH STEPSIZE IS CHOSEN SO THAT
ABSCLOCAL ERROR(D) .LE. RELERR*ABS(Y(I))+ABSERR FOR
EACH COMPONENT OF THE SOLUTION.
NEQN — THE NUMBER OF EQUATIONS
IFLAG — INDICATOR OF THE STATUS OF INTEGRATION
DTSGN/IOLD/INIT/KOP — THOSE QUANTITIES USED FOR
INTERNAL BOOKKEEPING BY STIFF WHICH ARE NEEDED FOR
CONTINUING INTEGRATION/ BUT WHICH ARE NORMALLY OF NO
INTEREST TO THE USER. THEY ARE PLACED IN THE COMMON
BLOCK TO AVOID LOCAL RETENTION BY STIFF BETWEEN CALLS.
SEE THE PROGRAM DOCUMENTATION FOR A DESCRIPTION OF THEIR
FUNCTIONING.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c'
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
11200
11300
11400
11500
11600
11700
11800
11900
12000
12100
12200
12300
12400
12500
12600
12700
12800
12900
13000
13100
13200
13300
13400
13500
13600
13700
13800
13900
14000
14100
14200
14300
14400
14500
14600
14700
14800
14900*
15000
15100
15200
15300
15400
15500
15600
15700
15800
15900
16000
16100
1 A?nn
I w£UU
16300
119
-------
C TO USE STIFF/ THE CALLING PROGRAM MUST
C
C 1. ALLOCATE SPACE FOR THE ARRAYS Y(NEQN)/ YP(NEQN)/
C IPS(NEQN)/ YYCNEQN/7)/ W(NEQN/NEQN)/ AND WORKCNEQN/9) .
C (Y AND YP MUST BE ARRAYS EVEN IF NEQN=1.)
C
C 2. DECLARE THE NAMES OF F AND FDER IN AN EXTERNAL
C STATEMENT/ AND PLACE THOSE NAMES IN THE CALL LIST.
C
C
C FIRST CALL TO STIFF
£»______________...._
»««_*__»___«•_•_*••»
C
C TO BEGIN INTEGRATION/ THE CALLING PROGRAM MUST
C
C 1. INITIALIZE THE FOLLOWING VARIABLES IN THE COMMON BLOCK
C SFDATA
C T — THE STARTING VALUE OF THE INDEPENDENT VARIABLE
C TOUT — THE VALUE OF THE INDEPENDENT VARIABLE AT WHICH
C OUTPUT IS DESIRED. (TOUT=T IS PERMITTED ON THE
C FIRST CALL ONLY/ IN WHICH CASE STIFF EVALUATES THE
C DERIVATIVES YP AND RETURNS WITH IFLAG=2.)
C RELERR — THE DESIRED RELATIVE ERROR TOLERANCE
C ABSERR — THE DESIRED ABSOLUTE ERROR TOLERANCE
C IFLAG — NORMALLY +1. IF THE INTEGRATION CANNOT
C PROCEED BEYOND TOUT/ SET IFLAG = -1.
C
C 2. INITIALIZE Y(*) TO CONTAIN THE STARTING VALUE OF THE
C SOLUTION.
C
C
C UPON SUCCESSFUL RETURN UFLAG=2)/ T=TOUT/ Y CONTAINS THE
C SOLUTION AT TOUT/ AND YP CONTAINS THE DERIVATIVES AT TOUT.
C
C OTHERWISE/ T IS THE LAST VALUE OF THE INDEPENDENT VARIABLE
C SUCCESSFULLY REACHED DURING THE INTEGRATION/ Y CONTAINS THE
C SOLUTION AT T/ YP CONTAINS THE DERIVATIVES AT T/ AND IFLAG
C INDICATES THE STATUS OF THE INTEGRATION AS DESCRIBED BELOW.
C
C
C SUBSEQUENT CALLS TO STIFF
C__ ______»
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C STIFF NORMALLY RETURNS WITH ALL INFORMATION NEEDED TO CONTINUEC
C INTEGRATION.
C
C IF THE INTEGRATION REACHED TOUT (IFLAG=2)/ THE USER NEED
C ONLY DEFINE A NEW TOUT AND CALL STIFF AGAIN. IF THE INTE-
C GRATION CANNOT CONTINUE INTERNALLY BEYOND THE NEW TOUT/
C THE USER MUST ALSO SET IFLAG=-2. IF THE NEW TOUT REPRESENTS
C A REVERSAL OF THE DIRECTION OF INTEGRATION/ IT IS NECESSARY
C THAT Y NOT HAVE BEEN ALTERED.
C
C
C
C
C
C
C
C
16400
16500
16600
16700
16800
16900
17000
17100
17200
17300
17400
17500
17600
17700
17800
17900
18000
18100
18200
18300
18400
18500
18600
18700
18800
18900
19000
19100
19200
19300
19400
19500
19600
19700
19800
19900
20000
20100
20200
20300
20400
20500
20600
20700
20800
20900
21000
21100
21200
21300
21400
21500
120
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
IF THE INTEGRATION WAS INTERRUPTED WITH ABS(IFLAG) = 3/ 4/ 5/
6/ OR 7/ INTEGRATION MAY BE CONTINUED SIMPLY BY CALLING STIFF
AGAIN. HOWEVER/ SUCH A RETURN INDICATES THAT THE
INTEGRATION IS NOT PROCEEDING SO SMOOTHLY AS MIGHT BE
EXPECTED/ SO THE USER SHOULD EXERCISE SOME CAUTION.
IF THE INTEGRATION WAS INTERRUPTED WITH IFLAG = 8/ THE USER
MUST GIVE ABSERR A NONZERO POSITIVE VALUE AND SET IFLAG = 2
(OR IFLAG=-2 IF IT IS NECESSARY TO STOP AT TOUT) IN ORDER TO
CONTINUE.
IF STIFF RETURNED WITH IFLAG = 9/ THE USER MUST CORRECT THE
INVALID INPUT AND SET IFLAG = 2 OR -2 (OR IFLAG =1 OR -1,
IF INTEGRATION HAS NOT YET STARTED) IN ORDER TO CONTINUE.
UPON RETURN/ THE VALUE OF IFLAG INDICATES THE STATUS OF THE
INTEGRATION AS FOLLOWS/
IFLAG = 2 — T=TOUT/ SO THAT INTEGRATION WAS SUCCESSFUL TO
THE OUTPUT POINT.
IFLAG = 3 — INTEGRATION WAS INTERRUPTED BECAUSE THE
REQUESTED ERROR TOLERANCES WERE TOO SMALL FOR
THE MACHINE PRECISION. THEY HAVE BEEN INCREASED
TO MINIMUM ACCEPTABLE VALUES.
IFLAG = 4 — INTEGRATION WAS INTERRUPTED BECAUSE THE
REQUESTED ERROR TOLERANCES COULD NOT BE
ATTAINED AT THE SMALLEST ALLOWABLE STEPSIZE.
THE TOLERANCES HAVE BEEN INCREASED TO VALUES
WHICH SHOULD ALLOW SUCCESSFUL CONTINUATION.
IFLAG = 5 — INTEGRATION WAS INTERRUPTED AFTER THE COUNTER
OF CALLS TO F EXCEEDED MAXNFE. THE COUNTER HAS
BEEN REDUCED BY MAXNFE.
IFLAG = 6 — INTEGRATION WAS INTERRUPTED BECAUSE THE
ERROR TOLERANCES WERE TOO SMALL TO PERMIT
THE ITERATION IN GEAR TO CONVERGE FOR THE
SMALLEST ALLOWABLE STEPSIZE. THE TOLERANCES
HAVE BEEN INCREASED TO VALUES WHICH WILL
ALLOW CONVERGENCE.
IFLAG = 7 — INTEGRATION WAS NOT CONTINUED BECAUSE THE
SPACING OF OUTPUT POINTS HAS BEEN SO SMALL AS TO
RESTRICT THE STEPSIZE FOR 50 CALLS TO GEAR/
INDICATING THAT THE CODE IS BEING USED IN A WAY
WHICH GREATLY IMPAIRS ITS EFFICIENCY. THE
MONITOR OF THIS EFFECT WAS RESET TO ALLOW
CONTINUATION IF DESIRED.
IFLAG = 8 — INTEGRATION WAS INTERRUPTED WHEN A COMPONENT OF
THE SOLUTION VANISHED WHILE A PURE RELATIVE
ERROR TEST WAS SPECIFIED. THE USER MUST SUPPLY
A NONZERO ABSERR TO CONTINUE INTEGRATION.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
21600
21700
21800
21900
22000
22100
22200
22300
22400
22500
22600
22700
22800
22900
23000
23100
oionn
OcUU
23300
23400
23500
23600
23700
23800
23900
24000
24100
24200
24300
24400
24500
24600
24700
24800
24900
25000
25100
25200
25300
25400
25500
25600
25700
25800
25900
26000
26100
26200
26300
26400
26500
26600
26700
121
-------
c
c
c
c
c
c
c
IFLAG = 9 — INVALID INPUT WAS PASSED TO STIFF. THE USER
MUST LOCATE AND CORRECT THE ERRO.R IN ORDER TO
CONTINUE INTEGRATION.
IFLAG =-3/ -4/ -5/ -6/ OR -7 — SAME AS IFLAG = 3/ 4/ 5/ 6,
OR 1, RESPECTIVELY/ EXCEPT THAT STIFF WAS CALLED
WITH IFLAG NEGATIVE.
C
C
C
C
C
C
C
c***********************************************************
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
THIS SUBROUTINE IS WRITTEN IN SINGLE PRECISION.
WITH MOST COMPILERS/ A DOUBLE PRECISION VERSION/ FOR USE WITH
DOUBLE PRECISION VERSIONS OF SUBROUTINES GEAR AND GINTRP/ IS
OBTAINED SIMPLY BY SETTING THE APPROPRIATE VALUE FOR FOURU
IN THE DATA STATEMENT BELOW AND BY REMOVING THE C FROM
COLUMN 1 OF THE NEXT SEVEN CARDS.
DOUBLE PRECISION ABS/ABSDT/ABSEPS/ABSERR/AMAX1/DABS/DMAX1/
1 DSIGN/DT/DTSGN/D1 /D2/EPS/FOURU/H/HOLD/
2 RELEPS/RELERR/RVAR/SIGN/T/TEND/TOUT/
3 W/WORK/X/Y/YP/YY
ABSCD1) =DABS(D1)
SIGN(D1/D2) = DSIGN(D1/D2)
AMAXKD1/D2) = DMAX1CD1/D2)
LOGICAL START
COMMON /SFDATA/ T/TOUT/RELERR/ABSERR/DTSGN/NEQN/IFLAG/
1 IOLD/INIT/KOP
COMMON /GEDATA/ X/HOLD/H/EPS/RVAR(8)/NEQ/K/NFE/
1 ICRASH/KOLD/ITST/IWEVAL/START
DIMENSION Y(NEQN)/YP(NEQN)/IPS(NEQN)/YY(NEQN/7)/
1 W(NEQN/NEQN)/WORK
EXTERNAL F/FDER
IN THE FOLLOWING DATA STATEMENT/ FOURU MUST BE SET TO A
VALUE APPROPRIATE TO THE MACHINE BEING USED.
FOURU IS 4.*U/ WHERE THE COMPUTER UNIT ROUNDOFF ERROR U IS
THE SMALLEST POSITIVE VALUE REPRESENTABLE IN THE MACHINE
SUCH THAT (1.+U).GT.1.
SOME REPRESENTATIVE VALUES ARE
U = 9.54E-7 FOURU = 3.82E-6 IBM 360/370 S.P.
U = 5.96E-8 FOURU = 2.39E-7 PDP-11 S.P.
U = 1.49E-8 FOURU = 5.97E-8 UNIVAC 1108
U = 7.45E-9 FOURU = 2.99E-8 PDP-10
U = 7.11E-15 FOURU = 2.85E-14 CDC 6000 SERIES
U = 2.22D-16 FOURU = 8.89D-16 IBM 360/370 D.P.
U = 1.39D-17 FOURU = 5.56D-17 PDP-11 D.P.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
26800
26900
27000
27100
27200
27300
27400
27500
27600
27700
27800
27900
28000
28100
28200
28300
28400
28500
28600
28700
28800
28900
29000
29100
29200
29300
29400
29500
29600
29700
29800
29900
30000
30100
30200
30300
30400
30500
30600
30700
30800
30900
31000
31100
31200
31300
31400
31500
31600
31700
31800
31900
122
-------
C
c
C THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER
C OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE.
C AS SET, THIS CORRESPONDS TO ABOUT 1500 STEPS.
DATA FOURU/2.39E-7/
THE
OF
AS
C
DATA MAXNFE/4000/
C
c****************************************************************c
c
C TEST FOR IMPROPER PARAMETERS
C
MFLAG=IABS(IFLAG)
IF UMFLAG.GT.7).OR.(MFLAG.EQ.O).OR.(NEQN.LT.1)) GO TO 200
IF ((RELERR.LT.O.O).OR.(ABSERR.LT.O.O)) GO TO 200
EPS=AMAX1(RELERR/ABSERR)
IF (EPS.LE.0.0) GO TO 200
C
NEQ=NEQN
DT=TOUT-T
ABSDT=ABS(DT)
TEND=T + 10.*DT
IFUFLAG.LT.O) TEND=TOUT
RELEPS=RELERR/EPS
ABSEPS=ABSERR/EPS
C
C IF THIS IS THE FIRST CALL/ GO TO INITIALIZATION SECTION.
C
IF (MFLAG.EQ.1) GO TO 150
C
C SINCE THIS IS NOT THE FIRST CALL/ T=TOUT IS NOT VALID
C INPUT.
C
IF (T.EQ.TOUT) GO TO 200
C
C IF INITIALIZATION WAS NOT COMPLETED ON THE FIRST CALL/ DO
C IT NOW.
C
IF (INIT.EQ.O) GO TO 180
C
C IF THE LAST STEP WAS MADE WITH IFLAG NEGATIVE/ OR IF THE
C DIRECTION OF INTEGRATION IS CHANGED/ RESTART THE INTEGRATION. C
C
IF (IOLD .LT. 0) GO TO 160
IF (DTSGN*DT .LT. 0.0) GO TO 160
C
C IF ALREADY PAST OUTPUT POINT/ INTERPOLATE AND RETURN.
C
50 IFCABS(X-T) .LT. ABSDT)GO TO 60
CALL GINTRPCTOUT/Y/YP/YY)
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
32000
32100
32200
32300
32400
32500
32600
32700
32800
32900
33000
33100
33200
33300
33400
33500
33600
33700
33800
33900
34000
34100
34200
34300
34400
34500
34600
34700
34800
34900
35000
35100
35200
35300
35400
35500
35600
35700
35800
35900
36000
36100
36200
36300
36400
36500
36600
36700
36800
36900
37000
37100
123
-------
IOLD=IFLAG 37200
IFLAG=2 37300
T=TOUT 37400
RETURN 37500
C C 37600
C IF CANNOT GO PAST OUTPUT POINT/ TEST WHETHER TOO CLOSE TO C 37700
C INTEGRATE WITH GEAR. C 37800
C C 37900
60 IF «IFLAG.LT.O).AND.(ABS(TOUT-X).LT.FOURU*ABS(X))) GO TO 210 38000
C C 38100
C TEST FOR TOO MUCH WORK. C 38200
C C 38300
IF(NFE .GT. MAXNFE) GO TO 230 38400
C C 38500
C LIMIT STEPSIZE/ SET WEIGHT VECTOR/ AND TAKE A STEP. C 38600
C C 38700
IF (ABS(H).LT.ABS(TEND-X)) GO TO 100 38800
KOP=KOP+1 38900
IF (KOP.GT.50) GO TO 290 39000
H=TEND-X 39100
C C 39200
C THE WEIGHT VECTOR IS PLACED IN Y<*). C 39300
C C 39400
100 DO 110 L=1/NEQN ' 39500
Y(L)=RELEPS*ABS(YY(L/1))+ABSEPS 39600
IF (Y(L).EQ.O.O) GO TO 270 39700
110 CONTINUE 39800
C C 39900
CALL GEAR(F/FDER/YY/W/IPS/Y/YP/WORK(1/1)/WORK(1/2)/ 40000
1 WORK(1/3)) 40100
C C 40200
C TEST WHETHER STEP WAS SUCCESSFUL. C 40300
C C 40400
IF (ICRASH.NE.O) GO TO 250 40500
C C 40600
C STEP SUCCEEDED. CONTINUE WITH INTEGRATION. C 40700
C C 40800
GO TO 50 40900
C C 41000
C****************************************************************c 41100
C C 41200
C INITIALIZATION SECTION — C 41300
C C 41400
C C 41500
C ON START AND RESTART/ SET WORK VARIABLES X AND YY/ STORE C 41600
C THE DIRECTION OF INTEGRATION/ AND INITIALIZE THE STEPSIZE. C 41700
C C 41800
150 NFE=0 41900
160 START=.TRUE. 42000
INIT=0 42100
KOP=0 42200
X=T 42300
124
-------
DO 170 L=1,NEQN
170 YY(L,1)=Y(L)
C
C IF T=TOUT ON FIRST CALL/ EVALUATE DERIVATIVES AND RETURN.
C
IF (T.NE.TOUT) GO TO 180
CALL F(X,YY,YP)
NFE=NFE-f1
IFLAG=2
RETURN
C
C SET STEPSIZE AND TRANSFER TO BEGIN INTEGRATION.
C
180 INIT=1
DTSGN=DT
H=SIGN(AMAX1(ABS(TOUT-X),FOURU*ABS(X))/TOUT-X)
GO TO 50
C
C
C
C
C
C
C
C
c****************************************************************c
C
C THE REMAINDER OF THE PROGRAM CONTAINS STATEMENTS TO HANDLE
C VARIOUS EXCEPTIONAL OCCURRENCES.
C
C
C INVALID INPUT WAS PASSED TO STIFF.
C
200 IFLAG=9
RETURN
C
C CANNOT PASS TOUT* BUT TOO CLOSE TO CALL GEAR. EXTRAPOLATE
C BY EULER METHOD AND RETURN.
C
210 H=TOUT-X
CALL F
NFE=NFE+1
IOLD=-1
RETURN
C
C NUMBER OF FUNCTION CALLS EXCEEDED MAXNFE. RESET COUNTER NFE.
C
230 IFLAG=ISIGN(5*IFLAG>
NFE=NFE-MAXNFE
GO TO 300
C
C ERROR TOLERANCES WERE TOO SMALL TO PERMIT A SUCCESSFUL STEP.
C INCREASE RELERR AND ABSERR, AND SET IFLAG TO INDICATE THE
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
42400
42500
42600
42700
42800
42900
43000
43100
43200
43300
43400
43500
43600
43700
43800
43900
44000
44100
44200
44300
44400
44500
44600
44700
44800
44900
45000
45100
45200
45300
45400
45500
45600
45700
45800
45900
46000
46100
46200
46300
46400
46500
46600
46700
46800
46900
47000
47100
47200
47300
47400
47500
125
-------
c
c
CAUSE OF THE DIFFICULTY.
250 IFLAG=ISIGN(4/IFLAG)
IF (ICRASH.EQ.2) IFLAG=ISIGN(3/IFLAG)
IF (ICRASH.EQ.3) IFLAG=ISIGN(6/IFLAG)
RELERR=EPS*RELEPS
ABSERR=EPS*ABSEPS
GO TO 300
A COMPONENT OF THE SOLUTION VANISHED WHEN A PURE RELATIVE
ERROR TEST WAS SPECIFIED. USER MUST SET ABSERR POSITIVE
IN ORDER TO PROCEED.
270 IFLAG=8
GO TO 300
THE OUTPUT INTERVAL HAS BEEN TOO SMALL TO ALLOW EFFICIENT
OPERATION OF STIFF.
290 KOP=0
IFLAG=ISIGN(7/IFLAG)
FOR ERROR RETURN/ SET T/ Y(*)/ AND YP(*) TO THE MOST RECENT
VALUES REACHED IN THE INTEGRATION.
300 T=X
IOLD=1
DO 310 L=1/NEQN
Y(L)=YY(L/1)
310 YP(L)=YY(L/2)/HOLD
RETURN
C
C*************************************************************
C
C END OF STIFF
C
END
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
47600
47700
47800
47900
48000
48100
48200
48300
48400
48500
48600
48700
48800
48900
49000
49100
49200
49300
49400
49500
49600
49700
49800
49900
50000
50100
50200
50300
50400
50500
50600
50700
50800
50900
51000
51100
51200
SUBROUTINE GEAR(FCT/FDER/Y,W,IPS,WT/YP/C,ERROR/SAVE)
c
c
c
c
c
c
c
c
c
c
SUBROUTINE GEAR USES A BACKWARD DIFFERENTIATION METHOD (OF
ORDER .LE. 6) TO PERFORM ONE STEP IN THE INTEGRATION OF A
SYSTEM OF NEQN FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS
DY/DT = F(T/Y) / Y = (Yd )/Y(2>/.. ./Y(NEQN))
TO SPECIFY THE DIFFERENTIAL EQUATIONS TO BE SOLVED, THE USER
MUST SUPPLY TWO SUBROUTINES. THE FIRST/ DENOTED HERE BY
FCT/ CALCULATES THE NEQN DERIVATIVES DY(I)/DT. THE SECOND/
C
C
C
C
C
C
C
C
C
C
10000
10100
10200
10300
10400
10500
10600
10700
10800
10900
11000
11100
11200
126
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
DENOTED HERE BY FDER/ CALCULATES THE NEQN**2 ELEMENTS OF
THE JACOBIAN MATRIX OF F(T/Y).
FCT MUST BE OF THE FORM FCT(T/Y/YP)/ WHERE Y AND YP ARE
ARRAYS EACH OF DIMENSION NEQN. FCT MUST CALCULATE THE
VALUES OF THE DERIVATIVES DY(I)/DT AT (T/Y) AND PLACE
THEM IN YP. T AND Y MUST NOT BE ALTERED BY FCT.
FDER MUST BE OF THE FORM FDER(T/Y,W/YP/DER)/ WHERE Y AND
YP ARE ARRAYS EACH OF DIMENSION NEQN/ AND W IS AN ARRAY
OF DIMENSIONS (NEQN/NEQN). DER IS A SUBROUTINE NAME.
WHEN FDER IS CALLED/ DER IS THE NAME WHICH THE USER HAS
SUPPLIED FOR FCT. FDER MUST CALCULATE THE VALUES OF ALL
THE PARTIAL DERIVATIVES OF F AT (T/Y). THE PARTIAL
DERIVATIVE OF THE M-TH COMPONENT OF F WITH RESPECT TO Y(N)
IS TO BE PLACED IN W(M/N). WHEN FDER IS CALLED/ YP
CONTAINS THE VALUES OF DY/DT AT (T/Y). FOR SOME PROBLEMS/
THESE VALUES MAY BE USEFUL IN SIMPLIFYING THE EVALUATION
OF THE JACOBIAN. T/ Y/ AND YP MUST NOT BE ALTERED BY FDER.
FOR FDER/ THE USER MAY SUPPLY THE DEFAULT SUBROUTINE
QJACOB/ WHICH USES FCT TO COMPUTE THE JACOBIAN BY
NUMERICAL DIFFERENCING. SINCE THIS PROCEDURE INVOLVES
NEQN CALLS TO FCT/ ITS USE IS NOT RECOMMENDED UNLESS THE
JACOBIAN IS DIFFICULT TO COMPUTE DIRECTLY.
SUBROUTINE GEAR IS INTENDED SPECIFICALLY FOR THE SOLUTION OF
STIFF PROBLEMS. FOR NONSTIFF PROBLEMS/ GEAR WILL BE LESS
EFFICIENT THAN A CODE BASED ON ADAMS METHODS.
INFORMATION IS EXCHANGED WITH GEAR BOTH THROUGH THE ARGUMENT
LIST AND THROUGH THE COMMON BLOCK GEDATA.
THE PARAMETERS IN THE ARGUMENT LIST ARE
FCT — THE NAME OF THE USER-SUPPLIED SUBROUTINE TO
CALCULATE THE DERIVATIVES
FDER — THE NAME OF THE USER-SUPPLIED SUBROUTINE TO
CALCULATE THE JACOBIAN MATRIX
Y — AN ARRAY OF DIMENSIONS (NEQN/7) WHOSE VALUE MUST NOT
BE ALTERED BETWEEN CALLS TO GEAR. THE FIRST NEQN
LOCATIONS/ Y(1/1)/Y(2/1)/.../Y(NEQN/1) CONTAIN THE
SOLUTION VALUES
W/IPS __ RESPECTIVELY/ AN ARRAY OF DIMENSIONS (NEQN/NEQN)
AND AN INTEGER ARRAY OF DIMENSION NEQN/ WHOSE VALUES
MUST NOT BE ALTERED BETWEEN CALLS TO GEAR
WT — AN ARRAY OF DIMENSION NEQN USED IN SPECIFYING THE
ERROR TEST.
YP/C/ERROR/SAVE — RESPECTIVELY/ ARRAYS OF DIMENSIONS NEQN/
NEQN/ NEQN/ AND (NEQN/7)/ WHICH ARE USED FOR WORKING
STORAGE BY GEAR. IT IS NOT NECESSARY TO PRESERVE
THEIR VALUES BETWEEN STEPS
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
11300
11400
11500
11600
11700
11800
11900
12000
12100
12200
12300
12400
12500
12600
12700
12800
12900
13000
13100
13200
13300
13400
13500
13600
13700
13800
13900
14000
14100
14200
14300
14400
14500
14600
14700
14800
14900
15000
15100
15200
15300
15400
15500
15600
15700
15800
15900
16000
16100
16200
16300
16400
127
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
THE PARAMETERS IN THE COMMON BLOCK 6EOATA ARE
T -- THE INDEPENDENT VARIABLE
HOLD — THE STEPSIZE USED FOR THE LAST SUCCESSFUL STEP
H — THE STEPSIZE TO BE USED IN ATTEMPTING THE NEXT STEP
EPS — THE DESIRED TOLERANCE FOR THE ERROR TEST
NEQN — THE NUMBER OF EQUATIONS BEING INTEGRATED
K — THE ORDER OF THE METHOD TO BE USED FOR THE NEXT STEP
ICRASH — AN INDICATOR OF ERROR CONDITIONS UPON RETURN
NFE — A COUNTER WHICH IS INCREMENTED BY 1 AFTER EACH
CALL TO FCT
START — A LOGICAL VARIABLE WHICH MUST BE SET .TRUE. ON
THE FIRST CALL TO GEAR
AHINV/A(7)/KOLD/ITST/IWEVAL — VARIABLES USED FOR INTERNAL
BOOKKEEPING BY GEAR WHOSE VALUES ARE NEEDED FOR
CONTINUING INTEGRATION. THEIR VALUES ARE NORMALLY OF
NO INTEREST TO THE USER/ AND THEY ARE PLACED IN THE
COMMON BLOCK ONLY TO AVOID LOCAL RETENTION BETWEEN
CALLS.
EXCEPT FOR H/ EPS/ AND NFE/ THE VALUES IN THE COMMON BLOCK
MUST NOT BE ALTERED BETWEEN CALLS TO GEAR.
THE COUNTER NFE IS PROVIDED SO THAT THE USER CAN MONITOR THE
AMOUNT OF WORK BEING DONE. ITS VALUE HAS NO MEANING TO GEAR.
GEAR AUTOMATICALLY ADJUSTS THE STEPSIZE AND THE ORDER OF THE
METHOD SO THAT THE MAX NORM OF THE VECTOR WITH COMPONENTS
LOCAL ERROR(I)/WT(I) IS LESS THAN EPS FOR EACH STEP. THE
ARRAY WT ALLOWS THE USER TO SPECIFY AN ERROR TEST WHICH IS
APPROPRIATE FOR HIS PROBLEM.
FOR EXAMPLE/
'
WTU) = 1.0 SPECIFIES ABSOLUTE ERROR
WTU) = ABSnn
IVcUU
19300
19400
19500
19600
19700
19800
19900
20000
20100
20200
20300
20400
20500
20600
20700
20800
20900
21000
21100
21200
21300
21400
21500
21600
128
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
C_
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
ALL COMPONENTS OF WT MUST BE NONZERO.
UPON RETURN, THE VALUE OF ICRASH INDICATES THE ACTION WHICH
WAS TAKEN.
ICRASH = 0 — A STEP OF LENGTH HOLD WAS SUCCESSFULLY
TAKEN TO THE VALUE NOW IN T. THE SOLUTION AT T
IS IN Y(*,1), AND Y(*/2) CONTAINS THE DERIVATIVES
MULTIPLIED BY HOLD.
ICRASH = 1 — NO STEP WAS TAKEN BECAUSE H WAS SMALLER THAN
PERMITTED BY THE MACHINE PRECISION. H WAS INCREASED
TO AN ACCEPTABLE VALUE.
ICRASH = 2 — NO STEP WAS TAKEN BECAUSE EPS WAS TOO SMALL
FOR THE MACHINE PRECISION. EPS WAS INCREASED TO AN
ACCEPTABLE VALUE.
ICRASH = 3 — NO STEP WAS TAKEN BECAUSE CORRECTOR
CONVERGENCE COULD NOT BE ACHIEVED AT THE SMALLEST
ALLOWED STEPSIZE. EPS IS INCREASED TO A VALUE WHICH
ALLOWS THE CONVERGENCE TEST TO BE PASSED.
ICRASH = 4 — NO STEP WAS TAKEN BECAUSE THE ERROR TEST
COULD NOT BE SATISFIED AT THE SMALLEST ALLOWED STEPSIZE.
EPS IS INCREASED TO A VALUE WHICH IS ACHIEVABLE.
FIRST CALL TO GEAR
__________________
TO BEGIN THE SOLUTION/ THE CALLING PROGRAM MUST
1. PROVIDE STORAGE FOR THE ARRAYS YCNEQN/7), W(NEQN,NEQN),
IPS(NEQN), WT(NEQN), YP(NEQN), C(NEQN)/ ERROR(NEQN)/
AND SAVECNEQN/7)
2. PLACE THE INITIAL VALUES OF THE SOLUTION IN
Y(1,1>, Y<2,1)/ .../ Y(NEQN,1>.
3. SET WT(1)/ ...* WT(NEQN) TO NONZERO VALUES TO PROVIDE
THE ERROR TEST DESIRED.
4. INITIALIZE THE FOLLOWING VARIABLES IN THE COMMON BLOCK
GEDATA
- __ THE STARTING VALUE OF THE INDEPENDENT VARIABLE
H — A NOMINAL STEPSIZE FOR THE FIRST STEP/
INDICATING THE DIRECTION OF INTEGRATION AND THE
MAXIMUM STEPSIZE. H WILL BE REDUCED AS NECESSARY
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
21700
21800
O 4 OOrt
21900
22000
22100
22200
22300
22400
22500
22600
22700
22800
22900
23000
23100
23200
23300
23400
23500
23600
23700
23800
23900
24000
24100
24200
24300
24400
24500
pAAnn
_HOUU
24700
24800
24900
25000
25100
25200
25300
25400
25500
25600
25700
25800
25900
26000
26100
26200
26300
26400
26500
26600
26700
26800
129
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c*
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
TO MEET THE ERROR CRITERION.
EPS — THE DESIRED ERROR TOLERANCE
NEQN — THE NUMBER OF EQUATIONS
NFE — SET TO 0 TO BEGIN A COUNT OF FUNCTION CALLS
START — SET .TRUE.
4. PLACE THE ACTUAL NAMES OF THE USER SUPPLIED SUBROUTINES
FCT AND FDER IN THE CALL LIST/ AND DECLARE THOSE NAMES
IN AN EXTERNAL STATEMENT.
SUBSEQUENT CALLS TO GEAR
GEAR RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE THE
INTEGRATION. TO MAINTAIN RELATIVE ERROR TESTS AS DESCRIBED
ABOVE, WT(*) MUST BE UPDATED AFTER EACH STEP. CALLING GEAR
AGAIN THEN ADVANCES THE SOLUTION ANOTHER STEP. TO CONTINUE
INTEGRATION USING THE REVISED STEPSIZE OR TOLERANCE
PROVIDED WITH AN ERROR RETURN (ICRASH.NE.O), IT IS ONLY
NECESSARY TO CALL GEAR AGAIN.
NORMALLY, THE INTEGRATION IS CONTINUED TO THE FIRST STEP
BEYOND THE DESIRED OUTPUT POINT, AND THE SOLUTION IS INTER-
POLATED TO THE OUTPUT POINT USING THE SUBROUTINE GINTRP.
IF IT IS IMPOSSIBLE TO INTEGRATE BEYOND THE OUTPUT POINT, H
MAY BE REDUCED TO HIT THE OUTPUT POINT. OTHERWISE, H SHOULD
NOT BE ALTERED, SINCE CHANGING H IMPAIRS THE EFFICIENCY OF
THE CODE. H MUST NEVER BE INCREASED. EPS MAY BE CHANGED
BETWEEN CALLS, BUT LARGE DECREASES SHOULD BE AVOIDED, SINCE
THIS WILL INVALIDATE THE STEP SELECTION MECHANISM AND LEAD
TO STEP FAILURES.
************************Vf*************************************^
THIS SUBROUTINE IS WRITTEN IN SINGLE PRECISION.
WITH MOST COMPILERS, A DOUBLE PRECISION VERSION MAY BE
OBTAINED SIMPLY BY PROVIDING A DOUBLE PRECISION VERSION OF
DECOMP, BY SETTING THE APPROPRIATE VALUE FOR FOURU IN THE
DATA STATEMENT BELOW, AND BY REMOVING THE C FROM COLUMN 1 OF
THE NEXT TEN CARDS.
DOUBLE PRECISION A,ABS,ABSH,ADEN,AHINV,AMAX1,ANUM,BND,C,
1 CNORM,D,DABS,DMAX1 ,DSIGN,DSQRT,DUM1 ,DUM2,
2 D1,EPS,EPSMIN,ERDWN,ERROR,ERTST,ERUP,
3 FOURU,GFACT,H,HMIN,HOLD,R,RATIO,RDWN,RUP,
4 R1 ,SAVE,SIGN,SQRT,SUM,T,TOLD,W,WT, Y,
5 YNORM,YP,YPNORM
ABS(DUM1)=DABS(DUM1)
AMAX1 (DUM1,DUM2)=DMAX1 (DUM1,DUM2)
SIGN(DUM1,DUM2)=DSIGN(DUM1,DUM2)
SQRT(DUM1 )=DSQRT(DUM1 )
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
c
c
c
c
c
*c
c
c
c
c
c
c
c
c
26900
27000
27100
27200
27300
27400
27500
27600
27700
27800
27900
28000
oomn
co IUU
28200
28300
28400
28500
28600
28700
28800
28900
29000
29100
29200
29300
29400
29500
29600
29700
29800
29900
30000
30100
30200
30300
30400
30500
30600
30700
30800
30900
31000
31100
31200
31300
31400
31500
31600
31700
31800
31900
32000
130
-------
LOGICAL START
C
EXTERNAL FCT
C
COMMON /GEDATA/ T/HOLD/H/EPS/AHINV/A(7)/NEQN/K/NFE/
1 ICRASH/KOLD/ITST/IWEVAL/START
C
DIMENSION Y(NEQN/7)/SAVE(NEQN/7)/YP(NEQN)/WT(NEQN)/
1 ERROR(NEQN) /C (NEQN) / IPS (NEQN) /WCNEQN/NEQN) /
2 ERTST(6)/ERUP(5)/ERDWN(6)/ANUM(27)/ADEN(6)
C
DATA ERTST/2. 0/4. 5/7. 333/10. 42/13. 7/1 7. 15/
DATA ERUP/3. 0/6. 0/9. 167/12. 5/15. 98/
DATA ERDWN/1 .0/1 .0/0. 5/0. 1667/0. 04133/0. 008267/
C
DATA ANUM/1./1./2./3./1./6./11./6./1./24./50./35./10./
1 1./120./274./225./85./15./1./720./1764./
2 1624./735./175./21./1./
DATA ADEN/-1 ./-3./-11 ./-50./-274./-1764./
C
C
C IN THE FOLLOWING DATA STATEMENT/ FOURU MUST BE SET TO A
C VALUE APPROPRIATE TO THE MACHINE BEING USED.
C
C FOURU IS 4.*U/ WHERE THE COMPUTER UNIT ROUNDOFF ERROR U IS
C THE SMALLEST POSITIVE VALUE REPRESENTABLE IN THE MACHINE
C SUCH THAT (1.+U).GT.1.
C
C SOME REPRESENTATIVE VALUES ARE
C U = 9.54E-7 FOURU = 3.82E-6 IBM 360/370 S.P.
C U = 5.96E-8 FOURU = 2.39E-7 PDP-11 S.P.
C U = 1.49E-8 FOURU = 5.97E-8 UNIVAC 1108
C U = 7.45E-9 FOURU = 2.99E-8 PDP-10
C U = 7.11E-15 FOURU = 2.85E-14 CDC 6000 SERIES
C U = 2.22D-16 FOURU = 8.89D-16 IBM 360/370 D.P.
C U = 1.39D-17 FOURU = 5.56D-17 PDP-11 D.P.
C
DATA FOURU/2.39E-7/
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c****************************************************************c
c ******
C BEGIN BLOCK 0 ******
c ******
c***********************
c
C IF STEPSIZE IS TOO SMALL FOR MACHINE PRECISION/ INCREASE IT
C TO AN ACCEPTABLE VALUE AND RETURN TO CALLING PROGRAM.
C
HMIN=FOURU*ABS(T)
IF (ABS(H).GE.HMIN) GO TO 5
H=SIGN(HMIN/H)
ICRASH=1
C
C
c
c
c
c
c
c
32100
32200
32300
32400
32500
32600
32700
32800
32900
33000
33100
33200
33300
33400
33500
33600
33700
33800
33900
34000
34100
34200
34300
34400
34500
34.600
34700
34800
34900
35000
35100
35200
35300
35400
35500
35600
35700
35800
35900
36000
36100
36200
36300
36400
36500
36600
36700
36800
36900
37000
37100
37200
131
-------
C
c
c
c
c
c
c
c
c
c
c
c
c
c
C
RETURN
IF THE REQUESTED ERROR TOLERANCE IS TOO SMALL FOR THE MACHINE
PRECISION/ INCREASE IT TO AN ACCEPTABLE VALUE AND RETURN TO
THE CALLING PROGRAM.
5 YNORM=0.0
DO 10 I=1/NEQN
10 YNORM=AMAX1(YNORM,ABS(YU,1)/WTU)))
EPSMIN=FOURU*YNORM
IF (EPS.GE.EPSMIN) GO TO 15
EPS=EPSMIN*(1.+FOURU)
ICRASH=2
RETURN
IF THIS IS
INITIALIZE
THE FIRST
Y(*/2).
STEP/ CHOOSE A STARTING STEPSIZE AND
50
15 IF (.NOT.START) GO TO
CALL FCT(T/Y/YP)
NFE=NFE+1
YPNORM=0.0
DO 20 I=1/NEQN
20 YPNORM=AMAX1 (YPNORM/ABS(YP(I)/WTU)))
ABSH=ABS(H)
IF
-------
c***********************
c
C IF ORDER HAS BEEN CHANGED/ RECOMPUTE COEFFICIENTS A(*)
C
100 IF (K.EQ.KOLD) GO TO 120
KP1=K+1
INDEX=(K*KP1)/2 - 1
DO 110 J=1/KP1
J1=INDEX+J '
110 A(J)=ANUM(J1)/ADEN(K)
ITST=KP1
IWEVAL=1
KOLD=K
C
C
C
C
C
C
C
IF THE STEPSIZE HAS BEEN CHANGED/ REVISE Y TO CONFORM TO THE
NEW STEPSIZE.
120 IF (H.EQ.HOLD) GO TO 130
RATIO=H/HOLD
R1=1.0
DO 125 J=2/KP1
R1=R1*RATIO
DO 125 I=1/NEQN
125 Y(I/J)=SAVE(I/J)*R1
ITST=KP1
IWEVAL=1
SET LIMIT FOR TEST OF CONVERGENCE
130 BND=EPS/FLOAT(2*NEQN*(K+2))
C
C***********************
c ******
C END BLOCK 1 ******
C ******
c****************************************************************c
c ******
C BEGIN BLOCK 2 ******
C ******
C***********************
C
C SOLVE FOR Y AT T+H BY USING NEWTON ITERATION
C
C
C
C
C
C
200 T=T+H
OBTAIN STARTING VALUE FOR ITERATION BY MULTIPLYING Y BY
PASCAL TRIANGLE MATRIX/ AND INITIALIZE ERRORC*) AND THE
ITERATION COUNTER ITER.
DO 205 J=2/KP1
DO 205 J1=J/KP1
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
42500
42600
42700
42800
42900
43000
43100
43200
43300
43400
43500
43600
43700
43800
43900
44000
44100
44200
44300
44400
44500
44600
44700
44800
44900
45000
45100
45200
45300
45400
45500
45600
45700
45800
45900
46000
46100
46200
46300
46400
46500
46600
46700
46800
46900
47000
47100
47200
47300
47400
47500
47600
133
-------
C
c
c
c
c
c
c
c
c
c
c
c
J2=K+J-J1
DO 205 I=1,NEQN
205 YCI,J2)=Y(I,J2)+Y(I,J2+1)
DO 210 I=1,NEQN
210 ERRORU>=0.0
ITER=0
C
C
C
C
C
C
W(*,*) IS REEVALUATED ONLY IF CONVERGENCE HAS ALREADY FAILED,
OR IF THERE HAS BEEN A CHANGE OF ORDER OR STEPSIZE. W IS
NOT REEVALUATED DURING THE ITERATION. (NEQN=1 IS TREATED
AS A SPECIAL CASE.)
IF (IWEVAL.LT.1) GO TO 225
CALL FCTCT/Y/.YP)
NFE=NFE+1
CALL FDER(T,Y,W,YP,FCT)
IWEVAL=-1
IF (NEQN.EQ.1) GO TO 220
AHINV=1./(A(1)*H)
DO 215 I=1,NEQN
215 W(I,I)=W(I,I>+AHINV
C
C
C
C
C
C
DECOMP PERFORMS AN LU FACTORIZATION OF W, WITH PARTIAL
PIVOTING. IPS(*) IS USED TO HOLD INFORMATION ON THE ROW
INTERCHANGES. C IS USED FOR WORKING STORAGE, AND J1 IS SET
TO -1 IF W IS FOUND TO BE SINGULAR.
CALL DECOMPCNEQN/W/.IPS/C/J1)
IF (J1.LT.O) GO TO 267
GO TO 230
NEQN=1
220 W(1,1)=1.+W(1,1)*A(1)*H
IF (WdxD.EQ.O.O) GO TO 267
GO TO 255
UP TO 3 ITERATIONS ARE TAKEN. THE ITERATION IS STOPPED WHEN
THE NORM OF THE CORRECTION C(*> IS LESS THAN BND. THE
SUCCESSIVE CORRECTIONS ARE ACCUMULATED IN ERROR(*).
225 CALL FCT(T/Y,YP)
NFE=NFE+1
IF (NEQN.EQ.1) GO TO 255
SOLVE FOR C(*) USING THE FACTORIZED FORM OF W RETURNED BY
DECOMP.
230 IP=IPS(1)
CC1) = (YUP,2)-YPUP)*H)*AHINV
DO 240 I=2/NEQN
IP=IPS(I)
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
47700
47800
47900
48000
48100
48200
48300
48400
48500
48600
48700
48800
48900
49000
49100
49200
49300
49400
49500
49600
49700
49800
49900
50000
50100
50200
50300
50400
50500
50600
50700
50800
50900
51000
51100
51200
51300
51400
51500
51600
51700
51800
51900
52000
52100
52200
52300
52400
52500
52600
52700
52800
134
-------
C
C
C
C
C
C
C
C
C
SUPNO.O
DO 235 J=1/IM1
235 SUM=SUM+W(IP/J)*C(J)
240 CU) = (YUP/2)-YP(IP)*H)*AHINV-SUM
IP=IPS(NEQN)
C(NEQN)=C(NEQN)/WdP/NEQN)
DO 250 IBACK=2/NEQN
I=NEQN-IBACK+1
IP=IPS(I)
IP1=I-H
SUM=0.0
DO 245 J=IP1/NEQN
245 SUM=SUM+W(IP/J)*C(J)
250 C(I)=(C(I)-SUM)/W(IP/I)
GO TO 260
C
C
C
C
C
C
C
C
NEQN=1
255 C(1)==YU/1)+AC1)*C
-------
C
c
c
c
c
c
c
270 Y(Is1)=SAVE(Is1>
GO TO 120
RESTORE Y
275 DO 280 J=1rKP1
DO 280 I=1xNEQN
280 YU,J)=SAVE(I/J)
IF THIS IS THE FIRST CONVERGENCE FAILURE, TRY AGAIN WITH W
REEVALUATED.
IWEVAL=1
IF (ABS(H).GE.HMIN) GO TO 200
C
C CONVERGENCE COULD NOT BE ACHIEVED AT THE SMALLEST POSSIBLE
C KSTEPSIZE. INCREASE EPS TO A VALUE WHICH WILL ALLOW
C CONVERGENCE AND RETURN TO CALLING PROGRAM.
C
GFACT=CNORM/BND
IF (GFACT.LT.2.) GFACT=2.
EPS=EPS*GFACT*(1.+FOURU)
H=SIGN(HMIN,H)
ICRASH=3
RETURN
C
C
C
C
C
CONVERGENCE WAS ACHIEVED.
THE ERROR TOLERANCE.
TEST WHETHER THE SOLUTION MEETS
285 IWEVAL=0
D=0.
DO 290 I=1/NEQN
290 D=AMAX1(D,ABS
-------
C
c
c
c
c
c
c
c
T=TOLD
DO 310 I=1/NEQN
310 YU,1)=SAVE(I,1)
CHECK POSSIBILITIES FOR CONTINUATION.
IFAIL=IFAIL+1
IF (IFAIL.GT.3) GO TO 320
IF <(IFAIL.LT.3).AND.(K.GT.1)) GO TO 325
THIS IS THE THIRD FAILURE, OR ELSE K=1 ON THE FIRST OR SECOND
FAILURE. H IS HALVED/ Y(*,2) IS EVALUATED, AND THE ORDER IS
SET TO 1 FOR THE NEXT ATTEMPT.
IFAIL=3
H=0.5*H
IF (ABS(H).LT.HMIN) H=SIGN(HMIN,H)
CALL FCT(T,Y,YP>
NFE=NFE-H
DO 315 I=1,NEQN
SAVE(I,2)=YPCI)*H
315 Y(I,2)=SAVEU,2)
HOLD=H
K=1
KP1=2
KOLD=1
ITST=2
IWEVAL=1
GO TO 200
C
C IFAIL.GT.3 AND K=1. CHOOSE OPTIMUM STEPSIZE FOR NEXT ATTEMPT.C
C
320 H=H/SQRT(2.*D)
IF (ABS(H).LT.HMIN) GO TO 340
GO TO 120
C
C FIRST OR SECOND FAILURE AND K.GE.2. TEST WHETHER ORDER
C SHOULD BE LOWERED, AND CHOOSE STEPSIZE FOR NEXT ATTEMPT.
C
325 R=1.2*(2.*D)**(1./FLOAT(KPD)
D1=0.
DO 330 I=1,NEQN
330 D1=AMAX1(D1,ABS(Y
-------
C
c
c
c
c
c
c
c
H=H/R
IF (ABS(H).LT.HMIN) GO TO 340
GO TO 120
THE ORDER IS LOWERED. THE STEP IS DECREASED BY AT LEAST 9
PERCENT/ BUT THE DECREASE IS LIMITED TO A FACTOR OF 1/10.
335 K=K-1
IF (RDWN.LT.1.1) RDWN=1.1
IF (RDWN.GT.10.) RDWN=10.
H=H/R
IF (ABS(H).GE.HMIN) GO TO 100
KP1=K+1
ERROR TOLERANCE CANNOT BE MET AT SMALLEST ALLOWED STEPSIZE.
INCREASE TOLERANCE/ RESTORE Y/ AND RETURN TO CALLING PROGRAM.
340 GFACT=(HMIN/ABS(H))**KP1
IF (GFACT.LT.2.) GFACT=2.
EPS=EPS*GFACT*(1.+FOURU)
H=SIGN(HMIN/H)
DO 345 J=2/KP1
DO 345 I=1/NEQN
345 Y(I/J)=SAVE(I/J)
ICRASH=4
RETURN
C***********************
C ******
C END BLOCK 3 ******
C ******
c****************************************************************c
c ******
C BEGIN BLOCK 4 ******
C ******
c***********************
c
THE STEP WAS SUCCESSFUL. STORE THE STEPSIZE USED AND CORRECT
THE REMAINING COMPONENTS OF Y.
C
C
C
C
C
C
400 HOLD=H
IF (KP1.LT.3) GO TO 410
DO 405 J=3/KP1
DO 405 I=1/NEQN
405 YU/J)=YCI/JHA(J)*ERRORCI>
410 ITST=ITST-1
REVISION OF ORDER AND STEPSIZE IS CONSIDERED ONLY IF ITST=0.
IF (ITST.EQ.O) GO TO 420
IF <(ITST.GT.1).OR.(K.EQ.6)> RETURN
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
68500
68600
68700
68800
68900
69000
69100
69200
69300
69400
69500
69600
69700
69800
69900
70000
70100
70200
70300
70400
70500
70600
70700
70800
70900
71000
71100
71200
71300
71400
71500
71600
71700
71800
71900
72000
72100
72200
72300
72400
72500
72600
72700
72800
72900
73000
73100
73200
73300
73400
73500
73600
138
-------
c
C ITST=1 AND K.NE.6 — SAVE ERRORC*) IN Y<*/7). IF NEXT STEP
C SUCCEEDS/ IT WILL BE USED IN DECIDING WHETHER TO RAISE ORDER.
C
DO 415 I=1/NEQN
415 Y(I/7)=ERROR(I)
RETURN
C
C DETERMINE APPROPRIATE ORDER AND STEPSIZE FOR NEXT STEP.
C
420 KD=0
R=1.2*(2.*D)**(1./FLOAT(KPD)
C
C IF K=1/ CANNOT LOWER ORDER. OTHERWISE/ CHECK WHETHER ORDER
C K-1 ALLOWS LARGER STEPSIZE.
C
IF (K.EQ.1) GO TO 430
D1=0.
DO 425 I=1/NEQN
425 D1 =AMAX1 (D1 /ABS (Y (I/KP1 ) /WT (I) ) )
RDWN=1.3*(2.*D1/(ERDWN(K)*EPS»**(1./FLOAT(K))
IF (RDWN.GT.R) GO TO 430
C
C ORDER K-1 IS BETTER THAN ORDER K.
C
R=RDWN
KD=-1
C
C IF K=6/ CANNOT RAISE ORDER. OTHERWISE/ CHECK WHETHER ORDER
C K+1 ALLOWS LARGER STEPSIZE.
C
C
430 IF (K.EQ.6) GO TO 445
D1=0.
DO 435 I=1/NEQN
435 D1=AMAX1 (D1/ABS((ERROR(I)-Y(I/7))/WT(I)))
RUP=1.4*(2.*D1/(ERUP(K)*EPS))**(1./FLOAT(K+2))
IF (RUP.GT.R) GO TO 445
C
C ORDER WILL BE RAISED/ SO APPEND A NEW COLUMN TO Y.
C
KD=1
R=RUP
D1=A(KP1)/FLOAT(KP1)
DO 440 I=1/NEQN
440 Y(I/KPH-1)=ERROR(I)*D1
C
C IF ORDER IS UNCHANGED AND THE STEPSIZE CHANGE WOULD BE LESS
C THAN 12 PERCENT INCREASE/ OR LESS THAN 1 PERCENT DECREASE/ NO
C CHANGE IS MADE. OTHERWISE/ THE NEW ORDER AND STEPSIZE ARE
C SET. ANOTHER TEST WILL BE MADE AFTER K+1 STEPS. STEP
C INCREASE IS LIMITED TO A FACTOR OF 2.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
c
c
73700
73800
73900
74000
74100
74200
74300
74400
74500
74600
74700
74800
74900
75000
75100
75200
75300
75400
75500
75600
75700
75800
75900
76000
76100
76200
76300
76400
76500
76600
76700
76800
76900
77000
77100
77200
77300
77400
77500
77600
77700
77800
77900
78000
78100
78200
78300
78400
78500
78600
78700
78800
139
-------
C C 78900
445 ITST=KP1 79000
IF ((KD.EQ.O).AND.(ABS(R-0.95).LT.0.06)) RETURN 79100
K=K+KD 79200
IF (R.LT.0.5) R=0.5 79300
H=SIGN(AMAX1(ABS(H/R),FOURU*ABS(T))/H) 79400
RETURN ' 79500
C C 79600
C*********************** C 79700
C ****** C 79800
C END BLOCK 4 ****** C 79900
C ****** C 80000
C****************************************************************C 80100
C C 80200
C END OF GEAR C 80300
C C 80400
END 80500
SUBROUTINE DECOMP(N/A/IPS/SCALES/JS) 10000
C C 10100
C****************************************************************C 10200
C C 10300
C THE NUMERICAL SOLUTION OF THE SYSTEM OF LINEAR EQUATIONS C 10400
C AX = B IS VERY SIMPLE ONCE THE N BY N MATRIX A HAS BEEN C 10500
C FACTORED INTO A PRODUCT OF LOWER AND UPPER TRIANGULAR C 10600
C MATRICES. C 10700
C C 10800
C DECOMP USES GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING AND C 10900
C SCALING TO FIND A TRANSPOSITION MATRIX P AND LOWER AND UPPER C 11000
C TRIANGULAR MATRICES L AND U SUCH THAT ALL THE DIAGONAL C 11100
C ELEMENTS OF L ARE 1 AND PLU = A. C 11200
C C 11300
C L AND U ARE STORED IN PLACE OF A AS THE MATRIX P(L-I)+PU C 11400
C = P(L-I+U) C 11500
C C 11600
C THE PERMUTATION MATRIX IS SPECIFIED BY IPS(*). ALL THE C 11700
C ELEMENTS OF THE K-TH COLUMN OF P ARE ZERO/ EXCEPT FOR THE C 11800
C IPS(K)-TH/ WHICH IS 1. (NOTE THAT THE TRANSPOSE OF P IS C 11900
C ITS INVERSE.) C 12000
C C 12100
C SCALES(*) IS USED FOR WORKING STORAGE. C 12200
C C 12300
C DECOMP RETURNS WITH JS=1 IF THE FACTORIZATION IS SUCCESSFUL/ C 12400
C OR WITH JS=-1 IF A IS FOUND TO BE SINGULAR. C 12500
C C 12600
C****************************************************************C 12700
C C 12800
C THIS SUBROUTINE IS WRITTEN IN SINGLE PRECISION. C 12900
C TO OBTAIN A DOUBLE PRECISION VERSION/ REMOVE THE C FROM C 13000
C COLUMN 1 OF THE NEXT FOUR CARDS. C 13100
C DOUBLE PRECISION A/ABS/AMAX1/BIG/DABS/DMAX1/D1/D2/PIVOT/ 13200
140
-------
C 1 ROWNRM/SCALES/SIZE 13300
C ABS(D1)=DABS(D1) 13400
C AMAX1(D1/D2)=DMAX1(D1/D2) 1350Q
C C 13600
DIMENSION A(N,N),IPS(N),SCALES(N) 13700
C C 13800
JS=-1 13900
C C 14000
DO 200 1=1/N 14100
IPS(I)=I 142QO
ROWNRM=0.0 14300
DO 100 J=1xN 14400
100 ROWNRM=AMAX1(ROWNRM/ABS(AU,J)» 14500
IF(ROWNRM.EQ.O.O) RETURN 14600
200 SCALES(I)=1.0/ROWNRM 14700
C C 14800
NM1=N-1 14900
DO 400 1=1,NM1 15000
BIG=0.0 15100
DO 300 J=1,N 15200
JP=IPS(J) 15300
SIZE=ABS(A(JP,I))*SCALES(JP) 15400
IF(SIZE.LE.BIG) GO TO 300 15500
BIG=SIZE 15600
IDXPIV=J 15700
300 CONTINUE 15800
IFCBIG.EQ.0.0) RETURN 15900
IF(IDXPIV.EQ.I) GO TO 350 16000
K=IPS(I) 16100
IPS(I)=IPS(IDXPIV) 16200
IPS(IDXPIV)=K 16300
350 IP=IPS(I) 16400
PIVOT=A(IP/I) 16500
IP1=I+1 16600
DO 400 J=IP1/N 16700
JP=IPS(J) 16800
A(JP/I)=A(JPrI)/PIVOT 16900
DO 400 K=IP1,N 17000
400 A(JP,K>=A(JP,K)-A(JP/I>*A(IP,K) 17100
C C 17200
KP=IPS(N) 17300
IF(A(KP/N).EQ.O.O) RETURN 17400
C C 17500
JS=1 17600
C C 17700
RETURN 17800
C C 17900
END 18000
SUBROUTINE GINTRP(TOUT/YOUT/YPOUT/.Y) 10000
C C 10100
141
-------
c****************************************************************c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c****************************************************************c
c
c
c
c
c
c
c
SUBROUTINE 6INTRP MAY BE USED IN CONJUNCTION WITH SUBROUTINE
GEAR TO INTERPOLATE THE SOLUTION TO THE NONMESH POINT TOUT.
GINTRP SHOULD BE CALLED WHEN THE INTEGRATION HAS PROCEEDED TO
THE FIRST MESH POINT BEYOND TOUT.
SUBROUTINE GEAR APPROXIMATES THE SOLUTION FOR X NEAR T BY A
K DEGREE POLYNOMIAL P(X). THE COEFFICIENT OF THE S**L TERM
OF P(X) IS STORED IN Y(*,L+1), L=0,1,.. .,K, WHERE
S = (X-T)/HOLD
SUBROUTINE GINTRP APPROXIMATES THE SOLUTION AND ITS
DERIVATIVE AT TOUT BY YOUT(*)=P(TOUT) AND YPOUT(*)=P'(TOUT),
RESPECTIVELY.
THE PARAMETERS IN THE ARGUMENT LIST ARE
TOUT — THE VALUE OF THE INDEPENDENT VARIABLE FOR WHICH
THE SOLUTION IS TO BE FOUND.
YOUT — AN ARRAY OF DIMENSION NEQN TO HOLD THE VALUE OF
THE SOLUTION AT TOUT
YPOUT — AN ARRAY OF DIMENSION NEQN TO HOLD THE VALUE OF
THE DERIVATIVE AT TOUT.
Y — THE ARRAY OF DIMENSIONS (NEQN,7) RETURNED AT T BY
SUBROUTINE GEAR.
THE VALUES OF T, HOLD, NEQN, AND K ARE PASSED FROM SUBROUTINE
GEAR THROUGH THE COMMON BLOCK GEDATA.
THIS PROGRAM IS WRITTEN IN SINGLE PRECISION. TO OBTAIN A
DOUBLE PRECISION VERSION, REMOVE THE C FROM COLUMN 1 OF THE
NEXT TWO CARDS.
DOUBLE PRECISION T,HOLD,RVAR,TOUT,YOUT,YPOUT,Y,S,FACT,
1 TEMPI,TEMP2
LOGICAL LVAR
DIMENSION YOUT(NEQN),YPOUT(NEQN),Y(NEQN,7),FACT(7)
COMMON /GEDATA/ T,HOLD,RVAR(10),NEQN,K,IVAR(5),LVAR
S=(TOUT-T)/HOLD
KP1=K+1
KP2=K+2
DO 100 J=2,KP1
100 FACT(J)=FLOAT(J-1)/HOLD
DO 300 1=1,NEQN
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c •
10200
10300
10400
10500
10600
10700
10800
10900
11000
11100
11200
11300
11400
11500
11600
11700
11800
11900
12000
12100
12200
12300
12400
12500
12600
12700
12800
12900
13000
13100
13200
13300
13400
13500
13600
13700
13800
13900
14000
14100
14200
14300
14400
14500
14600
14700
14800
14900
15000
15100
15200
15300
142
-------
TEMPI=0.0
TEMP2=0.0
DO 200 JBACK=1/K
J=KP2-JBACK
TEMP1=TEMP1*S+YU/J)
200 TEMP2=TEMP2*S+Y(I/J)*FACT(J)
YOUT(I)=TEMP1*S+YU/1)
300 YPOUT(I)=TEMP2
RETURN
C
C
C
END OF GINTRP
END
C
C
C
15400
15500
15600
15700
15800
15900
16000
16100
16200
16300
16400
16500
16600
16700
SUBROUTINE QJACOBCT/Y/W/YP/FCT)
C
c****************************************************************c
C
C THIS SUBROUTINE MAY BE USED WITH GEAR TO CALCULATE THE
C JACOBIAN MATRIX OF F(T/Y) BY NUMERICAL DIFFERENCING.
C NORMALLY/ QJACOB SHOULD BE USED ONLY IF THE JACOBIAN IS SO
C COMPLEX AS TO MAKE A DIRECT CALCULATION INCONVENIENT.
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
THIS PROGRAM IS WRITTEN IN SINGLE PRECISION.
TO OBTAIN A DOUBLE PRECISION VERSION/ REMOVE THE C FROM
COLUMN 1 OF THE NEXT FOUR CARDS.
DOUBLE PRECISION T/Y/W/YP/RV1/RV2/EPS/D1/D2/ABS/DABS/
1 AMAX1/DMAX1/DEL/YY/FOURU
ABS(D1)=DABS(D1)
AMAX1(D1/D2)=DMAX1(D1/D2)
LOGICAL LV
COMMON /GEDATA/ RV1(3)/EPS/RV2(8)/NEQN/K/NFE/IV(4)/LV
DIMENSION Y(NEQN)/YP(NEQN)/W(NEQN/NEQN)
IN THE FOLLOWING DATA STATEMENT/ FOURU MUST BE SET TO A
VALUE APPROPRIATE TO THE MACHINE BEING USED.
FOURU IS 4.*U/ WHERE THE COMPUTER UNIT ROUNDOFF ERROR U IS
THE SMALLEST POSITIVE VALUE REPRESENTABLE IN THE MACHINE
SUCH THAT <1.+U).GT.1.
SOME REPRESENTATIVE VALUES ARE
U = 9.54E-7 FOURU = 3.82E-6
U = 5.96E-8 FOURU = 2.39E-7
U = 1.49E-8 FOURU = 5.97E-8
IBM 360/370 S.P.
PDP-11 S.P.
UNIVAC 1108
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
10000
10100
10200
10300
10400
10500
10600
10700
10800
10900
11000
11100
11200
11300
11400
11500
11600
11700
11800
11900
12000
12100
12200
12300
12400
12500
12600
12700
12800
12900
13000
13100
13200
13300
1%400
13500
143
-------
C U = 7.45E-9 FOURU = 2.99E-8 PDP-10 C 13600
C U = 7.11E-15 FOURU = 2.85E-14 CDC 6000 SERIES C 13700
C U = 2.22D-16 FOURU = 8.89D-16 IBM 360/370 D.P. C 13800
C U = 1.39D-17 FOURU = 5.56D-17 PDP-11 D.P. C 13900
C C 14000
DATA FOURU/2.39E-7/ 14100
C C 14200
C C 14300
DO 200 I=1/>NEQN 14400
YY=Y(I) 14500
DEL=AMAX1(EPS**2,AMAX1(EPS,FOURU)*ABS(YY» 14600
Y(I)=Y(I)+DEL 14700
CALL FCT(T,Y,W(1,D) 14800
DO 100 J=1xNEQN 14900
100 W(J,I)=(WCJ,I)-YP(J))/DEL 15000
200 YU)=YY 15100
C C 15200
NFE=NFE+NEQN 15300
C C 15400
RETURN 15500
C C 15600
END 15700
144
-------
APPENDIX D
SOLVING PARTIAL DIFFERENTIAL EQUATIONS
Models describing the transport of substances through continuous media
are naturally represented as initial value problems for systems of partial -
differential equations. A simple example illustrating the essential features
of such a model is
3t
|f = au - bv (b) (D-l)
u(x,t0) = f(x) v(x,t0) = g(x) (c)
Here, u(x,t) may be regarded as the concentration of a chemical species
dissolved in the water of a (one dimensional) stream which flows with
velocity V. v(x,t) is the concentration of the same species bound to the
sediments of the stream bed. The first term on the right hand side of (a)
describes diffusion through the water, the second term accounts for the
motion of the water, and the remaining terms model the exchange of the
species between water and sediment. A more complex and realistic model of
this sort is the mercury model described in Lassiter( 1976) .
A general approach to solving (D-l) is to approximate it by a system of
first order ordinary differential equations. This is done by a process of
semidiscretization :
(1) The spatial dependence of u(x,t) and v(x,t) is represented by u^
and v.(t), which approximate u(x.,t) and v(x^,t) on a set of grid
points x., i=l,2,...,N.
(2) The partial derivatives appearing in (a) are approximated by finite
differences of the u^'s.
The resulting approximation is
= Z>A2(i) - 7AJ(i) - SU + 2nr
vf(t) = au. - bvL (D-2)
u.(to) = f(x.) v.(to) = g(x )
1 x x x i - 1, 2, ..., N
145
-------
where A*(i) and A2(i) represent the finite difference approximations of
3u/3x and 32u/3x2 at (x£,t). (D-2) is a system of 2N coupled first order
ordinary differential equations which may be solved numerically using the
codes in this report, once the forms of A1 and A2 are specified. To keep
considerations simple, suppose that the grid points are chosen to be equally
spaced:
x. = XD + ih
Then several possibilities for the finite differences immediately come to
mind,
for 3u/3x:
u. - u.
A*(i) = Ai(i) = -ii^ (forward difference)
u. - u.
A*(i) = Aa(i) = r (backward difference)
u. - u.
A!/-\ Al/'N 1+1 1-1 f *. -ij'jr.e \
A (i) = A|(i) = ST (central difference)
2n
for 32u/3x2:
u. , - 2u. + u.
A2(i) = A?(i) =
• , «•**• • *•* • -
2/-N , 2,-x 1+1 1 1-1
A2(i) -
A u; -
Considerable care is necessary in choosing which approximation to use. A
particularly bad choice is to use A| for 32u/3x2. The unsuitability of this
approximation can best be appreciated by considering the special case V=Q.
Then the even- and odd-numbered grid points are completely uncoupled, so that
the problem being solved is actually two independent problems, each of which
uses only half of the available information. Since any errors occurring in
the even (odd) points are propagated only among the even (odd) points, a
common manifestation in the solution is the appearance of an unphysical
oscillatory behavior in the spatial variable. (Such oscillations are visible
in the computed solutions shown in the reference cited above.) A similar,
though less easily visualizable, difficulty occurs with the use of AS for
3u/3x, because As(i) is independent of the value of u^. In general, it is
safest to use A? in combination with either Ai or Al, so that (D-l) might be
approximated by either
V f \ Vt \
u. = r-5-(u. -2u.+u. .) - —(u. -u.) - au. + bv.
i h* i+l ii-l n i+l i i i
v' = au. - bv. (D-3a)
111
u.(t0) = f(x.) v.(t0) = g(x )
i = 1, 2, ..., N
146
-------
or
v:=au^- bv± (D-3b)
ui(t0) = f(x£) v (t0) = g(x.)
1 i = 1, 2 ..... N
It is important to understand that the problem presented to the
integration routines is (D-3), which is only an approximation of (D-l).
While requesting increasingly accurate error tolerances will produce increas-
ingly accurate solutions of (D-3), the integration routines cannot reduce the
errors incurred in going from (D-l) to (D-3). Unless care is taken with this
initial approximation, it may be found that the solution of (D-3) has proper-
ties quite different from those expected for the solution of (D-l).
Besides the choice of A1 and A2, the most important factor affecting the
accuracy of the semidiscretization process is the spacing h of the spatial
grid. If u(x,t) or v(x,t) varies rapidly on a scale smaller than ft, then the
finite difference approximations of 3u/3x and 92u/3x2 cannot be accurate, no
matter what form is chosen. In particular, h must be small enough that the
initial conditions f(x) and g(x) are accurately represented. The quest for
increased accuracy can quickly lead to very large systems of equations, since
the number of equations doubles each time h is halved. Computing time and
storage requirements will often limit the accuracy which it is feasible to
achieve.
Additionally, the grid must be finite in extent. In effect, one cannot
approximate the initial value problem (D-l), but only the initial-boundary
value problem which is (D-l) with the additional condition that u(xj,t) =
v(xj,t) = u(xN,t) = v(xN,t) = 0. Thus, the approximation cannot be accurate
unless u(x,t) and v(x,t) are of negligible magnitude near Xj and XN over the
entire interval [t0,T] on which the integration is performed.
Given the model (D-3a) or (D-3b), it is tempting to use it to follow the
fate of a "point spill" of pollutant at x by using the initial conditions
u, (t0) = A, u.(t0) = 0 for i t k
k x (D-4)
v.(t0) = 0 for all i
If this is done for (D-3a), one immediately finds at t0,
D V
uk+l
147
-------
If V is sufficiently large, uiT_i will be negative, so that u^.j will
immediately assume a (physically meaningless) negative value. (The situation
would be even worse if the model had incorporated the differences AS and Af ,
which will always produce a negative value for "k-j, regardless of the value
of V.) For (D-3b), however, the derivatives at t0 are
. D V
u , = TjA + -rA
k+ 1 h2- h
and the problem of negative values does not occur.
The real difficulty here is that the grid spacing h is not sufficiently
small to accurately represent the initial condition. Indeed, the smooth
function f(x) = u(x,t0) which (D-4) is supposed to approximate has not even
been specified. Nevertheless, it may be permissible to ignore the origins of
(D-3b), and to use it and (D-4) as a sufficiently accurate model of the
system being studied. If this is done, it must be borne in mind that the
grid spacing h becomes an essential feature of the model: changing h will
produce a quantitatively different solution.
This discussion has been intended only to illuminate some of the
pitfalls which may be encountered in solving problems similar to (D-l). A
thorough exploration of the numerical solution of partial differential
equations may be found in the text of Ames (1969).
148
-------
/PI TECHNICAL REPORT DATA
(fiease read Instructions on the reverse before completing)
1. REPORT NO.
EPA-600/3-80-037
3. RECIPIENT'S ACCESSION-NO.
4. TITLE AND SUBTITLE
Efficient Algorithms for Solving Systems of Ordinary
Differential Equations for Ecosystem Modeling
5. REPORT DATE
March 1980 issuing date
6. PERFORMING ORGANIZATION CODE
I. AUTHOR(S)
John Malanchuk, John Otis, and Hubert Bouver
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Institute for Man and Environment
State University of New York
Plattsburgh, New York 12901
10. PROGRAM ELEMENT NO.
A28B1A
11. CONTRACT/GRANT NO.
R805452
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Research Laboratory—Athens GA
Office of Research and Development
U.S. Environmental Protection Agency
Athens, Georgia 30605
13. TYPE OF REPORT AND PERIOD COVERED
Final, 7/77-8/78
14. SPONSORING AGENCY CODE
EPA/600/01
15. SUPPLEMENTARY NOTES
16. ABSTRACT
This report presents three packages of subroutines for the numerical solution
of systems of first order ordinary differential equations. The three integration
methods were chosen to provide an efficient means of solving the full range of initial
value problems encountered in basin ecosystem modeling. The subroutines have been
designed to handle all aspects of the numerical integration as automatically as
possible so that their successful use does not require a sophisticated knowledge of
the techniques of numerical analysis.
After a general introduction to the principles and concepts of numerical inte-
gration, the particular methods used in the codes are discussed in some detail. Full
instructions for use are given, and some simple examples are presented. The complete
code listings, in ANSI FORTRAN IV, are contained in appendices.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS
COSATI Field/Group
Systems analysis
Mathematical models
Environments
62B
68D
72E
8. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (ThisReport)
UNCLASSIFIED
1. NO. OF PAGES
155
20. SECURITY CLASS (This page)
UNCLASSIFIED
22. PRICE
EPA Form 2220-1 (9-73)
149
ft US, GOVERNMENT PRINTING OFFICE: 1980-657-146/5637
------- |