# zvode: Solver for Ordinary Differential Equations (ODE) for COMPLEX... In deSolve: Solvers for Initial Value Problems of Differential Equations ('ODE', 'DAE', 'DDE')

## Description

Solves the initial value problem for stiff or nonstiff systems of ordinary differential equations (ODE) in the form:

dy/dt = f(t,y)

where dy and y are complex variables.

The R function `zvode` provides an interface to the FORTRAN ODE solver of the same name, written by Peter N. Brown, Alan C. Hindmarsh and George D. Byrne.

## Usage

 ```1 2 3 4 5 6 7``` ```zvode(y, times, func, parms, rtol = 1e-6, atol = 1e-6, jacfunc = NULL, jactype = "fullint", mf = NULL, verbose = FALSE, tcrit = NULL, hmin = 0, hmax = NULL, hini = 0, ynames = TRUE, maxord = NULL, bandup = NULL, banddown = NULL, maxsteps = 5000, dllname = NULL, initfunc = dllname, initpar = parms, rpar = NULL, ipar = NULL, nout = 0, outnames = NULL, forcings = NULL, initforc = NULL, fcontrol = NULL, ...) ```

## Arguments

 `y ` the initial (state) values for the ODE system. If `y` has a name attribute, the names will be used to label the output matrix. y has to be complex `times ` time sequence for which output is wanted; the first value of `times` must be the initial time; if only one step is to be taken; set `times = NULL`. `func ` either an R-function that computes the values of the derivatives in the ODE system (the model definition) at time `t`, or a character string giving the name of a compiled function in a dynamically loaded shared library. If `func` is an R-function, it must be defined as: `func <- function(t, y, parms, ...)`. `t` is the current time point in the integration, `y` is the current estimate of the variables in the ODE system. If the initial values `y` has a `names` attribute, the names will be available inside `func`. `parms` is a vector or list of parameters; ... (optional) are any other arguments passed to the function. The return value of `func` should be a list, whose first element is a vector containing the derivatives of `y` with respect to `time`, and whose next elements are global values that are required at each point in `times`. The derivatives must be specified in the same order as the state variables `y`. They should be complex numbers. If `func` is a string, then `dllname` must give the name of the shared library (without extension) which must be loaded before `zvode()` is called. See package vignette `"compiledCode"` for more details. `parms ` vector or list of parameters used in `func` or `jacfunc`. `rtol ` relative error tolerance, either a scalar or an array as long as `y`. See details. `atol ` absolute error tolerance, either a scalar or an array as long as `y`. See details. `jacfunc ` if not `NULL`, an R function that computes the Jacobian of the system of differential equations dydot(i)/dy(j), or a string giving the name of a function or subroutine in ‘dllname’ that computes the Jacobian (see vignette `"compiledCode"` for more about this option). In some circumstances, supplying `jacfunc` can speed up the computations, if the system is stiff. The R calling sequence for `jacfunc` is identical to that of `func`. If the Jacobian is a full matrix, `jacfunc` should return a matrix dydot/dy, where the ith row contains the derivative of dy_i/dt with respect to y_j, or a vector containing the matrix elements by columns (the way R and FORTRAN store matrices). Its elements should be complex numbers. If the Jacobian is banded, `jacfunc` should return a matrix containing only the nonzero bands of the Jacobian, rotated row-wise. See first example of `lsode`. `jactype ` the structure of the Jacobian, one of `"fullint"`, `"fullusr"`, `"bandusr"` or `"bandint"` - either full or banded and estimated internally or by user; overruled if `mf` is not `NULL`. `mf ` the "method flag" passed to function `zvode` - overrules `jactype` - provides more options than `jactype` - see details. `verbose ` if TRUE: full output to the screen, e.g. will print the `diagnostiscs` of the integration - see details. `tcrit ` if not `NULL`, then `zvode` cannot integrate past `tcrit`. The FORTRAN routine `dvode` overshoots its targets (times points in the vector `times`), and interpolates values for the desired time points. If there is a time beyond which integration should not proceed (perhaps because of a singularity), that should be provided in `tcrit`. `hmin ` an optional minimum value of the integration stepsize. In special situations this parameter may speed up computations with the cost of precision. Don't use hmin if you don't know why! `hmax ` an optional maximum value of the integration stepsize. If not specified, hmax is set to the largest difference in `times`, to avoid that the simulation possibly ignores short-term events. If 0, no maximal size is specified. `hini ` initial step size to be attempted; if 0, the initial step size is determined by the solver. `ynames ` logical; if `FALSE`: names of state variables are not passed to function `func` ; this may speed up the simulation especially for multi-D models. `maxord ` the maximum order to be allowed. `NULL` uses the default, i.e. order 12 if implicit Adams method (`meth = 1`), order 5 if BDF method (`meth = 2`). Reduce maxord to save storage space. `bandup ` number of non-zero bands above the diagonal, in case the Jacobian is banded. `banddown ` number of non-zero bands below the diagonal, in case the Jacobian is banded. `maxsteps ` maximal number of steps per output interval taken by the solver. `dllname ` a string giving the name of the shared library (without extension) that contains all the compiled function or subroutine definitions refered to in `func` and `jacfunc`. See package vignette `"compiledCode"`. `initfunc ` if not `NULL`, the name of the initialisation function (which initialises values of parameters), as provided in ‘dllname’. See package vignette `"compiledCode"`. `initpar ` only when ‘dllname’ is specified and an initialisation function `initfunc` is in the dll: the parameters passed to the initialiser, to initialise the common blocks (FORTRAN) or global variables (C, C++). `rpar ` only when ‘dllname’ is specified: a vector with double precision values passed to the DLL-functions whose names are specified by `func` and `jacfunc`. `ipar ` only when ‘dllname’ is specified: a vector with integer values passed to the dll-functions whose names are specified by `func` and `jacfunc`. `nout ` only used if `dllname` is specified and the model is defined in compiled code: the number of output variables calculated in the compiled function `func`, present in the shared library. Note: it is not automatically checked whether this is indeed the number of output variables calculated in the DLL - you have to perform this check in the code - See package vignette `"compiledCode"`. `outnames ` only used if ‘dllname’ is specified and `nout` > 0: the names of output variables calculated in the compiled function `func`, present in the shared library. These names will be used to label the output matrix. `forcings ` only used if ‘dllname’ is specified: a list with the forcing function data sets, each present as a two-columned matrix, with (time,value); interpolation outside the interval [min(`times`), max(`times`)] is done by taking the value at the closest data extreme. See forcings or package vignette `"compiledCode"`. `initforc ` if not `NULL`, the name of the forcing function initialisation function, as provided in ‘dllname’. It MUST be present if `forcings` has been given a value. See forcings or package vignette `"compiledCode"`. `fcontrol ` A list of control parameters for the forcing functions. forcings or package vignette `"compiledCode"` `... ` additional arguments passed to `func` and `jacfunc` allowing this to be a generic function.

