Description Usage Arguments Details Value Author(s) References See Also Examples
Solves the initial value problem for stiff or nonstiff systems of either:
a system of ordinary differential equations (ODE) of the form
y' = f(t,y,...)
or
a system of linearly implicit DAES in the form
M y' = f(t,y)
The R function bimd
provides an interface to the Fortran DAE
solver bimd, written by Cecilia Magherini and Luigi Bugnano.
It implements a Blended Implicit Methods of order 4681012 with step size control and continuous output.
The system of DAE's is written as an R function or can be defined in compiled code that has been dynamically loaded.
1 2 3 4 5 6 7 8 9  bimd(y, times, func, parms, nind = c(length(y), 0, 0),
rtol = 1e6, atol = 1e6, jacfunc = NULL, jactype = "fullint",
mass = NULL, massup = NULL, massdown = NULL, verbose = FALSE,
hmax = NULL, hini = 0, ynames = TRUE, minord = NULL,
maxord = NULL, bandup = NULL, banddown = NULL,
maxsteps = 1e4, maxnewtit = c(10, 12, 14, 16, 18), wrkpars = NULL,
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 DAE or 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 DAE or 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 
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. 
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 
verbose 
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 
minord 
the minimum order to be allowed, >= 3 and <= 9.

maxord 
the maximum order to be allowed, >= 
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 
maximal number of steps taken by the solver, for the entire integration. This is different from the settings of this argument in the solvers from package deSolve! 
maxnewtit 
A fivevalued integer vector, with the maximal number of splittingNewton iterations for the solution of the iplicit system in each step for order 4, 6, 8, 10 and 12 respectively. The default is c(10, 12, 14, 16, 18) 
wrkpars 
A 12valued real vector, with extra input parameters,
put in the work vector work, at position work[3:14] in the fortran
code  see details in fortran code. 
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 
... 
additional arguments passed to 
The work is done by the FORTRAN 77 subroutine bimd
, whose
documentation should be consulted for details.
There are four standard choices for the jacobian which can be specified with
jactype
.
The options for jactype are
a full Jacobian, calculated internally by the solver.
a full Jacobian, specified by user
function jacfunc
.
a banded Jacobian, specified by user
function jacfunc
; the size of the bands specified by
bandup
and banddown
.
a banded Jacobian, calculated by bimd;
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 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
).
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 ‘bimd’
returns with an unrecoverable error. If y
has a names
attribute, it will be used to label the columns of the output value.
Karline Soetaert <[email protected].nl>
Francesca Mazzia
L.BRUGNANO, C.MAGHERINI, F.MUGNAI. Blended Implicit Methods for the Numerical Solution of DAE problems. Jour. CAM 189 (2006) 3450.
L.BRUGNANO, C.MAGHERINI The BiM code for the numerical solution of ODEs Jour. CAM 164165 (2004) 145158.
L.BRUGNANO, C.MAGHERINI Some Linear Algebra issues concerning the implementation of Blended Implicit Methods Numer. Linear Alg. Appl. 12 (2005) 305314.
L.BRUGNANO, C.MAGHERINI Economical Error Estimates for Block Implicit Methods for ODEs via Deferred Correction. Appl. Numer. Math. 56 (2006) 608617.
L.BRUGNANO, C.MAGHERINI Blended Implementation of Block Implicit Methods for ODEs Appl. Numer. Math. 42 (2002) 2945.
gamd
another DAE solver from package deTestSet
,
mebdfi
another DAE solver from package deTestSet
,
daspk
another DAE solver from package deSolve
,
ode
for a general interface to most of the ODE solvers
from package deSolve
,
ode.1D
for integrating 1D models,
ode.2D
for integrating 2D models,
ode.3D
for integrating 3D models,
mebdfi
for integrating DAE models,
dopri853
for the DormandPrince RungeKutta method of order 8(53)
diagnostics
to print diagnostic messages.
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 137 138 139 140  ## =======================================================================
## 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 < bimd(yini, times, f1, parms = 0, jactype = "fullint")
plot(out)
## stiff method, usergenerated full Jacobian
out2 < bimd(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 < bimd(yini, times, f1, parms = 0, jactype = "bandint",
bandup = 1, banddown = 1)
## stiff method, usergenerated banded Jacobian
out4 < bimd(yini, times, f1, parms = 0, jactype = "bandusr",
jacfunc = bandjac, bandup = 1, banddown = 1)
## =======================================================================
## Example 2:
## 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 < bimd(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 < bimd(y = yini, mass = Mass, times = times, func = caraxisfun,
parms = NULL, nind = index)
plot(out, which = 1:4, type = "l", lwd = 2)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.