Implicit RungeKutta 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 3stage
RADAU IIA method.
It implements the implicit RungeKutta 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
1 2 3 4 5 6 7 8 9  radau(y, times, func, parms, nind = c(length(y), 0, 0),
rtol = 1e6, atol = 1e6, 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
y 
the initial (state) values for the ODE system. If 
times 
time sequence for which output is wanted; the first
value of 
func 
either an Rfunction that computes the values of the derivatives in the ODE system (the model definition) at time t, or the righthand side of the equation M dy/dt = f(t,y) if a DAE. (if
If
The return value of If 
parms 
vector or list of parameters used in 
nind 
if a DAE system: a threevalued vector with the number of variables of index 1, 2, 3 respectively. The equations must be defined such that the index 1 variables precede the index 2 variables which in turn precede the index 3 variables. The sum of the variables of different index should equal N, the total number of variables. This has implications on the scaling of the variables, i.e. index 2 variables are scaled by 1/h, index 3 variables are scaled by 1/h^2. 
rtol 
relative error tolerance, either a
scalar or an array as long as 
atol 
absolute error tolerance, either a scalar or an array as
long as 
jacfunc 
if not In some circumstances, supplying
If the Jacobian is a full matrix,

jactype 
the structure of the Jacobian, one of

mass 
the mass matrix.
If not If 
massup 
number of nonzero bands above the diagonal of the 
massdown 
number of nonzero bands below the diagonal of the 
rootfunc 
if not 
verbose 
if 
nroot 
only used if ‘dllname’ is specified: the number of
constraint functions whose roots are desired during the integration;
if 
hmax 
an optional maximum value of the integration stepsize. If
not specified, 
hini 
initial step size to be attempted; if 0, the initial step size is set equal to 1e6. Usually 1e3 to 1e5 is good for stiff equations 
ynames 
logical, if 
bandup 
number of nonzero bands above the diagonal, in case the Jacobian is banded. 
banddown 
number of nonzero bands below the diagonal, in case the Jacobian is banded. 
maxsteps 
average maximal number of steps per output interval
taken by the solver. This argument is defined such as to ensure
compatibility with the Livermoresolvers. RADAU only accepts the maximal
number of steps for the entire integration, and this is calculated
as 
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 dllfunctions whose names are
specified by 
ipar 
only when ‘dllname’ is specified: a vector with
integer values passed to the dllfunctions 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 twocolumned 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 
events 
A list that specifies events, i.e. when the value of a state variable is suddenly changed. See events for more information. 
lags 
A list that specifies timelags, i.e. the number of steps that has to be kept. To be used for delay differential equations. See timelags, dede for more information. 
... 
additional arguments passed to 
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 bybandup
andbanddown
. jactype = "bandint"
a banded Jacobian, calculated by radau; the size of the bands specified by
bandup
andbanddown
.
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 Rfunction. 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 Differentialalgebraic problems. Springer series in computational mathematics 14, SpringerVerlag, second edition.
See Also

ode
for a general interface to most of the ODE solvers , 
ode.1D
for integrating 1D models, 
ode.2D
for integrating 2D models, 
ode.3D
for integrating 3D models, 
daspk
for integrating DAE models up to index 1
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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136  ## =======================================================================
## 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, usergenerated full Jacobian
out2 < radau(yini, times, f1, parms = 0, jactype = "fullusr",
jacfunc = fulljac)
## stiff method, internallygenerated 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, usergenerated 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 at http://www.dm.uniba.it/~testset.
## =============================================================================
## Problem is written as M*y' = f(t,y,p).
## caraxisfun implements the righthand 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 < (L0Ll) * xl/Ll + 2 * lam2 * (xlxr) + lam1*xb
dvl < (L0Ll) * yl/Ll + 2 * lam2 * (ylyr) + lam1*yb  k * g
dur < (L0Lr) * (xrxb)/Lr  2 * lam2 * (xlxr)
dvr < (L0Lr) * (yryb)/Lr  2 * lam2 * (ylyr)  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)