## Details

see `vode`, the double precision version, for details.

## Value

A matrix of class `deSolve` with up to as many rows as elements in `times` and as many columns as elements in `y` plus the number of "global" values returned in the next elements of the return from `func`, plus and additional column for the time value. There will be a row for each element in `times` unless the FORTRAN routine ‘zvode’ returns with an unrecoverable error. If `y` has a names attribute, it will be used to label the columns of the output value.

## Note

From version 1.10.4, the default of atol was changed from 1e-8 to 1e-6, to be consistent with the other solvers.

The following text is adapted from the zvode.f source code:

When using `zvode` for a stiff system, it should only be used for the case in which the function f is analytic, that is, when each f(i) is an analytic function of each y(j). Analyticity means that the partial derivative df(i)/dy(j) is a unique complex number, and this fact is critical in the way `zvode` solves the dense or banded linear systems that arise in the stiff case. For a complex stiff ODE system in which f is not analytic, `zvode` is likely to have convergence failures, and for this problem one should instead use `ode` on the equivalent real system (in the real and imaginary parts of y).

## Author(s)

Karline Soetaert <[email protected]>

## References

P. N. Brown, G. D. Byrne, and A. C. Hindmarsh, 1989. VODE: A Variable Coefficient ODE Solver, SIAM J. Sci. Stat. Comput., 10, pp. 1038-1051.
Also, LLNL Report UCRL-98412, June 1988.

G. D. Byrne and A. C. Hindmarsh, 1975. A Polyalgorithm for the Numerical Solution of Ordinary Differential Equations. ACM Trans. Math. Software, 1, pp. 71-96.

A. C. Hindmarsh and G. D. Byrne, 1977. EPISODE: An Effective Package for the Integration of Systems of Ordinary Differential Equations. LLNL Report UCID-30112, Rev. 1.

G. D. Byrne and A. C. Hindmarsh, 1976. EPISODEB: An Experimental Package for the Integration of Systems of Ordinary Differential Equations with Banded Jacobians. LLNL Report UCID-30132, April 1976.

A. C. Hindmarsh, 1983. ODEPACK, a Systematized Collection of ODE Solvers. in Scientific Computing, R. S. Stepleman et al., eds., North-Holland, Amsterdam, pp. 55-64.

K. R. Jackson and R. Sacks-Davis, 1980. An Alternative Implementation of Variable Step-Size Multistep Formulas for Stiff ODEs. ACM Trans. Math. Software, 6, pp. 295-318.

Netlib: http://www.netlib.org

`vode` for the double precision version
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59``` ```## ======================================================================= ## Example 1 - very simple example ## df/dt = 1i*f, where 1i is the imaginary unit ## The initial value is f(0) = 1 = 1+0i ## ======================================================================= ZODE <- function(Time, f, Pars) { df <- 1i*f return(list(df)) } pars <- NULL yini <- c(f = 1+0i) times <- seq(0, 2*pi, length = 100) out <- zvode(func = ZODE, y = yini, parms = pars, times = times, atol = 1e-10, rtol = 1e-10) # The analytical solution to this ODE is the exp-function: # f(t) = exp(1i*t) # = cos(t)+1i*sin(t) (due to Euler's equation) analytical.solution <- exp(1i * times) ## compare numerical and analytical solution tail(cbind(out[,2], analytical.solution)) ## ======================================================================= ## Example 2 - example in "zvode.f", ## df/dt = 1i*f (same as above ODE) ## dg/dt = -1i*g*g*f (an additional ODE depending on f) ## ## Initial values are ## g(0) = 1/2.1 and ## z(0) = 1 ## ======================================================================= ZODE2<-function(Time,State,Pars) { with(as.list(State), { df <- 1i * f dg <- -1i * g*g * f return(list(c(df, dg))) }) } yini <- c(f = 1 + 0i, g = 1/2.1 + 0i) times <- seq(0, 2*pi, length = 100) out <- zvode(func = ZODE2, y = yini, parms = NULL, times = times, atol = 1e-10, rtol = 1e-10) ## The analytical solution is ## f(t) = exp(1i*t) (same as above) ## g(t) = 1/(f(t) + 1.1) analytical <- cbind(f = exp(1i * times), g = 1/(exp(1i * times) + 1.1)) ## compare numerical solution and the two analytical ones: tail(cbind(out[,2], analytical[,1])) ```