C ALGORITHM 687, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 17, NO. 1, PP. 1-10. MARCH, 1991.
*Initial Value Ordinary Differential Equations
1 Do you require an introduction to the initial value problem
tree?
2* The aim of this decision tree is to assist users in the
selection of a suitable code or codes to solve initial value
ordinary differential equations (IVODES), for example, equations
of the form y'=f(t,y), given y at an initial value of t, y(t0) = y0,
say.
Several approaches have been developed for the numerical
solution of initial value problems, of which three have proved
extremely useful, namely those based upon
- Adams formulae (Adams)
- Runge-Kutta formulae (RK)
- Backward differentiation formulae (BDF)
The majority of today's useful IVODE FORTRAN software comprise
instances of the methods represented by each of these formulae.
The decision tree isolates a generic approach that is most
appropriate for a given problem, and then determines which of
these codes, based on this approach, are suitable.
The codes recommended in this tree are selected mainly from the
the HARWELL (edition 6), IMSL (version 9.2), NAG (mark 12) and
SLATEC (version 3) libraries. It should be noted that the code LSODE,
although not available in a library, is considered to be of
library standard. Where few codes of library standard are
available, reliable research codes are suggested where appropriate.
Information on the availability of the codes is included in the tree.
There are a few situations in which currently no reliable software is
available, and seeking expert advice is recommended in these cases.
In situations in which several codes are specified the codes are
given in order of preference, separated by blank lines. In some
cases two or more codes are considered of approximately equal
preference. These are grouped together, with the library codes
given first, in alphabetical order, followed by the research
codes.
The current version of the tree is based mainly on
functionality. It is important to be aware that it is assumed
that the unit roundoff is less than 10**(-12). If this is not
the case then roundoff may dominate and questions in the
decision making process concerning high accuracy may be
inappropriate. If unit roundoff is greater than this value then
you may still use the system to obtain a recommended code.
However, you should be aware that the effects of a large unit
roundoff can be quite serious especially on some of the
recommended codes for "stiff" problems and that the effect may
never have been tested for some of the research codes.
Ease-of-use has been considered when making recommendations.
When an easy-to-use code is not available, we have attempted to
decide which parameter selection, which program outline
(pseudo-code), and which routine to help with intermediate
output (interpolation driver code) will best assist the user.
As the tree is to accommodate users with varying degrees of relevant
IVODE experience "help" information, in the form of explanations of
terminology and guidance on how to answer the questions, is provided
where it seems most appropriate.
The majority of codes for initial value problems are highly
sophisticated and may seem complicated for a new user to apply
correctly. The example calling programs, typically supplied with a
library, e.g. with the NAG library, are a useful way in which to get
accustomed to these codes. It is recommended that a new user should
first attempt to solve the problem (or a substitute "simple" problem)
with the simpler, straightforward codes, such as
- IMSL DVERK with the default options
- NAG D02BAF
- SLATEC DERKF (or RKF45) with default options.
This enables the user to become familiar with library codes and often
it provides the desired solution, as there is a large class of
problems that can be handled successfully by the simpler codes. Even
in situations in which the simpler codes fail, information concerning
the problem may be acquired that will be of assistance when returning
to the tree to seek a more appropriate code.
There are several situations in which the user may need to terminate
an IVODE NIT session at an intermediate point along a decision path.
For example, some reformulation or analysis of the problem may be
required in order to continue further. NIT has a "jump" facility
which can be used to enable the user to start a session from
questions other than the first, so the user can temporarily
abandon a session and then continue later from the "left-off" point.
These points are indicated in the tree by the specification of
the question/box number. The user should use Q to quit the
current session and then use J at the start of the new session
to jump to the relevant question.
It should be emphasized that the tree is dynamic in nature. Current
codes are continually being reviewed and new codes developed. The
authors would hope to review the tree on a regular basis. As more
experience is obtained and as new codes become available, the tree can
be modified to reflect the changes. For example, it is hoped that more
codes will be available to deal with implicit and differential
algebraic problems.
3* Is your problem an explicit system of the form
y' = f(t,y), y(t0)=y0,
or can it be converted easily to this form? (Question/box no. 3)
4 Is your problem stiff? (Question/box no. 4)
5* Is your problem implicit OR differential algebraic OR a
combination of these? For example, is it like one of the following
(i) My' = f(t,y),
(ii) y' = f(t,y,z), 0 = g(t,y,z) ? (Question/box no. 5)
6* Is there a possibility that the Jacobian matrix (df/dy), has
eigenvalues near or on the imaginary axis? For example, this
could occur if y" = c**2 y. (Question/box no. 6)
7* Is the function f(t,y), very expensive to evaluate OR do you
require high accuracy results?
8* Try
- NAG D02NGF, D02NHF or D02NJF in the cases where the
Jacobian (df/dy), of your differential system is full,
banded or sparse respectively.
- DASSL (a research code).
If neither is available try
- LSODI (a research code) in the linearly implicit case,
e.g. My' = f(t,y).
- HARWELL DC03 in the explicit mixed differential case,
e.g. y' = f(t,y,z), 0 = g(t,y,z).
Do you wish to see
9 Seek advice.
10 Is the system small, i.e. less than 20 equations?
11 Is the system small, i.e. less than 20 equations?
12* Does the definition of the differential equation depend on a
"status switch" (or a "status vector") which can change during
the integration, and are all circumstances where the status
switch can change values known explicitly? (These circumstances
determine collectively the zeros of an associated switching
function which is assumed to be available.) (Question/box
no. 12)
13* Does the derivative function f(t,y) have any potential
discontinuities? (Question/box no. 13)
14* Does the definition of the differential equation depend on a
"status switch" (or a "status vector") which can change during
the integration, and are all circumstances where the status
switch can change values known explicitly? (These circumstances
determine collectively the zeros of an associated switching
function which is assumed to be available.) (Question\box
no. 14)
15* Can the full Jacobian matrix fit in main memory?
(Question/box no. 15)
16* Does the definition of the differential equation depend on a
"status switch" (or a "status vector") which can change during
the integration, and are all circumstances where the status
switch can change values known explicitly? (These circumstances
determine collectively the zeros of an associated switching
function which is assumed to be available.) (Question/box
no. 16)
17 Can the full Jacobian matrix fit in main memory? (Question/box
no. 17)
18* Does the definition of the differential equation depend on a
"status switch" (or a "status vector") which can change during
the integration, and are all circumstances where the status
switch can change values known explicitly? (These circumstances
determine collectively the zeros of an associated switching
function which is assumed to be available.) (Question/box
no. 18)
19* Do you wish to find the points at which a function g(t,y,y') is
zero, for example when one of the components of y takes a known
value?
20 Is the Jacobian matrix dense, i.e. more than 5% of the entries
non-zero? (Question/box n. 20)
21* Is the Jacobian matrix banded and will about 2mn storage
locations fit in main memory, where m is the bandwidth and n is
the number of equations in the system? (Question/box no. 21)
22 Is the Jacobian matrix dense, i.e. more than 5% of the entries
non-zero? (Question/box no. 22)
23* Is the Jacobian matrix banded and will about 2mn storage
locations fit in main memory, where m is the bandwidth and n is
the number of equations in the system? (Question/box no. 23)
24 Are there many discontinuities in low order derivatives?
(Question/box no. 24)
25* Do you wish to find the points at which a function g(t,y,y') is
zero, for example when one of the components of y takes a known
value?
26 Does the Jacobian matrix have a banded or staircase structure?
(Question/box no. 26)
27* Does the definition of the differential equation depend on a
"status switch" (or a "status vector") which can change during
the integration, and are all circumstances where the status
switch can change values known explicitly? (These circumstances
determine collectively the zeros of an associated switching
function which is assumed to be available.) (Question/box
no. 27)
28* Is the number of non-zero entries of the Jacobian
sufficiently small that about 5 times this number easily
fits in main memory? (Question/box no. 28)
29 Does the Jacobian matrix have a banded or staircase structure?
30* Does the definition of the differential equation depend on a
"status switch" (or a "status vector") which can change during
the integration, and are all circumstances where the status
switch can change values known explicitly? (These circumstances
determine collectively the zeros of an associated switching
function which is assumed to be available.) (Question/box
no. 30)
31* Is the number of non-zero entries of the Jacobian
sufficiently small that about 5 times this number easily
fits in main memory? (Question/box no. 31)
32* Try a Runge-Kutta code with IVODE pseudocode 3
- IMSL DVERK with interpolant given in Enright et al.
"Effective solution of discontinuous IVPs using a Runge-
Kutta formula pair with interpolants", University of
Manchester Numerical Analysis Report No. 113, 1986.
- NAG D02PAF with interpolation routine D02XAF or D02XBF.
Do you wish to see
33* Try an Adams code
- RDEAM (a research code)
or with IVODE pseudocode 3 try
- SLATEC DEABM with IVODE driver 1 for interpolation.
- LSODE with MF = 10 and IVODE driver 3 for interpolation.
- NAG D02QAF with interpolation routine D02XGF or D02XHF.
Do you wish to see
34* If there are frequent discontinuities the following codes may
not be suitable. In this case back up the decision path and
follow the path for inexpensive function evaluations, or seek
advice.
Try an Adams code
- RDEAM (a research code).
- NAG D02CHF if only the first such point is required.
- LSODAR (a research code).
or with IVODE pseudocode 1
- SLATEC DEABM with IVODE driver 1 for interpolation.
- LSODE with MF = 10 and IVODE driver 3 for interpolation.
- NAG D02QAF with interpolation routine D02XGF or D02XHF.
Do you wish to see
35* Try an Adams code
- SLATEC DEABM.
- NAG D02CBF.
- LSODE with MF = 10.
- IMSL DGEAR with METH = 1 and MITER = 0.
Do you wish to see
36* Try a Runge-Kutta code
- NAG D02BHF if only one such point is required.
or with IVODE pseudocode 1
- IMSL DVERK with interpolant given in Enright et al's
"Interpolants for Runge-Kutta Formulas", ACM T.O.M.S. 12,
3, p. 193-218.
- NAG D02PAF with interpolation routine D02XAF or D02XBF and
root finder C05AZF.
Do you wish to see
37* Try a Runge-Kutta code
- NAG D02BBF.
- IMSL DVERK with interpolant given in Enright et al's
"Interpolants for Runge-Kutta Formulas", ACM T.O.M.S. 12,
3, p. 193-218.
- NAG D02PAF with interpolation routine D02XAF or D02XBF.
Do you wish to see
38* Does the definition of the differential equation depend on a
"status switch" (or a "status vector") which can change during
the integration, and are all circumstances where the status
switch can change values known explicitly? (These circumstances
determine collectively the zeros of an associated switching
function which is assumed to be available.) (Question/box
no. 38)
39* See an expert about using out-of-core linear algebra, if
available, or a method based on an iterative solution of the
algebraic equations. Otherwise, try the non-stiff decision path
(if the problem is only mildly stiff this may be satisfactory).
40* Does the definition of the differential equation depend on a
"status switch" (or a "status vector") which can change during
the integration, and are all circumstances where the status
switch can change values known explicitly? (These circumstances
determine collectively the zeros of an associated switching
function which is assumed to be available.) (Question/box
no. 40)
41* See an expert about using out-of-core linear algebra, if
available, or a method based on an iterative solution of the
algebraic equations. Otherwise, try the non-stiff decision path
(if the problem is only mildly stiff this may be satisfactory).
42* Try a Runge-Kutta code with IVODE pseudocode 2
- IMSL DVERK with interpolant given in Enright et al.
"Effective solution of discontinuous IVPs using a
Runge-Kutta formula pair with interpolants", University of
Manchester Numerical Analysis Report No. 113, 1986.
- NAG D02PAF.
Do you wish to see
43* With IVODE pseudocode 3 try
- NAG D02NBF with setup routine D02NWF to obtain the blended
formulae, setting MORDER = 4, and with interpolation
routine D02XJF.
- SECDER (a research code).
or try a BDF code with maximum order 2
- HARWELL DC03 with KMAXX = 2 in common block DC03Z.
- NAG D02NBF with MORDER = 2 in setup routine D02NVF, and
with interpolation routine D02XKF.
- SLATEC DEBDF with MAXORD = 2 in common block DEBDF1, and
IVODE driver 2 for interpolation.
- LSODE with MF = 21 or 22, IWORK(5) = 2, and IVODE driver 3
for interpolation.
Do you wish to see
44* Do you wish to find the points at which a function g(t,y,y')
is zero, for example when a component of y takes a known value?
45* With IVODE pseudocode 3 try
- NAG D02NCF with setup routine D02NWF to obtain the blended
formulae, setting MORDER = 4, and with interpolation
routine D02XJF.
or try a banded BDF code with maximum order 2
- HARWELL DC03 with KMAXX = 2 in common block DC03Z, and
INFO(6) = 1.
- NAG D02NCF with MORDER = 2 in setup routine D02NVF, and
with interpolation routine D02XKF.
- SLATEC DEBDF with INFO(6) = 1, MAXORD = 2 in common block
DEBDF1, and IVODE driver 2 for interpolation.
- LSODE with MF = 24 or 25, IWORK(5) = 2, and IVODE driver 3
for interpolation.
Do you wish to see
46* Do you wish to find the points at which a function g(t,y,y')
is zero, for example when a component of y takes a known value?
47* With IVODE pseudocode 3 try
- NAG D02NDF with setup routine D02NWF to obtain the blended
formulae, setting MORDER = 4, and with interpolation
routine D02XJF.
or try a sparse BDF code with maximum order 2
- HARWELL DC03 with KMAXX = 2 in common block DC03Z, and
INFO(6) = 2 or 3.
- NAG D02NDF with MORDER = 2 in setup routine D02NVF, and
with interpolation routine D02XKF.
- LSODES (a research code) with IWORK(5) = 2, and IVODE
driver 3 for interpolation.
Do you wish to see
48* Do you wish to find the points at which a function g(t,y,y')
is zero, for example when a component of y takes a known value?
49* Try a BDF code with IVODE pseudocode 3
- HARWELL DC03
- NAG D02NBF with setup routine D02NVF and with interpolation
routine D02XKF.
- SLATEC DEBDF with IVODE driver 2 for interpolation.
- LSODE with MF = 21 or 22 and IVODE driver 3 for
interpolation.
Do you wish to see
50* Do you wish to find the points at which a function g(t,y,y')
is zero, for example when a component of y takes a known value?
51* Try a banded BDF code with IVODE pseudocode 3
- HARWELL DC03 with INFO(6) = 1.
- NAG D02NCF with setup routine D02NVF and with interpolation
routine D02XKF.
- SLATEC DEBDF with INFO(6) = 1 and IVODE driver 2 for
interpolation.
- LSODE with MF = 24 or 25 and IVODE driver 3 for
interpolation.
Do you wish to see
52* Do you wish to find the points at which a function g(t,y,y')
is zero, for example when a component of y takes a known value?
53* Try a sparse BDF code with IVODE pseudocode 3
- HARWELL DC03 with INFO(6) = 2 or 3.
- NAG D02NDF with setup routine D02NVF and with interpolation
routine D02XKF.
- LSODES (a research code) with IVODE driver 3 for
interpolation.
Do you wish to see
54* Do you wish to find the points at which a function g(t,y,y')
is zero, for example when a component of y takes a known value?
55* Try
- NAG D02NBF with setup routine D02NWF to obtain the blended
formulae, setting MORDER = 4, with interpolation routine
D02XJF and IVODE pseudocode 1.
- SECDER (a research code) with IVODE pseudocode 1
or try a BDF code with maximum order 2
- RDEBD (a research code) with MAXORD = 2 in common block
RDEBD1.
- LSODAR (a research code) with JT = 1 or 2, IWORK(9) = 2.
or with IVODE pseudocode 1
- HARWELL DC03 with KMAXX = 2 in common block DC03Z.
- NAG D02NBF with MORDER = 2 in setup routine D02NVF, and
with interpolation routine D02XKF.
- SLATEC DEBDF with MAXORD = 2 in common block DEBDF1 and
IVODE driver 2 for interpolation.
- LSODE with MF = 21 or 22, IWORK(5) = 2, and IVODE driver 3
for interpolation.
Do you wish to see
56* Try
- NAG D02NBF with setup routine D02NWF to obtain the blended
formulae, setting MORDER = 4.
- SECDER or STRIDE (research codes)
or try a BDF code with maximum order 2
- HARWELL DC03 with KMAXX = 2 in common block DC03Z.
- NAG D02NBF with MORDER = 2 in setup routine D02NVF.
- SLATEC DEBDF with MAXORD = 2 in common block DEBDF1.
- LSODE with MF = 21 or 22, IWORK(5) = 2.
Do you wish to see
57* Try
- NAG D02NCF with setup routine D02NWF to obtain the blended
formulae, setting MORDER = 4, with interpolation routine
D02XJF and IVODE pseudocode 1.
or try a banded BDF code with maximum order 2
- RDEBD (a research code) with INFO(6) = 1, and MAXORD = 2 in
common block RDEBD1.
- LSODAR (a research code) with JT = 4 or 5, IWORK(9) = 2.
or with IVODE pseudocode 1
- HARWELL DC03 with KMAXX = 2 in common block DC03Z and
INFO(6) = 1.
- NAG D02NCF with MORDER = 2 in setup routine D02NVF, and
with interpolation routine D02XKF.
- SLATEC DEBDF with INFO(6) = 1, MAXORD = 2 in common block
DEBDF1 and IVODE driver 2 for interpolation.
- LSODE with MF = 24 or 25, IWORK(5) = 2, and IVODE driver 3
for interpolation.
Do you wish to see
58* Try
- NAG D02NCF with setup routine D02NWF to obtain the blended
formulae, setting MORDER = 4.
or try a banded BDF code with maximum order 2
- HARWELL DC03 with KMAXX = 2 in common block DC03Z and
INFO(6) = 1.
- NAG D02NCF with MORDER = 2 in setup routine D02NVF.
- SLATEC DEBDF with INFO(6) = 1, MAXORD = 2 in common block
DEBDF1.
- LSODE with MF = 24 or 25, IWORK(5) = 2.
Do you wish to see
59* With IVODE pseudocode 1 try
- NAG D02NDF with setup routine D02NWF to obtain the blended
formulae, setting MORDER = 4, with interpolation routine
D02XJF.
or try a sparse BDF code with maximum order 2
- HARWELL DC03 with KMAXX = 2 in common block DC03Z and
INFO(6) = 2 or 3.
- NAG D02NDF with MORDER = 2 in setup routine D02NVF, and
with interpolation routine D02XKF.
- LSODES (a research code) with IWORK(5) = 2 and IVODE driver
3 for interpolation.
Do you wish to see
60* Try
- NAG D02NDF with setup routine D02NWF to obtain the blended
formulae, setting MORDER = 4.
or try a sparse BDF code with maximum order 2
- HARWELL DC03 with KMAXX = 2 in common block DC03Z and
INFO(6) = 2 or 3.
- NAG D02NDF with MORDER = 2 in setup routine D02NVF.
- LSODES (a research code) with IWORK(5) = 2.
Do you wish to see
61* Try a BDF code
- NAG D02EJF if only the first such point is required.
- RDEBD (a research code).
- LSODAR (a research code) with JT = 1 or 2.
or with IVODE pseudocode 1
- HARWELL DC03.
- NAG D02NBF with setup routine D02NVF, and with
interpolation routine D02XKF.
- SLATEC DEBDF with IVODE driver 2 for interpolation.
- LSODE with MF = 21 or 22, and IVODE driver 3 for
interpolation.
Do you wish to see
62* Try a BDF code
- NAG D02EJF.
- HARWELL DC03.
- NAG D02NBF with setup routine D02NVF.
- SLATEC DEBDF.
- LSODE with MF = 21 or 22.
- IMSL DGEAR with METH = 2 and MITER = 1, 2 or 3.
Do you wish to see
63* Try a banded BDF code
- RDEBD (a research code) with INFO(6) = 1.
- LSODAR (a research code) with JT = 4 or 5.
or with IVODE pseudocode 1
- HARWELL DC03 with INFO(6) = 1.
- NAG D02NCF with setup routine D02NVF, and with
interpolation routine D02XKF.
- SLATEC DEBDF with INFO(6) = 1, and IVODE driver 2 for
interpolation.
- LSODE with MF = 24 or 25, and IVODE driver 3 for
interpolation.
Do you wish to see
64* Try a banded BDF code
- HARWELL DC03 with INFO(6) = 1.
- NAG D02NCF with setup routine D02NVF.
- SLATEC DEBDF with INFO(6) = 1.
- LSODE with MF = 24 or 25.
- IMSL DGEAR with METH = 2, MITER = -1 or -2.
Do you wish to see
65* Try a sparse BDF code with IVODE pseudocode 1
- HARWELL DC03 with INFO(6) = 2 or 3.
- NAG D02NDF with setup routine D02NVF, and with
interpolation routine D02XKF.
- LSODES (a research code) with IVODE driver 3 for
interpolation.
Do you wish to see
66* Try a sparse BDF code
- HARWELL DC03 with INFO(6) = 2 or 3.
- NAG D02NDF with setup routine D02NVF.
- LSODES (a research code).
Do you wish to see
801* PSEUDO-CODE 1.
Pseudo-code to locate the values of t at which
the function g(t,y,y') is zero.
*
DRIVER
*
set initial conditions of differential equation
*
determine initial sign of g(t,y,y')
*
set options of routine to return after every attempted step
*
determine initial step h
*
repeat until t = tend or integration completed or error
condition detected
*
invoke routine attempting a step from t to t + h
returns with step = h-bar <= h taken, y, y' and
recommended next step h-star
*
if g(t,y,y') has changed sign then
*
using interpolant z(t), find s such that
g(s,z(s),z'(s)) = 0
*
if integration to be terminated then
*
set "integration completed" indicator
*
else
*
take care with initial setting of sign of g and value
of y'
*
set h = h-star and signal a restart
*
end if
*
else
*
h = h-star
*
end if
*
end repeat
*
end DRIVER
*
802* PSEUDO-CODE 2.
*
Pseudo-code to handle discontinuities without the
use of an associated "switching" function.
*
DRIVER
*
set initial conditions of differential equation
*
set options of routine to return after every step
*
repeat until t = tend or integration completed or error
condition detected
*
invoke routine attempting a step from t to t + h
returns with h-bar <= h taken, y, y', and recommended next
step h-star
*
if h-bar < h and h-star >= 2h-bar then
*
using interpolant z(t), search for a sufficiently
accurate bracketing (TL,TH), of the discontinuity by
sampling z'(t) - f(t,z(t))
*
if search successful then
*
set t = TH, y = z(TH)
*
if integration to be continued then
*
set the next attempted step to be h and signal a
restart
*
else
*
set "integration completed" indicator
*
end if
*
end if
*
end if
*
end repeat
*
end DRIVER
*
803* PSEUDO-CODE 3.
*
Pseudo-code to handle discontinuities of the
"switching" type.
*
DRIVER
*
set initial conditions of differential equation
*
set initial status vector
*
set options of routine to return after every step
*
repeat until t = tend or integration completed or error
condition detected
*
invoke routine attempting a step from t to t + h
*
returns with h-bar <= h taken, y, y' and recommended next
step h-star
if g(t,y) has changed sign then
*
using interpolant z(t), find TL and TH such that
s is in (TL,TH) and g(s,z(s)) = 0, and (TH - TL)
is sufficiently small
*
if integration to be continued then
*
set t = TH, y = z(TH) and signal restart
*
set status vector to redefine differential equation
reset definition of g(t,y) if necessary
*
else
*
set "integration completed" indicator
*
end if
*
end if
*
end repeat
*
end DRIVER
*
804* Driver 1: Interpolation driver for DEABM
*
SUBROUTINE INTRPD(NEQ, OUTX, OUTY, OUTYP, RWORK, LRW,
+ IWORK, LIW)
*
* INTRPD is a driver for interpolation for the DEABM code.
* It returns approximate values of y(x) and y'(x) for x in
* the current interval.
*
* Parameters:
*
* On entry:
*
* NEQ -- the number of equations
* OUTX -- the value of the independent variable x, where
* the interpolation is required
* RWORK -- the real work area RWORK used in DEABM
* LRW -- the declared length of RWORK (in user's
* dimension)
* IWORK -- the integer work area used in DEABM
* LIW -- the declared length of IWORK (in user's
* dimension)
*
* On return:
*
* OUTY -- the interpolated values y (x), i=1,2,..,NEQ
* i
* OUTYP -- the interpolated derivative values
* y '(x), i=1,2,...,NEQ
* i
*
DIMENSION OUTY(NEQ), OUTYP(NEQ), RWORK(LRW), IWORK(LIW)
*
IYY = 2*NEQ + 22
IPHI = 3*NEQ + IYY
IPSI = 16*NEQ + 24 + IPHI
*
CALL INTRP(RWORK(13), RWORK(IYY), OUTX, OUTY, OUTYP,
+ NEQ, IWORK(28), RWORK(IPHI), RWORK(IPSI))
*
RETURN
*
END
*
805* Driver 2: Interpolation driver for DEBDF.
*
*
SUBROUTINE INTYDD(NEQ, OUTX, OUTY, OUTYP, RWORK, LRW,
+ IFLAG)
*
* INTYDD is a driver for interpolation for the DEBDF code.
* It returns the approximate values of y(x) and y'(x) for
* values of x in the current interval.
*
* Parameters:
*
* On entry:
*
* NEQ -- the number of equations
* OUTX -- the value of the independent variable x where
* interpolation is required
* RWORK -- the real work area used in DEBDF
* LRW -- declared length of RWORK (in user's
* dimension)
*
* On return:
*
* OUTY -- the interpolated values y (x), i = 1,2,...,NEQ
* i
* OUTYP -- the interpolated derivative values
* y '(x), i=1,2,..,NEQ
* i
* IFLAG -- an error flag:
*
* on return IFLAG = 0 if no error occurred
* IFLAG = -1 if x is invalid, i.e. outside
* the current interval
*
DIMENSION OUTY(NEQ), OUTYP(NEQ), RWORK(LRW)
*
IYH = NEQ + 240
*
CALL INTYD(OUTX, 0, RWORK(IYH), NEQ, OUTY, IFLAG)
*
IF (IFLAG .EQ. 0) THEN
*
CALL INTYD(OUTX, 1, RWORK(IYH), NEQ, OUTYP, IFLAG)
*
ELSE
*
* OUTX is not valid
*
IFLAG = -1
*
END IF
*
RETURN
*
END
*
806* Driver 3: Interpolation driver for LSODE or LSODES
*
SUBROUTINE INTDYD(NEQ, OUTX, OUTY, OUTYP, RWORK, LRW,
+ IFLAG)
*
* INTDYD is a driver for interpolation for the LSODE family
* of codes. It returns the approximate values of y(x) and
* y'(x) for values of x in the current interval.
*
* Parameters:
*
* On entry:
*
* NEQ -- the number of equations
* OUTX -- the value of the independent variable x where
* interpolation is required
* RWORK -- the real work area used in LSODE
* LRW -- declared length of RWORK (in user's
* dimension)
*
* On return:
*
* OUTY -- the interpolated values y (x), i = 1,2,...,NEQ
* i
* OUTYP -- the interpolated derivative values
* y'(x), i = 1,2,...,NEQ
* i
* IFLAG -- an error flag:
*
* on return IFLAG = 0 if no error occurred
* IFLAG = -1 if x is invalid, i.e. outside the
* current interval
*
DIMENSION OUTY(NEQ), OUTYP(NEQ), RWORK(LRW)
*
CALL INTDY(OUTX, 0, RWORK(21), NEQ, OUTY, IFLAG)
*
IF (IFLAG .EQ. 0) THEN
*
CALL INTDY(OUTX, 1, RWORK(21), NEQ, OUTYP, IFLAG)
*
ELSE
*
* OUTX is not valid
*
IFLAG = -1
*
END IF
*
RETURN
*
END
*
903* Most initial value routines require the equations to be in
explicit first-order form y' = f(t,y). A more general explicit
system may contain derivatives of higher order,
(k) (k-1)
C y = f(t,y,y',y'',...,y ),
*
but it can then be converted to first-order by introducing new
variables. E.g. consider the third-order explicit equation
*
z''' + (1-z)z'' + (zz')**2 = 0.
*
Let y1 = z, y2 = z', y3 = z'', then the following first-order
system is obtained:
*
y1' = y2,
y2' = y3,
y3' = -(1-y1)y3 - (y1y2)**2.
*
For most problems, converting to an equivalent first-order
system does not lead to significant extra work, but in some
special cases a direct approach can have advantages. Examples of
such non-standard problems would be higher-order differential
equations where f depends only on t and y. All the codes
recommended here require conversion to a first order system.
*
Complex equations to be integrated along the real line may
be converted into systems of real equations by splitting each
equation into its real and imaginary parts and solving for real
and imaginary parts of each component of the solution. Note that
this process doubles the size of the system and may not always
be appropriate, but we recommend no codes for integrating
complex equations directly.
*
If your system is not explicit but can be transformed to
explicit form inexpensively, then it may be treated as an
explicit problem. For example a problem of the form
*
My' = f(t,y)
*
where M is constant and diagonal (or has some other simple
structure) may be treated by calling a routine for explicit
equations and solving for y' in the user defined subroutine
for evaluating y'.
*
904* Stiffness is a characteristic present in differential
equations arising from models of real systems where there are
interactions taking place on more than one time scale. One
example would be in a chemical kinetics model where some
transient reactions are occurring on a time scale of a few
microseconds or less, while a slower steady-state reaction
takes place on a time scale of a second or more. Also ordinary
differential equations that arise from the method of lines
applied to parabolic problems are often mildly stiff.
*
These problems are mathematically well-conditioned in that their
solutions are generally smooth and small changes in the initial
conditions will lead to small changes in the solution.
Unfortunately, some numerical techniques, such as explicit
Runge-Kutta or Adams, are inappropriate for this class of
problem since, in order to assure numerical stability, the
step-size of these methods must be severely constrained.
*
One should suspect stiffness when a problem is very expensive to
integrate with standard methods and the cost of integration (or
average step-size required) seems to be relatively insensitive
to the accuracy requested. Stiffness is inevitable if the
eigenvalues of the Jacobian matrix (df/dy) of the differential
equation have negative real parts and
*
1/|Re(lambda)| < < the range of integration,
*
where lambda is the eigenvalue with largest negative real part.
*
Methods designed for stiff systems allow much larger step-sizes
to be used outside the transient region but require more work
per step. The cost of these methods is often dominated by the
cost of setting-up and solving linear algebraic equations at
each step. For this reason, when solving large systems of stiff
equations, it is important that the linear algebraic techniques
exploit any structure that may exist in the problem.
*
If the stiffness properties of a problem are not known, one of
the following codes could be used to help determine whether or
not the problem is stiff.
*
- NAG D02BDF with STIFF not equal to 0.0 for stiffness check
- SLATEC DERKF, checking for stiffness using IDID.
*
If a problem is stiff, but was not originally an explicit
problem, treat it in its original form for the purposes of
solution unless transformation to explicit form is very
inexpensive.
*
905* Implicit systems of the form
*
F(t,y,y') = 0,
*
are considered to be differential algebraic rather than merely
implicit if the matrix (dF/dy') is structurally singular, i.e.
singular for every choice of y and y'.
*
906* For instance, this happens frequently in the following
applications:
*
- the problem arises from the method of lines space
discretisation of a hyperbolic partial differential equation
e.g. a wave equation;
*
- the problem arises from structural dynamics,
*
e.g. My'' + Cy' + Ky = f,
y(a) = ya, y'(a) = ypa;
*
- there is a likelihood of high frequency oscillations.
*
The poor performance of high order BDF algorithms on such
problems has been acknowledged, and now this is an active
research area with several approaches and codes under
development.
*
If the whereabouts of the eigenvalues is not known, assume it is
close to the imaginary axis, or attempt to compute the
eigenvalues for some fixed t.
*
907* As a guide, the derivatives can be considered expensive if
the number of arithmetic operations (or equivalents e.g.
elementary mathematical functions, data accesses) per equation
is greater than 50.
*
In this context, high accuracy means k > 6 significant figures
are desired in the solution components and the unit roundoff,
uround <= 10**(-(k+5)).
*
908* LSODI handles linearly implicit equations of the form
*
M(t,y)y' = g(t,y).
*
HARWELL DC03 handles explicit mixed differential algebraic
equations of the form
*
0 = f (t,y ,...,y ), i = 1,...,NALGB
i 1 NEQ
*
y' = f (t,y ,...,y ), i = NALGB+1,...,NEQ.
i i 1 NEQ
*
912* A differential equation which can be presented in this form
has a potential discontinuity at these switching points. Such
discontinuites can be handled effectively provided one exploits
the fact that a "switching function", is known.
*
An example of such a problem where the switch depends on the
dependent variable is:
*
{ y for SW = 1
y' = {
{ -y/2 for SW = 0
*
where SW = 1 initially and is reset to 0 when y - 2 >= 0
*
and when SW = 0 it is reset to 1 whenever y - 1 <= 0.
*
If SW is reset depending only on the values of the independent
variable t then the values of the switching points should be
determined a priori and the integration restarted at these
predetermined points. When SW depends on the dependent variable
the switching points should be determined as the integration
proceeds and the integration restarted at each point. (See
C.W. Gear "Root finding in simulation and ODE solution",
Dept. of Computer Science Report No. UIUCDCS-R-82-1116, University
of Illinois at Urbana-Champaign, 1982.)
*
Sometimes it is inconvenient or even impossible to present problems
with potential discontinuities in this form (the associated
switching function may be defined only implicitly or it may have
a large number of zeros). In these cases a negative response at
this stage would be appropriate.
*
921* A matrix J, is banded if the non-zero elements only appear
along diagonal bands about the main diagonal from top left to
bottom right. The bandwidth m is the smallest integer such that
*
|i-j| > m => J(i,j) = 0.
*
Sometimes the bands contain zeros and the matrix appears to have
a staircase structure. In these situations it is important that
the bandwidth, as defined above, be determined correctly for use
with banded routines.
*
924* Many discontinuities of low order would imply frequent
integration restarts and hence result in a severe impact on
efficiency.
*
943* SLATEC DEBDF or RDEBD (a research code): to set the maximum
order to 2, the following common block should be declared in the
user's calling routine
*
COMMON/DEBDF1/RARR(218),IARR1(26),MAXORD,IARR2(6) for DEBDF
*
COMMON/RDEBD1/RARR(218),IARR1(26),MAXORD,IARR2(6) for RDEBD1
*
DEBDF1 or RDEBD should be called in one step mode, and MAXORD
assigned the value 2 after the first step. No other variables in
the common block should be changed.
*
HARWELL DC03: to set the maximum order to 2, the common
block DC03Z, given in the library documentation, should be
declared in the user's calling routine, and KMAXX assigned the
value 2 before the first call to DC03. No other variables in the
common block should be changed.
*
947* HARWELL DC03: to set the maximum order to 2, the common
block DC03Z, given in the library documentation, should be
declared in the user's calling routine, and KMAXX assigned the
value 2 before the first call to DC03. No other variables in the
common block should be changed.
*
999* Code References
*
* HARWELL Subroutine Library, Computing Division, Harwell
Laboratory, Oxfordshire, England.
*
HARWELL DC03
*
It is possible to obtain DC03 without obtaining the entire
library.
*
* IMSL Library, 7500 Bellaire Boulevard, Houston, Texas,
77036-5085.
*
IMSL DVERK, DGEAR
*
* NAG Library, 256 Banbury Road, Oxford OX2 7DE, England.
*
NAG D02BAF, D02BBF, D02BDF, D02BHF, D02CBF, D02CHF, D02EJF,
D02NBF, D02NCF, D02NDF, D02NGF, D02NHF, D02NJF, D02PAF, D02QAF
*
* SLATEC Library, Los Alamos Scientific Laboratory, Los Alamos,
Mew Mexico 87545.
*
SLATEC DEABM, DEBDF, DERKF
*
The SLATEC codes are available through NETLIB.
*
* DASSL - L. R. Petzold, Sandia National Laboratories,
Livermore, CA 94550.
*
Reference: A description of DASSL: A differential/algebraic
system solver, Proc. IMACS World Congress,
Montreal, Canada, August 1982.
*
* LSODAR - L. R. Petzold, Sandia National Laboratories,
Livermore, CA 94550.
- A. C. Hindmarsh, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
*
Reference: Automatic selection of methods for solving stiff
and non-stiff systems of ordinary differential
equations, Sandia National Laboratories Report
SAND80-8230, September 1980.
*
ODEPACK, a systemised collection of ODE solvers, in
Scientific Computing, R.S. Stepleman et al. (eds.),
North-Holland, Amsterdam, 1983 (vol 1 of IMACS
Trans. on Scientific Computation), pp. 55-64.
*
* LSODE - A. C. Hindmarsh, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
*
Reference: LSODE and LSODI: Two new initial value ordinary
differential equation solvers, ACM Signum
Newsletter, vol 15, no. 4 (1980) pp. 10-11.
*
ODEPACK, a systemised collection of ODE solvers, in
Scientific Computing, R.S. Stepleman et al. (eds.),
North-Holland, Amsterdam, 1983 (vol 1 of IMACS
Trans. on Scientific Computation), pp. 55-64.
*
* LSODES - A. C. Hindmarsh, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
- A. H. Sherman, J.S. Nolen & Associates, Houston,
Texas.
*
Reference: GEARS: A package for the solution of sparse stiff
ordinary differential equations, in A.M. Erisman et
al. (eds.), Electrical Power Problems: The
Mathematical Challenge, SIAM, Philadelphia, 1980,
pp. 190-200.
*
ODEPACK, a systemised collection of ODE solvers, in
Scientific Computing, R.S. Stepleman et al. (eds.),
North-Holland, Amsterdam, 1983 (vol 1 of IMACS
Trans. on Scientific Computation), pp. 55-64.
*
* LSODI - J. F. Painter, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
- A. C. Hindmarsh, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
*
Reference: LSODE and LSODI: Two new initial value ordinary
differential equation solvers, ACM Signum
Newsletter, vol 15, no. 4 (1980) pp. 10-11.
*
ODEPACK, a systemised collection of ODE solvers, in
Scientific Computing, R.S. Stepleman et al. (eds.),
North-Holland, Amsterdam, 1983 (vol 1 of IMACS
Trans. on Scientific Computation), pp. 55-64.
*
* RDEAM - H. A. Watts, Sandia National Laboratories,
Albuquerque, NM 87185.
*
Reference: RDEAM - An Adams ODE Code with Root Solving
Capability, Sandia Report SAND85-1595, December
1985.
*
* RDEBD - H. A. Watts, Sandia National Laboratories,
Albuquerque, NM 87185.
*
Reference: Backward Differentiation Formulae Revisited:
Improvements in DEBDF and a New Root Solving Code
RDEBD, Sandia Report SAND85-2676, July 1986.
*
* RKF45 - H. A. Watts, Sandia National Laboratories,
Albuquerque, NM 87185.
- L. F. Shampine, Southern Methodist University,
Dallas, TX 75275.
*
Reference: Computer Methods for Mathematical Computation,
G. E. Forsythe, M. A. Malcolm and C. B. Moler,
Prentice-Hall, Englewood, Cliffs, NJ (1977).
*
* SECDER - C. A. Addison, CMI, Fantoftvegen 38, N-5036
Fantoft, Bergen, Norway.
*
Reference: Implementing a Stiff Method based upon the Second
Derivative Formulas, Dept. of Computer Science
Technical Report No. 130/79, University of Toronto,
Toronto, Ontario M5S 1A7, Canada.
*
* STRIDE - J. C. Butcher, University of Auckland, Private Bag,
Auckland, New Zealand.
- K. Burrage, University of Auckland, Private Bag,
Auckland, New Zealand.
- F. H. Chipman, Acadia University, Wolfville, NS.
B0P 1X0 Canada.
*
Reference: STRIDE - Stable Runge-Kutta integrator for
differential equations, Dept. of Mathematics.
Report Series 150, University of Auckland,
Auckland, New Zealand.
*
*
1002 yes
1003 no
3004 yes
3005 no
4006 yes
4007 no
5008 yes
5009 no
60010 yes
60011 no
70012 yes
70013 no
80999 the code references?
100014 yes
100015 no
110016 yes
110017 no
120024 yes
120025 no
130018 yes
130019 no
140043 yes
140044 no
150020 yes
150021 no
160049 yes
160050 no
170022 yes
170023 no
180032 yes
180042 no
190036 yes
190037 no
200014 yes
200026 no
210027 yes
210028 no
220016 yes
220029 no
230030 yes
230031 no
240032 yes
240033 no
250034 yes
250035 no
260027 yes
260038 no
270045 yes
270046 no
280038 yes
280039 no
290030 yes
290040 no
300051 yes
300052 no
310040 yes
310041 no
320803 IVODE pseudocode 3?
320999 the code references?
330803 IVODE pseudocode 3?
330804 IVODE driver 1 for interpolation?
330806 IVODE driver 3 for interpolation?
330999 the code references?
340801 IVODE pseudocode 1?
340804 IVODE driver 1 for interpolation?
340806 IVODE driver 3 for interpolation?
340999 the code references?
350999 the code references?
360801 IVODE pseudocode 1?
360999 the code references?
370999 the code references?
380047 yes
380048 no
400053 yes
400054 no
420802 IVODE pseudocode 2?
420999 the code references?
430803 IVODE pseudocode 3?
430805 IVODE driver 2 for interpolation?
430806 IVODE driver 3 for interpolation?
430999 the code references?
440055 yes
440056 no
450803 IVODE pseudocode 3?
450805 IVODE driver 2 for interpolation?
450806 IVODE driver 3 for interpolation?
450999 the code references?
460057 yes
460058 no
470803 IVODE pseudocode 3?
470806 IVODE driver 3 for interpolation?
470999 the code references?
480059 yes
480060 no
490803 IVODE pseudocode 3?
490805 IVODE driver 2 for interpolation?
490806 IVODE driver 3 for interpolation?
490999 the code references?
500061 yes
500062 no
510803 IVODE pseudocode 3?
510805 IVODE driver 2 for interpolation?
510806 IVODE driver 3 for interpolation?
510999 the code references?
520063 yes
520064 no
530803 IVODE pseudocode 3?
530806 IVODE driver 3 for interpolation?
530999 the code references?
540065 yes
540066 no
550801 IVODE pseudocode 1?
550805 IVODE driver 2 for interpolation?
550806 IVODE driver 3 for interpolation?
550999 the code references?
560999 the code references?
570801 IVODE pseudocode 1?
570805 IVODE driver 2 for interpolation?
570806 IVODE driver 3 for interpolation?
570999 the code references?
580999 the code references?
590801 IVODE pseudocode 1?
590806 IVODE driver 3 for interpolation?
590999 the code references?
600999 the code references?
610801 IVODE pseudocode 1?
610805 IVODE driver 2 for interpolation?
610806 IVODE driver 3 for interpolation?
610999 the code references?
620999 the code references?
630801 IVODE pseudocode 1?
630805 IVODE driver 2 for interpolation?
630806 IVODE driver 3 for interpolation?
630999 the code references?
640999 the code references?
650801 IVODE pseudocode 1?
650806 IVODE driver 3 for interpolation?
650999 the code references?
660999 the code references?
30903?
40904?
50905?
60906?
70907?
80908?
120912?
140912?
160912?
180912?
210921?
230921?
240924?
260921?
270912?
290921?
300912?
380912?
400912?
430943?
450943?
470947?
550943?
560943?
570943?
580943?
590947?
600947?
*
1 1002 2*
1 1003 3*
3* 3004 4*
3* 3005 5*
4* 4006 6*
4* 4007 7*
5* 5008 8*
5* 5009 9
6* 60010 10
6* 60011 11
7* 70012 12*
7* 70013 13
8* 80999 999*
10 100014 14*
10 100015 15*
11 110016 16*
11 110017 17
12* 120024 24*
12* 120025 25*
13 130018 18*
13 130019 19*
14* 140043 43*
14* 140044 44*
15* 150020 20*
15* 150021 21*
16* 160049 49*
16* 160050 50*
17 170022 22
17 170023 23*
18* 180032 32*
18* 180042 42*
19* 190036 36*
19* 190037 37*
20* 200014 14*
20* 200026 26*
21* 210027 27*
21* 210028 28*
22 220016 16*
22 220029 29*
23* 230030 30*
23* 230031 31*
24* 240032 32*
24* 240033 33*
25* 250034 34*
25* 250035 35*
26* 260027 27*
26* 260038 38*
27* 270045 45*
27* 270046 46*
28* 280038 38*
28* 280039 39*
29* 290030 30*
29* 290040 40*
30* 300051 51*
30* 300052 52*
31* 310040 40*
31* 310041 41*
32* 320803 803*
32* 320999 999*
33* 330803 803*
33* 330804 804*
33* 330806 806*
33* 330999 999*
34* 340801 801*
34* 340804 804*
34* 340806 806*
34* 340999 999*
35* 350999 999*
36* 360801 801*
36* 360999 999*
37* 370999 999*
38* 380047 47*
38* 380048 48*
40* 400053 53*
40* 400054 54*
42* 420802 802*
42* 420999 999*
43* 430803 803*
43* 430805 805*
43* 430806 806*
43* 430999 999*
44* 440055 55*
44* 440056 56*
45* 450803 803*
45* 450805 805*
45* 450806 806*
45* 450999 999*
46* 460057 57*
46* 460058 58*
47* 470803 803*
47* 470806 806*
47* 470999 999*
48* 480059 59*
48* 480060 60*
49* 490803 803*
49* 490805 805*
49* 490806 806*
49* 490999 999*
50* 500061 61*
50* 500062 62*
51* 510803 803*
51* 510805 805*
51* 510806 806*
51* 510999 999*
52* 520063 63*
52* 520064 64*
53* 530803 803*
53* 530806 806*
53* 530999 999*
54* 540065 65*
54* 540066 66*
55* 550801 801*
55* 550805 805*
55* 550806 806*
55* 550999 999*
56* 560999 999*
57* 570801 801*
57* 570805 805*
57* 570806 806*
57* 570999 999*
58* 580999 999*
59* 590801 801*
59* 590806 806*
59* 590999 999*
60* 600999 999*
61* 610801 801*
61* 610805 805*
61* 610806 806*
61* 610999 999*
62* 620999 999*
63* 630801 801*
63* 630805 805*
63* 630806 806*
63* 630999 999*
64* 640999 999*
65* 650801 801*
65* 650806 806*
65* 650999 999*
66* 660999 999*
3* 30903 903*
4* 40904 904*
5* 50905 905*
6* 60906 906*
7* 70907 907*
8 80908 908*
12* 120912 912*
14* 140912 912*
16* 160912 912*
18* 180912 912*
21* 210921 921*
23* 230921 921*
24* 240924 924*
26* 260921 921*
27* 270912 912*
29* 290921 921*
30* 300912 912*
38* 380912 912*
40* 400912 912*
43* 430943 943*
45* 450943 943*
47* 470947 947*
55* 550943 943*
56* 560943 943*
57* 570943 943*
58* 580943 943*
59* 590947 947*
60* 600947 947*
S 1
/
PSEUDO-CODE 1.
*
Pseudo-code to locate the values of t at which
the function g(t,y,y') is zero.
*
DRIVER
*
set initial conditions of differential equation
*
determine initial sign of g(t,y,y')
*
set options of routine to return after every attempted step
*
determine initial step h
*
repeat until t = tend or integration completed or error
condition detected
*
invoke routine attempting a step from t to t + h
returns with step = h-bar <= h taken, y, y' and
recommended next step h-star
*
if g(t,y,y') has changed sign then
*
using interpolant z(t), find s such that
g(s,z(s),z'(s)) = 0
*
if integration to be terminated then
*
set ``integration completed'' indicator
*
else
*
take care with initial setting of sign of g and value
of y'
*
set h = h-star and signal a restart
*
end if
*
else
*
h = h-star
*
end if
*
end repeat
*
end DRIVER
*
*****************************************************************
*
PSEUDO-CODE 2.
*
Pseudo-code to handle discontinuities without the
use of an associated "switching" function.
*
DRIVER
*
set initial conditions of differential equation
*
set options of routine to return after every step
*
repeat until t = tend or integration completed or error
condition detected
*
invoke routine attempting a step from t to t + h
returns with h-bar <= h taken, y, y', and recommended next
step h-star
*
if h-bar < h and h-star >= 2h-bar then
*
using interpolant z(t), search for a sufficiently
accurate bracketing (TL,TH), of the discontinuity by
sampling z'(t) - f(t,z(t))
*
if search successful then
*
set t = TH, y = z(TH)
*
if integration to be continued then
*
set the next attempted step to be h and signal a
restart
*
else
*
set ``integration completed'' indicator
*
end if
*
end if
*
end if
*
end repeat
*
end DRIVER
*
*****************************************************************
*
PSEUDO-CODE 3.
*
Pseudo-code to handle discontinuities of the
``switching'' type.
*
DRIVER
*
set initial conditions of differential equation
*
set initial status vector
*
set options of routine to return after every step
*
repeat until t = tend or integration completed or error
condition detected
*
invoke routine attempting a step from t to t + h
*
returns with h-bar <= h taken, y, y' and recommended next
step h-star
if g(t,y) has changed sign then
*
using interpolant z(t), find TL and TH such that
s is in (TL,TH) and g(s,z(s)) = 0, and (TH - TL)
is sufficiently small
*
if integration to be continued then
*
set t = TH, y = z(TH) and signal restart
*
set status vector to redefine differential equation
reset definition of g(t,y) if necessary
*
else
*
set ``integration completed'' indicator
*
end if
*
end if
*
end repeat
*
end DRIVER
*
*****************************************************************
*
Driver 1: Interpolation driver for DEABM
*
*
*
SUBROUTINE INTRPD(NEQ, OUTX, OUTY, OUTYP, RWORK, LRW,
+ IWORK, LIW)
*
* INTRPD is a driver for interpolation for the DEABM code.
* It returns approximate values of y(x) and y'(x) for x in
* the current interval.
*
* Parameters:
*
* On entry:
*
* NEQ -- the number of equations
* OUTX -- the value of the independent variable x, where
* the interpolation is required
* RWORK -- the real work area RWORK used in DEABM
* LRW -- the declared length of RWORK (in user's
* dimension)
* IWORK -- the integer work area used in DEABM
* LIW -- the declared length of IWORK (in user's
* dimension)
*
* On return:
*
* OUTY -- the interpolated values y (x), i=1,2,..,NEQ
* i
* OUTYP -- the interpolated derivative values
* y '(x), i=1,2,...,NEQ
* i
*
DIMENSION OUTY(NEQ), OUTYP(NEQ), RWORK(LRW), IWORK(LIW)
*
IYY = 2*NEQ + 22
IPHI = 3*NEQ + IYY
IPSI = 16*NEQ + 24 + IPHI
*
CALL INTRP(RWORK(13), RWORK(IYY), OUTX, OUTY, OUTYP,
+ NEQ, IWORK(28), RWORK(IPHI), RWORK(IPSI))
*
RETURN
*
END
*
*****************************************************************
*
* Driver 2: Interpolation driver for DEBDF.
*
*
SUBROUTINE INTYDD(NEQ, OUTX, OUTY, OUTYP, RWORK, LRW,
+ IFLAG)
*
* INTYDD is a driver for interpolation for the DEBDF code.
* It returns the approximate values of y(x) and y'(x) for
* values of x in the current interval.
*
* Parameters:
*
* On entry:
*
* NEQ -- the number of equations
* OUTX -- the value of the independent variable x where
* interpolation is required
* RWORK -- the real work area used in DEBDF
* LRW -- declared length of RWORK (in user's
* dimension)
*
* On return:
*
* OUTY -- the interpolated values y (x), i = 1,2,...,NEQ
* i
* OUTYP -- the interpolated derivative values
* y '(x), i=1,2,..,NEQ
* i
* IFLAG -- an error flag:
*
* on return IFLAG = 0 if no error occurred
* IFLAG = -1 if x is invalid, i.e. outside
* the current interval
*
DIMENSION OUTY(NEQ), OUTYP(NEQ), RWORK(LRW)
*
IYH = NEQ + 240
*
CALL INTYD(OUTX, 0, RWORK(IYH), NEQ, OUTY, IFLAG)
*
IF (IFLAG .EQ. 0) THEN
*
CALL INTYD(OUTX, 1, RWORK(IYH), NEQ, OUTYP, IFLAG)
*
ELSE
*
* OUTX is not valid
*
IFLAG = -1
*
END IF
*
RETURN
*
END
*
*****************************************************************
*
* Driver 3: Interpolation driver for LSODE or LSODES
*
SUBROUTINE INTDYD(NEQ, OUTX, OUTY, OUTYP, RWORK, LRW,
+ IFLAG)
*
* INTDYD is a driver for interpolation for the LSODE family
* of codes. It returns the approximate values of y(x) and
* y'(x) for values of x in the current interval.
*
* Parameters:
*
* On entry:
*
* NEQ -- the number of equations
* OUTX -- the value of the independent variable x where
* interpolation is required
* RWORK -- the real work area used in LSODE
* LRW -- declared length of RWORK (in user's
* dimension)
*
* On return:
*
* OUTY -- the interpolated values y (x), i = 1,2,...,NEQ
* i
* OUTYP -- the interpolated derivative values
* y'(x), i = 1,2,...,NEQ
* i
* IFLAG -- an error flag:
*
* on return IFLAG = 0 if no error occurred
* IFLAG = -1 if x is invalid, i.e. outside the
* current interval
*
DIMENSION OUTY(NEQ), OUTYP(NEQ), RWORK(LRW)
*
CALL INTDY(OUTX, 0, RWORK(21), NEQ, OUTY, IFLAG)
*
IF (IFLAG .EQ. 0) THEN
*
CALL INTDY(OUTX, 1, RWORK(21), NEQ, OUTYP, IFLAG)
*
ELSE
*
* OUTX is not valid
*
IFLAG = -1
*
END IF
*
RETURN
*
END
*
*****************************************************************
*
Code References
*
* HARWELL Subroutine Library, Computing Division, Harwell
Laboratory, Oxfordshire, UK.
*
HARWELL DC03
*
It is possible to obtain DC03 without obtaining the entire
library.
*
* IMSL Library, 7500 Bellaire Boulevard, Houston, TX 77036-5085.
*
IMSL DVERK, DGEAR
*
* NAG Library, The NAG Office, Wilkinson House, Jordanhill Rd.,
Oxford OX2 8DR, UK.
*
NAG D02BAF, D02BBF, D02BDF, D02BHF, D02CBF, D02CHF, D02EJF,
D02NBF, D02NCF, D02NDF, D02NGF, D02NHF, D02NJF, D02PAF, D02QAF
*
* SLATEC Library, National Energy Software Center, Argonne National
Laboratory, 9700 South Cass Avenue, IL 60439
*
SLATEC DEABM, DEBDF, DERKF
*
*
* DASSL - L. R. Petzold, Sandia National Laboratories,
Livermore, CA 94550.
*
Reference: A description of DASSL: A differential/algebraic
system solver, Proc. IMACS World Congress,
Montreal, Canada, August 1982.
*
* LSODAR - L. R. Petzold, Sandia National Laboratories,
Livermore, CA 94550.
- A. C. Hindmarsh, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
*
Reference: Automatic selection of methods for solving stiff
and non-stiff systems of ordinary differential
equations, Sandia National Laboratories Report
SAND80-8230, September 1980.
*
ODEPACK, a systemised collection of ODE solvers, in
Scientific Computing, R.S. Stepleman et al. (eds.),
North-Holland, Amsterdam, 1983 (vol 1 of IMACS
Trans. on Scientific Computation), pp. 55-64.
*
* LSODE - A. C. Hindmarsh, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
*
Reference: LSODE and LSODI: Two new initial value ordinary
differential equation solvers, ACM Signum
Newsletter, vol 15, no. 4 (1980) pp. 10-11.
*
ODEPACK, a systemised collection of ODE solvers, in
Scientific Computing, R.S. Stepleman et al. (eds.),
North-Holland, Amsterdam, 1983 (vol 1 of IMACS
Trans. on Scientific Computation), pp. 55-64.
*
* LSODES - A. C. Hindmarsh, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
- A. H. Sherman, J.S. Nolen & Associates, Houston,
Texas.
*
Reference: GEARS: A package for the solution of sparse stiff
ordinary differential equations, in A.M. Erisman et
al. (eds.), Electrical Power Problems: The
Mathematical Challenge, SIAM, Philadelphia, 1980,
pp. 190-200.
*
ODEPACK, a systemised collection of ODE solvers, in
Scientific Computing, R.S. Stepleman et al. (eds.),
North-Holland, Amsterdam, 1983 (vol 1 of IMACS
Trans. on Scientific Computation), pp. 55-64.
*
* LSODI - J. F. Painter, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
- A. C. Hindmarsh, Lawrence Livermore National
Laboratory, Livermore, CA 94550.
*
Reference: LSODE and LSODI: Two new initial value ordinary
differential equation solvers, ACM Signum
Newsletter, vol 15, no. 4 (1980) pp. 10-11.
*
ODEPACK, a systemised collection of ODE solvers, in
Scientific Computing, R.S. Stepleman et al. (eds.),
North-Holland, Amsterdam, 1983 (vol 1 of IMACS
Trans. on Scientific Computation), pp. 55-64.
*
* RDEAM - H. A. Watts, Sandia National Laboratories,
Albuquerque, NM 87185.
*
Reference: RDEAM - An Adams ODE Code with Root Solving
Capability, Sandia Report SAND85-1595, December
1985.
*
* RDEBD - H. A. Watts, Sandia National Laboratories,
Albuquerque, NM 87185.
*
Reference: Backward Differentiation Formulae Revisited:
Improvements in DEBDF and a New Root Solving Code
RDEBD, Sandia Report SAND85-2676, July 1986.
*
* RKF45 - H. A. Watts, Sandia National Laboratories,
Albuquerque, NM 87185.
- L. F. Shampine, Southern Methodist University,
Dallas, TX 75275.
*
Reference: Computer Methods for Mathematical Computation,
G. E. Forsythe, M. A. Malcolm and C. B. Moler,
Prentice-Hall, Englewood, Cliffs, NJ (1977).
*
* SECDER - C. A. Addison, CMI, Fantoftvegen 38, N-5036
Fantoft, Bergen, Norway.
*
Reference: Implementing a Stiff Method based upon the Second
Derivative Formulas, Dept. of Computer Science
Technical Report No. 130/79, University of Toronto,
Toronto, Ontario M5S 1A7, Canada.
*
* STRIDE - J. C. Butcher, University of Auckland, Private Bag,
Auckland, New Zealand.
- K. Burrage, University of Auckland, Private Bag,
Auckland, New Zealand.
- F. H. Chipman, Acadia University, Wolfville, NS.
B0P 1X0 Canada.
*
Reference: STRIDE - Stable Runge-Kutta integrator for
differential equations, Dept. of Mathematics.
*