# lsode: Solver for Ordinary Differential Equations (ODE) 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)

.

The R function lsode provides an interface to the FORTRAN ODE solver of the same name, written by Alan C. Hindmarsh and Andrew H. Sherman.

It combines parts of the code lsodar and can thus find the root of at least one of a set of constraint functions g(i) of the independent and dependent variables. This can be used to stop the simulation or to trigger events, i.e. a sudden change in one of the state variables.

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

In contrast to lsoda, the user has to specify whether or not the problem is stiff and choose the appropriate solution method.

lsode is very similar to vode, but uses a fixed-step-interpolate method rather than the variable-coefficient method in vode. In addition, in vode it is possible to choose whether or not a copy of the Jacobian is saved for reuse in the corrector iteration algorithm; In lsode, a copy is not kept.

## Usage

 1 2 3 4 5 6 7 8 lsode(y, times, func, parms, rtol = 1e-6, atol = 1e-6, jacfunc = NULL, jactype = "fullint", mf = NULL, rootfunc = NULL, verbose = FALSE, nroot = 0, 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, events=NULL, lags = NULL,...)

## Details

The work is done by the FORTRAN subroutine lsode, whose documentation should be consulted for details (it is included as comments in the source file ‘src/opkdmain.f’). The implementation is based on the November, 2003 version of lsode, from Netlib.

Before using the integrator lsode, the user has to decide whether or not the problem is stiff.

If the problem is nonstiff, use method flag mf = 10, which selects a nonstiff (Adams) method, no Jacobian used.
If the problem is stiff, there are four standard choices which can be specified with jactype or mf.

The options for jactype are

jactype = "fullint"

a full Jacobian, calculated internally by lsode, corresponds to mf = 22,

jactype = "fullusr"

a full Jacobian, specified by user function jacfunc, corresponds to mf = 21,

jactype = "bandusr"

a banded Jacobian, specified by user function jacfunc; the size of the bands specified by bandup and banddown, corresponds to mf = 24,

jactype = "bandint"

a banded Jacobian, calculated by lsode; the size of the bands specified by bandup and banddown, corresponds to mf = 25.

More options are available when specifying mf directly.
The legal values of mf are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25.
mf is a positive two-digit integer, mf = (10*METH + MITER), where

METH

indicates the basic linear multistep method: METH = 1 means the implicit Adams method. METH = 2 means the method based on backward differentiation formulas (BDF-s).

MITER

indicates the corrector iteration method: MITER = 0 means functional iteration (no Jacobian matrix is involved). MITER = 1 means chord iteration with a user-supplied full (NEQ by NEQ) Jacobian. MITER = 2 means chord iteration with an internally generated (difference quotient) full Jacobian (using NEQ extra calls to func per df/dy value). MITER = 3 means chord iteration with an internally generated diagonal Jacobian approximation (using 1 extra call to func per df/dy evaluation). MITER = 4 means chord iteration with a user-supplied banded Jacobian. MITER = 5 means chord iteration with an internally generated banded Jacobian (using ML+MU+1 extra calls to func per df/dy evaluation).

If MITER = 1 or 4, the user must supply a subroutine jacfunc.

Inspection of the example below shows how to specify both a banded and full Jacobian.

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

The diagnostics of the integration can be printed to screen by calling diagnostics. If verbose = TRUE, the diagnostics will written to the screen at the end of the integration.

See vignette("deSolve") for an explanation of each element in the vectors containing the diagnostic properties and how to directly access them.

Models may be defined in compiled C or FORTRAN code, as well as in an R-function. See package vignette "compiledCode" for details.

Examples in both C and FORTRAN are in the ‘dynload’ subdirectory of the deSolve package directory.

lsode can find the root of at least one of a set of constraint functions rootfunc of the independent and dependent variables. It then returns the solution at the root if that occurs sooner than the specified stop condition, and otherwise returns the solution according the specified stop condition.

Caution: Because of numerical errors in the function rootfun due to roundoff and integration error, lsode may return false roots, or return the same root at two or more nearly equal values of time.

## 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 ‘lsode’ returns with an unrecoverable error. If y has a names attribute, it will be used to label the columns of the output value.

## Author(s)

Karline Soetaert <[email protected]>

## References

Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE Solvers," in Scientific Computing, R. S. Stepleman, et al., Eds. (North-Holland, Amsterdam, 1983), pp. 55-64.

• rk,

• rk4 and euler for Runge-Kutta integrators.

• lsoda, lsodes, lsodar, vode, daspk for other solvers of the Livermore family,

