bvpcol {bvpSolve}  R Documentation 
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.
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.
## ============================================================================= ## 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)