Version:	1.40
Title:	Solvers for Initial Value Problems of Differential Equations ('ODE', 'DAE', 'DDE')
Author:	Karline Soetaert [aut], Thomas Petzoldt [aut, cre], R. Woodrow Setzer [aut], Peter N. Brown [ctb] (files ddaspk.f, dvode.f, zvode.f), George D. Byrne [ctb] (files dvode.f, zvode.f), Ernst Hairer [ctb] (files radau5.f, radau5a), Alan C. Hindmarsh [ctb] (files ddaspk.f, dlsode.f, dvode.f, zvode.f, opdkmain.f, opdka1.f), Cleve Moler [ctb] (file dlinpck.f), Linda R. Petzold [ctb] (files ddaspk.f, dlsoda.f), Youcef Saad [ctb] (file dsparsk.f), Clement W. Ulrich [ctb] (file ddaspk.f)
Maintainer:	Thomas Petzoldt <thomas.petzoldt@tu-dresden.de>
Depends:	R (≥ 3.3.0)
Imports:	methods, graphics, grDevices, stats
Suggests:	scatterplot3d, FME
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 of 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.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
URL:	http://desolve.r-forge.r-project.org/
LazyData:	yes
NeedsCompilation:	yes
Packaged:	2023-11-27 22:26:41 UTC; thpe
Repository:	CRAN
Date/Publication:	2023-11-27 23:30:02 UTC

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:r-sig-dynamic-models@r-project.org

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

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 <karline.soetaert@nioz.nl>

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

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 <setzer.woodrow@epa.gov>

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

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:

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

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

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 <karline.soetaert@nioz.nl>

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

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 <karline.soetaert@nioz.nl>

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)

Internal deSolve Functions

Description

Internal deSolve functions, these are not to be called by the user.

Usage

timestep(prev = TRUE)

Arguments

Details

Function timestep is intended to return the current or next timestep of the integration. It works only under specific circumstances and should not be used by the end user.

Instead of this, please see the example below for a pure R solution.

See Also

diagnostics for information about the time steps used,
lagvalue and lagderiv that can be used for DDEs.

Examples

###################################################
### This example shows how to retrieve information 
### about the used time steps.
###################################################

## a function closure  ('lexical scoping')
modelClosure <- function(t0) {
  t.old <- t.act <- t0
  function(t, y, parms) {
    t.old  <<- t.act
    t.act  <<- t
    cat(t, "\t", t - t.old, "\n")
    with (as.list(c(y, parms)), {
      dP <- a * P      -  b * P * K
      dK <- b * P * K  -  c * K
      list(c(dP, dK))
    })
  }
}

model <- modelClosure(0) # initialization

parms <- c(a = 0.1, b = 0.1, c = 0.1)

y <- c(P = 1, K = 2)

out <- ode(y = y, func = model, times = c(0, 2),
 parms = parms, method = "lsoda")

ls() # prove that t.old and t.new are local within 'model' 

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

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

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

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:

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

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

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:

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

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

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 <setzer.woodrow@epa.gov>

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

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 <karline.soetaert@nioz.nl>

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")
  

Solver for Ordinary Differential Equations (ODE)

Description

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

dy/dt = f(t,y)

.

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

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

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

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

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

Usage

lsode(y, times, func, parms, rtol = 1e-6, atol = 1e-6,  
  jacfunc = NULL, jactype = "fullint", mf = NULL, rootfunc = NULL,
  verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0, hmax = NULL, 
  hini = 0, ynames = TRUE, maxord = NULL, bandup = NULL, banddown = NULL,
  maxsteps = 5000, dllname = NULL, initfunc = dllname,
  initpar = parms, rpar = NULL, ipar = NULL, nout = 0,
  outnames = NULL, forcings=NULL, initforc = NULL, 
  fcontrol=NULL, events=NULL, lags = NULL,...)

Arguments

Details

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

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

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

The options for jactype are

jactype = "fullint"

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

jactype = "fullusr"

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

jactype = "bandusr"

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

jactype = "bandint"

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

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

METH

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

MITER

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

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

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

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

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

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

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

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.

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

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

Value

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

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

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.

See Also

• rk,

• rk4 and euler for Runge-Kutta integrators.

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

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

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

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

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

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

diagnostics to print diagnostic messages.

Examples

## =======================================================================
## Example 1:
##   Various ways to solve the same model.
## =======================================================================

## the model, 5 state variables
f1 <- function  (t, y, parms) {
  ydot <- vector(len = 5)

  ydot[1] <-  0.1*y[1] -0.2*y[2]
  ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3]
  ydot[3] <-           -0.3*y[2] +0.1*y[3] -0.2*y[4]
  ydot[4] <-                     -0.3*y[3] +0.1*y[4] -0.2*y[5]
  ydot[5] <-                               -0.3*y[4] +0.1*y[5]

  return(list(ydot))
}

## the Jacobian, written as a full matrix
fulljac <- function  (t, y, parms) {
  jac <- matrix(nrow = 5, ncol = 5, byrow = TRUE,
                data = c(0.1, -0.2,  0  ,  0  ,  0  ,
                        -0.3,  0.1, -0.2,  0  ,  0  ,
                         0  , -0.3,  0.1, -0.2,  0  ,
                         0  ,  0  , -0.3,  0.1, -0.2,
                         0  ,  0  ,  0  , -0.3,  0.1))
  return(jac)
}

## the Jacobian, written in banded form
bandjac <- function  (t, y, parms) {
  jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE,
                data = c( 0  , -0.2, -0.2, -0.2, -0.2,
                          0.1,  0.1,  0.1,  0.1,  0.1,
                         -0.3, -0.3, -0.3, -0.3,    0))
  return(jac)
}

## initial conditions and output times
yini  <- 1:5
times <- 1:20

## default: stiff method, internally generated, full Jacobian
out   <- lsode(yini, times, f1, parms = 0, jactype = "fullint")

## stiff method, user-generated full Jacobian
out2  <- lsode(yini, times, f1, parms = 0, jactype = "fullusr",
              jacfunc = fulljac)

## stiff method, internally-generated banded Jacobian
## one nonzero band above (up) and below(down) the diagonal
out3  <- lsode(yini, times, f1, parms = 0, jactype = "bandint",
                              bandup = 1, banddown = 1)

## stiff method, user-generated banded Jacobian
out4  <- lsode(yini, times, f1, parms = 0, jactype = "bandusr",
              jacfunc = bandjac, bandup = 1, banddown = 1)

## non-stiff method
out5  <- lsode(yini, times, f1, parms = 0, mf = 10)

## =======================================================================
## Example 2:
##   diffusion on a 2-D grid
##   partially specified Jacobian
## =======================================================================

diffusion2D <- function(t, Y, par) {
   y <- matrix(nrow = n, ncol = n, data = Y)
   dY   <- r*y     # production

   ## diffusion in X-direction; boundaries = 0-concentration
   Flux <- -Dx * rbind(y[1,],(y[2:n,]-y[1:(n-1),]),-y[n,])/dx
   dY   <- dY - (Flux[2:(n+1),]-Flux[1:n,])/dx

   ## diffusion in Y-direction
   Flux <- -Dy * cbind(y[,1],(y[,2:n]-y[,1:(n-1)]),-y[,n])/dy
   dY    <- dY - (Flux[,2:(n+1)]-Flux[,1:n])/dy

   return(list(as.vector(dY)))
}

## parameters
dy    <- dx <- 1   # grid size
Dy    <- Dx <- 1   # diffusion coeff, X- and Y-direction
r     <- 0.025     # production rate
times <- c(0, 1)

n  <- 50
y  <- matrix(nrow = n, ncol = n, 0)

pa <- par(ask = FALSE)

## initial condition
for (i in 1:n) {
  for (j in 1:n) {
    dst <- (i - n/2)^2 + (j - n/2)^2
    y[i, j] <- max(0, 1 - 1/(n*n) * (dst - n)^2)
  }
}
filled.contour(y, color.palette = terrain.colors)

## =======================================================================
##   jacfunc need not be estimated exactly
##   a crude approximation, with a smaller bandwidth will do.
##   Here the half-bandwidth 1 is used, whereas the true
##   half-bandwidths are equal to n.
##   This corresponds to ignoring the y-direction coupling in the ODEs.
## =======================================================================

print(system.time(
  for (i in 1:20) {
    out  <-  lsode(func = diffusion2D, y = as.vector(y), times = times,
              parms = NULL, jactype = "bandint", bandup = 1, banddown = 1)

    filled.contour(matrix(nrow = n, ncol = n, out[2,-1]), zlim = c(0,1),
                  color.palette = terrain.colors, main = i)

    y <- out[2, -1]
  }
))
par(ask = pa)

Solver for Ordinary Differential Equations (ODE) With Sparse Jacobian

Description

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

dy/dt = f(t,y)

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

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

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

Usage