• ode for a general interface to most of the ODE solvers,

• ode.band for solving models with a banded Jacobian,

• ode.1D for integrating 1-D models,

• ode.2D for integrating 2-D models,

• ode.3D for integrating 3-D models,

diagnostics to print diagnostic messages.

## Examples

 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 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 ## ======================================================================= ## Example 1: ## Various ways to solve the same model. ## ======================================================================= ## the model, 5 state variables f1 <- function (t, y, parms) { ydot <- vector(len = 5) ydot[1] <- 0.1*y[1] -0.2*y[2] ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3] ydot[3] <- -0.3*y[2] +0.1*y[3] -0.2*y[4] ydot[4] <- -0.3*y[3] +0.1*y[4] -0.2*y[5] ydot[5] <- -0.3*y[4] +0.1*y[5] return(list(ydot)) } ## the Jacobian, written as a full matrix fulljac <- function (t, y, parms) { jac <- matrix(nrow = 5, ncol = 5, byrow = TRUE, data = c(0.1, -0.2, 0 , 0 , 0 , -0.3, 0.1, -0.2, 0 , 0 , 0 , -0.3, 0.1, -0.2, 0 , 0 , 0 , -0.3, 0.1, -0.2, 0 , 0 , 0 , -0.3, 0.1)) return(jac) } ## the Jacobian, written in banded form bandjac <- function (t, y, parms) { jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE, data = c( 0 , -0.2, -0.2, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, -0.3, -0.3, -0.3, -0.3, 0)) return(jac) } ## initial conditions and output times yini <- 1:5 times <- 1:20 ## default: stiff method, internally generated, full Jacobian out <- lsode(yini, times, f1, parms = 0, jactype = "fullint") ## stiff method, user-generated full Jacobian out2 <- lsode(yini, times, f1, parms = 0, jactype = "fullusr", jacfunc = fulljac) ## stiff method, internally-generated banded Jacobian ## one nonzero band above (up) and below(down) the diagonal out3 <- lsode(yini, times, f1, parms = 0, jactype = "bandint", bandup = 1, banddown = 1) ## stiff method, user-generated banded Jacobian out4 <- lsode(yini, times, f1, parms = 0, jactype = "bandusr", jacfunc = bandjac, bandup = 1, banddown = 1) ## non-stiff method out5 <- lsode(yini, times, f1, parms = 0, mf = 10) ## ======================================================================= ## Example 2: ## diffusion on a 2-D grid ## partially specified Jacobian ## ======================================================================= diffusion2D <- function(t, Y, par) { y <- matrix(nrow = n, ncol = n, data = Y) dY <- r*y # production ## diffusion in X-direction; boundaries = 0-concentration Flux <- -Dx * rbind(y[1,],(y[2:n,]-y[1:(n-1),]),-y[n,])/dx dY <- dY - (Flux[2:(n+1),]-Flux[1:n,])/dx ## diffusion in Y-direction Flux <- -Dy * cbind(y[,1],(y[,2:n]-y[,1:(n-1)]),-y[,n])/dy dY <- dY - (Flux[,2:(n+1)]-Flux[,1:n])/dy return(list(as.vector(dY))) } ## parameters dy <- dx <- 1 # grid size Dy <- Dx <- 1 # diffusion coeff, X- and Y-direction r <- 0.025 # production rate times <- c(0, 1) n <- 50 y <- matrix(nrow = n, ncol = n, 0) pa <- par(ask = FALSE) ## initial condition for (i in 1:n) { for (j in 1:n) { dst <- (i - n/2)^2 + (j - n/2)^2 y[i, j] <- max(0, 1 - 1/(n*n) * (dst - n)^2) } } filled.contour(y, color.palette = terrain.colors) ## ======================================================================= ## jacfunc need not be estimated exactly ## a crude approximation, with a smaller bandwidth will do. ## Here the half-bandwidth 1 is used, whereas the true ## half-bandwidths are equal to n. ## This corresponds to ignoring the y-direction coupling in the ODEs. ## ======================================================================= print(system.time( for (i in 1:20) { out <- lsode(func = diffusion2D, y = as.vector(y), times = times, parms = NULL, jactype = "bandint", bandup = 1, banddown = 1) filled.contour(matrix(nrow = n, ncol = n, out[2,-1]), zlim = c(0,1), color.palette = terrain.colors, main = i) y <- out[2, -1] } )) par(ask = pa)

### Example output

user  system elapsed
0.747   0.004   0.757

deSolve documentation built on July 17, 2017, 3 a.m.