rk4 | R Documentation |
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.
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, ...)
y |
the initial (state) values for the ODE system. If |
times |
times at which explicit estimates for |
func |
either an R-function that computes the values of the derivatives in the ODE system (the model definition) at time t, or a character string giving the name of a compiled function in a dynamically loaded shared library. If The return value of If |
parms |
vector or list of parameters used in |
nspec |
for 1D models only: the number of species (components)
in the model. If |
dimens |
for 1D models only: the number of boxes in the
model. If |
names |
for 1D models only: the names of the components; used for plotting. |
verbose |
a logical value that, when |
ynames |
if |
dllname |
a string giving the name of the shared library
(without extension) that contains all the compiled function or
subroutine definitions refered to in |
initfunc |
if not |
initpar |
only when ‘dllname’ is specified and an
initialisation function |
rpar |
only when ‘dllname’ is specified: a vector with
double precision values passed to the DLL-functions whose names are
specified by |
ipar |
only when ‘dllname’ is specified: a vector with
integer values passed to the dll-functions whose names are specified
by |
nout |
only used if |
outnames |
only used if ‘dllname’ is specified and
|
forcings |
only used if ‘dllname’ is specified: a list with
the forcing function data sets, each present as a two-columned matrix,
with (time, value); interpolation outside the interval
[min( See forcings or package vignette |
initforc |
if not |
fcontrol |
A list of control parameters for the forcing functions.
See forcings or vignette |
... |
additional arguments passed to |
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.
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.
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.
Thomas Petzoldt thomas.petzoldt@tu-dresden.de
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.
## =======================================================================
## 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"))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.