lsodes(y, times, func, parms, rtol = 1e-6, atol = 1e-6, 
  jacvec = NULL, sparsetype = "sparseint", nnz = NULL,
  inz = NULL,  rootfunc = NULL,
  verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0,
  hmax = NULL, hini = 0, ynames = TRUE, maxord = NULL,
  maxsteps = 5000, lrw = NULL, liw = NULL, 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

 
DLSODES- RWORK length is insufficient to proceed.                               
  Length needed is .ge. LENRW (=I1), exceeds LRW (=I2)                    
  In above message,  I1 =     27627   I2 =     25932 

Details

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

lsodes is applied for stiff problems, where the Jacobian has a sparse structure.

There are several choices depending on whether jacvec is specified and depending on the setting of sparsetype.

If function jacvec is present, then it should return the j-th column of the Jacobian matrix.

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

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

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

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

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

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

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

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

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

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 ‘doc/examples/dynload’ subdirectory of the deSolve package directory.

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

Value

A matrix of class deSolve with up to as many rows as elements in times and as many columns as elements in y plus the number of "global" values returned in the next elements of the return from func, plus and additional column for the time value. There will be a row for each element in times unless the FORTRAN routine ‘lsodes’ 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 <karline.soetaert@nioz.nl>

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.

S. C. Eisenstat, M. C. Gursky, M. H. Schultz, and A. H. Sherman, Yale Sparse Matrix Package: I. The Symmetric Codes, Int. J. Num. Meth. Eng., 18 (1982), pp. 1145-1151.

S. C. Eisenstat, M. C. Gursky, M. H. Schultz, and A. H. Sherman, Yale Sparse Matrix Package: II. The Nonsymmetric Codes, Research Report No. 114, Dept. of Computer Sciences, Yale University, 1977.

See Also

• rk,

• rk4 and euler for Runge-Kutta integrators.

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

## Various ways to solve the same model.

## =======================================================================
## The example from lsodes source code
## A chemical model
## =======================================================================

n  <- 12
y  <- rep(1, n)
dy <- rep(0, n)

times <- c(0, 0.1*(10^(0:4)))

rtol <- 1.0e-4
atol <- 1.0e-6

parms <- c(rk1  = 0.1,   rk2 = 10.0, rk3 = 50.0,  rk4 = 2.5,  rk5 = 0.1,
           rk6  = 10.0,  rk7 = 50.0, rk8 = 2.5,   rk9 = 50.0, rk10 = 5.0,
           rk11 = 50.0, rk12 = 50.0,rk13 = 50.0, rk14 = 30.0,
           rk15 = 100.0,rk16 = 2.5, rk17 = 100.0,rk18 = 2.5,
           rk19 = 50.0, rk20 = 50.0)

#
chemistry <- function (time, Y, pars) {
  with (as.list(pars), {
    dy[1] <- -rk1 *Y[1]
    dy[2] <-  rk1 *Y[1]        + rk11*rk14*Y[4]  + rk19*rk14*Y[5]  -
              rk3 *Y[2]*Y[3]   - rk15*Y[2]*Y[12] - rk2*Y[2]
    dy[3] <-  rk2 *Y[2]        - rk5 *Y[3]       - rk3*Y[2]*Y[3]   -
              rk7*Y[10]*Y[3]   + rk11*rk14*Y[4]   + rk12*rk14*Y[6]
    dy[4] <-  rk3 *Y[2]*Y[3]   - rk11*rk14*Y[4]  - rk4*Y[4]
    dy[5] <-  rk15*Y[2]*Y[12]  - rk19*rk14*Y[5]  - rk16*Y[5]
    dy[6] <-  rk7 *Y[10]*Y[3]  - rk12*rk14*Y[6]  - rk8*Y[6]
    dy[7] <-  rk17*Y[10]*Y[12] - rk20*rk14*Y[7]  - rk18*Y[7]
    dy[8] <-  rk9 *Y[10]       - rk13*rk14*Y[8]  - rk10*Y[8]
    dy[9] <-  rk4 *Y[4]        + rk16*Y[5]       + rk8*Y[6]         +
              rk18*Y[7]
    dy[10] <- rk5 *Y[3]        + rk12*rk14*Y[6]  + rk20*rk14*Y[7]   +
              rk13*rk14*Y[8]   - rk7 *Y[10]*Y[3] - rk17*Y[10]*Y[12] -
              rk6 *Y[10]       - rk9*Y[10]
    dy[11] <- rk10*Y[8]
    dy[12] <- rk6 *Y[10]       + rk19*rk14*Y[5]  + rk20*rk14*Y[7]   -
              rk15*Y[2]*Y[12]  - rk17*Y[10]*Y[12]
    return(list(dy))
  })
}

## =======================================================================
## application 1. lsodes estimates the structure of the Jacobian
##                and calculates the Jacobian by differences
## =======================================================================
out <- lsodes(func = chemistry, y = y, parms = parms, times = times,
              atol = atol, rtol = rtol, verbose = TRUE)

## =======================================================================
## application 2. the structure of the Jacobian is input
##                lsodes calculates the Jacobian by differences
##                this is not so efficient...
## =======================================================================

## elements of Jacobian that are not zero
nonzero <-  matrix(nc = 2, byrow = TRUE, data = c(
  1, 1,   2, 1,    # influence of sp1 on rate of change of others
  2, 2,   3, 2,   4, 2,   5, 2,  12, 2,
  2, 3,   3, 3,   4, 3,   6, 3,  10, 3,
  2, 4,   3, 4,   4, 4,   9, 4,  # d (dyi)/dy4
  2, 5,   5, 5,   9, 5,  12, 5,
  3, 6,   6, 6,   9, 6,  10, 6,
  7, 7,   9, 7,  10, 7,  12, 7,
  8, 8,  10, 8,  11, 8,
  3,10,   6,10,   7,10,   8,10,  10,10,  12,10,
  2,12,   5,12,   7,12,  10,12,  12,12)
)

## when run, the default length of rwork is too small
## lsodes will tell the length actually needed
# out2 <- lsodes(func = chemistry, y = y, parms = parms, times = times,
#              inz = nonzero, atol = atol,rtol = rtol)  #gives warning
out2 <- lsodes(func = chemistry, y = y, parms = parms, times = times, 
            sparsetype = "sparseusr", inz = nonzero,   
             atol = atol, rtol = rtol, verbose = TRUE, lrw = 353)
                            
## =======================================================================
## application 3. lsodes estimates the structure of the Jacobian
##                the Jacobian (vector) function is input
## =======================================================================
chemjac <- function (time, Y, j, pars) {
  with (as.list(pars), {
    PDJ <- rep(0,n)

    if (j == 1){
       PDJ[1] <- -rk1
       PDJ[2] <- rk1
    } else if (j == 2) {
       PDJ[2] <- -rk3*Y[3] - rk15*Y[12] - rk2
       PDJ[3] <- rk2 - rk3*Y[3]
       PDJ[4] <- rk3*Y[3]
       PDJ[5] <- rk15*Y[12]
       PDJ[12] <- -rk15*Y[12]
    } else if (j == 3) {
       PDJ[2] <- -rk3*Y[2]
       PDJ[3] <- -rk5 - rk3*Y[2] - rk7*Y[10]
       PDJ[4] <- rk3*Y[2]
       PDJ[6] <- rk7*Y[10]
       PDJ[10] <- rk5 - rk7*Y[10]
    } else if (j == 4) {
       PDJ[2] <- rk11*rk14
       PDJ[3] <- rk11*rk14
       PDJ[4] <- -rk11*rk14 - rk4
       PDJ[9] <- rk4
    } else if (j == 5) {
       PDJ[2] <- rk19*rk14
       PDJ[5] <- -rk19*rk14 - rk16
       PDJ[9] <- rk16
       PDJ[12] <- rk19*rk14
    } else if (j == 6) {
       PDJ[3] <- rk12*rk14
       PDJ[6] <- -rk12*rk14 - rk8
       PDJ[9] <- rk8
       PDJ[10] <- rk12*rk14
    } else if (j == 7) {
       PDJ[7] <- -rk20*rk14 - rk18
       PDJ[9] <- rk18
       PDJ[10] <- rk20*rk14
       PDJ[12] <- rk20*rk14
    } else if (j == 8) {
       PDJ[8] <- -rk13*rk14 - rk10
       PDJ[10] <- rk13*rk14
       PDJ[11] <- rk10
    } else if (j == 10) {
       PDJ[3] <- -rk7*Y[3]
       PDJ[6] <- rk7*Y[3]
       PDJ[7] <- rk17*Y[12]
       PDJ[8] <- rk9
       PDJ[10] <- -rk7*Y[3] - rk17*Y[12] - rk6 - rk9
       PDJ[12] <- rk6 - rk17*Y[12]
    } else if (j == 12) {
       PDJ[2] <- -rk15*Y[2]
       PDJ[5] <- rk15*Y[2]
       PDJ[7] <- rk17*Y[10]
       PDJ[10] <- -rk17*Y[10]
       PDJ[12] <- -rk15*Y[2] - rk17*Y[10]
    }
    return(PDJ)
  })
} 

out3 <- lsodes(func = chemistry, y = y, parms = parms, times = times,
              jacvec = chemjac, atol = atol, rtol = rtol)

## =======================================================================
## application 4. The structure of the Jacobian (nonzero elements) AND
##                the Jacobian (vector) function is input
## =======================================================================
out4 <- lsodes(func = chemistry, y = y, parms = parms, times = times,
               lrw = 351, sparsetype = "sparseusr", inz = nonzero,
               jacvec = chemjac, atol = atol, rtol = rtol,
               verbose = TRUE)
               
# The sparsejan variant 
# note: errors in inz may cause R to break, so this is not without danger...
# out5 <- lsodes(func = chemistry, y = y, parms = parms, times = times,
#               jacvec = chemjac, atol = atol, rtol = rtol, sparsetype = "sparsejan",
#               inz = c(1,3,8,13,17,21,25,29,32,32,38,38,43,                   # ian
#               1,2, 2,3,4,5,12, 2,3,4,6,10, 2,3,4,9, 2,5,9,12, 3,6,9,10,      # jan 
#               7,9,10,12, 8,10,11, 3,6,7,8,10,12, 2,5,7,10,12), lrw = 343) 

General Solver for Ordinary Differential Equations

Description

Solves a system of ordinary differential equations; a wrapper around the implemented ODE solvers

Usage

ode(y, times, func, parms, 
method = c("lsoda", "lsode", "lsodes", "lsodar", "vode", "daspk",
           "euler", "rk4", "ode23", "ode45", "radau", 
           "bdf", "bdf_d", "adams", "impAdams", "impAdams_d", "iteration"), ...)

## S3 method for class 'deSolve'
print(x, ...)
## S3 method for class 'deSolve'
summary(object, select = NULL, which = select, 
                 subset = NULL, ...)

Arguments

Details

This is simply a wrapper around the various ode solvers.

See package vignette for information about specifying the model in compiled code.

See the selected integrator for the additional options.

The default integrator used is lsoda.

The option method = "bdf" provdes a handle to the backward differentiation formula (it is equal to using method = "lsode"). It is best suited to solve stiff (systems of) equations.

The option method = "bdf_d" selects the backward differentiation formula that uses Jacobi-Newton iteration (neglecting the off-diagonal elements of the Jacobian (it is equal to using method = "lsode", mf = 23). It is best suited to solve stiff (systems of) equations.

method = "adams" triggers the Adams method that uses functional iteration (no Jacobian used); (equal to method = "lsode", mf = 10. It is often the best choice for solving non-stiff (systems of) equations. Note: when functional iteration is used, the method is often said to be explicit, although it is in fact implicit.

method = "impAdams" selects the implicit Adams method that uses Newton- Raphson iteration (equal to method = "lsode", mf = 12.

method = "impAdams_d" selects the implicit Adams method that uses Jacobi- Newton iteration, i.e. neglecting all off-diagonal elements (equal to method = "lsode", mf = 13.

For very stiff systems, method = "daspk" may outperform method = "bdf".

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 <karline.soetaert@nioz.nl>

See Also

• plot.deSolve for plotting the outputs,

• dede general solver for delay differential equations

• 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,

• aquaphy, ccl4model, where ode is used,

• lsoda, lsode, lsodes, lsodar, vode, daspk, radau,

• rk, rkMethod for additional Runge-Kutta methods,

• forcings and events,

• diagnostics to print diagnostic messages.

Examples

## =======================================================================
## Example1: Predator-Prey Lotka-Volterra model (with logistic prey)
## =======================================================================

LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion    <- rIng  * Prey * Predator
    GrowthPrey   <- rGrow * Prey * (1 - Prey/K)
    MortPredator <- rMort * Predator

    dPrey        <- GrowthPrey - Ingestion
    dPredator    <- Ingestion * assEff - MortPredator

    return(list(c(dPrey, dPredator)))
  })
}

pars  <- c(rIng   = 0.2,    # /day, rate of ingestion
           rGrow  = 1.0,    # /day, growth rate of prey
           rMort  = 0.2 ,   # /day, mortality rate of predator
           assEff = 0.5,    # -, assimilation efficiency
           K      = 10)     # mmol/m3, carrying capacity

yini  <- c(Prey = 1, Predator = 2)
times <- seq(0, 200, by = 1)
out   <- ode(yini, times, LVmod, pars)
summary(out)

## Default plot method
plot(out)

## User specified plotting
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "time", ylab = "Conc",
        main = "Lotka-Volterra", lwd = 2)
legend("topright", c("prey", "predator"), col = 1:2, lty = 1:2)

## =======================================================================
## Example2: Substrate-Producer-Consumer Lotka-Volterra model
## =======================================================================

## Note:
## Function sigimp passed as an argument (input) to model
##   (see also lsoda and rk examples)

SPCmod <- function(t, x, parms, input)  {
  with(as.list(c(parms, x)), {
    import <- input(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)
  })
}

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

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

## external signal with rectangle impulse
signal <- data.frame(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
xstart <- c(S = 1, P = 1, C = 1)

## Solve model
out <- ode(y = xstart, times = times,
           func = SPCmod, parms = parms, input = sigimp)

## Default plot method
plot(out)

## User specified plotting
mf <- par(mfrow = c(1, 2))
matplot(out[,1], out[,2:4], type = "l", xlab = "time", ylab = "state")
legend("topright", col = 1:3, lty = 1:3, legend = c("S", "P", "C"))
plot(out[,"P"], out[,"C"], type = "l", lwd = 2, xlab = "producer",
  ylab = "consumer")
par(mfrow = mf)

## =======================================================================
## Example3: Discrete time model - using method = "iteration"
##           The host-parasitoid model from Soetaert and Herman, 2009, 
##           Springer - p. 284.
## =======================================================================

Parasite <- function(t, y, ks) {
  P <- y[1]
  H <- y[2]
  f    <- A * P / (ks + H)
  Pnew <- H * (1 - exp(-f))
  Hnew <- H * exp(rH * (1 - H) - f)
  
  list (c(Pnew, Hnew))
}
rH <- 2.82 # rate of increase
A  <- 100  # attack rate
ks <- 15   # half-saturation density

out <- ode(func = Parasite, y = c(P = 0.5, H = 0.5), times = 0:50, parms = ks,
           method = "iteration")
            
out2<- ode(func = Parasite, y = c(P = 0.5, H = 0.5), times = 0:50, parms = 25,
           method = "iteration")

out3<- ode(func = Parasite, y = c(P = 0.5, H = 0.5), times = 0:50, parms = 35,
           method = "iteration")

## Plot all 3 scenarios in one figure
plot(out, out2, out3, lty = 1, lwd = 2)

## Same like "out", but *output* every two steps
## hini = 1 ensures that the same *internal* timestep of 1 is used
outb <- ode(func = Parasite, y = c(P = 0.5, H = 0.5),
            times = seq(0, 50, 2), hini = 1, parms = ks,
            method = "iteration")
plot(out, outb, type = c("l", "p"))

## Not run: 
## =======================================================================
## Example4: Playing with the Jacobian options - see e.g. lsoda help page
##
## IMPORTANT: The following example is temporarily broken because of 
##            incompatibility with R 3.0 on some systems.
##            A fix is on the way.
## =======================================================================

## a stiff equation, exponential decay, run 500 times
stiff <- function(t, y, p) {   # y and r are a 500-valued vector
  list(- r * y)
}

N    <- 500
r    <- runif(N, 15, 20)
yini <- runif(N, 1, 40)

times <- 0:10

## Using the default
print(system.time(
  out <- ode(y = yini, parms = NULL, times = times, func = stiff)
))
# diagnostics(out) shows that the method used = bdf (2), so it it stiff

## Specify that the Jacobian is banded, with nonzero values on the 
## diagonal, i.e. the bandwidth up and down = 0 

print(system.time(
  out2 <- ode(y = yini, parms = NULL, times = times, func = stiff,
              jactype = "bandint", bandup = 0, banddown = 0)
))

## Now we also specify the Jacobian function

jacob <- function(t, y, p) -r

print(system.time(
  out3 <- ode(y = yini, parms = NULL, times = times, func = stiff, 
              jacfunc = jacob, jactype = "bandusr", 
              bandup = 0, banddown = 0)
))
## The larger the value of N, the larger the time gain...

## End(Not run)

Solver For Multicomponent 1-D Ordinary Differential Equations

Description

Solves a system of ordinary differential equations resulting from 1-Dimensional partial differential equations that have been converted to ODEs by numerical differencing.

Usage

ode.1D(y, times, func, parms, nspec = NULL, dimens = NULL, 
   method= c("lsoda", "lsode", "lsodes", "lsodar", "vode", "daspk",
   "euler", "rk4", "ode23", "ode45", "radau", "bdf", "adams", "impAdams",
   "iteration"),
   names = NULL, bandwidth = 1, restructure = FALSE, ...)

Arguments

Details

This is the method of choice for multi-species 1-dimensional models, that are only subjected to transport between adjacent layers.

More specifically, this method is to be used if the state variables are arranged per species:

A[1], A[2], A[3],.... B[1], B[2], B[3],.... (for species A, B))

Two methods are implemented.

• The default method rearranges the state variables as A[1], B[1], ... A[2], B[2], ... A[3], B[3], .... This reformulation leads to a banded Jacobian with (upper and lower) half bandwidth = number of species.

  Then the selected integrator solves the banded problem.

• The second method uses lsodes. Based on the dimension of the problem, the method first calculates the sparsity pattern of the Jacobian, under the assumption that transport is only occurring between adjacent layers. Then lsodes is called to solve the problem.

  As lsodes is used to integrate, it may be necessary to specify the length of the real work array, lrw.

  Although a reasonable guess of lrw is made, it is possible that this will be too low. In this case, ode.1D will return with an error message telling the size of the work array actually needed. In the second try then, set lrw equal to this number.

  For instance, if you get the error:

     
  DLSODES- RWORK length is insufficient to proceed.                               
    Length needed is .ge. LENRW (=I1), exceeds LRW (=I2)                    
    In above message,  I1 =     27627   I2 =     25932 
  

  set lrw equal to 27627 or a higher value

If the model is specified in compiled code (in a DLL), then option 2, based on lsodes is the only solution method.

For single-species 1-D models, you may also use ode.band.

See the selected integrator for the additional options.

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.

The output will have the attributes istate, and rstate, two vectors with several useful elements. The first element of istate returns the conditions under which the last call to the integrator returned. Normal is istate = 2. If verbose = TRUE, the settings of istate and rstate will be written to the screen. See the help for the selected integrator for details.

Note

It is advisable though not mandatory to specify both nspec and dimens. In this case, the solver can check whether the input makes sense (i.e. if nspec * dimens == length(y)).

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

See Also

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

• ode.band for integrating models with a banded Jacobian

• ode.2D for integrating 2-D models

• ode.3D for integrating 3-D models

• lsodes,lsode, lsoda, lsodar,vode for the integration options.

diagnostics to print diagnostic messages.

Examples

## =======================================================================
## example 1
## a predator and its prey diffusing on a flat surface
## in concentric circles
## 1-D model with using cylindrical coordinates
## Lotka-Volterra type biology
## =======================================================================

## ================
## Model equations
## ================

lvmod <- function (time, state, parms, N, rr, ri, dr, dri) {
  with (as.list(parms), {
    PREY <- state[1:N]
    PRED <- state[(N+1):(2*N)]

    ## Fluxes due to diffusion
    ## at internal and external boundaries: zero gradient
    FluxPrey <- -Da * diff(c(PREY[1], PREY, PREY[N]))/dri
    FluxPred <- -Da * diff(c(PRED[1], PRED, PRED[N]))/dri

    ## Biology: Lotka-Volterra model
    Ingestion     <- rIng  * PREY * PRED
    GrowthPrey    <- rGrow * PREY * (1-PREY/cap)
    MortPredator  <- rMort * PRED

    ## Rate of change = Flux gradient + Biology
    dPREY    <- -diff(ri * FluxPrey)/rr/dr   +
                GrowthPrey - Ingestion
    dPRED    <- -diff(ri * FluxPred)/rr/dr   +
                Ingestion * assEff - MortPredator

    return (list(c(dPREY, dPRED)))
  })
}

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

## model parameters:

R  <- 20                        # total radius of surface, m
N  <- 100                       # 100 concentric circles
dr <- R/N                       # thickness of each layer
r  <- seq(dr/2,by = dr,len = N) # distance of center to mid-layer
ri <- seq(0,by = dr,len = N+1)  # distance to layer interface
dri <- dr                       # dispersion distances

parms <- c(Da     = 0.05,       # m2/d, dispersion coefficient
           rIng   = 0.2,        # /day, rate of ingestion
           rGrow  = 1.0,        # /day, growth rate of prey
           rMort  = 0.2 ,       # /day, mortality rate of pred
           assEff = 0.5,        # -, assimilation efficiency
           cap    = 10)         # density, carrying capacity

## Initial conditions: both present in central circle (box 1) only
state    <- rep(0, 2 * N)
state[1] <- state[N + 1] <- 10
                
## RUNNING the model:
times  <- seq(0, 200, by = 1)   # output wanted at these time intervals

## the model is solved by the two implemented methods:
## 1. Default: banded reformulation
print(system.time(
  out <- ode.1D(y = state, times = times, func = lvmod, parms = parms,
                nspec = 2, names = c("PREY", "PRED"),
                N = N, rr = r, ri = ri, dr = dr, dri = dri)
))

## 2. Using sparse method
print(system.time(
  out2 <- ode.1D(y = state, times = times, func = lvmod, parms = parms,
                 nspec = 2, names = c("PREY","PRED"), 
                 N = N, rr = r, ri = ri, dr = dr, dri = dri,
                 method = "lsodes")
))

## ================
## Plotting output
## ================
# the data in 'out' consist of: 1st col times, 2-N+1: the prey
# N+2:2*N+1: predators

PREY   <- out[, 2:(N + 1)]

filled.contour(x = times, y = r, PREY, color = topo.colors,
               xlab = "time, days", ylab = "Distance, m",
               main = "Prey density")
# similar:
image(out, which = "PREY", grid = r, xlab = "time, days", 
      legend = TRUE, ylab = "Distance, m", main = "Prey density")

image(out2, grid = r)

# summaries of 1-D variables
summary(out)

# 1-D plots:
matplot.1D(out, type = "l", subset = time == 10)
matplot.1D(out, type = "l", subset = time > 10 & time < 20)

## =======================================================================
## Example 2.
## Biochemical Oxygen Demand (BOD) and oxygen (O2) dynamics
## in a river
## =======================================================================

## ================
## Model equations
## ================
O2BOD <- function(t, state, pars) {
  BOD <- state[1:N]
  O2  <- state[(N+1):(2*N)]

  ## BOD dynamics
  FluxBOD <- v * c(BOD_0, BOD)   # fluxes due to water transport
  FluxO2  <- v * c(O2_0, O2)
  
  BODrate <- r * BOD             # 1-st order consumption

  ## rate of change = flux gradient  - consumption + reaeration (O2)
  dBOD         <- -diff(FluxBOD)/dx - BODrate
  dO2          <- -diff(FluxO2)/dx  - BODrate      +  p * (O2sat-O2)

  return(list(c(dBOD = dBOD, dO2 = dO2)))
}
 
 
## ==================
## Model application
## ==================
## parameters
dx      <- 25        # grid size of 25 meters
v       <- 1e3       # velocity, m/day
x       <- seq(dx/2, 5000, by = dx)  # m, distance from river
N       <- length(x)
r       <- 0.05      # /day, first-order decay of BOD
p       <- 0.5       # /day, air-sea exchange rate 
O2sat   <- 300       # mmol/m3 saturated oxygen conc
O2_0    <- 200       # mmol/m3 riverine oxygen conc
BOD_0   <- 1000      # mmol/m3 riverine BOD concentration

## initial conditions:
state <- c(rep(200, N), rep(200, N))
times <- seq(0, 20, by = 0.1)

## running the model
##  step 1  : model spinup
out <- ode.1D(y = state, times, O2BOD, parms = NULL, 
              nspec = 2, names = c("BOD", "O2"))

## ================
## Plotting output
## ================
## select oxygen (first column of out:time, then BOD, then O2
O2   <- out[, (N + 2):(2 * N + 1)]
color = topo.colors

filled.contour(x = times, y = x, O2, color = color, nlevels = 50,
               xlab = "time, days", ylab = "Distance from river, m",
               main = "Oxygen")
               
## or quicker plotting:
image(out, grid = x,  xlab = "time, days", ylab = "Distance from river, m")               

Solver for 2-Dimensional Ordinary Differential Equations

Description

Solves a system of ordinary differential equations resulting from 2-Dimensional partial differential equations that have been converted to ODEs by numerical differencing.

Usage

ode.2D(y, times, func, parms, nspec = NULL, dimens,
  method= c("lsodes", "euler", "rk4", "ode23", "ode45", "adams", "iteration"),
  names = NULL, cyclicBnd = NULL, ...)

Arguments

Details

This is the method of choice for 2-dimensional models, that are only subjected to transport between adjacent layers.

Based on the dimension of the problem, and if lsodes is used as the integrator, the method first calculates the sparsity pattern of the Jacobian, under the assumption that transport is only occurring between adjacent layers. Then lsodes is called to solve the problem.

If the model is not stiff, then it is more efficient to use one of the explicit integration routines

In some cases, a cyclic boundary condition exists. This is when the first boxes in x-or y-direction interact with the last boxes. In this case, there will be extra non-zero fringes in the Jacobian which need to be taken into account. The occurrence of cyclic boundaries can be toggled on by specifying argument cyclicBnd. For innstance, cyclicBnd = 1 indicates that a cyclic boundary is required only for the x-direction, whereas cyclicBnd = c(1,2) imposes a cyclic boundary for both x- and y-direction. The default is no cyclic boundaries.

If lsodes is used to integrate, it will probably be necessary to specify the length of the real work array, lrw.

Although a reasonable guess of lrw is made, it is likely that this will be too low. In this case, ode.2D will return with an error message telling the size of the work array actually needed. In the second try then, set lrw equal to this number.

For instance, if you get the error:

   
DLSODES- RWORK length is insufficient to proceed.                               
  Length needed is .ge. LENRW (=I1), exceeds LRW (=I2)                    
  In above message,  I1 =     27627   I2 =     25932 

set lrw equal to 27627 or a higher value.

See lsodes for the additional options.

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.

The output will have the attributes istate, and rstate, two vectors with several useful elements. The first element of istate returns the conditions under which the last call to the integrator returned. Normal is istate = 2. If verbose = TRUE, the settings of istate and rstate will be written to the screen. See the help for the selected integrator for details.

Note

It is advisable though not mandatory to specify both nspec and dimens. In this case, the solver can check whether the input makes sense (as nspec * dimens[1] * dimens[2] == length(y)).

Do not use this method for problems that are not 2D!

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

See Also

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

• ode.band for integrating models with a banded Jacobian

• ode.1D for integrating 1-D models

• ode.3D for integrating 3-D models

• lsodes for the integration options.

diagnostics to print diagnostic messages.

Examples

## =======================================================================
## A Lotka-Volterra predator-prey model with predator and prey
## dispersing in 2 dimensions
## =======================================================================

## ==================
## Model definitions
## ==================

lvmod2D <- function (time, state, pars, N, Da, dx) {
  NN <- N*N
  Prey <- matrix(nrow = N, ncol = N,state[1:NN])
  Pred <- matrix(nrow = N, ncol = N,state[(NN+1):(2*NN)])

  with (as.list(pars), {
    ## Biology
    dPrey <- rGrow * Prey * (1- Prey/K) - rIng  * Prey * Pred
    dPred <- rIng  * Prey * Pred*assEff - rMort * Pred

    zero <- rep(0, N)

    ## 1. Fluxes in x-direction; zero fluxes near boundaries
    FluxPrey <- -Da * rbind(zero,(Prey[2:N,] - Prey[1:(N-1),]), zero)/dx
    FluxPred <- -Da * rbind(zero,(Pred[2:N,] - Pred[1:(N-1),]), zero)/dx

    ## Add flux gradient to rate of change
    dPrey    <- dPrey - (FluxPrey[2:(N+1),] - FluxPrey[1:N,])/dx
    dPred    <- dPred - (FluxPred[2:(N+1),] - FluxPred[1:N,])/dx

    ## 2. Fluxes in y-direction; zero fluxes near boundaries
    FluxPrey <- -Da * cbind(zero,(Prey[,2:N] - Prey[,1:(N-1)]), zero)/dx
    FluxPred <- -Da * cbind(zero,(Pred[,2:N] - Pred[,1:(N-1)]), zero)/dx

    ## Add flux gradient to rate of change
    dPrey    <- dPrey - (FluxPrey[,2:(N+1)] - FluxPrey[,1:N])/dx
    dPred    <- dPred - (FluxPred[,2:(N+1)] - FluxPred[,1:N])/dx

    return(list(c(as.vector(dPrey), as.vector(dPred))))
 })
}


## ===================
## Model applications
## ===================

pars    <- c(rIng   = 0.2,    # /day, rate of ingestion
             rGrow  = 1.0,    # /day, growth rate of prey
             rMort  = 0.2 ,   # /day, mortality rate of predator
             assEff = 0.5,    # -, assimilation efficiency
             K      = 5  )    # mmol/m3, carrying capacity

R  <- 20                      # total length of surface, m
N  <- 50                      # number of boxes in one direction
dx <- R/N                     # thickness of each layer
Da <- 0.05                    # m2/d, dispersion coefficient

NN <- N*N                     # total number of boxes

## initial conditions
yini    <- rep(0, 2*N*N)
cc      <- c((NN/2):(NN/2+1)+N/2, (NN/2):(NN/2+1)-N/2)
yini[cc] <- yini[NN+cc] <- 1

## solve model (5000 state variables...  use Cash-Karp Runge-Kutta method
times   <- seq(0, 50, by = 1)
out <- ode.2D(y = yini, times = times, func = lvmod2D, parms = pars,
              dimens = c(N, N), names = c("Prey", "Pred"),
              N = N, dx = dx, Da = Da, method = rkMethod("rk45ck"))

diagnostics(out)
summary(out)

# Mean of prey concentration at each time step
Prey <- subset(out, select = "Prey", arr = TRUE)
dim(Prey)
MeanPrey <- apply(Prey, MARGIN = 3, FUN = mean)
plot(times, MeanPrey)

## Not run: 
## plot results
Col <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                          "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

for (i in seq(1, length(times), by = 1))
   image(Prey[ , ,i],
   col = Col(100), xlab = , zlim = range(out[,2:(NN+1)]))

## similar, plotting both and adding a margin text with times:
image(out, xlab = "x", ylab = "y", mtext = paste("time = ", times))

## End(Not run)

select <- c(1, 40)
image(out, xlab = "x", ylab = "y", mtext = "Lotka-Volterra in 2-D",
      subset = select, mfrow = c(2,2), legend =  TRUE)


# plot prey and pred at t = 10; first use subset to select data
prey10 <-  matrix (nrow = N, ncol = N, 
     data = subset(out, select = "Prey", subset = (time == 10)))
pred10 <-  matrix (nrow = N, ncol = N, 
     data = subset(out, select = "Pred", subset = (time == 10)))

mf <- par(mfrow = c(1, 2))
image(prey10)
image(pred10)
par (mfrow = mf)

# same, using deSolve's image:
image(out, subset = (time == 10))

## =======================================================================
## An example with a cyclic boundary condition.
## Diffusion in 2-D; extra flux on 2 boundaries,
## cyclic boundary in y
## =======================================================================


diffusion2D <- function(t, Y, par) {
  y    <- matrix(nrow = nx, ncol = ny, data = Y)  # vector to 2-D matrix
  dY   <- -r * y       # consumption
  BNDx <- rep(1, nx)   # boundary concentration
  BNDy <- rep(1, ny)   # boundary concentration

  ## diffusion in X-direction; boundaries=imposed concentration
  Flux <- -Dx * rbind(y[1,] - BNDy, (y[2:nx,] - y[1:(nx-1),]), BNDy - y[nx,])/dx
  dY   <- dY - (Flux[2:(nx+1),] - Flux[1:nx,])/dx

  ## diffusion in Y-direction
  Flux <- -Dy * cbind(y[,1] - BNDx, (y[,2:ny]-y[,1:(ny-1)]), BNDx - y[,ny])/dy
  dY   <- dY - (Flux[,2:(ny+1)] - Flux[,1:ny])/dy

  ## extra flux on two sides
  dY[,1] <- dY[,1] + 10
  dY[1,] <- dY[1,] + 10

  ## and exchange between sides on y-direction
  dY[,ny] <- dY[,ny] + (y[,1] - y[,ny]) * 10
  return(list(as.vector(dY)))
}

## parameters
dy <- dx <- 1  # grid size
Dy <- Dx <- 1  # diffusion coeff, X- and Y-direction
r  <- 0.05     # consumption rate

nx <- 50
ny <- 100
y  <- matrix(nrow = nx, ncol = ny, 1)

## model most efficiently solved with lsodes - need to specify lrw

print(system.time(
  ST3 <- ode.2D(y, times = 1:100, func = diffusion2D, parms = NULL,
                dimens = c(nx, ny), verbose = TRUE, names = "Y",
                lrw = 400000, atol = 1e-10, rtol = 1e-10, cyclicBnd = 2)
))

# summary of 2-D variable
summary(ST3)

# plot output at t = 10
t10 <-  matrix (nrow = nx, ncol = ny, 
     data = subset(ST3, select = "Y", subset = (time == 10)))

persp(t10, theta = 30, border = NA, phi = 70, 
     col = "lightblue", shade = 0.5, box = FALSE)

# image plot, using deSolve's image function
image(ST3, subset = time == 10, method = "persp",
      theta = 30, border = NA, phi = 70, main = "", 
      col = "lightblue", shade = 0.5, box = FALSE)

## Not run: 
zlim <- range(ST3[, -1])
for (i in 2:nrow(ST3)) {
  y <- matrix(nrow = nx, ncol = ny, data = ST3[i, -1])
  filled.contour(y, zlim = zlim, main = i)
}

# same
image(ST3, method = "filled.contour")

## End(Not run)

Solver for 3-Dimensional Ordinary Differential Equations

Description

Solves a system of ordinary differential equations resulting from 3-Dimensional partial differential equations that have been converted to ODEs by numerical differencing.

Usage

ode.3D(y, times, func, parms, nspec = NULL, dimens, 
  method = c("lsodes", "euler", "rk4", "ode23", "ode45", "adams", "iteration"),
  names = NULL, cyclicBnd = NULL, ...)

Arguments

Details

This is the method of choice for 3-dimensional models, that are only subjected to transport between adjacent layers.

Based on the dimension of the problem, the method first calculates the sparsity pattern of the Jacobian, under the assumption that transport is only occurring between adjacent layers. Then lsodes is called to solve the problem.

As lsodes is used to integrate, it will probably be necessary to specify the length of the real work array, lrw.

Although a reasonable guess of lrw is made, it is likely that this will be too low.

In this case, ode.2D will return with an error message telling the size of the work array actually needed. In the second try then, set lrw equal to this number.

For instance, if you get the error:

   
DLSODES- RWORK length is insufficient to proceed.                               
  Length needed is .ge. LENRW (=I1), exceeds LRW (=I2)                    
  In above message,  I1 =     27627   I2 =     25932 

set lrw equal to 27627 or a higher value.

See lsodes for the additional options.

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.

The output will have the attributes istate, and rstate, two vectors with several useful elements. The first element of istate returns the conditions under which the last call to the integrator returned. Normal is istate = 2. If verbose = TRUE, the settings of istate and rstate will be written to the screen. See the help for the selected integrator for details.

Note

It is advisable though not mandatory to specify both nspec and dimens. In this case, the solver can check whether the input makes sense (as nspec*dimens[1]*dimens[2]*dimens[3] == length(y)).

Do not use this method for problems that are not 3D!

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

See Also

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

• ode.band for integrating models with a banded Jacobian

• ode.1D for integrating 1-D models

• ode.2D for integrating 2-D models

• lsodes for the integration options.

diagnostics to print diagnostic messages.

Examples

## =======================================================================
## Diffusion in 3-D; imposed boundary conditions
## =======================================================================
diffusion3D <- function(t, Y, par) {

  ## function to bind two matrices to an array
  mbind <- function (Mat1, Array, Mat2, along = 1) {
    dimens <- dim(Array) + c(0, 0, 2)
    if (along == 3)
      array(dim = dimens, data = c(Mat1, Array, Mat2))
    else if (along == 1)
      aperm(array(dim = dimens,
        data=c(Mat1, aperm(Array, c(3, 2, 1)), Mat2)), c(3, 2, 1))
    else if (along == 2)
      aperm(array(dim = dimens,
        data = c(Mat1, aperm(Array, c(1, 3, 2)), Mat2)), c(1, 3, 2))
  }

  yy  <- array(dim=c(n, n, n), data = Y)      # vector to 3-D array
  dY  <- -r*yy                                # consumption
  BND <- matrix(nrow = n, ncol = n, data = 1) # boundary concentration

  ## diffusion in x-direction
  ## new array including boundary concentrations in X-direction
  BNDx <- mbind(BND, yy, BND, along = 1)
  ## diffusive Flux
  Flux <- -Dx * (BNDx[2:(n+2),,] - BNDx[1:(n+1),,])/dx
  ## rate of change = - flux gradient
  dY[] <- dY[] - (Flux[2:(n+1),,] - Flux[1:n,,])/dx

  ## diffusion in y-direction
  BNDy <- mbind(BND, yy, BND, along = 2)
  Flux <- -Dy * (BNDy[,2:(n+2),] - BNDy[,1:(n+1),])/dy
  dY[] <- dY[] - (Flux[,2:(n+1),] - Flux[,1:n,])/dy

  ## diffusion in z-direction
  BNDz <- mbind(BND, yy, BND, along = 3)
  Flux <- -Dz * (BNDz[,,2:(n+2)] - BNDz[,,1:(n+1)])/dz
  dY[] <- dY[] - (Flux[,,2:(n+1)] - Flux[,,1:n])/dz

  return(list(as.vector(dY)))
}

## parameters
dy <- dx <- dz <-1   # grid size
Dy <- Dx <- Dz <-1   # diffusion coeff, X- and Y-direction
r  <- 0.025          # consumption rate

n  <- 10
y  <- array(dim=c(n,n,n),data=10.)

## use lsodes, the default (for n>20, Runge-Kutta more efficient)
print(system.time(
  RES <- ode.3D(y, func = diffusion3D, parms = NULL, dimens = c(n, n, n),
                times = 1:20, lrw = 120000, atol = 1e-10,
                rtol = 1e-10, verbose = TRUE)
))

y <- array(dim = c(n, n, n), data = RES[nrow(RES), -1])
filled.contour(y[, , n/2], color.palette = terrain.colors)

summary(RES)

## Not run: 
for (i in 2:nrow(RES)) {
  y <- array(dim=c(n,n,n),data=RES[i,-1])
  filled.contour(y[,,n/2],main=i,color.palette=terrain.colors)
}

## End(Not run)

Solver for Ordinary Differential Equations; Assumes a Banded Jacobian

Description

Solves a system of ordinary differential equations.

Assumes a banded Jacobian matrix, but does not rearrange the state variables (in contrast to ode.1D). Suitable for 1-D models that include transport only between adjacent layers and that model only one species.

Usage

ode.band(y, times, func, parms, nspec = NULL, dimens = NULL, 
  bandup = nspec, banddown = nspec, method = "lsode", names = NULL, 
  ...)

Arguments

Details

This is the method of choice for single-species 1-D reactive transport models.

For multi-species 1-D models, this method can only be used if the state variables are arranged per box, per species (e.g. A[1], B[1], A[2], B[2], A[3], B[3], ... for species A, B). By default, the model function will have the species arranged as A[1], A[2], A[3], ... B[1], B[2], B[3], ... in this case, use ode.1D.

See the selected integrator for the additional options.

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.

The output will have the attributes istate and rstate, two vectors with several elements. See the help for the selected integrator for details. the first element of istate returns the conditions under which the last call to the integrator returned. Normal is istate = 2. If verbose = TRUE, the settings of istate and rstate will be written to the screen.

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

See Also

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

• ode.1D for integrating 1-D models

• ode.2D for integrating 2-D models

• ode.3D for integrating 3-D models

• lsode, lsoda, lsodar, vode for the integration options.

diagnostics to print diagnostic messages.

Examples

## =======================================================================
## The Aphid model from Soetaert and Herman, 2009.
## A practical guide to ecological modelling.
## Using R as a simulation platform. Springer.
## =======================================================================

## 1-D diffusion model

## ================
## Model equations
## ================
Aphid <- function(t, APHIDS, parameters) {
  deltax  <- c (0.5*delx, rep(delx, numboxes-1), 0.5*delx)
  Flux    <- -D*diff(c(0, APHIDS, 0))/deltax
  dAPHIDS <- -diff(Flux)/delx + APHIDS*r

  list(dAPHIDS)   # the output
}
  
## ==================
## Model application
## ==================

## the model parameters:
D         <- 0.3    # m2/day  diffusion rate
r         <- 0.01   # /day    net growth rate
delx      <- 1      # m       thickness of boxes
numboxes  <- 60 

## distance of boxes on plant, m, 1 m intervals
Distance  <- seq(from = 0.5, by = delx, length.out = numboxes)

## Initial conditions, ind/m2
## aphids present only on two central boxes
APHIDS        <- rep(0, times = numboxes)
APHIDS[30:31] <- 1
state         <- c(APHIDS = APHIDS)      # initialise state variables 
                  
## RUNNING the model:
times <- seq(0, 200, by = 1)   # output wanted at these time intervals
out   <- ode.band(state, times, Aphid, parms = 0, 
                  nspec = 1, names = "Aphid")

## ================
## Plotting output
## ================
image(out, grid = Distance, method = "filled.contour", 
      xlab = "time, days", ylab = "Distance on plant, m",
      main = "Aphid density on a row of plants")

matplot.1D(out, grid = Distance, type = "l", 
   subset = time %in% seq(0, 200, by = 10))

# add an observed dataset to 1-D plot (make sure to use correct name):
data <- cbind(dist  = c(0,10, 20,  30,  40, 50, 60), 
              Aphid = c(0,0.1,0.25,0.5,0.25,0.1,0))

matplot.1D(out, grid = Distance, type = "l", 
   subset = time %in% seq(0, 200, by = 10), 
   obs = data, obspar = list(pch = 18, cex = 2, col="red"))
## Not run: 
plot.1D(out, grid = Distance, type = "l")


## End(Not run)

Plot, Image and Histogram Method for deSolve Objects

Description

Plot the output of numeric integration routines.

Usage

## S3 method for class 'deSolve'
plot(x, ..., select = NULL, which = select, ask = NULL,
              obs = NULL, obspar = list(), subset = NULL)
## S3 method for class 'deSolve'
hist(x, select = 1:(ncol(x)-1), which  = select, ask = NULL,
              subset = NULL, ...)
## S3 method for class 'deSolve'
image(x, select = NULL, which = select, ask = NULL,
              add.contour = FALSE, grid = NULL,
              method = "image", legend = FALSE, subset = NULL, ...)
## S3 method for class 'deSolve'
subset(x, subset = NULL, select = NULL,
             which = select, arr = FALSE, ...)

plot.1D (x, ..., select = NULL, which = select, ask = NULL,
         obs = NULL, obspar = list(), grid = NULL,
         xyswap = FALSE, delay = 0, vertical = FALSE, subset = NULL)

matplot.0D(x, ..., select = NULL, which = select,
           obs = NULL, obspar = list(), subset = NULL,
           legend = list(x = "topright"))

matplot.1D(x, select = NULL, which = select, ask = NULL,
         obs = NULL, obspar = list(), grid = NULL,
         xyswap = FALSE, vertical = FALSE, subset = NULL, ...)

Arguments

Details

The number of panels per page is automatically determined up to 3 x 3 (par(mfrow = c(3, 3))). This default can be overwritten by specifying user-defined settings for mfrow or mfcol. Set mfrow equal to NULL to avoid the plotting function to change user-defined mfrow or mfcol settings.

Other graphical parameters can be passed as well. Parameters are vectorized, either according to the number of plots (xlab, ylab, main, sub, xlim, ylim, log, asp, ann, axes, frame.plot, panel.first, panel.last, cex.lab, cex.axis, cex.main) or according to the number of lines within one plot (other parameters e.g. col, lty, lwd etc.) so it is possible to assign specific axis labels to individual plots, resp. different plotting style. Plotting parameter ylim, or xlim can also be a list to assign different axis limits to individual plots.

Similarly, the graphical parameters for observed data, as passed by obspar can be vectorized, according to the number of observed data sets.

Image plots will only work for 1-D and 2-D variables, as solved with ode.1D and ode.2D. In the first case, an image with times as x- and the grid as y-axis will be created. In the second case, an x-y plot will be created, for all times. Unless ask = FALSE, the user will be asked to confirm page changes. Via argument mtext, it is possible to label each page in case of 2D output.

For images, it is possible to pass an argument method which can take the values "image" (default), "filled.contour", "contour" or "persp", in order to use the respective plotting method.

plot and matplot.0D will always have times on the x-axis. For problems solved with ode.1D, it may be more useful to use plot.1D or matplot.1D which will plot how spatial variables change with time. These plots will have the grid on the x-axis.

Value

Function subset called with arr = FALSE will return a matrix with up to as many rows as selected by subset and as many columns as selected variables.

When arr = TRUE then an array will be outputted with dimensions equal to the dimension of the selected variable, augmented with the number of rows selected by subset. This means that the last dimension points to times.

Function subset also has an attribute that contains the times selected.

See Also

deSolve, ode, print.deSolve,

hist image matplot, plot.default for the underlying functions from package graphics,

ode.2D, for an example of using subset with arr = TRUE.

Examples

## =======================================================================
## Example 1. A Predator-Prey model with 4 species in matrix formulation
## =======================================================================

LVmatrix <- function(t, n, parms) {
  with(parms, {
    dn <- r * n + n * (A %*% n)
    return(list(c(dn)))
  })
}
parms <- list(
  r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
  A = matrix(c(0.0, 0.0, -0.2, 0.01,      # prey 1
               0.0, 0.0, 0.02, -0.1,      # prey 2
               0.2, 0.02, 0.0, 0.0,       # predator 1; prefers prey 1
               0.01, 0.1, 0.0, 0.0),      # predator 2; prefers prey 2
               nrow = 4, ncol = 4, byrow=TRUE)
)
times <- seq(from = 0, to = 500, by = 0.1)
y     <- c(prey1 = 1, prey2 = 1, pred1 = 2, pred2 = 2)

out <- ode(y, times, LVmatrix, parms)

## Basic line plot
plot(out, type = "l")

## User-specified axis labels
plot(out, type = "l", ylab = c("Prey 1", "Prey 2", "Pred 1", "Pred 2"),
  xlab = "Time (d)", main = "Time Series")

## Set user-defined mfrow
pm <- par (mfrow = c(2, 2))

## "mfrow=NULL" keeps user-defined mfrow
plot(out, which = c("prey1", "pred2"), mfrow = NULL, type = "l", lwd = 2)

plot(out[,"prey1"], out[,"pred1"], xlab="prey1",
  ylab = "pred1", type = "l", lwd = 2)
plot(out[,"prey2"], out[,"pred2"], xlab = "prey2",
  ylab = "pred2", type = "l",lwd = 2)

## restore graphics parameters
par ("mfrow" = pm)

## Plot all in one figure, using matplot
matplot.0D(out, lwd = 2)

## Split y-variables in two groups
matplot.0D(out, which = list(c(1,3), c(2,4)),
           lty = c(1,2,1,2), col=c(4,4,5,5),
           ylab = c("prey1,pred1", "prey2,pred2"))

## =======================================================================
## Example 2. Add second and third output, and observations
## =======================================================================

# New runs with different parameter settings
parms2      <- parms
parms2$r[1] <- 0.2
out2 <- ode(y, times, LVmatrix, parms2)

# New runs with different parameter settings
parms3      <- parms
parms3$r[1] <- 0.05
out3 <- ode(y, times, LVmatrix, parms3)

# plot all three outputs
plot(out, out2, out3, type = "l",
     ylab = c("Prey 1", "Prey 2", "Pred 1", "Pred 2"),
     xlab = "Time (d)", main = c("Prey 1", "Prey 2", "Pred 1", "Pred 2"),
     col = c("red", "blue", "darkred"))


## 'observed' data
obs <- as.data.frame(out[out[,1] %in% seq(10, 500, by = 30), ])

plot(out, which = "prey1", type = "l", obs = obs,
     obspar = list(pch = 18, cex = 2))

plot(out, type = "l", obs = obs, col = "red")

matplot.0D(out, which = c("prey1", "pred1"), type = "l", obs = obs)

## second set of 'observed' data and two outputs
obs2 <- as.data.frame(out2[out2[,1] %in% seq(10, 500, by = 50), ])

## manual xlim, log
plot(out, out2, type = "l", obs = list(obs, obs2), col = c("red", "blue"),
      obspar = list(pch = 18:19, cex = 2, col = c("red", "blue")),
      log = c("y", ""), which = c("prey1", "prey1"),
      xlim = list(c(100, 500), c(0, 400)))

## data in 'long' format
OBS <- data.frame(name = c(rep("prey1", 3), rep("prey2", 2)),
                  time = c(10, 100, 250, 10, 400),
                  value = c(0.05, 0.04, 0.7, 0.5, 1))
OBS
plot(out, obs = OBS, obspar = c(pch = 18, cex = 2))

# a subset only:
plot(out, subset = prey1 < 0.5, type = "p")

# Simple histogram
hist(out, col = "darkblue", breaks = 50)

hist(out, col = "darkblue", breaks = 50, subset = prey1<1 & prey2 < 1)

# different parameters per plot
hist(out, col = c("darkblue", "red", "orange", "black"),
     breaks = c(10,50))

## =======================================================================
## The Aphid model from Soetaert and Herman, 2009.
## A practical guide to ecological modelling.
## Using R as a simulation platform. Springer.
## =======================================================================

## 1-D diffusion model

## ================
## Model equations
## ================
Aphid <- function(t, APHIDS, parameters) {
  deltax  <- c (0.5*delx, rep(delx, numboxes - 1), 0.5*delx)
  Flux    <- -D * diff(c(0, APHIDS, 0))/deltax
  dAPHIDS <- -diff(Flux)/delx + APHIDS * r
  list(dAPHIDS, Flux = Flux)
}

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

## the model parameters:
D         <- 0.3    # m2/day  diffusion rate
r         <- 0.01   # /day    net growth rate
delx      <- 1      # m       thickness of boxes
numboxes  <- 60

## distance of boxes on plant, m, 1 m intervals
Distance  <- seq(from = 0.5, by = delx, length.out = numboxes)

## Initial conditions, ind/m2
## aphids present only on two central boxes
APHIDS        <- rep(0, times = numboxes)
APHIDS[30:31] <- 1
state         <- c(APHIDS = APHIDS)      # initialise state variables

## RUNNING the model:
times <- seq(0, 200, by = 1)   # output wanted at these time intervals
out   <- ode.1D(state, times, Aphid, parms = 0, nspec = 1, names = "Aphid")

image(out, grid = Distance, main = "Aphid model", ylab = "distance, m",
  legend = TRUE)

## restricting time
image(out, grid = Distance, main = "Aphid model", ylab = "distance, m",
  legend = TRUE, subset = time < 100)

image(out, grid = Distance, main = "Aphid model", ylab = "distance, m",
  method = "persp", border = NA, theta = 30)

FluxAphid <- subset(out, select = "Flux", subset = time < 50)

matplot.1D(out, type = "l", lwd = 2, xyswap = TRUE, lty = 1)

matplot.1D(out, type = "l", lwd = 2, xyswap = TRUE, lty = 1,
           subset = time < 50)

matplot.1D(out, type = "l", lwd = 2, xyswap = TRUE, lty = 1,
           subset = time %in% seq(0, 200, by = 10), col = "grey")

## Not run: 
  plot(out, ask = FALSE, mfrow = c(1, 1))
  plot.1D(out, ask = FALSE, type = "l", lwd = 2, xyswap = TRUE)

## End(Not run)

## see help file for ode.2D for images of 2D variables

Implicit Runge-Kutta RADAU IIA

Description

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

dy/dt = f(t,y)

or linearly implicit differential algebraic equations in the form:

M dy/dt = f(t,y)

.

The R function radau provides an interface to the Fortran solver RADAU5, written by Ernst Hairer and G. Wanner, which implements the 3-stage RADAU IIA method. It implements the implicit Runge-Kutta method of order 5 with step size control and continuous output. The system of ODEs or DAEs is written as an R function or can be defined in compiled code that has been dynamically loaded.

Usage

radau(y, times, func, parms, nind = c(length(y), 0, 0), 
  rtol = 1e-6, atol = 1e-6, jacfunc = NULL, jactype = "fullint", 
  mass = NULL, massup = NULL, massdown = NULL, rootfunc = NULL,
  verbose = FALSE, nroot = 0, hmax = NULL, hini = 0, ynames = TRUE,
  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

Details

The work is done by the FORTRAN subroutine RADAU5, whose documentation should be consulted for details. The implementation is based on the Fortran 77 version from January 18, 2002.

There are four standard choices for the Jacobian which can be specified with jactype.

The options for jactype are

jactype = "fullint"

a full Jacobian, calculated internally by the solver.

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"

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

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

The input parameters rtol, and atol determine the error control performed by the solver, which roughly keeps the local error of y(i) below rtol(i)*abs(y(i))+atol(i).

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

See vignette("deSolve") from the deSolve package 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" from package deSolve for details.

Information about linking forcing functions to compiled code is in forcings (from package deSolve).

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

Value

A matrix of class deSolve with up to as many rows as elements in times and as many columns as elements in y plus the number of "global" values returned in the next elements of the return from func, plus and additional column for the time value. There will be a row for each element in times unless the FORTRAN routine 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

References

E. Hairer and G. Wanner, 1996. Solving Ordinary Differential Equations II. Stiff and Differential-algebraic problems. Springer series in computational mathematics 14, Springer-Verlag, second edition.

See Also

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

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

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

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

• daspk for integrating DAE models up to index 1

diagnostics to print diagnostic messages.

Examples

## =======================================================================
## Example 1: ODE
##   Various ways to solve the same model.
## =======================================================================

## the model, 5 state variables
f1 <- function  (t, y, parms) {
  ydot <- vector(len = 5)

  ydot[1] <-  0.1*y[1] -0.2*y[2]
  ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3]
  ydot[3] <-           -0.3*y[2] +0.1*y[3] -0.2*y[4]
  ydot[4] <-                     -0.3*y[3] +0.1*y[4] -0.2*y[5]
  ydot[5] <-                               -0.3*y[4] +0.1*y[5]

  return(list(ydot))
}

## the Jacobian, written as a full matrix
fulljac <- function  (t, y, parms) {
  jac <- matrix(nrow = 5, ncol = 5, byrow = TRUE,
                data = c(0.1, -0.2,  0  ,  0  ,  0  ,
                        -0.3,  0.1, -0.2,  0  ,  0  ,
                         0  , -0.3,  0.1, -0.2,  0  ,
                         0  ,  0  , -0.3,  0.1, -0.2,
                         0  ,  0  ,  0  , -0.3,  0.1))
  return(jac)
}

## the Jacobian, written in banded form
bandjac <- function  (t, y, parms) {
  jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE,
                data = c( 0  , -0.2, -0.2, -0.2, -0.2,
                          0.1,  0.1,  0.1,  0.1,  0.1,
                         -0.3, -0.3, -0.3, -0.3,    0))
  return(jac)
}

