65 Differential Equation Solvers
1 $SUBROUTINES: New Differential Equation Solving Method: LSODA (ADVAN13)
As of NM7, a differential equation solver has been introduced, called LSODA, and is accessed using ADVAN=13 or ADVAN13. This routine is useful for stiff and non-stiff equations. This is similar to the LSODI routine used by ADVAN9, except that ADVAN13 can at times execute more quickly than ADVAN9. The ADVAN13 differential equation solver has been shown to solve problems more quickly with the new estimation methods, whereas for classical NONMEM methods, selecting ADVAN6 or 9 may still be of greater advantage.
Example:
$SUBROUTINES ADVAN13 TRANS1 TOL=5Where TOL is the number of digits accuracy desired to integrate the differential equations (accuracy to within \(10^{-\text{TOL}}\)). The code to the differential equation solver is found in ..\source\LSODA.f90. On occasion, coded errors will be displayed if the algorithm is having trouble integrating the equations. These errors may usually be ignored, unless the error shows up frequently, and ultimately results in failure for the problem to complete. Typically the remedy is to increase or decrease TOL, but for those who desire to understand what the error codes mean, there are well documented comments on these at the beginning of LSODA.f90. They are printed here for convenience:
! ISTATE=An index used for input and output to specify the the state of the calculation.
!
! On input,the values of istate are as follows.
! 1 Means this is the first call for the problem (initializations will be done).
! See note below.
! 2 Means this is not the first call,and the calculation is to continue
! normally, with no change in any input parameters except possibly TOUT
! and ITASK. (If ITOL,RTOL,and/or ATOL are changed between calls with
! ISTATE=2,the new values will be used but not tested for legality.)
! 3 Means this is not the first call,and the calculation is to continue
! normally,but with a change in input parameters other than TOUT and ITASK.
! changes are allowed in NEQ,ITOL,RTOL,ATOL,IOPT,LRW,LIW,JT,ML,MU and any
! optional inputs except H0,MXORDN,AND MXORDS.
! (see IWORK description for ML and MU.)
! Note: A preliminary call with TOUT=T is not counted as a first call here,as
! no initialization or checking of input is done. (Such a call is sometimes
! useful for the purpose of outputting the initial conditions.) Thus the first
! call for which TOUT /= T requires ISTATE=1 on input.
!
! On output,istate has the following values and meanings.
! 1 Means nothing was done; TOUT=T and ISTATE=1 on input.
! 2 Means the integration was performed successfully.
! -1 Means an excessive amount of work (more than MXSTEP steps) was done on
! this call,before completing the requested task,but the integration was
! otherwise successful as far as T. (MXSTEP is an optional input and is
! normally 500.) TO continue,the user may simply reset ISTATE to a value > 1
! and call again (the excess work step counter will be reset to 0).
! In addition,the user may increase MXSTEP to avoid this error return
! (see below on optional inputs).
! -2 Means too much accuracy was requested for the precision of the machine
! being used. This was detected before completing the requested task,but
! the integration was successful as far as T. To continue,the tolerance
! parameters must be reset,and ISTATE must be set to 3. The optional output
! TOLSF may be used for this purpose. (Note: If this condition is detected
! before taking any steps,then an illegal input return (ISTATE=-3) occurs
! instead.)
! -3 Means illegal input was detected,before taking any integration steps.
! See written message for details.
! Note: If the solver detects an infinite loop of calls to the solver with
! illegal input,it will cause the run to stop.
! -4 Means there were repeated error test failures on one attempted step,before
! completing the requested task,but the integration was successful as far as T.
! The problem may have a singularity,or the input may be inappropriate.
! -5 Means there were repeated convergence test failures on one attempted step,
! before completing the requested task,but the integration was successful as
! far as T. This may be caused by an inaccurate jacobian matrix, if one is
! being used.
! -6 Means EWT(I) became zero for some I during the integration. Pure relative
! error control (ATOL(I)=0.0) was requested on a variable which has now
! vanished. The integration was successful as far as T.
! -7 Means the length of RWORK and/or IWORK was too small to proceed,but the
! integration was successful as far as T. This happens when DLSODA chooses
! to switch methods but LRW and/or LIw is too small for the new method.
!
! Note: Since the normal output value of ISTATE is 2, it does not need to be
! reset for normal continuation. Also,since a negative input value of ISTATE
! will be regarded as illegal, a negative output value requires the user to
! change it, and possibly other inputs,before calling the solver again.1.1 ATOL (NM72)
A $EST option when using ADVAN13 is the absolute tolerance. The ATOL for ADVAN13 by default is 12 (that is, precision is \(10^{-12}\)). Usually the problem runs quickly when using ADVAN13 with this setting. On occasion, however, you may want to reduce ATOL (usually set it equal to that of TOL), and improve speeds of up to 3 to 4 fold. ATOL may be set at the $EST or $COV command. The absolute tolerance is set to the same ATOL for all compartments.
As of NM73, ATOL also acts on ADVAN9’s differential equation solver, where by default absolute significant digits accuracy (absolute tolerance) is 12. As of NM74, for ADVAN9, ATOL=99 specifies that variable, calculated ATOL values are to be derived, in accordance with an algorithm that is present in ADVAN9. Do not use ATOL=99 for ADVAN13, 14, and 15, as it does not work properly for these algorithms.
The relative tolerance is still set by TOL by the $SUBROUTINES, $COV, or $TOL record. As of NM74, ATOL acts on ADVAN14, ADVAN15, ADVAN16, ADVAN17, and ADVAN18 as well.
1.2 MXSTEP (NM73)
Additional control may be obtained by setting the maximum number of integration steps (default is 10000 for ADVAN13, 14, 15, 16, and 17 and 2147483647 for ADVAN9):
$PK
MXSTEP=5000Suitable for ADVAN9, ADVAN13, ADVAN14, ADVAN15, ADVAN16, ADVAN17, and ADVAN18 maximum integration steps can be controlled by this variable.
2 $SUBROUTINES TOL, ATOL, SSTOL, and SSATOL: Additional control of relative and absolute tolerance (NM74)
Additional TOL type options may be set at $SUBROUTINES, which will allow these settings to be in effect throughout the $PROB, including during simulation.
2.1 ATOL=n
The absolute tolerance can be set at the $SUBROUTINES record, similar to TOL, the relative tolerance. ATOL is used for ADVAN9, ADVAN13, ADVAN14, ADVAN15, ADVAN16, ADVAN17, and ADVAN18. If not set, ATOL is set to 1.0E-12 by default. Note an ATOL set at the $EST record will override that set at the $SUB record, during the estimation. An ATOL set at the $COV record will override that set at the $SUBR record. Similarly, a relative TOL set at the $COV record will override a TOL set at the $SUBR record.
2.2 SSTOL=n
The relative tolerance for steady state evaluations can be set at the $SUBROUTINES record. If not specified, the SSTOL is set to TOL.
2.3 SSATOL=n
The absolute tolerance for steady state evaluations can be set at the $SUBROUTINES record. If not specified, the SSATOL is set to ATOL.
2.4 TOLC=n
The relative tolerance for the FOCE/LAPLACE $COV step can be set at the $SUBROUTINES record. TOLC is used for ADVAN9, ADVAN13, ADVAN14, ADVAN15, ADVAN16, ADVAN17, and ADVAN18. If not set TOLC defaults to the value of TOL. Note a TOL set at the $COV record will override TOLC set at the $SUB record, during the estimation.
2.5 ATOLC=n
The absolute tolerance for the FOCE/LAPLACE $COV step can be set at the $SUBROUTINES record. ATOLC is used for ADVAN9, ADVAN13, ADVAN14, ADVAN15, ADVAN16, ADVAN17, and ADVAN18. If not set, ATOLC is set to ATOL. Note an ATOL set at the $COV record will override the ATOLC set at the $SUB record.
2.6 SSTOLC=n
The relative tolerance for steady state evaluations during the FOCE/Laplace $COV step can be set at the $SUBROUTINES record. If not specified, the SSTOLC is set to SSTOL.
2.7 SSATOLC=n
The absolute tolerance for steady state evaluations during the FOCE/Laplace $COV step can be set at the $SUBROUTINES record. If not specified, the SSATOLC is set to SSATOL.
2.8 $TOL
An alternative method to setting TOL, ATOL, SSTOL, and SSATOL, TOLC, ATOLC, SSTOLC, and SSATOLC at the $SUBR record is by using the $TOL record, and specific ones for each compartment may be entered. For example, following the $TOL record, the following may be entered, one on each line:
NRD(0) = 4 ; Set SSTOL
NRD(1) = 5 ; Set TOL for compartment 1
NRD(2) = 6 ; Set TOL for compartment 2
...
ANRD(0) = 4 ; Set SSATOL
ANRD(1) = 7 ; Set ATOL for compartment 1
ANRD(2) = 6 ; Set ATOL for compartment 2
...
NRDC(0) = 8 ; Set SSTOLC
NRDC(1) = 9 ; Set TOLC for compartment 1
NRDC(2) = 10 ; Set TOLC for compartment 2
...
ANRDC(0) = 3 ; Set SSATOLC
ANRDC(1) = 4 ; Set ATOLC for compartment 1
ANRDC(2) = 5 ; Set ATOLC for compartment 2Note that if NRD(0) and ANRD(0) are not defined, the values for NRD(1) and ANRD(1) are used, respectively. Similarly, if NRDC(0) and ANRDC(0) are not defined, the values for NRDC(1) and ANRDC(1) are used, respectively. For additional compartments not specified, the tolerance of the last compartment specified will be used. So in the above example, NRD(3)=NRD(2), NRD(4)=NRD(2), etc, and ANRD(3)=ANRD(2), ANRD(4)=ANRD(2), etc.
ATOL set at the $EST record and TOL and ATOL set at the $COV record will temporarily override these values set in $TOL, for all compartments (except the SS tolerances). These overrides at the $EST and $COV records are a legacy of pre-nm74 versions, and should no longer be used, now that $SUBR tolerance options and $TOL statements can allow the user to provide detailed settings to various compartment, steady state, estimation, and covariance tolerances.
As of nm75, you may specify compartment by its name defined in $MODEL. For example:
$MODEL COMP=(DEPOT,INITIALOFF,DEFDOSE) COMP=(CENTRAL,DEFOBS,NOOFF)
...
$TOL
NRD(DEPOT)=5
NRD(central)=6Finally, you may supply a TOL routine that assigns values of NRD and ANRD specifically for each NONMEM step (estimation, covariance, simulation, table/scatter step, simulation, initial parameters estimate, nonparametric). For example, create a toluser.f90 file, as
SUBROUTINE TOL(NRD,ANRD,NRDC,ANRDC)
USE NMPRD_INT, ONLY: IPROB
USE NM_BAYES_INT, ONLY: NM_STEP,BASE_STEP,EST_STEP,COV_STEP, &
TABLE_STEP,SIML_STEP,INE_STEP, NONP_STEP
INTEGER :: NRD(0:*), ANRD(0:*), NRDC(0:*), ANRDC(0:*)
IF(NM_STEP==EST_STEP) THEN
NRD(1)=6
ANRD(1)=10
ELSE IF (NM_STEP==COV_STEP) THEN
NRD(1)=7
ANRD(1)=8
ELSE IF (NM_STEP==TABLE_STEP) THEN
NRD(1)=8
ANRD(1)=7
ELSE
NRD(1)=9
ANRD(1)=12
ENDIF
IF(IPROB>1) THEN
NRD(1)=NRD(1)+1
ANRD(1)=ANRD(1)+1
ENDIF
RETURN
ENDand incorporate using $SUBR:
$SUBROUTINES ADVAN13 TRANS1 TOL=toluser.f90You can have tolerances settings defined for each NONMEM step. Notice that NRD and ANRD defined with the XXX_STEP condition takes the place of having to define NRDC and ANRDC. You may even define tolerances for specific problems, etc., as shown above. More of the USE declarations for relevant parameters are listed in ..\util\nonmem_reserved_general. Compartment names defined in $MODEL are not available in user-defined TOL subroutines.
3 $SUBROUTINES: Yet Another New Differential Equation Solving Method: CVODES (ADVAN14) (NM74)
The CVODES ordinary differential equation (ODE) solver system allows advanced control of how the ODE’s are solved. This system is for users with large sets of ODE’s to solve, and who have some understanding of ODE solving methods. The CVODES is a subset of the Sundials system from Lawrence Livermore National Laboratory, by Alan Hindmarch and Radu Serban, and is a descendent of the LSODA (ADVAN13) system. For many problems, the LSODA (ADVAN13) method works just as efficiently as CVODES (ADVAN14), but there may be cases of very complicated models where the user may benefit from using ADVAN14.
ADVAN14 can be used with its default settings, by simply setting
$SUBROUTINES ADVAN14 ...If the user wishes to change the settings of the CVODES system, he should make a copy of ..\pr\CVODEU.f90 to the run directory (or other location), modify its contents, preferably rename it, such as CVODEU2.f90, so you know it is modified, then the $SUBROUTINES statement in the control stream file should include the OTHER statement:
$SUBROUTINES ADVAN14 .. OTHER=CVODEU2.f90Include a path if your modified CVODEU file is not in your run path, e.g.:
$SUBROUTINES ADVAN14 .. OTHER=\my_favorate_settings\CVODEU2.f90While there are some comments in the CVODEU.f90 routine to identify the options, a thorough understanding requires consulting the SUNDIALS manual ..\guides\cvs_guide.pdf. Furthermore, the user may develop their own Jacobian evaluation and other routines, etc. Using ADVAN14 in a customized manner is not for the faint-hearted.
Options ATOL and MXSTEP are available to ADVAN14, as with ADVAN13.
ADVAN14 also has a root-finding algorithm. See the comments in CVODEU.f90 for its use. See also ..\examples\cmax14.ctl and ..\examples\cmax14u.f90 for using the root-finder to the find the exact tmax and cmax in an absorption-central compartment model.
For the c routines to have been compiled, the proper c compiler (variable ccc) and compiler options (variable ccop) needed to have been set in the SETUP74 script at the time of building NONMEM. If this was not done at the time of building, you can modify the compiler and options in the appropriate cvode_build* script located in the ..\pr directory, and execute it from a terminal window.
4 $SUBROUTINES: Yet Another New Differential Equation Solving Method: IDAS (ADVAN15) (NM74)
The IDAS differential-algebraic equation (DAE) solver system allows advanced control of how the ODE’s with equilibrium systems (algebraic equations) are solved. This system is for users with large sets of ODE’s to solve, and who have some understanding of ODE solving methods. The IDAS is a subset of the Sundials system from Lawrence Livermore National Laboratory, by Alan Hindmarch and Radu Serban, and is a descendent of the LSODI1 (ADVAN9) system. For the few tests performed, the IDA system can work more efficiently than LSODI1 (ADVAN9) depending on the settings of parameters in IDAU.f90 (see below). There may also be cases of very complicated models where the user may benefit from using ADVAN15.
ADVAN15 can be used with its default settings, by simply setting
$SUBROUTINES ADVAN15 ...If the user wishes to change the parameters of the IDAS system, he should make a copy of ..\pr\IDAU.f90 to the run directory (or other location), modify its contents, preferably rename it, such as IDAU2.f90, so you know it is modified, then the $SUBROUTINES statement in the control stream file should include the OTHER statement:
$SUBROUTINES ADVAN15 .. OTHER=IDAU2.f90Include a path if your modified IDAU file is not in your run path, e.g.:
$SUBROUTINES ADVAN15 .. OTHER=\my_favoriate_settings\IDAU2.f90While there are some comments in the IDAU.f90 routine to identify the options, a thorough understanding requires consulting the SUNDIALS manual ..\guides\idas_guide.pdf. Furthermore, the user may develop their own Jacobian evaluation and other routines, etc. Using ADVAN15 in a customized manner is not for the faint-hearted.
Options ATOL and MXSTEP are available to ADVAN15, as with ADVAN13.
ADVAN15 also has a root-finding algorithm. See the comments in IDAU.f90 for its use.
For the c routines to have been compiled, the proper c compiler (variable ccc) and compiler options (variable ccop) needed to have been set in the SETUP74 script at the time of building NONMEM. If this was not done at the time of building, you can modify the compiler and options in the appropriate ida_build* script located in the ..\pr directory, and execute it from a terminal window.
5 ITASK_, STOP_TIME: Avoiding overshoot in ADVAN9, ADVAN13, ADVAN14, and ADVAN15
These LSODA based routines use an algorithm for integration that overshoots the integration interval during calls to DES, but still accurately evaluates at the end of the integration interval when all calculations are completed. However, you may wish to capture a maximal or minimal value during $DES, and the overshoot should be turned off for this purpose. This is readily done by setting ITASK_=4 in $PK:
$PK
ITASK_=4You may also specify a STOP_TIME (Tcrit) past which it should not integrate, if it is different from the end of the normal integration interval:
IF(TIME==4.0) STOP_TIME=5.0To set back to default (end of normal integration interval),
STOP_TIME=-1.0d+300ITASK_=3 will also work for ADVAN14 and ADVAN15. Do not use ITASK_=2 for any of the routines, or ITASK_=3 for ADVAN9 or ADVAN13, as inaccurate evaluations of the state variables will result. You can read about ITASK by searching for the word in ..\pr\LSODI1.f90 (ADVAN9), ..\pr\LSODA.f90 (ADVAN13), (..\pr\CVODE.f90, ..\pr\cvode\fcvodes.c, ..\pr\cvode\cvodes.c, ..\guides\cv*.pdf) = ADVAN14, and (..\pr\IDA.f90, ..\pr\cvode\fidas.c, ..\pr\cvode\idas.c, ..\guides\ida*.pdf) = ADVAN15.
6 HMIN_, HINIT_, HMAX_: Avoiding missed short events in ADVAN9, ADVAN13, ADVAN14, and ADVAN15 (nm760)
These reserved variables are used to modify the integration step size behavior in ADVAN9, ADVAN13, ADVAN14, and ADVAN15. When a modeled short, large rate input occurs at some time during the simulation, that is not a scheduled bolus event or is not easily specified with MTIME(), the LSODA type ODE solvers can miss this event. For example, a gamma function shaped input with a large N>20 parameter, can create a stiff, pseudo-bolus input that may be missed by one of these ADVAN’s. One way to avoid missing this input is to set the integration step size to an upper limit:
$PK
HMAX_=0.01One may also set the initial step size (HINIT_), or minimal step size (HMIN_) within the $PK block in a similar manner. To revert to default action, set HINIT_, HMIN_, or HMAX_ to 0.
The usual technique for handling short-term rate inputs is to set the RTOL and ATOL to very high levels, such as between 13 and 15, but this may not work for every problem, and one may need to resort to explicitly limiting the HMAX_ or HINIT_ in the manner shown above.
7 $SUBROUTINES: Yet Another New Differential Equation Solving Method: Delay Differential Equations (ADVAN16, ADVAN17, ADVAN18) (NM75)
For use of these delay differential equation (DDE) solvers, please see Using the Delay differential equation Solvers with the doexpand program for Discrete Delay Problems (nm75) in I.87 ddexpand Utility Program for Modeling Discrete Time Delays by method of Steps (NM74).
Delay Differential Equation solvers are available as ADVAN16 and ADVAN18 (no Equilibrium Compartments), and ADVAN17 (with Equilibrium Compartments).