Package 'deSolve'

General Solvers for Initial Value Problems of Ordinary Differential Equations (ODE), Partial Differential Equations (PDE), Differential Algebraic Equations (DAE) and delay differential equations (DDE).

Description

Functions that solve initial value problems of a system of first-order ordinary differential equations (ODE), of partial differential equations (PDE), of differential algebraic equations (DAE) and delay differential equations.

The functions provide an interface to the FORTRAN functions lsoda, lsodar, lsode, lsodes of the ODEPACK collection, to the FORTRAN functions dvode, zvode and daspk, and a C-implementation of solvers of the Runge-Kutta family with fixed or variable time steps.

The package contains routines designed for solving ODEs resulting from 1-D, 2-D and 3-D partial differential equations (PDE) that have been converted to ODEs by numerical differencing. It includes root-finding (or event location) and provides access to lagged variables and derivatives.

The system of differential equations is written as an R function or defined in compiled code that has been dynamically loaded, see package vignette compiledCode for details. The solvers may be used as part of a modeling package for differential equations, or for parameter estimation using any appropriate modeling tool for non-linear models in R such as optim, nls, nlm or nlme, or FME.

Package Vignettes, Examples, Online Resources

• Solving Initial Value Differential Equations in R (pdf, R code)

• Writing Code in Compiled Languages (pdf, R code)

• Examples in R (code), and in Fortran or C (doc/dynload, doc/dynload-dede)

• deSolve homepage: https://desolve.r-forge.r-project.org (Papers, Books, PDFs)

• Mailing list: mailto:[email protected]

Author(s)

Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer

References

Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer (2010): Solving Differential Equations in R: Package deSolve Journal of Statistical Software, 33(9), 1–25. doi:10.18637/jss.v033.i09

Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer (2010): Solving differential equations in R. The R Journal 2(2), 5-15. doi:10.32614/RJ-2010-013

Karline Soetaert, Thomas Petzoldt (2011): Solving ODEs, DAEs, DDEs and PDEs in R. Journal of Numerical Analysis, Industrial and Applied Mathematics (JNAIAM) 6(1-2), 51-65.

Karline Soetaert, Jeff Cash, Francesca Mazzia, (2012): Solving Differential Equations in R. Springer, 248 pp.

Alan 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.

L. R. Petzold, (1983): A Description of DASSL: A Differential/Algebraic System Solver, in Scientific Computing, R. S. Stepleman et al. (Eds.), North-Holland, Amsterdam, pp. 65-68.

P. N. Brown, G. D. Byrne, A. C. Hindmarsh (1989): VODE: A Variable Coefficient ODE Solver, SIAM J. Sci. Stat. Comput., 10, pp. 1038-1051. doi:10.1137/0910062

See also the references given on the specific help pages of the different methods.

See Also

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

ode.band for solving models with a banded Jacobian,

ode.1D, ode.2D, ode.3D, for integrating 1-D, 2-D and 3-D models,

dede for a general interface to the delay differential equation solvers,

lsoda, lsode, lsodes, lsodar, vode, for ODE solvers of the Livermore family,

daspk, for a DAE solver up to index 1, of the Livermore family,

radau for integrating DAEs up to index 3 using an implicit Runge-Kutta,

rk, rkMethod, rk4, euler for Runge-Kutta solvers,

DLLfunc, DLLres, for testing model implementations in compiled code,

forcings, events, for how to implement forcing functions (external variables) and events (sudden changes in state variables),

lagvalue, lagderiv, for how to get access to lagged values of state variables and derivatives.

Examples

library(deSolve)

## Chaos in the atmosphere
Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <-  a * X + Y * Z
    dY <-  b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}

parameters <- c(a = -8/3, b = -10, c = 28)
state      <- c(X = 1, Y = 1, Z = 1)
times      <- seq(0, 100, by = 0.01)

out <- ode(y = state, times = times, func = Lorenz, parms = parameters)

plot(out)

## add a 3D figure if package scatterplot3D is available
if (require(scatterplot3d))
  scatterplot3d(out[,-1], type = "l")

---

A Physiological Model of Unbalanced Algal Growth

Description

A phytoplankton model with uncoupled carbon and nitrogen assimilation as a function of light and Dissolved Inorganic Nitrogen (DIN) concentration.

Algal biomass is described via 3 different state variables:

• low molecular weight carbohydrates (LMW), the product of photosynthesis,

• storage molecules (RESERVE) and

• the biosynthetic and photosynthetic apparatus (PROTEINS).

All algal state variables are expressed in $\rm mmol\, C\, m^{-3}$. Only proteins contain nitrogen and chlorophyll, with a fixed stoichiometric ratio. As the relative amount of proteins changes in the algae, so does the N:C and the Chl:C ratio.

An additional state variable, dissolved inorganic nitrogen (DIN) has units of $\rm mmol\, N\, m^{-3}$.

The algae grow in a dilution culture (chemostat): there is constant inflow of DIN and outflow of culture water, including DIN and algae, at the same rate.

Two versions of the model are included.

• In the default model, there is a day-night illumination regime, i.e. the light is switched on and off at fixed times (where the sum of illuminated + dark period = 24 hours).

• In another version, the light is imposed as a forcing function data set.

This model is written in FORTRAN.

Usage

aquaphy(times, y, parms, PAR = NULL, ...)

Arguments

times	time sequence for which output is wanted; the first value of times must be the initial time,
y	the initial (state) values ("DIN", "PROTEIN", "RESERVE", "LMW"), in that order,
parms	vector or list with the aquaphy model parameters; see the example for the order in which these have to be defined.
PAR	a data set of the photosynthetically active radiation (light intensity), if NULL, on-off PAR is used,
...	any other parameters passed to the integrator ode (which solves the model).

Details

The model is implemented primarily to demonstrate the linking of FORTRAN with R-code.

The source can be found in the ‘doc/examples/dynload’ subdirectory of the package.

Author(s)

Karline Soetaert <[email protected]>

References

Lancelot, C., Veth, C. and Mathot, S. (1991). Modelling ice-edge phytoplankton bloom in the Scotia-Weddel sea sector of the Southern Ocean during spring 1988. Journal of Marine Systems 2, 333–346.

Soetaert, K. and Herman, P. (2008). A practical guide to ecological modelling. Using R as a simulation platform. Springer.

See Also

ccl4model, the CCl4 inhalation model.

Examples

## ======================================================
##
## Example 1. PAR an on-off function
##
## ======================================================


## -----------------------------
## the model parameters:
## -----------------------------

parameters <- c(maxPhotoSynt   = 0.125,      # mol C/mol C/hr
                rMortPHY       = 0.001,      # /hr
                alpha          = -0.125/150, # uEinst/m2/s/hr
                pExudation     = 0.0,        # -
                maxProteinSynt = 0.136,      # mol C/mol C/hr
                ksDIN          = 1.0,        # mmol N/m3
                minpLMW        = 0.05,       # mol C/mol C
                maxpLMW        = 0.15,       # mol C/mol C
                minQuotum      = 0.075,      # mol C/mol C
                maxStorage     = 0.23,       # /h
                respirationRate= 0.0001,     # /h
                pResp          = 0.4,        # -
                catabolismRate = 0.06,       # /h
                dilutionRate   = 0.01,       # /h
                rNCProtein     = 0.2,        # mol N/mol C
                inputDIN       = 10.0,       # mmol N/m3
                rChlN          = 1,          # g Chl/mol N
                parMean        = 250.,       # umol Phot/m2/s
                dayLength      = 15.         # hours
                )

## -----------------------------
## The initial conditions
## -----------------------------

state <- c(DIN    = 6.,     # mmol N/m3
          PROTEIN = 20.0,   # mmol C/m3
          RESERVE = 5.0,    # mmol C/m3
          LMW     = 1.0)    # mmol C/m3

## -----------------------------
## Running the model
## -----------------------------

times <- seq(0, 24*20, 1)

out <- as.data.frame(aquaphy(times, state, parameters))

## -----------------------------
## Plotting model output
## -----------------------------

par(mfrow = c(2, 2), oma = c(0, 0, 3, 0))
col <- grey(0.9)
ii <- 1:length(out$PAR)              

plot(times[ii], out$Chlorophyll[ii], type = "l",
      main = "Chlorophyll", xlab = "time, hours",ylab = "ug/l")
polygon(times[ii], out$PAR[ii]-10, col = col, border = NA); box()
lines(times[ii], out$Chlorophyll[ii], lwd = 2 )


plot (times[ii], out$DIN[ii], type = "l", main = "DIN",
      xlab = "time, hours",ylab = "mmolN/m3")
polygon(times[ii], out$PAR[ii]-10, col = col, border = NA); box()
lines(times[ii], out$DIN[ii], lwd = 2 )


plot (times[ii], out$NCratio[ii], type = "n", main = "NCratio",
      xlab = "time, hours", ylab = "molN/molC")
polygon(times[ii], out$PAR[ii]-10, col = col, border = NA); box()
lines(times[ii], out$NCratio[ii], lwd = 2 )


plot (times[ii], out$PhotoSynthesis[ii],type = "l",
       main = "PhotoSynthesis", xlab = "time, hours",
       ylab = "mmolC/m3/hr")
polygon(times[ii], out$PAR[ii]-10, col = col, border = NA); box()
lines(times[ii], out$PhotoSynthesis[ii], lwd = 2 )

mtext(outer = TRUE, side = 3, "AQUAPHY, PAR= on-off", cex = 1.5)

## -----------------------------
## Summary model output
## -----------------------------
t(summary(out))

## ======================================================
##
## Example 2. PAR a forcing function data set
##
## ======================================================

times <- seq(0, 24*20, 1)

## -----------------------------
## create the forcing functions
## -----------------------------

ftime  <- seq(0,500,by=0.5)
parval <- pmax(0,250 + 350*sin(ftime*2*pi/24)+
   (runif(length(ftime))-0.5)*250)
Par    <- matrix(nc=2,c(ftime,parval))


state <- c(DIN     = 6.,     # mmol N/m3
           PROTEIN = 20.0,   # mmol C/m3
           RESERVE = 5.0,    # mmol C/m3
           LMW     = 1.0)    # mmol C/m3
              
out <- aquaphy(times, state, parameters, Par)