## initial conditions and output times
yini  <- 1:5
times <- 1:20

## default: stiff method, internally generated, full Jacobian
out   <- radau(yini, times, f1, parms = 0)
plot(out)

## stiff method, user-generated full Jacobian
out2  <- radau(yini, times, f1, parms = 0, jactype = "fullusr",
               jacfunc = fulljac)

## stiff method, internally-generated banded Jacobian
## one nonzero band above (up) and below(down) the diagonal
out3  <- radau(yini, times, f1, parms = 0, jactype = "bandint",
               bandup = 1, banddown = 1)

## stiff method, user-generated banded Jacobian
out4  <- radau(yini, times, f1, parms = 0, jactype = "bandusr",
               jacfunc = bandjac, bandup = 1, banddown = 1)


## =======================================================================
## Example 2: ODE
##   stiff problem from chemical kinetics
## =======================================================================
Chemistry <- function (t, y, p) {
  dy1 <- -.04*y[1] + 1.e4*y[2]*y[3]
  dy2 <- .04*y[1] - 1.e4*y[2]*y[3] - 3.e7*y[2]^2
  dy3 <- 3.e7*y[2]^2
  list(c(dy1, dy2, dy3))
}

times <- 10^(seq(0, 10, by = 0.1))
yini <- c(y1 = 1.0, y2 = 0, y3 = 0)

