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. *