plot(out, which = c("PAR", "Chlorophyll", "DIN", "NCratio"), 
     xlab = "time, hours", 
     ylab = c("uEinst/m2/s", "ug/l", "mmolN/m3", "molN/molC"))

mtext(outer = TRUE, side = 3, "AQUAPHY, PAR=forcing", cex = 1.5)

# Now all variables plotted in one figure...
plot(out, which = 1:9, type = "l")

par(mfrow = c(1, 1))

---

Closed Chamber Study of CCl4 Metabolism by Rats.

Description

The results of a closed chamber experiment to determine metabolic parameters for CCl4 (carbon tetrachloride) in rats.

Usage

data(ccl4data)

Format

This data frame contains the following columns:

time

the time (in hours after starting the experiment).

initconc

initial chamber concentration (ppm).

animal

this is a repeated measures design; this variable indicates which animal the observation pertains to.

ChamberConc

chamber concentration at time, in ppm.

Source

Evans, et al. 1994 Applications of sensitivity analysis to a physiologically based pharmacokinetic model for carbon tetrachloride in rats. Toxicology and Applied Pharmacology 128: 36 – 44.

Examples

plot(ChamberConc ~ time, data = ccl4data, xlab = "Time (hours)",
  xlim = range(c(0, ccl4data$time)),
  ylab = "Chamber Concentration (ppm)", log = "y")
ccl4data.avg <- aggregate(ccl4data$ChamberConc,
                  by = ccl4data[c("time", "initconc")], mean)
points(x ~ time, data = ccl4data.avg, pch = 16)

---

The CCl4 Inhalation Model

Description

The CCl4 inhalation model implemented in .Fortran

Usage

ccl4model(times, y, parms, ...)

Arguments

times	time sequence for which the model has to be integrated.
y	the initial values for the state variables ("AI", "AAM", "AT", "AF", "AL", "CLT" and "AM"), in that order.
parms	vector or list holding the ccl4 model parameters; see the example for the order in which these have to be defined.
...	any other parameters passed to the integrator ode (which solves the model).

Details

The model is implemented primarily to demonstrate the linking of FORTRAN with R-code.

The source can be found in the ‘/doc/examples/dynload’ subdirectory of the package.

Author(s)

R. Woodrow Setzer <[email protected]>

See Also

Try demo(CCL4model) for how this model has been fitted to the dataset ccl4data,

aquaphy, another FORTRAN model, describing growth in aquatic phytoplankton.

Examples

## =================
## Parameter values
## =================

Pm <- c(
  ## Physiological parameters
  BW = 0.182,   # Body weight (kg)
  QP = 4.0  ,   # Alveolar ventilation rate (hr^-1)
  QC = 4.0  ,   # Cardiac output (hr^-1)
  VFC = 0.08,   # Fraction fat tissue (kg/(kg/BW))
  VLC = 0.04,   # Fraction liver tissue (kg/(kg/BW))
  VMC = 0.74,   # Fraction of muscle tissue (kg/(kg/BW))
  QFC = 0.05,   # Fractional blood flow to fat ((hr^-1)/QC
  QLC = 0.15,   # Fractional blood flow to liver ((hr^-1)/QC)
  QMC = 0.32,   # Fractional blood flow to muscle ((hr^-1)/QC)

  ## Chemical specific parameters for chemical
  PLA = 16.17,  # Liver/air partition coefficient
  PFA = 281.48, # Fat/air partition coefficient
  PMA = 13.3,   # Muscle/air partition coefficient
  PTA = 16.17,  # Viscera/air partition coefficient
  PB = 5.487,   # Blood/air partition coefficient
  MW = 153.8,   # Molecular weight (g/mol)
  VMAX = 0.04321671, # Max. velocity of metabolism (mg/hr) -calibrated
  KM = 0.4027255,    # Michaelis-Menten constant (mg/l) -calibrated

  ## Parameters for simulated experiment
  CONC = 1000,  # Inhaled concentration
  KL = 0.02,    # Loss rate from empty chamber /hr
  RATS = 1.0,   # Number of rats enclosed in chamber
  VCHC = 3.8    # Volume of closed chamber (l)
)

## ================
## State variables
## ================
y <- c(
  AI  = 21, # total mass , mg
  AAM = 0,
  AT  = 0,
  AF  = 0,
  AL  = 0,
  CLT = 0,  # area under the conc.-time curve in the liver
  AM  = 0   # the amount metabolized (AM)
)

## ==================
## Model application
## ==================

times <- seq(0, 6, by = 0.1)

## initial inhaled concentration-calibrated
conc <- c(26.496, 90.197, 245.15, 951.46) 

plot(ChamberConc ~ time, data = ccl4data, xlab = "Time (hours)",
       xlim = range(c(0, ccl4data$time)),
       ylab = "Chamber Concentration (ppm)",
       log = "y", main = "ccl4model")

for (cc in conc) {
  Pm["CONC"] <- cc

  VCH <- Pm[["VCHC"]] - Pm[["RATS"]] * Pm[["BW"]]
  AI0 <- VCH * Pm[["CONC"]] * Pm[["MW"]]/24450 
  y["AI"] <- AI0
  
  ## run the model:
  out <- as.data.frame(ccl4model(times, y, Pm))
  lines(out$time, out$CP, lwd = 2)
}

legend("topright", lty = c(NA, 1), pch = c(1, NA), lwd = c(NA, 2),
       legend = c("data", "model"))

## ================================== 
## An example with tracer injection
## ==================================
## every day, a conc of 2 is added to AI.
## 1. implemented as a data.frame
eventdat <- data.frame(var = rep("AI", 6), time = 1:6 ,
  value = rep(1, 6), method = rep("add", 6))

eventdat

print(system.time(
  out <-ccl4model(times, y, Pm, events = list(data = eventdat))
))

plot(out, mfrow = c(3, 4), type = "l", lwd = 2)


# 2. implemented as a function in a DLL!
print(system.time(
out2 <-ccl4model(times, y, Pm, events = list(func = "eventfun", time = 1:6))
))

plot(out2, mfrow=c(3, 4), type = "l", lwd = 2)

---

Check shared library (DLL/.so) of a compiled model.

Description

Check shared library (DLL/.so) of a compiled model and create a list of symbols.

Usage

checkDLL(func, jacfunc, dllname, initfunc, verbose, nout, outnames, JT = 1)

Arguments

func	character: name of the derivative function.
jacfunc	an R function, that computes the Jacobian of the system of differential equations $\partial\dot{y}_i/\partial y_j$, or a string giving the name of a function or subroutine in ‘dllname’ that computes the Jacobian.
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.
initfunc	the name of the initialisation function (which initialises values of parameters), as provided in ‘dllname’. See package vignette "compiledCode".
verbose	reserved for future extensions.
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.
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.
JT	integer specifying the type of the Jacobian. The default value of 1 must be set to 2 for solver lsodes if 'func' is specified in a DLL or inline compiled and if a 'jacfunc' is provided.

Details

The function checkDLL is normally called internally by the solver functions. It can be used to avoid overhead, when a small compiled model with a low number of integration steps is repeatedly called. The feature is currently only available for the lsoda solver.

Value

List of class deSolve.symbols with:

ModelInit	pointer to the init function of the DLL (class "externalptr").
Func	pointer to the derivative function in the DLL (class "externalptr").
JacFunc	pointer to the Jacobi function in the DLL (class "externalptr").
Nglobal	number of output variables calculated in the compiled function.
Nmtot	list of names of derivatives and output variables.

See Also

lsoda

Examples

## Not run: 
symbols <- checkDLL(func = "derivs", jacfunc = NULL, dllname = "lorenzc",
             initfunc = "initmod", verbose = TRUE, nout = 0,
             outnames = NULL, JT = 1)

## End(Not run)

---

Find Nearest Event for Each Time Step and Clean Time Steps to Avoid Doubles

Description

These functions can be used for checking time steps and events used by ode solver functions. They are normally called internally within the solvers.

Usage

nearestEvent(times, eventtimes)
cleanEventTimes(times, eventtimes, eps = .Machine$double.eps * 10)

Arguments

times	the vector of output times,
eventtimes	a vector with the event times,
eps	relative tolerance value below which two numbers are assumed to be numerically equal.

Details

In floating point arithmetics, problems can occur if values have to be compared for 'equality' but are only close to each other and not exactly the same.

The utility functions can be used to add all eventtimes to the output times vector, but without including times that are very close to an event.

This means that all values of eventtimes are contained but only the subset of times that have no close neighbors in eventtimes.

These checks are normally performed internally by the integration solvers.

Value

nearestEvent returns a vector with the closest events for each time step and

cleanEventTimes returns a vector with the output times without all those that are 'very close' to an event.

Author(s)

Thomas Petzoldt

See Also

events

Examples

events <- sort(c(0, 2, 3, 4 + 1e-10, 5, 7 - 1e-10,
                 7 + 6e-15, 7.5, 9, 24.9999, 25, 80, 1001, 1e300))
times  <- sort(c(0, 1:7, 4.5, 6.75, 7.5, 9.2, 9.0001, 25, 879, 1e3, 1e300+5))

nearest <- nearestEvent(times, events)
data.frame(times=times, nearest = nearest)

## typical usage: include all events in times after removing values that
## are numerically close together, events have priority 
times
unique_times <- cleanEventTimes(times, events)
newtimes <- sort(c(unique_times, events))
newtimes

---

Solver for Differential Algebraic Equations (DAE)

Description

Solves either:

• a system of ordinary differential equations (ODE) of the form

  $y' = f(t, y, ...)$

  or