out <- radau(func = Chemistry, times = times, y = yini, parms = NULL)
plot(out, log = "x", type = "l", lwd = 2)

## =============================================================================
## Example 3: DAE
## Car axis problem, index 3 DAE, 8 differential, 2 algebraic equations
## from
## F. Mazzia and C. Magherini. Test Set for Initial Value Problem Solvers,
## release 2.4. Department
## of Mathematics, University of Bari and INdAM, Research Unit of Bari,
## February 2008.
## Available from https://archimede.uniba.it/~testset/
## =============================================================================

## Problem is written as M*y' = f(t,y,p).
## caraxisfun implements the right-hand side:

caraxisfun <- function(t, y, parms) {
  with(as.list(y), {
  
    yb <- r * sin(w * t)
    xb <- sqrt(L * L - yb * yb)
    Ll <- sqrt(xl^2 + yl^2)
    Lr <- sqrt((xr - xb)^2 + (yr - yb)^2)
        
    dxl <- ul; dyl <- vl; dxr <- ur; dyr <- vr
        
    dul  <- (L0-Ll) * xl/Ll      + 2 * lam2 * (xl-xr) + lam1*xb
    dvl  <- (L0-Ll) * yl/Ll      + 2 * lam2 * (yl-yr) + lam1*yb - k * g
               
    dur  <- (L0-Lr) * (xr-xb)/Lr - 2 * lam2 * (xl-xr)
    dvr  <- (L0-Lr) * (yr-yb)/Lr - 2 * lam2 * (yl-yr) - k * g
        
    c1   <- xb * xl + yb * yl
    c2   <- (xl - xr)^2 + (yl - yr)^2 - L * L
        
    list(c(dxl, dyl, dxr, dyr, dul, dvl, dur, dvr, c1, c2))
  })
}

