Solves Boundary Value Problems For Ordinary Differential Equations (ODE) or semiexplicit DifferentialAlgebraic Equations (DAE) with index at most 2.
It is possible to solve stiff ODE systems, by using an automatic continuation strategy
This is an implementation of the fortran codes colsys.f, colnew.f and coldae.f written by respectively U. Ascher, J. christiansen and R.D. Russell (colsys), U. Ascher and G. Bader (colnew) and U. Ascher and C. Spiteri.
The continuation strategy is an implementation of the fortran code colmod written by J.R. Cash, M.H. Wright and F. Mazzia.
1 2 3 4 5 6 7 8 9  bvpcol (yini = NULL, x, func, yend = NULL, parms = NULL,
order = NULL, ynames = NULL, xguess = NULL, yguess = NULL,
jacfunc = NULL, bound = NULL, jacbound = NULL,
leftbc = NULL, posbound = NULL, islin = FALSE, nmax = 1000,
ncomp = NULL, atol = 1e8, colp = NULL, bspline = FALSE,
fullOut = TRUE, dllname = NULL, initfunc = dllname,
rpar = NULL, ipar = NULL, nout = 0, outnames = NULL,
forcings = NULL, initforc = NULL, fcontrol = NULL,
verbose = FALSE, epsini = NULL, eps = epsini, dae = NULL, ...)

yini 
either a vector with the initial (state) variable values for
the ODE system, or If If If 
x 
sequence of the independent variable 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 point If The return value of If the problem is a DAE, then the algebraic equations should be the last. If 
yend 
either a vector with the final (state) variable values for the
ODE system, or if If If 
parms 
vector or a list with parameters passed to If 
epsini 
the initial value of the continuation parameter. If