• a system of differential algebraic equations (DAE) of the form

  $F(t,y,y') = 0$

  or

• a system of linearly implicit DAES in the form

  $M y' = f(t, y)$

using a combination of backward differentiation formula (BDF) and a direct linear system solution method (dense or banded).

The R function daspk provides an interface to the FORTRAN DAE solver of the same name, written by Linda R. Petzold, Peter N. Brown, Alan C. Hindmarsh and Clement W. Ulrich.

The system of DE's is written as an R function (which may, of course, use .C, .Fortran, .Call, etc., to call foreign code) or be defined in compiled code that has been dynamically loaded.

Usage

daspk(y, times, func = NULL, parms, nind = c(length(y), 0, 0), 
  dy = NULL, res = NULL, nalg = 0, 
  rtol = 1e-6, atol = 1e-6, jacfunc = NULL,
  jacres = NULL, jactype = "fullint", mass = NULL, estini = NULL,
  verbose = FALSE, tcrit = NULL, hmin = 0, hmax = NULL,
  hini = 0, ynames = TRUE, maxord = 5, 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, ...)

Arguments

y	the initial (state) values for the DE system. If y has a name attribute, the names will be used to label the output matrix.
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	to be used if the model is an ODE, or a DAE written in linearly implicit form (M y' = f(t, y)). func should be an R-function that computes the values of the derivatives in the ODE system (the model definition) at time t. func 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, unless ynames is FALSE. 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 should be specified in the same order as the specification of the state variables y. Note that it is not possible to define func as a compiled function in a dynamically loaded shared library. Use res instead.
parms	vector or list of parameters used in func, jacfunc, or res
nind	if a DAE system: a three-valued vector with the number of variables of index 1, 2, 3 respectively. The equations must be defined such that the index 1 variables precede the index 2 variables which in turn precede the index 3 variables. The sum of the variables of different index should equal N, the total number of variables. Note that this has been added for consistency with radau. If used, then the variables are weighed differently than in the original daspk code, i.e. index 2 variables are scaled with 1/h, index 3 variables are scaled with 1/h^2. In some cases this allows daspk to solve index 2 or index 3 problems.
dy	the initial derivatives of the state variables of the DE system. Ignored if an ODE.
res	if a DAE system: either an R-function that computes the residual function $F(t,y,y')$ of the DAE system (the model defininition) at time t, or a character string giving the name of a compiled function in a dynamically loaded shared library. If res is a user-supplied R-function, it must be defined as: res <- function(t, y, dy, parms, ...). Here t is the current time point in the integration, y is the current estimate of the variables in the ODE system, dy are the corresponding derivatives. If the initial y or dy have a names attribute, the names will be available inside res, unless ynames is FALSE. parms is a vector of parameters. The return value of res should be a list, whose first element is a vector containing the residuals of the DAE system, i.e. $\delta = F(t,y,y')$, and whose next elements contain output variables that are required at each point in times. If res is a string, then dllname must give the name of the shared library (without extension) which must be loaded before daspk() is called (see package vignette "compiledCode" for more information).
nalg	if a DAE system: the number of algebraic equations (equations not involving derivatives). Algebraic equations should always be the last, i.e. preceeded by the differential equations. Only used if estini = 1.
rtol	relative error tolerance, either a scalar or a vector, one value for each y,
atol	absolute error tolerance, either a scalar or a vector, one value for each y.
jacfunc	if not NULL, an R function that computes the Jacobian of the system of differential equations. Only used in case the system is an ODE ($y' = f(t, y)$), specified by func. The R calling sequence for jacfunc is identical to that of func. If the Jacobian is a full matrix, jacfunc should return a matrix $\partial\dot{y}/\partial y$, 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). 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.
jacres	jacres and not jacfunc should be used if the system is specified by the residual function $F(t, y, y')$, i.e. jacres is used in conjunction with res. If jacres is an R-function, the calling sequence for jacres is identical to that of res, but with extra parameter cj. Thus it should be called as: jacres = func(t, y, dy, parms, cj, ...). Here t is the current time point in the integration, y is the current estimate of the variables in the ODE system, $y'$ are the corresponding derivatives and cj is a scalar, which is normally proportional to the inverse of the stepsize. If the initial y or dy have a names attribute, the names will be available inside jacres, unless ynames is FALSE. parms is a vector of parameters (which may have a names attribute). If the Jacobian is a full matrix, jacres should return the matrix $dG/dy + c_j\cdot dG/dy'$, where the $i$th row is the sum of the derivatives of $G_i$ with respect to $y_j$ and the scaled derivatives of $G_i$ with respect to $y'_j$. If the Jacobian is banded, jacres should return only the nonzero bands of the Jacobian, rotated rowwise. See details for the calling sequence when jacres is a string.
jactype	the structure of the Jacobian, one of "fullint", "fullusr", "bandusr" or "bandint" - either full or banded and estimated internally or by the user.
mass	the mass matrix. If not NULL, the problem is a linearly implicit DAE and defined as $M\, dy/dt = f(t,y)$. The mass-matrix $M$ should be of dimension $n*n$ where $n$ is the number of $y$-values. If mass=NULL then the model is either an ODE or a DAE, specified with res
estini	only if a DAE system, and if initial values of y and dy are not consistent (i.e. $F(t,y,dy) \neq 0$), setting estini = 1 or 2, will solve for them. If estini = 1: dy and the algebraic variables are estimated from y; in this case, the number of algebraic equations must be given (nalg). If estini = 2: y will be estimated from dy.
verbose	if TRUE: full output to the screen, e.g. will print the diagnostiscs of the integration - see details.
tcrit	the FORTRAN routine daspk 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 large models.
maxord	the maximum order to be allowed. Reduce maxord to save storage space ( <= 5)
bandup	number of non-zero bands above the diagonal, in case the Jacobian is banded (and jactype one of "bandint", "bandusr")
banddown	number of non-zero bands below the diagonal, in case the Jacobian is banded (and jactype one of "bandint", "bandusr")
maxsteps	maximal number of steps per output interval taken by the solver; will be recalculated to be at least 500 and a multiple of 500; if verbose is TRUE the solver will give a warning if more than 500 steps are taken, but it will continue till maxsteps steps. (Note this warning was always given in deSolve versions < 1.10.3).
dllname	a string giving the name of the shared library (without extension) that contains all the compiled function or subroutine definitions referred to in res and jacres. 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 res and jacres.
ipar	only when ‘dllname’ is specified: a vector with integer values passed to the dll-functions whose names are specified by res and jacres.
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 res, 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 res, 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. See forcings or vignette compiledCode.
events	A matrix or data frame that specifies events, i.e. when the value of a state variable is suddenly changed. See events for more information.
lags	A list that specifies timelags, i.e. the number of steps that has to be kept. To be used for delay differential equations. See timelags, dede for more information.
...	additional arguments passed to func, jacfunc, res and jacres, allowing this to be a generic function.

Details

The daspk solver uses the backward differentiation formulas of orders one through five (specified with maxord) to solve either:

• an ODE system of the form

  $y' = f(t,y,...)$

  or

• a DAE system of the form

  $y' = M f(t,y,...)$

  or

• a DAE system of the form

  $F(t,y,y') = 0$

  . The index of the DAE should be preferable <= 1.

ODEs are specified using argument func, DAEs are specified using argument res.

If a DAE system, Values for y and y' (argument dy) at the initial time must be given as input. Ideally, these values should be consistent, that is, if t, y, y' are the given initial values, they should satisfy F(t,y,y') = 0.
However, if consistent values are not known, in many cases daspk can solve for them: when estini = 1, y' and algebraic variables (their number specified with nalg) will be estimated, when estini = 2, y will be estimated.

The form of the Jacobian can be specified by jactype. This is one of:

jactype = "fullint":

a full Jacobian, calculated internally by daspk, the default,

jactype = "fullusr":

a full Jacobian, specified by user function jacfunc or jacres,

jactype = "bandusr":

a banded Jacobian, specified by user function jacfunc or jacres; the size of the bands specified by bandup and banddown,

jactype = "bandint":

a banded Jacobian, calculated by daspk; the size of the bands specified by bandup and banddown.

If jactype = "fullusr" or "bandusr" then the user must supply a subroutine jacfunc.

If jactype = "fullusr" or "bandusr" then the user must supply a subroutine jacfunc or jacres.

The input parameters rtol, and atol determine the error control performed by the solver. If the request for precision exceeds the capabilities of the machine, daspk will return an error code. See lsoda for details.

When the index of the variables is specified (argument nind), and higher index variables are present, then the equations are scaled such that equations corresponding to index 2 variables are multiplied with 1/h, for index 3 they are multiplied with 1/h^2, where h is the time step. This is not in the standard DASPK code, but has been added for consistency with solver radau. Because of this, daspk can solve certain index 2 or index 3 problems.

res and jacres may be defined in compiled C or FORTRAN code, as well as in an R-function. See package vignette "compiledCode" for details. Examples in FORTRAN are in the ‘dynload’ subdirectory of the deSolve package directory.

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.

More information about models defined in compiled code is in the package vignette ("compiledCode"); information about linking forcing functions to compiled code is in forcings.

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

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 or res, plus an additional column (the first) for the time value. There will be one row for each element in times unless the FORTRAN routine ‘daspk’ returns with an unrecoverable error. If y has a names attribute, it will be used to label the columns of the output value.

Note

In this version, the Krylov method is not (yet) supported.

From deSolve version 1.10.4 and above, the following changes were made

1. the argument list to daspk now also includes nind, the index of each variable. This is used to scale the variables, such that daspk in R can also solve certain index 2 or index 3 problems, which the original Fortran version may not be able to solve.

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

3. the multiple warnings from daspk when the number of steps exceed 500 were toggled off unless verbose is TRUE

Author(s)

Karline Soetaert <[email protected]>

References

L. R. Petzold, A Description of DASSL: A Differential/Algebraic System Solver, in Scientific Computing, R. S. Stepleman et al. (Eds.), North-Holland, Amsterdam, 1983, pp. 65-68.

K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, Elsevier, New York, 1989.

P. N. Brown and A. C. Hindmarsh, Reduced Storage Matrix Methods in Stiff ODE Systems, J. Applied Mathematics and Computation, 31 (1989), pp. 40-91. doi:10.1016/0096-3003(89)90110-0

P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Using Krylov Methods in the Solution of Large-Scale Differential-Algebraic Systems, SIAM J. Sci. Comp., 15 (1994), pp. 1467-1488. doi:10.1137/0915088

P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Consistent Initial Condition Calculation for Differential-Algebraic Systems, LLNL Report UCRL-JC-122175, August 1995; submitted to SIAM J. Sci. Comp.

Netlib: https://netlib.org

See Also

• radau for integrating DAEs up to index 3,

• rk,

• rk4 and euler for Runge-Kutta integrators.

• lsoda, lsode, lsodes, lsodar, vode, 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

## =======================================================================
## Coupled chemical reactions including an equilibrium
## modeled as (1) an ODE and (2) as a DAE
##
## The model describes three chemical species A,B,D:
## subjected to equilibrium reaction D <- > A + B
## D is produced at a constant rate, prod
## B is consumed at 1s-t order rate, r
## Chemical problem formulation 1: ODE
## =======================================================================

## Dissociation constant
K <- 1 

## parameters
pars <- c(
        ka   = 1e6,     # forward rate
        r    = 1,
        prod = 0.1)


Fun_ODE <- function (t, y, pars)
{
  with (as.list(c(y, pars)), {
    ra  <- ka*D        # forward rate
    rb  <- ka/K *A*B   # backward rate

    ## rates of changes
    dD  <- -ra + rb + prod
    dA  <-  ra - rb
    dB  <-  ra - rb - r*B
    return(list(dy = c(dA, dB, dD),
                CONC = A+B+D))
  })
}

## =======================================================================
## Chemical problem formulation 2: DAE
## 1. get rid of the fast reactions ra and rb by taking
## linear combinations   : dD+dA = prod (res1) and
##                         dB-dA = -r*B (res2)
## 2. In addition, the equilibrium condition (eq) reads:
## as ra = rb : ka*D = ka/K*A*B = >      K*D = A*B
## =======================================================================

Res_DAE <- function (t, y, yprime, pars)
{
  with (as.list(c(y, yprime, pars)), {

    ## residuals of lumped rates of changes
    res1 <- -dD - dA + prod
    res2 <- -dB + dA - r*B
    
    ## and the equilibrium equation
    eq   <- K*D - A*B

    return(list(c(res1, res2, eq),
                CONC = A+B+D))
  })
}

## =======================================================================
## Chemical problem formulation 3: Mass * Func
## Based on the DAE formulation
## =======================================================================

Mass_FUN <- function (t, y, pars) {
  with (as.list(c(y, pars)), {

    ## as above, but without the 
    f1 <- prod
    f2 <- - r*B
    
    ## and the equilibrium equation
    f3   <- K*D - A*B

    return(list(c(f1, f2, f3),
                CONC = A+B+D))
  })
}
Mass <- matrix(nrow = 3, ncol = 3, byrow = TRUE, 
  data=c(1,  0, 1,         # dA + 0 + dB
        -1,  1, 0,         # -dA + dB +0
         0,  0, 0))        # algebraic
         
times <- seq(0, 100, by = 2)

## Initial conc; D is in equilibrium with A,B
y     <- c(A = 2, B = 3, D = 2*3/K)

## ODE model solved with daspk
ODE <- daspk(y = y, times = times, func = Fun_ODE,
                     parms = pars, atol = 1e-10, rtol = 1e-10)

## Initial rate of change
dy  <- c(dA = 0, dB = 0, dD = 0) 

## DAE model solved with daspk
DAE <- daspk(y = y, dy = dy, times = times,
         res = Res_DAE, parms = pars, atol = 1e-10, rtol = 1e-10)

MASS<- daspk(y=y, times=times, func = Mass_FUN, parms = pars, mass = Mass)

## ================
## plotting output
## ================

plot(ODE, DAE, xlab = "time", ylab = "conc", type = c("l", "p"),
     pch = c(NA, 1))

legend("bottomright", lty = c(1, NA), pch = c(NA, 1),
  col = c("black", "red"), legend = c("ODE", "DAE"))

# difference between both implementations:
max(abs(ODE-DAE))

## =======================================================================
## same DAE model, now with the Jacobian
## =======================================================================
jacres_DAE <- function (t, y, yprime, pars, cj)
{
    with (as.list(c(y, yprime, pars)), {
##    res1 = -dD - dA + prod
      PD[1,1] <- -1*cj      # d(res1)/d(A)-cj*d(res1)/d(dA)
      PD[1,2] <- 0          # d(res1)/d(B)-cj*d(res1)/d(dB)
      PD[1,3] <- -1*cj      # d(res1)/d(D)-cj*d(res1)/d(dD)
##     res2 = -dB + dA - r*B
      PD[2,1] <- 1*cj
      PD[2,2] <- -r -1*cj
      PD[2,3] <- 0
##    eq = K*D - A*B
      PD[3,1] <- -B
      PD[3,2] <- -A
      PD[3,3] <- K
      return(PD)
   })
}

PD <- matrix(ncol = 3, nrow = 3, 0)

DAE2 <- daspk(y = y, dy = dy, times = times,
          res = Res_DAE, jacres = jacres_DAE, jactype = "fullusr",
          parms = pars, atol = 1e-10, rtol = 1e-10)
         
max(abs(DAE-DAE2))

## See \dynload subdirectory for a FORTRAN implementation of this model

## =======================================================================
## The chemical model as a DLL, with production a forcing function
## =======================================================================
times <- seq(0, 100, by = 2)

pars <- c(K = 1, ka   = 1e6, r    = 1)

## Initial conc; D is in equilibrium with A,B
y     <- c(A = 2, B = 3, D = as.double(2*3/pars["K"]))

## Initial rate of change
dy  <- c(dA = 0, dB = 0, dD = 0)

# production increases with time
prod <- matrix(ncol = 2, 
               data = c(seq(0, 100, by = 10), 0.1*(1+runif(11)*1)))

ODE_dll <- daspk(y = y, dy = dy, times = times, res = "chemres",
          dllname = "deSolve", initfunc = "initparms",
          initforc = "initforcs", parms = pars, forcings = prod,
          atol = 1e-10, rtol = 1e-10, nout = 2, 
          outnames = c("CONC","Prod"))

plot(ODE_dll, which = c("Prod", "D"), xlab = "time",
     ylab = c("/day", "conc"), main = c("production rate","D"))

---

General Solver for Delay Differential Equations.

Description

Function dede is a general solver for delay differential equations, i.e. equations where the derivative depends on past values of the state variables or their derivatives.

Usage

dede(y, times, func=NULL, parms,
    method = c( "lsoda", "lsode", "lsodes", "lsodar", "vode",
       "daspk", "bdf", "adams", "impAdams", "radau"), control = NULL, ...)

Arguments

y	the initial (state) values for the DE system, a vector. If y has a name attribute, the names will be used to label the output matrix.
times	time sequence for which output is wanted; the first value of times must be the initial time.
func	an R-function that computes the values of the derivatives in the ODE system (the model definition) at time $t$. func 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 DE 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. If method "daspk" is used, then func can be NULL, in which case res should be used.
parms	parameters passed to func.
method	the integrator to use, either a string ("lsoda", "lsode", "lsodes", "lsodar", "vode", "daspk", "bdf", "adams", "impAdams", "radau") or a function that performs the integration. The default integrator used is lsoda.
control	a list that can supply (1) the size of the history array, as control$mxhist; the default is 1e4 and (2) how to interpolate, as control$interpol, where 1 is hermitian interpolation, 2 is variable order interpolation, using the Nordsieck history array. Only for the two Adams methods is the second option recommended.
...	additional arguments passed to the integrator.

Details

Functions lagvalue and lagderiv are to be used with dede as they provide access to past (lagged) values of state variables and derivatives. The number of past values that are to be stored in a history matrix, can be specified in control$mxhist. The default value (if unspecified) is 1e4.

Cubic Hermite interpolation is used by default to obtain an accurate interpolant at the requested lagged time. For methods adams, impAdams, a more accurate interpolation method can be triggered by setting control$interpol = 2.

dede does not deal explicitly with propagated derivative discontinuities, but relies on the integrator to control the stepsize in the region of a discontinuity.

dede does not include methods to deal with delays that are smaller than the stepsize, although in some cases it may be possible to solve such models.

For these reasons, it can only solve rather simple delay differential equations.

When used together with integrator lsodar, or lsode, dde can simultaneously locate a root, and trigger an event. See last example.

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 second element of the return from func, plus an additional column (the first) for the time value. There will be one row for each element in times unless the integrator 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]>

See Also

lagvalue, lagderiv,for how to specify lagged variables and derivatives.

Examples

## =============================================================================
## A simple delay differential equation
## dy(t) = -y(t-1) ; y(t<0)=1
## =============================================================================

##-----------------------------
## the derivative function
##-----------------------------
derivs <- function(t, y, parms) {
  if (t < 1)
    dy <- -1
  else
    dy <- - lagvalue(t - 1)
  list(c(dy))
}

##-----------------------------
## initial values and times
##-----------------------------
yinit <- 1
times <- seq(0, 30, 0.1)

##-----------------------------
## solve the model
##-----------------------------
yout <- dede(y = yinit, times = times, func = derivs, parms = NULL)

##-----------------------------
## display, plot results
##-----------------------------
plot(yout, type = "l", lwd = 2, main = "dy/dt = -y(t-1)")

## =============================================================================
## The infectuous disease model of Hairer; two lags.
## example 4 from Shampine and Thompson, 2000
## solving delay differential equations with dde23
## =============================================================================

##-----------------------------
## the derivative function
##-----------------------------
derivs <- function(t,y,parms) {
  if (t < 1)
    lag1 <- 0.1
  else
    lag1 <- lagvalue(t - 1,2)
  if (t < 10)
    lag10 <- 0.1
  else
    lag10 <- lagvalue(t - 10,2)

  dy1 <- -y[1] * lag1 + lag10
  dy2 <-  y[1] * lag1 - y[2]
  dy3 <-  y[2] - lag10
  list(c(dy1, dy2, dy3))
}

##-----------------------------
## initial values and times
##-----------------------------
yinit <- c(5, 0.1, 1)
times <- seq(0, 40, by = 0.1)

##-----------------------------
## solve the model
##-----------------------------
system.time(
  yout <- dede(y = yinit, times = times, func = derivs, parms = NULL)
)

##-----------------------------
## display, plot results
##-----------------------------
matplot(yout[,1], yout[,-1], type = "l", lwd = 2, lty = 1,
  main = "Infectuous disease - Hairer")

## =============================================================================
## time lags + EVENTS triggered by a root function
## The two-wheeled suitcase model
## example 8 from Shampine and Thompson, 2000
## solving delay differential equations with dde23
## =============================================================================

##-----------------------------
## the derivative function
##-----------------------------
derivs <- function(t, y, parms) {
  if (t < tau)
    lag <- 0
  else
    lag <- lagvalue(t - tau)

  dy1 <- y[2]
  dy2 <- -sign(y[1]) * gam * cos(y[1]) +
          sin(y[1]) - bet * lag[1] + A * sin(omega * t + mu)
  list(c(dy1, dy2))
}

## root and event function
root <- function(t,y,parms) ifelse(t>0, return(y), return(1))
event <- function(t,y,parms) return(c(y[1], y[2]*0.931))

gam = 0.248; bet = 1; tau = 0.1; A = 0.75
omega = 1.37; mu = asin(gam/A)

##-----------------------------
## initial values and times
##-----------------------------
yinit <- c(y = 0, dy = 0)
times <- seq(0, 12, len = 1000)

##-----------------------------
## solve the model
##-----------------------------
## Note: use a solver that supports both root finding and events,
##       e.g. lsodar, lsode, lsoda, adams, bdf
yout <- dede(y = yinit, times = times, func = derivs, parms = NULL,
  method = "lsodar", rootfun = root, events = list(func = event, root = TRUE))

##-----------------------------
## display, plot results
##-----------------------------

plot(yout, which = 1, type = "l", lwd = 2, main = "suitcase model", mfrow = c(1,2))
plot(yout[,2], yout[,3], xlab = "y", ylab = "dy", type = "l", lwd = 2)

---

Print Diagnostic Characteristics of Solvers

Description

Prints several diagnostics of the simulation to the screen, e.g. number of steps taken, the last step size, ...

Usage

diagnostics(obj, ...)
## Default S3 method:
diagnostics(obj, ...)

Arguments

obj	is an output data structure produced by one of the solver routines.
...	optional arguments allowing to extend diagnostics as a generic function.

Details

Detailed information obout the success of a simulation is printed, if a diagnostics function exists for a specific solver routine. A warning is printed, if no class-specific diagnostics exists.

Please consult the class-specific help page for details.

See Also

diagnostics.deSolve for diagnostics of differential equaton solvers.

---

Print Diagnostic Characteristics of ODE and DAE Solvers

Description

Prints several diagnostics of the simulation to the screen, e.g. number of steps taken, the last step size, ...

Usage

## S3 method for class 'deSolve'
diagnostics(obj, Full = FALSE, ...)

Arguments

obj	is the output matrix as produced by one of the integration routines.
Full	when TRUE then all messages will be printed, including the ones that are not relevant for the solver. If FALSE, then only the relevant messages will be printed.
...	optional arguments allowing to extend diagnostics as a generic function.

Details

When the integration output is saved as a data.frame, then the required attributes are lost and method diagnostics will not work anymore.

Value

The integer and real vector with diagnostic values; for function lsodar also the root information.

See tables 2 and 3 in vignette("deSolve") for what these vectors contain.

Note: the number of function evaluations are *without* the extra calls performed to generate the ordinary output variables (if present).

Examples

## The famous Lorenz equations: chaos in the earth's atmosphere
## Lorenz 1963. J. Atmos. Sci. 20, 130-141.

chaos <- function(t, state, parameters) {
  with(as.list(c(state)), {

    dx     <- -8/3 * x + y * z
    dy     <- -10 * (y - z)
    dz     <- -x * y + 28 * y - z

    list(c(dx, dy, dz))
  })
}

state <- c(x = 1, y = 1, z = 1)
times <- seq(0, 50, 0.01)
out   <- vode(state, times, chaos, 0)
pairs(out, pch = ".")
diagnostics(out)

---

Evaluates a Derivative Function Represented in a DLL

Description

Calls a function, defined in a compiled language as a DLL

Usage

DLLfunc(func, times, y, parms, dllname,
  initfunc = dllname, rpar = NULL, ipar = NULL, nout = 0,
  outnames = NULL, forcings = NULL, initforc = NULL,
  fcontrol = NULL)

Arguments

func	the name of the function in the dynamically loaded shared library,
times	first value = the time at which the function needs to be evaluated,
y	the values of the dependent variables for which the function needs to be evaluated,
parms	the parameters that are passed to the initialiser function,
dllname	a string giving the name of the shared library (without extension) that contains the compiled function or subroutine definitions referred to in func,
initfunc	if not NULL, the name of the initialisation function (which initialises values of parameters), as provided in ‘dllname’. See details.
rpar	a vector with double precision values passed to the DLL-function func and jacfunc present in the DLL, via argument rpar,
ipar	a vector with integer values passed to the dll-function func and jacfunc present in the DLL, via function argument ipar,
nout	the number of output variables.
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.
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 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 package vignette "compiledCode".
fcontrol	A list of control parameters for the forcing functions. See package vignette "compiledCode".

Details

This function is meant to help developing FORTRAN or C models that are to be used to solve ordinary differential equations (ODE) in packages deSolve and/or rootSolve.

Value

a list containing:

dy	the rate of change estimated by the function,
var	the ordinary output variables of the function.

Author(s)

Karline Soetaert <[email protected]>

See Also

ode for a general interface to most of the ODE solvers

Examples

## ==========================================================================
## ex. 1
## ccl4model
## ==========================================================================
## Parameter values and initial conditions
## see example(ccl4model) for a more comprehensive implementation

Parms <- c(0.182, 4.0, 4.0, 0.08, 0.04, 0.74, 0.05, 0.15, 0.32,
        16.17, 281.48, 13.3, 16.17, 5.487, 153.8, 0.04321671,
        0.4027255, 1000, 0.02, 1.0, 3.8)

yini <- c(AI = 21, AAM = 0, AT = 0, AF = 0, AL = 0, CLT = 0,  AM = 0)

## the rate of change
DLLfunc(y = yini, dllname = "deSolve", func = "derivsccl4",
        initfunc = "initccl4", parms = Parms, times = 1,
        nout = 3, outnames = c("DOSE", "MASS", "CP")  )

## ==========================================================================
## ex. 2
## SCOC model
## ==========================================================================

## Forcing function "data"
Flux  <- matrix(ncol = 2, byrow = TRUE, data = c(1, 0.654, 2, 0.167))
parms <- c(k = 0.01)
Yini  <- 60

DLLfunc(y=Yini, times=1, func = "scocder",
    parms = parms, dllname = "deSolve",
    initforc = "scocforc",  forcings = Flux,
    initfunc = "scocpar", nout = 2,
    outnames = c("Mineralisation","Depo"))
## correct value = dy = flux - k * y = 0.654 - 0.01 * 60

DLLfunc(y = Yini, times = 2, func = "scocder",
    parms = parms, dllname = "deSolve",
    initforc = "scocforc",  forcings = Flux,
    initfunc = "scocpar", nout = 2,
    outnames = c("Mineralisation", "Depo"))

---

Evaluates a Residual Derivative Function Represented in a DLL

Description

Calls a residual function, $F(t,y,y')$ of a DAE system (differential algebraic equations) defined in a compiled language as a DLL.

To be used for testing the implementation of DAE problems in compiled code

Usage

DLLres(res, times, y, dy, parms, dllname,
  initfunc = dllname, rpar = NULL, ipar = NULL, nout = 0,
  outnames = NULL, forcings = NULL, initforc = NULL,
  fcontrol = NULL)

Arguments

res	the name of the function in the dynamically loaded shared library,
times	first value = the time at which the function needs to be evaluated,
y	the values of the dependent variables for which the function needs to be evaluated,
dy	the derivative of the values of the dependent variables for which the function needs to be evaluated,
parms	the parameters that are passed to the initialiser function,
dllname	a string giving the name of the shared library (without extension) that contains the compiled function or subroutine definitions referred to in func,
initfunc	if not NULL, the name of the initialisation function (which initialises values of parameters), as provided in ‘dllname’. See details,
rpar	a vector with double precision values passed to the DLL-function func and jacfunc present in the DLL, via argument rpar,
ipar	a vector with integer values passed to the DLL-function func and jacfunc present in the DLL, via function argument ipar,
nout	the number of output variables.
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.
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 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 package vignette "compiledCode".
fcontrol	A list of control parameters for the forcing functions. See package vignette "compiledCode".

Details

This function is meant to help developing FORTRAN or C models that are to be used to solve differential algebraic equations (DAE) in package deSolve.

Value

a list containing:

res	the residual of derivative estimated by the function
var	the ordinary output variables of the function

Author(s)

Karline Soetaert <[email protected]>

See Also

daspk to solve DAE problems

Examples

## =========================================================================
## Residuals from the daspk chemical model, production a forcing function
## =========================================================================
## Parameter values and initial conditions
## see example(daspk) for a more comprehensive implementation

pars <- c(K = 1, ka = 1e6, r = 1)

## Initial conc; D is in equilibrium with A,B
y <- c(A = 2, B = 3, D = 2 * 3/pars["K"])

## Initial rate of change
dy <- c(dA = 0, dB = 0, dD = 0)

## production increases with time
prod <- matrix(ncol = 2,
          data = c(seq(0, 100, by = 10), seq(0.1, 0.5, len = 11)))

DLLres(y = y, dy = dy, times = 5, res = "chemres",
          dllname = "deSolve", initfunc = "initparms",
          initforc = "initforcs", parms = pars, forcings = prod,
          nout = 2, outnames = c("CONC", "Prod"))

---

Implementing Events and Roots in Differential Equation Models.

Description

An event occurs when the value of a state variable is suddenly changed, e.g. because a value is added, subtracted, or multiplied. The integration routines cannot deal easily with such state variable changes. Typically these events occur only at specific times. In deSolve, events can be imposed by means of an input data.frame, that specifies at which time and how a certain state variable is altered, or via an event function.

Roots occur when a root function becomes zero. By default when a root is found, the simulation either stops (no event), or triggers an event.

Details

The events are specified by means of argument events passed to the integration routines.

events should be a list that contains one of the following:

func:

an R-function or the name of a function in compiled code that specifies the event,

data:

a data.frame that specifies the state variables, times, values and types of the events. Note that the event times must also be part of the integration output times, else the event will not take place. As from version 1.9.1, this is checked by the solver, and a warning message is produced if event times are missing in times; see also cleanEventTimes for utility functions to check and solve such issues.

time:

when events are specified by an event function: the times at which the events take place. Note that these event times must also be part of the integration output times exactly, else the event would not take place. As from version 1.9.1 this is checked by the solver, and an error message produced if event times are missing in times; see also cleanEventTimes for utility functions to check and solve such issues.

root:

when events are specified by a function and triggered by a root, this logical should be set equal to TRUE

terminalroot:

when events are triggered by a root, the default is that the simulation continues after the event is executed. In terminalroot, we can specify which roots should terminate the simulation.

maxroot:

when root = TRUE, the maximal number of times at with a root is found and that are kept; defaults to 100. If the number of roots > maxroot, then only the first maxroot will be outputted.

ties:

if events, as specified by a data.frame are "ordered", set to "ordered", the default is "notordered". This will save some computational time.

In case the events are specified by means of an R function (argument events$func), it must be defined as: 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 events$func. parms is a vector or list of parameters; ... (optional) are any other arguments passed to the function via the call to the integration method. The event function should return the y-values (some of which modified), as a vector.

If events$func is a string, this indicates that the events are specified by a function in compiled code. This function has as arguments, the number of state variables, the time, and the state variable vector. See package vignette "compiledCode" for more details.

In case events are specified by an R-function, this requires either: input of the time of the events, a vector as defined in events$time OR the specification of a root function. In the latter case, the model must be solved with an integration routine with root-finding capability

The root function itself should be specified with argument rootfunc. In this case, the integrator is informed that the simulation it to be continued after a root is found by setting events$root equal to TRUE.

If the events are specified by a data frame (argument events$data), this should contain the following columns (and in that order):

var:

the state variable name or number that is affected by the event

time:

the time at which the event is to take place; the solvers will check if the time is embraced by the simulation time

value:

the value, magnitude of the event

method:

which event is to take place; should be one of ("replace", "add", "multiply"); also allowed is to specify the number (1 = replace, 2 = add, 3 = multiply)

For instance, the following line

"v1" 10 2 "add"

will cause the value 2 to be added to a state variable, called "v1" at time = 10.

From deSolve version 1.9.1 the following routines have root-finding capability: lsoda, lsode, lsodes, and radau. For the first 3 integration methods, the root finding algorithm is based on the algorithm in solver LSODAR, and is implemented in FORTRAN. For radau, the root solving algorithm is written in C-code, and it works slightly different. Thus, some problems involving roots may be more efficiently solved with either lsoda, lsode, or lsodes, while other problems are more efficiently solved with radau.

If a root function is defined, but not an event function, then by default the solver will stop at a root. If this is not desirable, e.g. because we want to record the position of many roots, then a dummy "event" function can be defined which returns the values of the state variables - unaltered.

If roots and events are combined, and roots are found, then the output will have attribute troot which will contain the times at which a root was found (and the event trigerred). There will be at most events$maxroot such values. The default is 100.

See two last examples; also see example of ccl4model.

Author(s)

Karline Soetaert

See Also

forcings, for how to implement forcing functions.

lsodar, for more examples of roots

Examples

## =============================================================================
## 1. EVENTS in a data.frame
## =============================================================================

## derivative function: derivatives set to 0
derivs <- function(t, var, parms) {
  list(dvar = rep(0, 2))
}

yini <- c(v1 = 1, v2 = 2)
times <- seq(0, 10, by = 0.1)

eventdat <- data.frame(var = c("v1", "v2", "v2", "v1"),
                       time = c(1, 1, 5, 9) ,
                       value = c(1, 2, 3, 4),
                       method = c("add", "mult", "rep", "add"))
eventdat

out <- vode(func = derivs, y = yini, times = times, parms = NULL,
            events = list(data = eventdat))
plot(out)

##
eventdat <- data.frame(var = c(rep("v1", 10), rep("v2", 10)),
                       time = c(1:10, 1:10),
                       value = runif(20),
                       method = rep("add", 20))
eventdat

out <- ode(func = derivs, y = yini, times = times, parms = NULL,
           events = list(data = eventdat))

plot(out)

## =============================================================================
## 2. EVENTS in a function
## =============================================================================

## derivative function: rate of change v1 = 0, v2 reduced at first-order rate
derivs <- function(t, var, parms) {
   list(c(0, -0.5 * var[2]))
}


# events: add 1 to v1, multiply v2 with random number
eventfun <- function(t, y, parms){
  with (as.list(y),{
    v1 <- v1 + 1
    v2 <- 5 * runif(1)
    return(c(v1, v2))
  })
}

yini <- c(v1 = 1, v2 = 2)
times <- seq(0, 10, by = 0.1)

out <- ode(func = derivs, y = yini, times = times, parms = NULL,
           events = list(func = eventfun, time = c(1:9, 2.2, 2.4)) )
plot(out, type = "l")

## =============================================================================
## 3. EVENTS triggered by a root function
## =============================================================================

## derivative: simple first-order decay
derivs <- function(t, y, pars) {
  return(list(-0.1 * y))
}

## event triggered if state variable = 0.5
rootfun <- function (t, y, pars) {
  return(y - 0.5)
}

## sets state variable = 1
eventfun <- function(t, y, pars) {
  return(y = 1)
}

yini <- 2
times <- seq(0, 100, 0.1)

## uses ode to solve; root = TRUE specifies that the event is
## triggered by a root.
out <- ode(times = times, y = yini, func = derivs, parms = NULL,
           events = list(func = eventfun, root = TRUE),
           rootfun = rootfun)

plot(out, type = "l")

## time of the root:
troot <- attributes(out)$troot
points(troot, rep(0.5, length(troot)))


## =============================================================================
## 4. More ROOT examples: Rotation function
## =============================================================================
Rotate <- function(t, x, p )
  list(c( x[2],
         -x[1]  ))

## Root = when second state variable = 0
rootfun <- function(t, x, p) x[2]

## "event" returns state variables unchanged
eventfun <- function(t, x, p) x
times <- seq(from = 0, to = 15, by = 0.1)

## 1. No event: stops at first root
out1 <- ode(func = Rotate, y = c(5, 5), parms = 0,
           times = times, rootfun = rootfun)
tail(out1)

## 2. Continues till end of times and records the roots
out <- ode(func = Rotate, y = c(5, 5), parms = 0,
           times = times, rootfun = rootfun,
           events = list(func = eventfun, root = TRUE) )

plot(out)
troot <- attributes(out)$troot  # time of roots
points(troot,rep(0, length (troot)))

## Multiple roots:  either one of the state variables = 0
root2 <- function(t, x, p) x

out2 <- ode(func = Rotate, y = c(5, 5), parms = 0,
           times = times, rootfun = root2,
           events = list(func = eventfun, root = TRUE) )

plot(out2, which = 2)
troot <- attributes(out2)$troot
indroot <- attributes(out2)$indroot  # which root was found
points(troot, rep(0, length (troot)), col = indroot, pch = 18, cex = 2)

## Multiple roots and stop at first time root 1.
out3 <- ode(func = Rotate, y = c(5, 5), parms = 0,
      times = times, rootfun = root2,
      events = list(func = eventfun, root = TRUE, terminalroot = 1))


## =============================================================================
## 5. Stop at 5th root - only works with radau.
## =============================================================================
Rotate <- function(t, x, p )
  list(c( x[2],
         -x[1],
         0  ))

## Root = when second state variable = 0
root3  <- function(t, x, p)  c(x[2], x[3] - 5)
event3 <- function (t, x, p) c(x[1:2], x[3]+1)
times <- seq(0, 15, 0.1)
out3 <- ode(func = Rotate, y = c(x1 = 5, x2 = 5, nroot = 0),
      parms = 0, method = "radau",
      times = times, rootfun = root3,
      events = list(func = event3, root = TRUE, terminalroot = 2))
plot(out3)
attributes(out3)[c("troot", "nroot", "indroot")]

## =============================================================================
## 6 Event in R-code, model function in compiled code - based on vode example
## =============================================================================

times <- 1:365
Flux <- cbind(times, sin(pi*times/365)^2) # forcing function

# run without events
out <- ode(y = c(C = 1), times, func = "scocder", parms = c(k=0.01),
   dllname = "deSolve", initforc = "scocforc", forcings = Flux,
   initfunc = "scocpar", nout = 2, outnames = c("Mineralisation", "Depo"))

# Event halves the concentration
EventMin <- function(t, y , p) y/2

out2 <- ode(y = c(C = 1), times, func = "scocder", parms = c(k=0.01),
   dllname = "deSolve", initforc = "scocforc", forcings = Flux,
   initfunc = "scocpar", nout = 2, outnames = c("Mineralisation", "Depo"),
   events = list (func = EventMin, time = c(50.1, 200, 210.5)))

plot(out, out2)

---

Passing Forcing Functions to Models Written in R or Compiled Code.

Description

A forcing function is an external variable that is essential to the model, but not explicitly modeled. Rather, it is imposed as a time-series. Thus, if a model uses forcing variables, their value at each time point needs to be estimated by interpolation of a data series.

Details

The forcing functions are imposed as a data series, that contains the values of the forcings at specified times.

Models may be defined in compiled C or FORTRAN code, as well as in R.

If the model is defined in R code, it is most efficient to:

1. define a function that performs the linear interpolation, using R's approxfun. It is generally recommended to use rule = 2, such as to allow extrapolation outside of the time interval, especially when using the Livermore solvers, as these may exceed the last time point.

2. call this function within the model's derivative function, to interpolate at the current timestep.

See first example.

If the models are defined in compiled C or FORTRAN code, it is possible to use deSolves forcing function update algorithm. This is the compiled-code equivalent of approxfun or approx.

In this case:
1. the forcing function data series is provided by means of argument forcings,

2. initforc is the name of the forcing function initialisation function, as provided in ‘dllname’, while

3. fcontrol is a list used to finetune how the forcing update should be performed.

The fcontrol argument is a list that can supply any of the following components (conform the definitions in the approxfun function):

method

specifies the interpolation method to be used. Choices are "linear" or "constant",

rule

an integer describing how interpolation is to take place outside the interval [min(times), max(times)]. If rule is 1 then an error will be triggered and the calculation will stop if times extends the interval of the forcing function data set. If it is 2, the default, the value at the closest data extreme is used, a warning will be printed if verbose is TRUE,

Note that the default differs from the approx default.

f

For method = "constant" a number between 0 and 1 inclusive, indicating a compromise between left- and right-continuous step functions. If y0 and y1 are the values to the left and right of the point then the value is y0 * (1 - f) + y1 * f so that f = 0 is right-continuous and f = 1 is left-continuous,

ties

Handling of tied times values. Either a function with a single vector argument returning a single number result or the string "ordered".

Note that the default is "ordered", hence the existence of ties will NOT be investigated; in the C code this will mean that -if ties exist, the first value will be used; if the dataset is not ordered, then nonsense will be produced.

Alternative values for ties are mean, min etc

The defaults are:

fcontrol = list(method = "linear", rule = 2, f = 0, ties = "ordered")

Note that only ONE specification is allowed, even if there is more than one forcing function data set.

More information about models defined in compiled code is in the package vignette ("compiledCode").

Note

How to write compiled code is described in package vignette "compiledCode", which should be referred to for details.

This vignette also contains examples on how to pass forcing functions.

Author(s)

Karline Soetaert,

Thomas Petzoldt,

R. Woodrow Setzer

See Also

approx or approxfun, the R function,

events for how to implement events.

Examples

## =============================================================================
## FORCING FUNCTION: The sediment oxygen consumption example - R-code:
## =============================================================================

## Forcing function data
Flux <- matrix(ncol=2,byrow=TRUE,data=c(
  1, 0.654, 11, 0.167,   21, 0.060, 41, 0.070, 73,0.277, 83,0.186,
  93,0.140,103, 0.255,  113, 0.231,123, 0.309,133,1.127,143,1.923,
  153,1.091,163,1.001,  173, 1.691,183, 1.404,194,1.226,204,0.767,
  214, 0.893,224,0.737, 234,0.772,244, 0.726,254,0.624,264,0.439,
  274,0.168,284 ,0.280, 294,0.202,304, 0.193,315,0.286,325,0.599,
  335, 1.889,345, 0.996,355,0.681,365,1.135))

parms <- c(k=0.01)

times <- 1:365

## the model
sediment <- function( t, O2, k) 
  list (c(Depo(t) - k * O2), depo = Depo(t))

# the forcing functions; rule = 2 avoids NaNs in interpolation
Depo <- approxfun(x = Flux[,1], y = Flux[,2], method = "linear", rule = 2)

Out <- ode(times = times, func = sediment, y = c(O2 = 63), parms = parms)
  
## same forcing functions, now constant interpolation
Depo <- approxfun(x = Flux[,1], y = Flux[,2], method = "constant",
  f = 0.5, rule = 2)

Out2 <- ode(times = times, func = sediment, y = c(O2 = 63), parms = parms)

mf <- par(mfrow = c(2, 1))
plot (Out, which = "depo", type = "l", lwd = 2, mfrow = NULL)
lines(Out2[,"time"], Out2[,"depo"], col = "red", lwd = 2)

plot (Out, which = "O2", type = "l", lwd = 2, mfrow = NULL)
lines(Out2[,"time"], Out2[,"O2"], col = "red", lwd = 2)

## =============================================================================
## SCOC is the same model, as implemented in FORTRAN
## =============================================================================

out<- SCOC(times, parms = parms, Flux = Flux)

plot(out[,"time"], out[,"Depo"], type = "l", col = "red")
lines(out[,"time"], out[,"Mineralisation"], col = "blue")

## Constant interpolation of forcing function - left side of interval
fcontrol <- list(method = "constant")

out2 <- SCOC(times, parms = parms, Flux = Flux, fcontrol = fcontrol)

plot(out2[,"time"], out2[,"Depo"], type = "l", col = "red")
lines(out2[,"time"], out2[,"Mineralisation"], col = "blue")


## Not run: 
## =============================================================================
## show examples (see respective help pages for details)
## =============================================================================

example(aquaphy)

## show package vignette with tutorial about how to use compiled models
## + source code of the vignette
## + directory with C and FORTRAN sources
vignette("compiledCode")
edit(vignette("compiledCode"))
browseURL(paste(system.file(package = "deSolve"), "/doc", sep = ""))

## End(Not run)

---

Solver for Ordinary Differential Equations (ODE), Switching Automatically Between Stiff and Non-stiff Methods

Description

Solving initial value problems for stiff or non-stiff systems of first-order ordinary differential equations (ODEs).

The R function lsoda provides an interface to the FORTRAN ODE solver of the same name, written by Linda R. Petzold and Alan C. Hindmarsh.

The system of ODE's is written as an R function (which may, of course, use .C, .Fortran, .Call, etc., to call foreign code) or be defined in compiled code that has been dynamically loaded. A vector of parameters is passed to the ODEs, so the solver may be used as part of a modeling package for ODEs, or for parameter estimation using any appropriate modeling tool for non-linear models in R such as optim, nls, nlm or nlme

lsoda differs from the other integrators (except lsodar) in that it switches automatically between stiff and nonstiff methods. This means that the user does not have to determine whether the problem is stiff or not, and the solver will automatically choose the appropriate method. It always starts with the nonstiff method.

Usage

lsoda(y, times, func, parms, rtol = 1e-6, atol = 1e-6,
  jacfunc = NULL, jactype = "fullint", rootfunc = NULL,
  verbose = FALSE, nroot = 0, tcrit = NULL,
  hmin = 0, hmax = NULL, hini = 0, ynames = TRUE,
  maxordn = 12, maxords = 5, 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,...)

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.
times	times at which explicit estimates for y are desired. The first value in times must be the initial time.
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, or a list of symbols returned by checkDLL. 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. If func is a string, then dllname must give the name of the shared library (without extension) which must be loaded before lsoda() is called. See package vignette "compiledCode" for more details. func can also be a list of symbols returned by checkDLL to avoid overhead from repeated internal calls to this function. It is an experimental feature for special situations, when a small compiled model with a low number of integration steps is repeatedly called. It is currently only available for the lsoda solver, see example.
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 $\partial\dot{y}_i/\partial y_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 $\partial\dot{y}/\partial y$, 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). 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.
rootfunc	if not NULL, an R function that computes the function whose root has to be estimated or a string giving the name of a function or subroutine in ‘dllname’ that computes the root function. The R calling sequence for rootfunc is identical to that of func. rootfunc should return a vector with the function values whose root is sought. When rootfunc is provided, then lsodar will be called.
verbose	if TRUE: full output to the screen, e.g. will print the diagnostiscs of the integration - see details.
nroot	only used if ‘dllname’ is specified: the number of constraint functions whose roots are desired during the integration; if rootfunc is an R-function, the solver estimates the number of roots.
tcrit	if not NULL, then lsoda cannot integrate past tcrit. The FORTRAN routine lsoda 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 large models.
maxordn	the maximum order to be allowed in case the method is non-stiff. Should be <= 12. Reduce maxord to save storage space.
maxords	the maximum order to be allowed in case the method is stiff. Should be <= 5. 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. See forcings or vignette compiledCode.
events	A matrix or data frame that specifies events, i.e. when the value of a state variable is suddenly changed. See events for more information.
lags	A list that specifies timelags, i.e. the number of steps that has to be kept. To be used for delay differential equations. See timelags, dede for more information.
...	additional arguments passed to func and jacfunc allowing this to be a generic function.

Details

All the hard work is done by the FORTRAN subroutine lsoda, 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 12 November 2003 version of lsoda, from Netlib.

lsoda switches automatically between stiff and nonstiff methods. This means that the user does not have to determine whether the problem is stiff or not, and the solver will automatically choose the appropriate method. It always starts with the nonstiff method.

The form of the Jacobian can be specified by jactype which can take the following values:

"fullint"

a full Jacobian, calculated internally by lsoda, the default,

"fullusr"

a full Jacobian, specified by user function jacfunc,

"bandusr"

a banded Jacobian, specified by user function jacfunc the size of the bands specified by bandup and banddown,

"bandint"

banded Jacobian, calculated by lsoda; the size of the bands specified by bandup and banddown.

If jactype = "fullusr" or "bandusr" then the user must supply a subroutine jacfunc.

The following description of error control is adapted from the documentation of the lsoda source code (input arguments rtol and atol, above):

The input parameters rtol, and atol determine the error control performed by the solver. The solver will control the vector e of estimated local errors in y, according to an inequality of the form max-norm of ( e/ewt ) $\leq$ 1, where ewt is a vector of positive error weights. The values of rtol and atol should all be non-negative. The form of ewt is:

$\mathbf{rtol} \times \mathrm{abs}(\mathbf{y}) + \mathbf{atol}$

where multiplication of two vectors is element-by-element.

If the request for precision exceeds the capabilities of the machine, the FORTRAN subroutine lsoda will return an error code; under some circumstances, the R function lsoda will attempt a reasonable reduction of precision in order to get an answer. It will write a warning if it does so.

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.

More information about models defined in compiled code is in the package vignette ("compiledCode"); information about linking forcing functions to compiled code is in forcings.

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

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

Note

The ‘demo’ directory contains some examples of using gnls to estimate parameters in a dynamic model.

Author(s)

R. Woodrow Setzer <[email protected]>

References

Hindmarsh, Alan C. (1983) ODEPACK, A Systematized Collection of ODE Solvers; in p.55–64 of Stepleman, R.W. et al.[ed.] (1983) Scientific Computing, North-Holland, Amsterdam.

Petzold, Linda R. (1983) Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations. Siam J. Sci. Stat. Comput. 4, 136–148. doi:10.1137/0904010

Netlib: https://netlib.org

See Also

• rk, rkMethod, rk4 and euler for Runge-Kutta integrators.

• lsode, which can also find a root

• 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

## =======================================================================
## Example 1:
##   A simple resource limited Lotka-Volterra-Model
##
## Note:
## 1. parameter and state variable names made
##    accessible via "with" function
## 2. function sigimp accessible through lexical scoping
##    (see also ode and rk examples)
## =======================================================================

SPCmod <- function(t, x, parms) {
  with(as.list(c(parms, x)), {
    import <- sigimp(t)
    dS <- import - b*S*P + g*C     #substrate
    dP <- c*S*P  - d*C*P           #producer
    dC <- e*P*C  - f*C             #consumer
    res <- c(dS, dP, dC)
    list(res)
  })
}

## Parameters
parms  <- c(b = 0.0, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)

## vector of timesteps
times  <- seq(0, 100, length = 101)

## external signal with rectangle impulse
signal <- as.data.frame(list(times = times,
                            import = rep(0,length(times))))

signal$import[signal$times >= 10 & signal$times <= 11] <- 0.2

sigimp <- approxfun(signal$times, signal$import, rule = 2)


## Start values for steady state
y <- xstart <- c(S = 1, P = 1, C = 1)

## Solving
out <-  lsoda(xstart, times, SPCmod, parms)

## Plotting
mf <- par("mfrow")
plot(out, main = c("substrate", "producer", "consumer"))
plot(out[,"P"], out[,"C"], type = "l", xlab = "producer", ylab = "consumer")
par(mfrow = mf)

## =======================================================================
## Example 2:
##  from lsoda source code
## =======================================================================

## names makes this easier to read, but may slow down execution.
parms   <- c(k1 = 0.04, k2 = 1e4, k3 = 3e7)
my.atol <- c(1e-6,  1e-10,  1e-6)
times   <- c(0,4 * 10^(-1:10))

lsexamp <- function(t, y, p) {
  yd1 <- -p["k1"] * y[1] + p["k2"] * y[2]*y[3]
  yd3 <- p["k3"] * y[2]^2
  list(c(yd1, -yd1-yd3, yd3), c(massbalance = sum(y)))
}

exampjac <- function(t, y, p) {
  matrix(c(-p["k1"],	   p["k1"],       0,

            p["k2"]*y[3],
          - p["k2"]*y[3] - 2*p["k3"]*y[2],
                           2*p["k3"]*y[2],

          p["k2"]*y[2],  -p["k2"]*y[2],  0
          ), 3, 3)
}


## measure speed (here and below)
system.time(
  out <- lsoda(c(1, 0, 0), times, lsexamp, parms, rtol = 1e-4,
           atol = my.atol, hmax = Inf)
)
out

## This is what the authors of lsoda got for the example:

## the output of this program (on a cdc-7600 in single precision)
## is as follows..
##
## at t =  4.0000e-01   y =  9.851712e-01  3.386380e-05  1.479493e-02
## at t =  4.0000e+00   y =  9.055333e-01  2.240655e-05  9.444430e-02
## at t =  4.0000e+01   y =  7.158403e-01  9.186334e-06  2.841505e-01
## at t =  4.0000e+02   y =  4.505250e-01  3.222964e-06  5.494717e-01
## at t =  4.0000e+03   y =  1.831975e-01  8.941774e-07  8.168016e-01
## at t =  4.0000e+04   y =  3.898730e-02  1.621940e-07  9.610125e-01
## at t =  4.0000e+05   y =  4.936363e-03  1.984221e-08  9.950636e-01
## at t =  4.0000e+06   y =  5.161831e-04  2.065786e-09  9.994838e-01
## at t =  4.0000e+07   y =  5.179817e-05  2.072032e-10  9.999482e-01
## at t =  4.0000e+08   y =  5.283401e-06  2.113371e-11  9.999947e-01
## at t =  4.0000e+09   y =  4.659031e-07  1.863613e-12  9.999995e-01
## at t =  4.0000e+10   y =  1.404280e-08  5.617126e-14  1.000000e+00

## Using the analytic Jacobian speeds up execution a little :

system.time(
  outJ <- lsoda(c(1, 0, 0), times, lsexamp, parms, rtol = 1e-4,
            atol = my.atol, jacfunc = exampjac, jactype = "fullusr", hmax = Inf)
)

all.equal(as.data.frame(out), as.data.frame(outJ)) # TRUE
diagnostics(out)
diagnostics(outJ) # shows what lsoda did internally

---

Solver for Ordinary Differential Equations (ODE), Switching Automatically Between Stiff and Non-stiff Methods and With Root Finding

Description

Solving initial value problems for stiff or non-stiff systems of first-order ordinary differential equations (ODEs) and including root-finding.

The R function lsodar provides an interface to the FORTRAN ODE solver of the same name, written by Alan C. Hindmarsh and Linda R. Petzold.

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

lsodar differs from lsode in two respects.

• It switches automatically between stiff and nonstiff methods (similar as lsoda).

• It finds the root of at least one of a set of constraint functions g(i) of the independent and dependent variables.

Two uses of lsodar are:

• To stop the simulation when a certain condition is met

• To trigger events, i.e. sudden changes in one of the state variables when a certain condition is met.

when a particular condition is met.

Usage

lsodar(y, times, func, parms, rtol = 1e-6, atol = 1e-6, 
  jacfunc = NULL, jactype = "fullint", rootfunc = NULL,
  verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0,
  hmax = NULL, hini = 0, ynames = TRUE, maxordn = 12,
  maxords = 5, 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, ...)

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.
times	times at which explicit estimates for y are desired. The first value in times must be the initial time.
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. If func is a string, then dllname must give the name of the shared library (without extension) which must be loaded before lsodar() 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 $\partial\dot{y}_i/\partial y_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 $\partial\dot{y}/\partial y$, 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). 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.
rootfunc	if not NULL, an R function that computes the function whose root has to be estimated or a string giving the name of a function or subroutine in ‘dllname’ that computes the root function. The R calling sequence for rootfunc is identical to that of func. rootfunc should return a vector with the function values whose root is sought.
verbose	a logical value that, when TRUE, will print the diagnostiscs of the integration - see details.
nroot	only used if ‘dllname’ is specified: the number of constraint functions whose roots are desired during the integration; if rootfunc is an R-function, the solver estimates the number of roots.
tcrit	if not NULL, then lsodar cannot integrate past tcrit. The FORTRAN routine lsodar 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 large models.
maxordn	the maximum order to be allowed in case the method is non-stiff. Should be <= 12. Reduce maxord to save storage space.
maxords	the maximum order to be allowed in case the method is stiff. Should be <= 5. 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. See forcings or vignette compiledCode.
events	A matrix or data frame that specifies events, i.e. when the value of a state variable is suddenly changed. See events for more information.
lags	A list that specifies timelags, i.e. the number of steps that has to be kept. To be used for delay differential equations. See timelags, dede for more information.
...	additional arguments passed to func and jacfunc allowing this to be a generic function.

Details

The work is done by the FORTRAN subroutine lsodar, 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 lsodar, from Netlib.

lsodar switches automatically between stiff and nonstiff methods (similar as lsoda). This means that the user does not have to determine whether the problem is stiff or not, and the solver will automatically choose the appropriate method. It always starts with the nonstiff method.

lsodar 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, lsodar may return false roots, or return the same root at two or more nearly equal values of time.

The form of the Jacobian can be specified by jactype which can take the following values:

jactype = "fullint":

a full Jacobian, calculated internally by lsodar, the default,

jactype = "fullusr":

a full Jacobian, specified by user function jacfunc,

jactype = "bandusr":

a banded Jacobian, specified by user function jacfunc; the size of the bands specified by bandup and banddown,

jactype = "bandint":

banded Jacobian, calculated by lsodar; the size of the bands specified by bandup and banddown.

If jactype = "fullusr" or "bandusr" then the user must supply a subroutine jacfunc.

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

The output will have the attribute iroot, if a root was found iroot is a vector, its length equal to the number of constraint functions it will have a value of 1 for the constraint function whose root that has been found and 0 otherwise.

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.

More information about models defined in compiled code is in the package vignette ("compiledCode"); information about linking forcing functions to compiled code is in forcings.

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

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

If a root has been found, the output will have the attribute iroot, an integer indicating which root has been found.

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.

Linda R. Petzold, Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations, Siam J. Sci. Stat. Comput. 4 (1983), pp. 136-148. doi:10.1137/0904010

Kathie L. Hiebert and Lawrence F. Shampine, Implicitly Defined Output Points for Solutions of ODEs, Sandia Report SAND80-0180, February 1980.

Netlib: https://netlib.org

See Also

• roots for more examples on roots and events

• rk, rkMethod, rk4 and euler for Runge-Kutta integrators.

• lsoda, lsode, lsodes, 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

## =======================================================================
## Example 1:
##   from lsodar source code
## =======================================================================

Fun <- function (t, y, parms) {
  ydot <- vector(len = 3)
  ydot[1] <- -.04*y[1] + 1.e4*y[2]*y[3]
  ydot[3] <- 3.e7*y[2]*y[2]
  ydot[2] <- -ydot[1] - ydot[3]
  return(list(ydot, ytot = sum(y)))
}

rootFun <- function (t, y, parms) {
  yroot <- vector(len = 2)
  yroot[1] <- y[1] - 1.e-4
  yroot[2] <- y[3] - 1.e-2
  return(yroot)
}

y     <- c(1, 0, 0)
times <- c(0, 0.4*10^(0:8))

out   <- lsodar(y = y, times = times, fun = Fun, rootfun = rootFun,
                rtol = 1e-4, atol = c(1e-6, 1e-10, 1e-6), parms = NULL)
print(paste("root is found for eqn", which(attributes(out)$iroot == 1)))
print(out[nrow(out),])

diagnostics(out)
  
## =======================================================================
## Example 2:
##   using lsodar to estimate steady-state conditions
## =======================================================================

## Bacteria (Bac) are growing on a substrate (Sub)
model <- function(t, state, pars) {
  with (as.list(c(state, pars)), {
    ##        substrate uptake             death     respiration
    dBact <-  gmax*eff*Sub/(Sub+ks)*Bact - dB*Bact - rB*Bact
    dSub  <- -gmax    *Sub/(Sub+ks)*Bact + dB*Bact            + input

    return(list(c(dBact,dSub)))
  })
}

## root is the condition where sum of |rates of change|
## is very small

rootfun <- function (t, state, pars) {
  dstate <- unlist(model(t, state, pars)) # rate of change vector
  return(sum(abs(dstate)) - 1e-10)
}

pars <- list(Bini = 0.1, Sini = 100, gmax = 0.5, eff = 0.5,
             ks = 0.5, rB = 0.01, dB = 0.01, input = 0.1)

tout    <- c(0, 1e10)
state   <- c(Bact = pars$Bini, Sub = pars$Sini)
out     <- lsodar(state, tout, model, pars, rootfun = rootfun)
print(out)


## =======================================================================
## Example 3:
##   using lsodar to trigger an event
## =======================================================================

## a state variable is decaying at a first-order rate. 
## when it reaches the value 0.1, a random amount is added.

derivfun <- function (t,y,parms)
  list (-0.05 * y)

rootfun <- function (t,y,parms)
  return(y - 0.1) 

eventfun <- function(t,y,parms)
  return(y + runif(1))  

yini <- 0.8
times <- 0:200

out <- lsodar(func=derivfun, y = yini, times=times, 
  rootfunc = rootfun, events = list(func=eventfun, root = TRUE))

plot(out, type = "l", lwd = 2, main = "lsodar with event")

---