eps <- 0.01; M <- 10; k <- M * eps^2/2; 
L <- 1; L0 <- 0.5; r <- 0.1; w <- 10; g <- 1

yini <- c(xl = 0, yl = L0, xr = L, yr = L0,
          ul = -L0/L, vl = 0,
          ur = -L0/L, vr = 0,
          lam1 = 0, lam2 = 0)

# the mass matrix
Mass      <- diag(nrow = 10, 1)
Mass[5,5] <- Mass[6,6] <- Mass[7,7] <- Mass[8,8] <- M * eps * eps/2
Mass[9,9] <- Mass[10,10] <- 0
Mass

# index of the variables: 4 of index 1, 4 of index 2, 2 of index 3
index <- c(4, 4, 2)

times <- seq(0, 3, by = 0.01)
out <- radau(y = yini, mass = Mass, times = times, func = caraxisfun,
        parms = NULL, nind = index)

plot(out, which = 1:4, type = "l", lwd = 2)

Explicit One-Step Solvers for Ordinary Differential Equations (ODE)

Description

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

The R function rk is a top-level function that provides interfaces to a collection of common explicit one-step solvers of the Runge-Kutta family with fixed or variable time steps.

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

Usage

rk(y, times, func, parms, rtol = 1e-6, atol = 1e-6,
  verbose = FALSE, tcrit = NULL, hmin = 0, hmax = NULL,
  hini = hmax, ynames = TRUE, method = rkMethod("rk45dp7", ... ),
  maxsteps = 5000, dllname = NULL, initfunc = dllname,
  initpar = parms, rpar = NULL, ipar = NULL,
  nout = 0, outnames = NULL, forcings = NULL,
  initforc = NULL, fcontrol = NULL, events = NULL, ...)

