Ordinary Differential Equations (ODEs) Solvers
Ordinary differential equations (ODEs) can be solved for Phoenix models by several methods. These methods are available under the option Max ODE in the Run Options tab.
Some models use ODE systems expressed by 
, where
is a time-dependent column vector
is a (Jacobian) matrix with no time-dependent elements, but may contain the infusion rate information, etc.
For such ODEs, there is no need to use numerical ODE solvers that contain time-step errors (e.g., DVERK, LSODA, LSODE, DOPRI5) as there is a simple solution of 
, where (0) is the initial value of at t = 0 and is known.
To get the solution 
, it is important to calculate 
, which is a matrix. The Matrix Exponent ODE solver is used to calculate , and can only be used to calculate ODEs where the solution is . This solver does not need to repeatedly evaluate the ODEs and, therefore, is faster than the other numerical ODE solvers available in Phoenix. Furthermore, it is not influenced by the stiffness of the ODEs. Matrix Exponent uses a general method called “scaling and squaring”[1] to calculate . It uses the property:
= 
where
s is the scaling parameter and a non-negative integer
is estimated by Padé approximation
The 
inside () is “scaling” 
, and the superscript 
on is “squaring” ().
Matrix exponential ODE solvers
Numerical ODE solvers
A correctly rounded mathematical library[2] is used for some elementary functions implementation to reduce the differences in results between platforms, CPUs, and operating systems.
The ODE solvers used by Phoenix have three accuracy controls that are available in the Advanced options section of the Run Options tab:
ODE Rel. Tol.: Relative tolerance (RTOL)
ODE Abs. Tol.: Absolute tolerance (ATOL)
ODE max step: Maximum number of allowable steps or function evaluations (depending on the solver) to achieve the above accuracies (MAXSTEP)
For more information, see the “RTOL, ATOL, and MAXSTEP ODE error controls” section.
This is the default option for solving ordinary differential equations for Phoenix models and is based on the DGPADM subroutine in the Expokit package developed by Roger B. Sidje [3]. As its name implies, the Matrix Exponent option uses the matrix exponential for solving ODEs. Hence, it requires that the system can be expressible as , where is a matrix (the Jacobian) with elements that are constant over a time interval. It uses “scaling and squaring” method combined with Padè approximation to calculate . It is faster and more accurate than the numerical ODE solvers and does not suffer from stiffness of the ODEs. Additional information on using the matrix exponential for solving ODEs is available in the “Matrix Exponent” section in the PML Reference Guide.
For the case where the Matrix Exponent option cannot be used for the given model, Phoenix automatically switches the ODE solver to the non-stiff solver DVERK except in the case where the model involves gammaDelay function (for this case, it automatically switches to the non-stiff solver DOPRI5). Note that the actual ODE solver requested/used can be found at the beginning of the Core Output tab (Text Output Results section).
LSODE (Livermore Solver for Ordinary Differential Equations)[4] is the ODE solver used by Phoenix to solve stiff systems of the form dy/dt=f. It treats the Jacobian matrix df/dy as either a dense (full) or a banded matrix, and as either supplied automatically by Phoenix or internally approximated by difference quotients (if Phoenix fails to generate it).
Stiff systems typically arise from models with processes that proceed at different rates, for example, slow equilibration of a peripheral compartment with a central compartment relative to the rate at which elimination occurs from the central compartment. Generally, a stiff solver outperforms, in terms of both time and accuracy, a non-stiff solver for such systems but the non-stiff solver can be superior for non-stiff systems.
Nonlinear systems can change from stiff to non-stiff during evaluation, which can result in incorrect estimates and/or much longer execution times if the incorrect ODE solver is selected. LSODA [4] is a variant version of the LSODE solver and it automatically selects between non-stiff and stiff methods. It uses Adams methods (predictor corrector) in the non-stiff case, and Backward Differentiation Formula (BDF) methods (the Gear methods) in the stiff case. LSODA is a robust adaptive stepwise solver that uses the non-stiff method initially, and dynamically monitors data to decide which method to use.
Based on the DVERK solver[5], this is an explicit Runge-Kutta method founded on Verner’s fifth and sixth order pair of formulas. Due to its fast and robust performance, it is one of the widely used methods for solving non-stiff problems and is also the default solver used by Phoenix when the matrix exponent solver is not applicable.
It is well-known that the explicit Runge-Kutta methods cannot compete with specifically designed stiff methods, except for very mildly stiff systems. For the stiff problems, integration errors are expected. Switch to the LSODE/LSODA solver.
DOPRI5 is another widely used solver for non-stiff problems. Based on the DOPRI5 ODE solver by E. Hairer and G. Wanner[6], it is an explicit Runge-Kutta method of order (4)5, based on Dormand-Prince pair. More specifically, it uses six function evaluations to calculate fourth- and fifth-order accurate solutions. Numerical experiments show that DOPRI5 may be more reliable and efficient than DVERK, in certain stiff problems.