deSolve (version 1.8.1)

rkMethod: 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

method

a string constant naming one of the pre-defined methods of the Runge-Kutta family of solvers. The most common methods are the fixed-step methods "euler", "rk2", "rk4" or the variable step methods "

...

specification of a user-defined solver, see Value and example below.

Value

• A list with the following elements:
• ID
  name of the method (character)
• varstep
  boolean value specifying if the method allows for variable time step (TRUE) or not (FALSE).
• FSAL
  (first same as last) optional boolean value specifying if the method allows re-use of the last function evaluation (TRUE) or not (FALSE or NULL).
• A
  coefficient matrix of the method. As link{rk} supports only explicit methods, this matrix must be lower triangular. A must be a vector for fixed step methods where only the subdiagonal values are different from zero.
• b1
  coefficients of the lower order Runge-Kutta pair.
• b2
  coefficients of the higher order Runge-Kutta pair (optional, for embedded methods that allow variable time step).
• c
  coefficients for calculating the intermediate time steps.
• d
  optional coefficients for built-in polynomial interpolation of the outputs from internal steps (dense output), currently only available for method rk45dp7 (Dormand-Prince).
• densetype
  optional integer value specifying the dense output formula; currently densetype = 1 for rk45dp7 (Dormand-Prince) and densetype = 2 for rk45ck (Cash-Karp). Undefined values (e.g., densetype = NULL) disable dense output.
• stage
  number of function evaluations needed (corresponds to number of rows in A).
• Qerr
  global error order of the method, important for automatic time-step adjustment.
• nknots
  integer value specifying the order of interpolation polynomials for methods without dense output. If nknots < 2 (the default) then internal interpolation is switched off and integration is performed step by step between external time steps.

  If nknots is between 3 and 8, Neville-Aitken polynomials are used, which need at least nknots + 1 internal time steps. Local interpolation may speed up integration but can lead to local errors higher than the tolerance, especially if external and internal time steps are very different.

• alpha
  optional tuning parameter for stepsize adjustment. If alpha is omitted, it is set to $1/Qerr - 0.75 beta$. The default value is $1/Qerr$ (for beta = 0).
• beta
  optional tuning parameter for stepsize adjustment. Typical values are $0$ (default) or $0.4/Qerr$

Details

This function supplies method settings for rk or ode. If called without arguments, the names of all 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: lll{ rkMethod | Description "euler" | Euler's Method "rk2" | 2nd order Runge-Kutta, fixed time step (Heun's method) "rk4" | classical 4th order Runge-Kutta, fixed time step "rk23" | Runge-Kutta, order 2(3), Octave: ode23 "rk23bs", "ode23" | Bogacki-Shampine, order 2(3), Matlab: ode23 "rk34f" | Runge-Kutta-Fehlberg, order 3(4) "rk45ck" | Runge-Kutta Cash-Karp, order 4(5) "rk45f" | Runge-Kutta-Fehlberg, order 4(5), Octave: ode45, pair=1 "rk45e" | Runge-Kutta-England, order 4(5) "rk45dp6" | Dormand-Prince, order 4(5), local order 6 "rk45dp7", "ode45" | Dormand-Prince 4(5), local order 7 | (also known as dopri5, MATLAB: ode45, Octave: ode45, pair=0) "rk78f" | Runge-Kutta-Fehlberg, order 7(8) "rk78dp" | Dormand-Prince, order 7(8) } 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 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).

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.

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. http://octave.sourceforge.net/odepkg

Prince, P. J. and Dormand, J. R. (1981) High order embedded Runge-Kutta formulae, J. Comput. Appl. Math. 7(1), 67--75.

Runge, C. (1895) Ueber die numerische Aufloesung von Differentialgleichungen, Math. Ann. 46, 167--178. MATLAB (R) is a registed property of The Mathworks Inc. http://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)

ode(x, times, func, parms, method = rkMethod("rk4"))

ode(x, times, func, parms, method = "ode45")

o0 <- ode(x, times, func, parms, method = "lsoda")
o1 <- ode(x, times, func, parms, method = rkMethod("rk45dp7"))

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

## show differences between methods
par(mfrow=c(3,1))
matplot(o1[,1], o1[,-1], type="l", xlab="Time", main="State Variables")
matplot(o1[,1], o2[,-1], type="p", pch=16, add=TRUE)
matplot(o0[,1], o1[,-1]-o0[,-1], type="l", lty="solid", xlab="Time",
  main="Difference: rk45dp7 - lsoda")
matplot(o1[,1], o2[,-1]-o1[,-1], type="l", lty="dashed", xlab="Time",
  main="difference: Neville-Aitken - Dense Output")
abline(h=0, col="grey")

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

## define and use a new rk method
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
  )
)

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