Arguments

Details

Function rk is a generalized implementation that can be used to evaluate different solvers of the Runge-Kutta family of explicit ODE solvers. A pre-defined set of common method parameters is in function rkMethod which also allows to supply user-defined Butcher tables.

The input parameters rtol, and atol determine the error control performed by the solver. The solver will control the vector 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.

Models can be defined in R as a user-supplied R-function, that must be called as: yprime = func(t, y, parms). t is the current time point in the integration, y is the current estimate of the variables in the ODE system.

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 second element contains output variables that are required at each point in time. Examples are given below.

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

Note

Arguments rpar and ipar are provided for compatibility with lsoda.

Starting with version 1.8 implicit Runge-Kutta methods are also supported by this general rk interface, however their implementation is still experimental. Instead of this you may consider radau for a specific full implementation of an implicit Runge-Kutta method.

Author(s)

Thomas Petzoldt thomas.petzoldt@tu-dresden.de

References

Butcher, J. C. (1987) The numerical analysis of ordinary differential equations, Runge-Kutta and general linear methods, Wiley, Chichester and New York.

Engeln-Muellges, G. and Reutter, F. (1996) Numerik Algorithmen: Entscheidungshilfe zur Auswahl und Nutzung. VDI Verlag, Duesseldorf.

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.

Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (2007) Numerical Recipes in C. Cambridge University Press.

See Also

For most practical cases, solvers of the Livermore family (i.e. the ODEPACK solvers, see below) are superior. Some of them are also suitable for stiff ODEs, differential algebraic equations (DAEs), or partial differential equations (PDEs).

• rkMethod for a list of available Runge-Kutta parameter sets,

• rk4 and euler for special versions without interpolation (and less overhead),

