dae: General Solver for Differential Algebraic Equations

Description Usage Arguments Details Value Author(s) See Also Examples

View source: R/dae.R

Description

Solves a system of differential algebraic equations; a wrapper around the implemented DAE solvers

Usage

1
2
dae(y, times, parms, dy, res = NULL, func = NULL,
    method = c("mebdfi", "daspk", "radau", "gamd", "bimd"), ...)

Arguments

y

the initial (state) values for the DAE system, a vector. If y has a name attribute, the names will be used to label the output matrix.

times

time sequence for which output is wanted; the first value of times must be the initial time.

parms

vector or list of parameters used in res

dy

the initial derivatives of the state variables of the DAE system.

func

to be used if the model is an ODE, or a DAE written in linearly implicit form (M y' = f(t, y)). func should be an R-function that computes the values of the derivatives in the ODE system (the model definition) at time t.

func must be defined as: func <- function(t, y, parms,...).
t is the current time point in the integration, y is the current estimate of the variables in the ODE system. If the initial values y has a names attribute, the names will be available inside func, unless ynames is FALSE. parms is a vector or list of parameters. ... (optional) are any other arguments passed to the function.

The return value of func should be a list, whose first element is a vector containing the derivatives of y with respect to time, and whose next elements are global values that are required at each point in times. The derivatives should be specified in the same order as the specification of the state variables y.

res

either an R-function that computes the residual function F(t,y,y') of the DAE system (the model defininition) at time t, or a character string giving the name of a compiled function in a dynamically loaded shared library.

If res is a user-supplied R-function, it must be defined as: res <- function(t, y, dy, parms, ...).

Here t is the current time point in the integration, y is the current estimate of the variables in the DAE system, dy are the corresponding derivatives. If the initial y or dy have a names attribute, the names will be available inside res, unless ynames is FALSE. parms is a vector of parameters.

The return value of res should be a list, whose first element is a vector containing the residuals of the DAE system, i.e. delta = F(t,y,y'), and whose next elements contain output variables that are required at each point in times.

If res is a string, then dllname must give the name of the shared library (without extension) which must be loaded before dae() is called (see deSolve package vignette "compiledCode" for more information).

method

the solver to use, either a string ("mebdfi", "daspk"), "radau", "gamd" or a function that performs integration.

...

additional arguments passed to the solvers.

Details

This is simply a wrapper around the various dae solvers.

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

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

Author(s)

Karline Soetaert <[email protected]>

See Also

diagnostics to print diagnostic messages.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
## =======================================================================
## Chemical problem
## =======================================================================

daefun <- function(t, y, dy, parms) {
  with (as.list(c(y, dy, parms)), {
    res1 <- dA + dAB + lambda * A
    res2 <- dAB + dB
    alg <- B * A - K * AB
  list(c(res1, res2, alg), sumA = A + AB)
  })
}

parms <- c(lambda = 0.1, K = 1e-4)

yini  <- with(as.list(parms),
     c(A = 1, AB = 1, B = K))
dyini <- c(dA = 0, dAB = 0, dB = 0)

times <- 0:100

print(system.time(
out <-dae (y=yini, dy = dyini, times = times, res = daefun,
           parms = parms, method = "daspk")
))

plot(out, ylab = "conc.", xlab = "time", lwd = 2)
mtext("IVP DAE", side = 3, outer = TRUE, line = -1)

deTestSet documentation built on May 29, 2017, 2:56 p.m.