# stodes: Steady-state solver for ordinary differential equations (ODE)... In rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations

## Description

Estimates the steady-state condition for a system of ordinary differential equations (ODE) in the form:

dy/dt = f(t,y)

and where the jacobian matrix df/dy has an arbitrary sparse structure.

Uses a newton-raphson method, implemented in Fortran.

The system of ODE's is written as an R function or defined in compiled code that has been dynamically loaded.

## Usage

 1 2 3 4 5 6 7 8 stodes(y, time = 0, func, parms = NULL, rtol = 1e-6, atol = 1e-8, ctol = 1e-8, sparsetype = "sparseint", verbose = FALSE, nnz = NULL, inz = NULL, lrw = NULL, ngp = NULL, positive = FALSE, maxiter = 100, ynames = TRUE, dllname = NULL, initfunc = dllname, initpar = parms, rpar = NULL, ipar = NULL, nout = 0, outnames = NULL, forcings = NULL, initforc = NULL, fcontrol = NULL, spmethod = "yale", control = NULL, times = time, ...)

## Details

The work is done by a Fortran 77 routine that implements the Newton-Raphson method.

stodes is to be used for problems, where the Jacobian has a sparse structure.

There are several choices for the sparsity specification, selected by argument sparsetype.

• sparsetype = "sparseint". The sparsity is estimated by the solver, based on numerical differences. In this case, it is advisable to provide an estimate of the number of non-zero elements in the Jacobian (nnz). This value can be approximate; upon return the number of nonzero elements actually required will be known (1st element of attribute dims). In this case, inz need not be specified.

• sparsetype = "sparseusr". The sparsity is determined by the user. In this case, inz should be a matrix, containing indices (row, column) to the nonzero elements in the Jacobian matrix. The number of nonzeros nnz will be set equal to the number of rows in inz.

• sparsetype = "sparsejan". The sparsity is also determined by the user. In this case, inz should be a vector, containting the ian and jan elements of the sparse storage format, as used in the sparse solver. Elements of ian should be the first n+1 elements of this vector, and contain the starting locations in jan of columns 1.. n. jan contains the row indices of the nonzero locations of the jacobian, reading in columnwise order. The number of nonzeros nnz will be set equal to the length of inz - (n+1).

• sparsetype = "1D", "2D", "3D". The sparsity is estimated by the solver, based on numerical differences. Assumes finite differences in a 1D, 2D or 3D regular grid - used by functions ode.1D, ode.2D, ode.3D. Similar are "2Dmap", and "3Dmap", which also include a mapping variable (passed in nnz).

The Jacobian itself is always generated by the solver (i.e. there is no provision to provide an analytic Jacobian).

This is done by perturbing simulataneously a combination of state variables that do not affect each other.

This significantly reduces computing time. The number of groups with independent state variables can be given by ngp.

The input parameters rtol, atol and ctol determine the error control performed by the solver. See help for stode for details.

Models may be defined in compiled C or Fortran code, as well as in R. See package vignette for details on how to write models in compiled code.

When the spmethod equals ilut or ilutp, a number of parameters can be specified in argument control. They are:

fillin, the fill-in parameter. Each row of L and each row of U will have a maximum of lfil elements (excluding the diagonal element). lfil must be >= 0.

droptol, sets the threshold for dropping small terms in the factorization.

When ilutp is chosen the following arguments can also be specified:

permtol = tolerance ratio used to determne whether or not to permute two columns. At step i columns i and j are permuted when abs(a(i,j))*permtol .gt. abs(a(i,i)) [0 –> never permute; good values 0.1 to 0.01]

lenplufac = sets the working array - increase its value if a warning.

## Value

A list containing

 y a vector with the state variable values from the last iteration during estimation of steady-state condition of the system of equations. If y has a names attribute, it will be used to label the output values. ... the number of "global" values returned.

The output will have the attribute steady, which returns TRUE, if steady-state has been reached and the attribute precis with an estimate of the precision attained during each iteration, the mean absolute rate of change (sum(abs(dy))/n).

## Author(s)

Karline Soetaert <[email protected]>

## References

For a description of the Newton-Raphson method, e.g.

Press, WH, Teukolsky, SA, Vetterling, WT, Flannery, BP, 1996. Numerical Recipes in FORTRAN. The Art of Scientific computing. 2nd edition. Cambridge University Press.

When spmethod = "yale" then the algorithm uses linear algebra routines from the Yale sparse matrix package:

Eisenstat, S.C., Gursky, M.C., Schultz, M.H., Sherman, A.H., 1982. Yale Sparse Matrix Package. i. The symmetric codes. Int. J. Num. meth. Eng. 18, 1145-1151.

else the functions ilut and ilutp from sparsekit package are used:

Yousef Saad, 1994. SPARSKIT: a basic tool kit for sparse matrix computations. VERSION 2

Yousef Saad, 2003. Iterative methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics.

steady.band, to find the steady-state of ODE models with a banded Jacobian