• lsoda, lsode, lsodes, lsodar, vode, daspk for 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: Resource-producer-consumer Lotka-Volterra model
## =======================================================================

## Notes:
## - Parameters are a list, names accessible via "with" function
## - Function sigimp passed as an argument (input) to model
##  (see also ode and lsoda examples)

SPCmod <- function(t, x, parms, input)  {
  with(as.list(c(parms, x)), {
    import <- input(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)
  })
}

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

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

## external signal with rectangle impulse
signal <- data.frame(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
xstart <- c(S = 1, P = 1, C = 1)

## Euler method
out1  <- rk(xstart, times, SPCmod, parms, hini = 0.1, 
            input = sigimp, method = "euler")

## classical Runge-Kutta 4th order
out2 <- rk(xstart, times, SPCmod, parms, hini = 1, 
           input = sigimp, method = "rk4")

## Dormand-Prince method of order 5(4)
out3 <- rk(xstart, times, SPCmod, parms, hmax = 1, 
           input = sigimp, method = "rk45dp7")

mf <- par("mfrow")
## deSolve plot method for comparing scenarios
plot(out1, out2, out3, which = c("S", "P", "C"), 
     main = c ("Substrate", "Producer", "Consumer"),
     col =c("black", "red", "green"), 
     lty = c("solid", "dotted", "dotted"), lwd = c(1, 2, 1))

## user-specified plot function
plot (out1[,"P"], out1[,"C"], type = "l", xlab = "Producer", ylab = "Consumer")
lines(out2[,"P"], out2[,"C"], col = "red",   lty = "dotted", lwd = 2)
lines(out3[,"P"], out3[,"C"], col = "green", lty = "dotted")

legend("center", legend = c("euler", "rk4", "rk45dp7"),
  lty = c(1, 3, 3), lwd = c(1, 2, 1),
  col = c("black", "red", "green"))
par(mfrow = mf)

Solve System of ODE (Ordinary Differential Equation)s by Euler's Method or Classical Runge-Kutta 4th Order Integration.

Description

Solving initial value problems for systems of first-order ordinary differential equations (ODEs) using Euler's method or the classical Runge-Kutta 4th order integration.

Usage

euler(y, times, func, parms, verbose = FALSE, ynames = TRUE,
  dllname = NULL, initfunc = dllname, initpar = parms,
  rpar = NULL, ipar = NULL, nout = 0, outnames = NULL,
  forcings = NULL, initforc = NULL, fcontrol = NULL, ...)

rk4(y, times, func, parms, verbose = FALSE, ynames = TRUE,
  dllname = NULL, initfunc = dllname, initpar = parms,
  rpar = NULL, ipar = NULL, nout = 0, outnames = NULL,
  forcings = NULL, initforc = NULL, fcontrol = NULL, ...)

euler.1D(y, times, func, parms, nspec = NULL, dimens = NULL,
  names = NULL, verbose = FALSE, ynames = TRUE,
  dllname = NULL, initfunc = dllname, initpar = parms,
  rpar = NULL,  ipar = NULL, nout = 0, outnames = NULL,
  forcings = NULL, initforc = NULL, fcontrol = NULL, ...)

Arguments

Details

rk4 and euler are special versions of the two fixed step solvers with less overhead and less functionality (e.g. no interpolation and no events) compared to the generic Runge-Kutta codes called by ode resp. rk.

If you need different internal and external time steps or want to use events, please use: rk(y, times, func, parms, method = "rk4") or rk(y, times, func, parms, method = "euler").

See help pages of rk and rkMethod for details.

Function euler.1D essentially calls function euler but contains additional code to support plotting of 1D models, see ode.1D and plot.1D for details.

Value

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

Note

For most practical cases, solvers with flexible timestep (e.g. rk(method = "ode45") and especially solvers of the Livermore family (ODEPACK, e.g. lsoda) are superior.

Author(s)

Thomas Petzoldt thomas.petzoldt@tu-dresden.de

See Also

• rkMethod for a list of available Runge-Kutta parameter sets,

• rk for the more general Runge-Code,

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

• dede for integrating models with delay differential equations,

diagnostics to print diagnostic messages.

Examples

## =======================================================================
## Example: Analytical and numerical solutions of logistic growth
## =======================================================================

## the derivative of the logistic
logist <- function(t, x, parms) {
  with(as.list(parms), {
    dx <- r * x[1] * (1 - x[1]/K)
    list(dx)
  })
}

time  <- 0:100
N0    <- 0.1; r <- 0.5; K <- 100
parms <- c(r = r, K = K)
x <- c(N = N0)

## analytical solution
plot(time, K/(1 + (K/N0-1) * exp(-r*time)), ylim = c(0, 120),
  type = "l", col = "red", lwd = 2)

## reasonable numerical solution with rk4
time <- seq(0, 100, 2)
out <- as.data.frame(rk4(x, time, logist, parms))
points(out$time, out$N, pch = 16, col = "blue", cex = 0.5)

## same time step with euler, systematic under-estimation
time <- seq(0, 100, 2)
out <- as.data.frame(euler(x, time, logist, parms))
points(out$time, out$N, pch = 1)

## unstable result
time <- seq(0, 100, 4)
out <- as.data.frame(euler(x, time, logist, parms))
points(out$time, out$N, pch = 8, cex = 0.5)

## method with automatic time step
out <- as.data.frame(lsoda(x, time, logist, parms))
points(out$time, out$N, pch = 1, col = "green")

legend("bottomright",
  c("analytical","rk4, h=2", "euler, h=2",
    "euler, h=4", "lsoda"),
  lty = c(1, NA, NA, NA, NA), lwd = c(2, 1, 1, 1, 1),
  pch = c(NA, 16, 1, 8, 1),
  col = c("red", "blue", "black", "black", "green"))

Collection of Parameter Sets (Butcher Arrays) for the Runge-Kutta Family of ODE Solvers

Description

This function returns a list specifying coefficients and properties of ODE solver methods from the Runge-Kutta family.

Usage

rkMethod(method = NULL, ...)

Arguments

Details

This function supplies method settings for rk or ode. If called without arguments, the names of all currently implemented solvers of the Runge-Kutta family are returned.

The following comparison gives an idea how the algorithms of deSolve are related to similar algorithms of other simulation languages:

Note that this table is based on the Runge-Kutta coefficients only, but the algorithms differ also in their implementation, in their stepsize adaption strategy and interpolation methods.

The table reflects the state at time of writing and it is of course possible that implementations change.

Methods "rk45dp7" (alias "ode45") and "rk45ck" contain specific and efficient built-in interpolation schemes (dense output).

As an alternative, Neville-Aitken polynomials can be used to interpolate between time steps. This is available for all RK methods and may be useful to speed up computation if no dense-output formula is available. Note however, that this can introduce considerable local error; it is disabled by default (see nknots below).

Value

A list with the following elements:

Note

• Adaptive stepsize Runge-Kuttas are preferred if the solution contains parts when the states change fast, and parts when not much happens. They will take small steps over bumpy ground and long steps over uninteresting terrain.

• As a suggestion, one may use "rk23" (alias "ode23") for simple problems and "rk45dp7" (alias "ode45") for rough problems. The default solver is "rk45dp7" (alias "ode45"), because of its relatively high order (4), re-use of the last intermediate steps (FSAL = first same as last) and built-in polynomial interpolation (dense output).

• Solver "rk23bs", that supports also FSAL, may be useful for slightly stiff systems if demands on precision are relatively low.

• Another good choice, assuring medium accuracy, is the Cash-Karp Runge-Kutta method, "rk45ck".

• Classical "rk4" is traditionally used in cases where an adequate stepsize is known a-priori or if external forcing data are provided for fixed time steps only and frequent interpolation of external data needs to be avoided.

• Method "rk45dp7" (alias "ode45") contains an efficient built-in interpolation scheme (dense output) based on intermediate function evaluations.

Starting with version 1.8 implicit Runge-Kutta (irk) methods are also supported by the general rk interface, however their implementation is still experimental. Instead of this you may consider radau for a specific full implementation of an implicit Runge-Kutta method.

Author(s)

Thomas Petzoldt thomas.petzoldt@tu-dresden.de

References

Bogacki, P. and Shampine L.F. (1989) A 3(2) pair of Runge-Kutta formulas, Appl. Math. Lett. 2, 1–9.

Butcher, J. C. (1987) The numerical analysis of ordinary differential equations, Runge-Kutta and general linear methods, Wiley, Chichester and New York.

Cash, J. R. and Karp A. H., 1990. A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides, ACM Transactions on Mathematical Software 16, 201–222. doi:10.1145/79505.79507

Dormand, J. R. and Prince, P. J. (1980) A family of embedded Runge-Kutta formulae, J. Comput. Appl. Math. 6(1), 19–26.

Engeln-Muellges, G. and Reutter, F. (1996) Numerik Algorithmen: Entscheidungshilfe zur Auswahl und Nutzung. VDI Verlag, Duesseldorf.

Fehlberg, E. (1967) Klassische Runge-Kutta-Formeln fuenfter and siebenter Ordnung mit Schrittweiten-Kontrolle, Computing (Arch. Elektron. Rechnen) 4, 93–106.

Kutta, W. (1901) Beitrag zur naeherungsweisen Integration totaler Differentialgleichungen, Z. Math. Phys. 46, 435–453.

Octave-Forge - Extra Packages for GNU Octave, Package OdePkg. https://octave.sourceforge.io

Prince, P. J. and Dormand, J. R. (1981) High order embedded Runge-Kutta formulae, J. Comput. Appl. Math. 7(1), 67–75. doi:10.1016/0771-050X(81)90010-3

Runge, C. (1895) Ueber die numerische Aufloesung von Differentialgleichungen, Math. Ann. 46, 167–178.

MATLAB (R) is a registed property of The Mathworks Inc. https://www.mathworks.com/

See Also

rk, ode

Examples

rkMethod()          # returns the names of all available methods
rkMethod("rk45dp7") # parameters of the Dormand-Prince 5(4) method
rkMethod("ode45")   # an alias for the same method

func <- function(t, x, parms) {
  with(as.list(c(parms, x)),{
    dP  <- a * P      - b * C * P
    dC  <- b * P * C  - c * C
    res <- c(dP, dC)
    list(res)
  })
}
times  <- seq(0, 200, length = 101)
parms  <- c(a = 0.1, b = 0.1, c = 0.1)
x <- c(P = 2, C = 1)

## rk using ode45 as the default method
out <- rk(x, times, func, parms)

## all methods can be called also from 'ode' by using rkMethod
out <- ode(x, times, func, parms, method = rkMethod("rk4"))

## 'ode' has aliases for the most common RK methods
out <- ode(x, times, func, parms, method = "ode45")

##===========================================================================
## Comparison of local error from different interpolation methods
##===========================================================================

## lsoda with lower tolerances (1e-10) used as reference
o0 <- ode(x, times, func, parms, method = "lsoda", atol = 1e-10, rtol = 1e-10)

## rk45dp7 with hmax = 10 > delta_t  = 2
o1 <- ode(x, times, func, parms, method = rkMethod("rk45dp7"), hmax = 10)

## disable dense-output interpolation
## and use only Neville-Aitken polynomials instead
o2 <- ode(x, times, func, parms,
  method = rkMethod("rk45dp7", densetype = NULL, nknots = 5), hmax = 10)

## stop and go: disable interpolation completely
## and integrate explicitly between external time steps
o3 <- ode(x, times, func, parms,
  method = rkMethod("rk45dp7", densetype = NULL, nknots = 0, hmax=10))

## compare different interpolation methods with lsoda
mf <- par("mfrow" = c(4, 1))
matplot(o1[,1], o1[,-1], type = "l", xlab = "Time", main = "State Variables",
  ylab = "P, C")
matplot(o0[,1], o0[,-1] - o1[,-1], type = "l", xlab = "Time", ylab = "Diff.",
  main="Difference between lsoda and ode45 with dense output")
abline(h = 0, col = "grey")
matplot(o0[,1], o0[,-1] - o2[,-1], type = "l", xlab = "Time", ylab = "Diff.",
  main="Difference between lsoda and ode45 with Neville-Aitken")
abline(h = 0, col = "grey")
matplot(o0[,1], o0[,-1] - o3[,-1], type = "l", xlab = "Time", ylab = "Diff.",
  main="Difference between lsoda and ode45 in 'stop and go' mode")
abline(h = 0, col = "grey")
par(mf)

##===========================================================================
## rkMethod allows to define user-specified Runge-Kutta methods
##===========================================================================
out <- ode(x, times, func, parms,
          method = rkMethod(ID = "midpoint",
            varstep = FALSE,
            A      = c(0, 1/2),
            b1      = c(0, 1),
            c       = c(0, 1/2),
            stage   = 2,
            Qerr    = 1
          )
        )
plot(out)

## compare method diagnostics
times  <- seq(0, 200, length = 10)
o1 <- ode(x, times, func, parms, method = rkMethod("rk45ck"))
o2 <- ode(x, times, func, parms, method = rkMethod("rk78dp"))
diagnostics(o1)
diagnostics(o2)

A Sediment Model of Oxygen Consumption

Description

A model that describes oxygen consumption in a marine sediment.

One state variable:

• sedimentary organic carbon,

Organic carbon settles on the sediment surface (forcing function Flux) and decays at a constant rate.

The equation is simple:

\frac{dC}{dt} = Flux - k C

This model is written in FORTRAN.

Usage

SCOC(times, y = NULL, parms, Flux, ...)

Arguments

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 <karline.soetaert@nioz.nl>

References

Soetaert, K. and P.M.J. Herman, 2009. A Practical Guide to Ecological Modelling. Using R as a Simulation Platform. Springer, 372 pp.

See Also

ccl4model, the CCl4 inhalation model.

aquaphy, the algal growth model.

Examples

## 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
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")

Time Lagged Values of State Variables and Derivatives.

Description

Functions lagvalue and lagderiv provide access to past (lagged) values of state variables and derivatives.

They are to be used with function dede, to solve delay differential equations.

Usage

lagvalue(t, nr)
lagderiv(t, nr)

Arguments

Details

The lagvalue and lagderiv can only be called during the integration, the lagged time should not be smaller than the initial simulation time, nor should it be larger than the current simulation time.

Cubic Hermite interpolation is used to obtain an accurate interpolant at the requested lagged time.

Value

a scalar (or vector) with the lagged value(s).

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

See Also

dede, for how to implement delay differential equations.

Examples

## =============================================================================
## exercise 6 from Shampine and Thompson, 2000
## solving delay differential equations with dde23
##
## two lag values
## =============================================================================

##-----------------------------
## the derivative function
##-----------------------------
derivs <- function(t, y, parms) { 
  History <- function(t) c(cos(t), sin(t))
  if (t < 1)
    lag1 <- History(t - 1)[1]    
  else 
    lag1 <- lagvalue(t - 1)[1] # returns a vector; select first element

  if (t < 2)
    lag2 <- History(t - 2)[2]
  else 
    lag2 <- lagvalue(t - 2,2) # faster than lagvalue(t - 2)[2]

  dy1 <- lag1 * lag2
  dy2 <- -y[1] * lag2
  
  list(c(dy1, dy2), lag1 = lag1, lag2 = lag2)
}

##-----------------------------
## parameters
##-----------------------------

r <- 3.5; m <- 19

##-----------------------------
## initial values and times
##-----------------------------

yinit <- c(y1 = 0, y2 = 0)
times <- seq(0, 20, by = 0.01)

##-----------------------------
## solve the model  
##-----------------------------

yout <- dede(y = yinit, times = times, func = derivs,
  parms = NULL, atol = 1e-9)

##-----------------------------
## plot results
##-----------------------------

plot(yout, type = "l", lwd = 2)

## =============================================================================
## The predator-prey model with time lags, from Hale
## problem 1 from Shampine and Thompson, 2000
## solving delay differential equations with dde23
##
## a vector with lag valuess
## =============================================================================

##-----------------------------
## the derivative function
##-----------------------------
predprey <- function(t, y, parms) {
  tlag <- t - 1
  if (tlag < 0)
    ylag <- c(80, 30)
  else 
    ylag <- lagvalue(tlag)  # returns a vector
  
  dy1 <- a * y[1] * (1 - y[1]/m) + b * y[1] * y[2]
  dy2 <- c * y[2] + d * ylag[1] * ylag[2]
  list(c(dy1, dy2))
}

##-----------------------------
## parameters
##-----------------------------

a <- 0.25; b <- -0.01; c <- -1 ; d <- 0.01; m <- 200

##-----------------------------
## initial values and times
##-----------------------------

yinit <- c(y1 = 80, y2 = 30)
times <- seq(0, 100, by = 0.01)

#-----------------------------
# solve the model  
#-----------------------------

yout <- dede(y = yinit, times = times, func = predprey, parms = NULL)

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

plot(yout, type = "l", lwd = 2, main = "Predator-prey model", mfrow = c(2, 2))
plot(yout[,2], yout[,3], xlab = "y1", ylab = "y2", type = "l", lwd = 2)

diagnostics(yout)

## =============================================================================
##
## A neutral delay differential equation (lagged derivative)   
##  y't = -y'(t-1), y(t) t < 0 = 1/t
##
## =============================================================================

##-----------------------------
## the derivative function
##-----------------------------
derivs <- function(t, y, parms) {
  tlag <- t - 1
  if (tlag < 0)
    dylag <- -1
  else
    dylag <- lagderiv(tlag)

  list(c(dy = -dylag), dylag = dylag)
}

##-----------------------------
## initial values and times
##-----------------------------

yinit <- 0
times <- seq(0, 4, 0.001)

##-----------------------------
## solve the model  
##-----------------------------

yout <- dede(y = yinit, times = times, func = derivs, parms = NULL)

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

plot(yout, type = "l", lwd = 2)

Solver for Ordinary Differential Equations (ODE)

Description

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

dy/dt = f(t,y)

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

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

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

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

Usage

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

Arguments

Details

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

If the problem is nonstiff, use method flag mf = 10, which selects a nonstiff (Adams) method, no Jacobian used.

If the problem is stiff, there are four standard choices which can be specified with jactype or mf.

The options for jactype are

jac = "fullint":

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

jac = "fullusr":

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

jac = "bandusr":

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

jac = "bandint":

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

More options are available when specifying mf directly.

The legal values of mf are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25, -11, -12, -14, -15, -21, -22, -24, -25.

mf is a signed two-digit integer, mf = JSV*(10*METH + MITER), where

JSV = SIGN(mf)

indicates the Jacobian-saving strategy: JSV = 1 means a copy of the Jacobian is saved for reuse in the corrector iteration algorithm. JSV = -1 means a copy of the Jacobian is not saved.

METH

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

MITER

indicates the corrector iteration method: MITER = 0 means functional iteration (no Jacobian matrix is involved).

MITER = 1 means chord iteration with a user-supplied full (NEQ by NEQ) Jacobian.

MITER = 2 means chord iteration with an internally generated (difference quotient) full Jacobian (using NEQ extra calls to func per df/dy value).

MITER = 3 means chord iteration with an internally generated diagonal Jacobian approximation (using 1 extra call to func per df/dy evaluation).

MITER = 4 means chord iteration with a user-supplied banded Jacobian.

MITER = 5 means chord iteration with an internally generated banded Jacobian (using ML+MU+1 extra calls to func per df/dy evaluation).

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

The example for integrator lsode demonstrates how to specify both a banded and full Jacobian.

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, vode will return an error code. See lsoda for details.

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

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

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

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

Note

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

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

References

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

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

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

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

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

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

Netlib: https://netlib.org

See Also

• rk,

• rk4 and euler for Runge-Kutta integrators.

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

## =======================================================================
## ex. 1
## 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, 100, 0.01)

out   <- vode(state, times, chaos, 0)

plot(out, type = "l")   # all versus time
plot(out[,"x"], out[,"y"], type = "l", main = "Lorenz butterfly",
  xlab = "x", ylab = "y")


## =======================================================================
## ex. 2
## SCOC model, in FORTRAN  - to see the FORTRAN code:
## browseURL(paste(system.file(package="deSolve"),
##                             "/doc/examples/dynload/scoc.f",sep=""))
## example from Soetaert and Herman, 2009, chapter 3. (simplified)
## =======================================================================

## 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)

meanDepo <- mean(approx(Flux[,1], Flux[,2], xout = seq(1, 365, by = 1))$y)

Yini <- c(y = as.double(meanDepo/parms))

times <- 1:365
out <- vode(Yini, times, func = "scocder",
    parms = parms, dllname = "deSolve",
    initforc = "scocforc", forcings = Flux,
    initfunc = "scocpar", nout = 2,
    outnames = c("Mineralisation", "Depo"))

matplot(out[,1], out[,c("Depo", "Mineralisation")], 
        type = "l", col = c("red", "blue"), xlab = "time", ylab = "Depo")

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

out2 <- vode(Yini, times, func = "scocder",
    parms = parms, dllname = "deSolve",
    initforc = "scocforc",  forcings = Flux, fcontrol = fcontrol,
    initfunc = "scocpar", nout = 2,
    outnames = c("Mineralisation", "Depo"))
matplot(out2[,1], out2[,c("Depo", "Mineralisation")], 
        type = "l", col = c("red", "blue"), xlab = "time", ylab = "Depo")

## Constant interpolation of forcing function - middle of interval
fcontrol <- list(method = "constant", f = 0.5)

out3 <- vode(Yini, times, func = "scocder",
    parms = parms, dllname = "deSolve",
    initforc = "scocforc",  forcings = Flux, fcontrol = fcontrol,
    initfunc = "scocpar", nout = 2,
    outnames = c("Mineralisation", "Depo"))

matplot(out3[,1], out3[,c("Depo", "Mineralisation")], 
        type = "l", col = c("red", "blue"), xlab = "time", ylab = "Depo")

plot(out, out2, out3)

Solver for Ordinary Differential Equations (ODE) for COMPLEX variables

Description

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

dy/dt = f(t,y)

where dy and y are complex variables.

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

Usage

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

Arguments

Details

see vode, the double precision version, for details.

Value

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

Note

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

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

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

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

References

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

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

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

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

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

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

Netlib: https://netlib.org

See Also

vode for the double precision version

Examples

## =======================================================================
## Example 1 - very simple example 
## df/dt = 1i*f, where 1i is the imaginary unit
## The initial value is f(0) = 1 = 1+0i
## =======================================================================

ZODE <- function(Time, f, Pars) {
  df <-  1i*f
  return(list(df))
}

pars    <- NULL
yini    <- c(f = 1+0i)
times   <- seq(0, 2*pi, length = 100)
out     <- zvode(func = ZODE, y = yini, parms = pars, times = times,
  atol = 1e-10, rtol = 1e-10)

# The analytical solution to this ODE is the exp-function:
# f(t) = exp(1i*t)
#      = cos(t)+1i*sin(t)  (due to Euler's equation)

analytical.solution  <- exp(1i * times) 

## compare numerical and analytical solution
tail(cbind(out[,2], analytical.solution))


## =======================================================================
## Example 2 - example in "zvode.f",  
## df/dt = 1i*f        (same as above ODE)
## dg/dt = -1i*g*g*f   (an additional ODE depending on f)
##
## Initial values are
## g(0) = 1/2.1 and
## z(0) = 1  
## =======================================================================

ZODE2<-function(Time,State,Pars) {
  with(as.list(State), {
    df <- 1i * f
    dg <- -1i * g*g * f
    return(list(c(df, dg)))
  })
}

yini    <- c(f = 1 + 0i, g = 1/2.1 + 0i)
times   <- seq(0, 2*pi, length = 100)
out     <- zvode(func = ZODE2, y = yini, parms = NULL, times = times,
  atol = 1e-10, rtol = 1e-10)


## The analytical solution is
## f(t) = exp(1i*t)   (same as above)
## g(t) = 1/(f(t) + 1.1)

analytical <- cbind(f = exp(1i * times), g = 1/(exp(1i * times) + 1.1))

## compare numerical solution and the two analytical ones:
tail(cbind(out[,2], analytical[,1]))