eps 
the desired value of precision for which the user would like
to solve the problem. 
ynames 
The names of the variables; used to label the output, and avaliable within the functions. If 
xguess 
Initial grid Supplying 
yguess 
First guess values of if the rows of It is also allowed to pass the output of a previous run for continuation.
This will use the information that is stored in the attributes
See example 3b. 
jacfunc 
jacobian (optional)  either an Rfunction that evaluates the
jacobian of If If 
bound 
boundary function (optional)  only if If 
jacbound 
jacobian of the boundary function (optional)  only if
If If 
leftbc 
only if 
posbound 
only used if 
islin 
set to 
nmax 
maximal number of subintervals during the calculation. 
order 
the order of each derivative in For higherorder derivatives, specifying the order can improve computational efficiency, but this interface is more complex. If 
ncomp 
used if the model is specified by compiled code, the number of
components (or equations). See package vignette Also to be used if the boundary conditions are specified by 
atol 
error tolerance, a scalar. 
colp 
number of collocation points per subinterval. 
bspline 
if 
fullOut 
if set to 
dllname 
a string giving the name of the shared library
(without extension) that contains all the compiled function or
subroutine definitions referred to in See package vignette 
initfunc 
if not 
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 function is specified in compiled code 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 package vignette 
initforc 
if not See package vignette 
fcontrol 
A list of control parameters for the forcing functions. See package vignette 
verbose 
when 
dae 
if the problem is a DAE, should be a list containing the
See example 5 
... 
additional arguments passed to the model functions. 
If eps
does not have a value and dae
= NULL,
then the method is based on an implementation
of the Collocation methods called "colnew" and
"colsys" to solve multipoint boundary value problems of ordinary
differential equations.
The ODEs and boundary conditions are made available through the
userprovided routines, func
and vectors yini
and yend
or (optionally) bound
. bvpcol
can also solve multipoint
boundary value problems (see one but last example).
The corresponding partial derivatives are optionally available through the
userprovided routines, jacfunc
and jacbound
. Default is that
they are automatically generated by R, using numerical differences.
The userrequested tolerance is provided through atol
.
If the function terminates because the maximum number of subintervals was exceeded, then it is recommended that 'the program be run again with a larger value for this maximum.'
If eps
does have a value, then the method is based on an implementation
of the Collocation methods called "colmod".
The type of problems which this is designed to solve typically
involve a small positive parameter 0 < eps << 1.
As eps becomes progressively smaller, the problem normally becomes
increasingly difficult to approximate numerically (for example, due
to the appearance of narrow transition layers in the profile of the
analytic solution).
The idea of continuation is to solve a chain of problems in which the parameter eps decreases monotonically towards some desired value. That is, a sequence of problems is attempted to be solved:
epsini > eps1 > eps2 > eps3 > ..... > eps > 0
where epsini
is a user provided starting value and eps
is a
user desired final value for the parameter.
If dae
is not NULL, then it is assumed that a DAE has to be solved.
In that case, dae
should contain give the index
of the DAE and the number of algebraic
equations (nalg
).
(this part comes from the comments in the code coldae). With respect to the dae, it should be noted that the code does not explicitly check the index of the problem, so if the index is > 2 then the code will not work well. The number of boundary conditions required is independent of the index. it is the user's responsibility to ensure that these conditions are consistent with the constraints. The conditions at the left end point must include a subset equivalent to specifying the index2 constraints there. For an index2 problem in hessenberg form, the projected collocation method of Ascher and Petzold [2] is used.
A matrix of class bvpSolve
, with up to as many rows as elements in
x
and as many columns
as elements in yini
plus the number of "global" values returned
in the second element of the return from func
, plus an additional
column (the first) for the x
value.
There will be one row for each element in x
unless the solver returns
with an unrecoverable error.
If ynames
is given, or yini
, yend
has a names attribute,
or yguess
has named rows, the names will be used to label the
columns of the output value.
The output will also have attributes
istate
and rstate
which contain the collocation output required e.g. for continuation of a
problem, unless fullOutput
is FALSE
colnew.f (Bader and Ascher, 1987), is a modification of the code
colsys.f (Ascher, Christiansen and Russell, 1981), which incorporates
a new basis representation replacing bsplines, and improvements for
the linear and nonlinear algebraic equation solvers. To toggle on/off
colsys, set bspline
= TRUE
/FALSE
colmod is a revised version of the package colnew by Bader and Ascher (1987), which in turn is a modification of the package colsys by Ascher, Christiansen and Russell (1981). Colmod has been adapted to allow an automatic continuation strategy to be used (Cash et al., 1995).
The mesh selection algorithm used in colmod differs from that used in colnew
Karline Soetaert <karline.soetaert@nioz.nl>
U. Ascher, J. Christiansen and R. D. Russell, (1981) collocation software for boundaryvalue odes, acm trans. math software 7, 209222.
G. Bader and U. Ascher, (1987) a new basis implementation for a mixed order boundary value ode solver, siam j. scient. stat. comput. 8, 487483.
U. Ascher, J. Christiansen and R.D. Russell, (1979) a collocation solver for mixed order systems of boundary value problems, math. comp. 33, 659679.
U. Ascher, J. Christiansen and R.D. Russell, (1979) colsys  a collocation code for boundary value problems, lecture notes comp.sc. 76, springer verlag, B. Childs et. al. (eds.), 164185.
J. R. Cash, G. Moore and R. W. Wright, (1995) an automatic continuation strategy for the solution of singularly perturbed linear twopoint boundary value problems, j. comp. phys. 122, 266279.
U. Ascher and R. Spiteri, 1994. collocation software for boundary value differentialalgebraic equations, siam j. scient. stat. comput. 15, 938952.
U. Ascher and L. Petzold, 1991. projected implicit rungekutta methods for differential algebraic equations, siam j. num. anal. 28 (1991), 10971120.
bvpshoot
for the shooting method
bvptwp for a MIRK formula
diagnostics.bvpSolve
, for a description of diagnostic messages
approx.bvpSolve
, for approximating solution in new values
plot.bvpSolve
, for a description of plotting the output of the
BVP solvers.
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 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316  ## =============================================================================
## Example 1: simple standard problem
## solve the BVP ODE:
## d2y/dt^2=3py/(p+t^2)^2
## y(t= 0.1)=0.1/sqrt(p+0.01)
## y(t= 0.1)= 0.1/sqrt(p+0.01)
## where p = 1e5
##
## analytical solution y(t) = t/sqrt(p + t^2).
##
## The problem is rewritten as a system of 2 ODEs:
## dy=y2
## dy2=3p*y/(p+t^2)^2
## =============================================================================
#
# Derivative function
#
fun < function(t, y, pars) {
dy1 < y[2]
dy2 <  3 * p * y[1] / (p+t*t)^2
return(list(c(dy1,
dy2))) }
# parameter value
p < 1e5
# initial and final condition; second conditions unknown
init < c(0.1 / sqrt(p+0.01), NA)
end < c( 0.1 / sqrt(p+0.01), NA)
# Solve bvp
sol < bvpcol(yini = init, yend = end,
x = seq(0.1, 0.1, by = 0.001), func = fun)
plot(sol, which = 1)
# add analytical solution
curve(x/sqrt(p+x*x), add = TRUE, type = "p")
diagnostics(sol)
zoom < approx(sol, xout = seq(0.005, 0.005, by = 0.0001))
plot(zoom, which = 1, main = "zoom in on [0.0005,0.0005]")
## =============================================================================
## Example 1b:
## Same problem, now solved as a secondorder equation
## and with different value of "p".
## =============================================================================
fun2 < function(t, y, pars)
{ dy <  3 * p * y[1] / (p+t*t)^2
list(dy)
}
p < 1e4
sol2 < bvpcol(yini = init, yend = end, order = 2,
x = seq(0.1, 0.1, by = 0.001), func = fun2)
# plot both runs at once:
plot(sol, sol2, which = 1)
## =============================================================================
## Example 1c: simple
## solve d2y/dx2 + 1/x*dy/dx + (11/(4x^2)y = sqrt(x)*cos(x),
## on the interval [1,6] and with boundary conditions:
## y(1)=1, y(6)=0.5
##
## Write as set of 2 odes
## dy/dx = y2
## dy2/dx =  1/x*dy/dx  (11/(4x^2)y + sqrt(x)*cos(x)
## =============================================================================
f2 < function(x, y, parms)
{
dy < y[2]
dy2 < 1/x * y[2] (11/(4*x^2))*y[1] + sqrt(x)*cos(x)
list(c(dy, dy2))
}
x < seq(1, 6, 0.1)
sol < bvpcol(yini = c(1, NA), yend = c(0.5, NA), bspline = TRUE,
x = x, func = f2)
plot(sol, which = 1)
# add the analytic solution
curve(0.0588713*cos(x)/sqrt(x) + 1/4*sqrt(x)*cos(x)+0.740071*sin(x)/sqrt(x)+
1/4*x^(3/2)*sin(x), add = TRUE, type = "l")
## =============================================================================
## Example 2. Uses continuation
## Test problem 24
## =============================================================================
Prob24< function(t, y, ks) { #eps is called ks here
A < 1+t*t
AA < 2*t
ga < 1.4
list(c(y[2],(((1+ga)/2 ks*AA)*y[1]*y[2]y[2]/y[1]
(AA/A)*(1(ga1)*y[1]^2/2))/(ks*A*y[1])))
}
ini < c(0.9129, NA)
end < c(0.375, NA)
xguess < c(0, 1)
yguess < matrix(nrow = 2, ncol = 2, 0.9 )
# bvpcol works with eps NOT too small, and good initial condition ...
sol < bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
xguess = xguess, yguess = yguess,
parms = 0.1, func = Prob24, verbose = FALSE)
# when continuation is used: does not need a good initial condition
sol2 < bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 0.05, func = Prob24,
eps = 0.05)
#zoom < approx(sol2, xout = seq(0.01, 0.02, by = 0.0001))
#plot(zoom, which = 1, main = "zoom in on [0.01, 0.02]")
sol3 < bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 0.01, func = Prob24 , eps = 0.01)
sol4 < bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 0.001, func = Prob24, eps = 0.001)
# This takes a long time
## Not run:
print(system.time(
sol5 < bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 1e4, func = Prob24, eps = 1e4)
))
## End(Not run)
plot(sol, sol2, sol3, sol4, which = 1, main = "test problem 24",
lwd = 2)
legend("topright", col = 1:4, lty = 1:4, lwd = 2,
legend = c("0.1", "0.05", "0.01", "0.001"), title = "eps")
## =============================================================================
## Example 3  solved with specification of boundary, and jacobians
## d4y/dx4 =R(dy/dx*d2y/dx2 y*dy3/dx3)
## y(0)=y'(0)=0
## y(1)=1, y'(1)=0
##
## dy/dx = y2
## dy2/dx = y3 (=d2y/dx2)
## dy3/dx = y4 (=d3y/dx3)
## dy4/dx = R*(y2*y3 y*y4)
## =============================================================================
# derivative function: 4 firstorder derivatives
f1st< function(x, y, S) {
list(c(y[2],
y[3],
y[4],
1/S*(y[2]*y[3]  y[1]*y[4]) ))
}
# jacobian of derivative function
df1st < function(x, y, S) {
matrix(nrow = 4, ncol = 4, byrow = TRUE, data = c(
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1,
1*y[4]/S, y[3]/S, y[2]/S, y[1]/S))
}
# boundary
g2 < function(i, y, S) {
if (i == 1) return (y[1])
if (i == 2) return (y[2])
if (i == 3) return (y[1]  1)
if (i == 4) return (y[2])
}
# jacobian of boundary
dg2 < function(i, y, S) {
if (i == 1) return(c(1, 0, 0, 0))
if (i == 2) return(c(0, 1, 0, 0))
if (i == 3) return(c(1, 0, 0, 0))
if (i == 4) return(c(0, 1, 0, 0))
}
# we use posbound to specify the position of boundary conditions
# we can also use leftbc = 2 rather than posbound = c(0,0,1,1)
S < 1/100
sol < bvpcol(x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f1st, parms = S, eps = S,
bound = g2, jacfunc = df1st, jacbound = dg2)
plot(sol)
## =============================================================================
## Example 3b  solved with specification of boundary, and jacobians
## and as a higherorder derivative
## d4y/dx4 =R(dy/dx*d2y/dx2 y*dy3/dx3)
## y(0)=y'(0)=0
## y(1)=1, y'(1)=0
## =============================================================================
# derivative function: one fourthorder derivative
f4th < function(x, y, S) {
list(1/S * (y[2]*y[3]  y[1]*y[4]))
}
# jacobian of derivative function
df4th < function(x, y, S) {
matrix(nrow = 1, ncol = 4, byrow = TRUE, data = c(
1*y[4]/S, y[3]/S, y[2]/S, y[1]/S))
}
# boundary function  same as previous example
# jacobian of boundary  same as previous
# order = 4 specifies the equation to be 4th order
# solve with bspline false
S < 1/100
sol < bvpcol (x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f4th, order = 4,
parms = S, eps = S, bound = g2, jacfunc = df4th,
jacbound = dg2 )
plot(sol)
# Use (manual) continuation to find solution of a more difficult example
# Previous solution collocation from sol passed ("guess = sol")
sol2 < bvpcol(x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f4th,
parms = 1e6, order = 4, eps = 1e6,
bound = g2, jacfunc = df4th, jacbound = dg2 )
# plot both at same time
plot(sol, sol2, lwd = 2)
legend("bottomright", leg = c(100, 10000), title = "R = ",
col = 1:2, lty = 1:2, lwd = 2)
## =============================================================================
## Example 4  a multipoint bvp
## dy1 = (y2  1)/2
## dy2 = (y1*y2  x)/mu
## over interval [0,1]
## y1(1) = 0; y2(0.5) = 1
## =============================================================================
multip < function (x, y, p) {
list(c((y[2]  1)/2,
(y[1]*y[2]  x)/mu))
}
bound < function (i, y, p) {
if (i == 1) y[2] 1 # at x=0.5: y2=1
else y[1] # at x= 1: y1=0
}
mu < 0.1
sol < bvpcol(func = multip, bound = bound,
x = seq(0, 1, 0.01), posbound = c(0.5, 1))
plot(sol)
# check boundary value
sol[sol[,1] == 0.5,]
## =============================================================================
## Example 5  a bvp DAE
## =============================================================================
bvpdae < function(t, x, ks, ...) {
p1 < p2 < sin(t)
dp1 < dp2 < cos(t)
dx1 < (ks + x[2]  p2)*x[4] + dp1
dx2 < dp2
dx3 < x[4]
res < (x[1]  p1)*(x[4]  exp(t))
list(c(dx1, dx2, dx3, res), res = res)
}
boundfun < function(i, x, par, ...) {
if (i == 1) return(x[1]  sin(0))
if (i == 2) return(x[3]  1)
if (i == 3) return(x[2]  sin(1))
if (i == 4) return((x[1]  sin(1))*(x[4]  exp(1))) # Not used here..
}
x < seq(0, 1, by = 0.01)
mass < diag(nrow = 4) ; mass[4, 4] < 0
# solved using boundfun
out < bvpcol (func = bvpdae, bound = boundfun, x = x,
parms = 1e4, ncomp = 4, leftbc = 2,
dae = list(index = 2, nalg = 1))
# solved using yini, yend
out1 < bvpcol (func = bvpdae, x = x, parms = 1e4,
yini = c(sin(0), NA, 1, NA),
yend = c(NA, sin(1), NA, NA),
dae = list(index = 2, nalg = 1))
# the analytic solution
ana < cbind(x, "1" = sin(x), "2" = sin(x), "3" = 1, "4" = 0, res = 0)
plot(out, out1, obs = ana)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.