bvpcol {bvpSolve} R Documentation

Solves multipoint boundary value problems of ordinary differential equations or differential algebraic equations, using a collocation method.

Description

Solves Boundary Value Problems For Ordinary Differential Equations (ODE) or semi-explicit Differential-Algebraic 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.

Usage

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 = 1e-8, 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, ...)

Details

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 multi-point boundary value problems of ordinary differential equations.

The ODEs and boundary conditions are made available through the user-provided 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 user-provided routines, jacfunc and jacbound. Default is that they are automatically generated by R, using numerical differences.

The user-requested 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 index-2 constraints there. For an index-2 problem in hessenberg form, the projected collocation method of Ascher and Petzold  is used.

Value

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

Note

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 b-splines, 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

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

References

U. Ascher, J. Christiansen and R. D. Russell, (1981) collocation software for boundary-value odes, acm trans. math software 7, 209-222.

G. Bader and U. Ascher, (1987) a new basis implementation for a mixed order boundary value ode solver, siam j. scient. stat. comput. 8, 487-483.

U. Ascher, J. Christiansen and R.D. Russell, (1979) a collocation solver for mixed order systems of boundary value problems, math. comp. 33, 659-679.

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.), 164-185.

J. R. Cash, G. Moore and R. W. Wright, (1995) an automatic continuation strategy for the solution of singularly perturbed linear two-point boundary value problems, j. comp. phys. 122, 266-279.

U. Ascher and R. Spiteri, 1994. collocation software for boundary value differential-algebraic equations, siam j. scient. stat. comput. 15, 938-952.

U. Ascher and L. Petzold, 1991. projected implicit runge-kutta methods for differential- algebraic equations, siam j. num. anal. 28 (1991), 1097-1120.

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.

Examples

## =============================================================================
## 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 = 1e-5
##
## 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
dy2 <- - 3 * p * y / (p+t*t)^2
return(list(c(dy1,
dy2))) }

# parameter value
p    <- 1e-5

# 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)

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 second-order equation
## and with different value of "p".
## =============================================================================

fun2 <- function(t, y, pars)
{ dy <- - 3 * p * y / (p+t*t)^2
list(dy)
}

p <- 1e-4
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 + (1-1/(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 - (1-1/(4x^2)y + sqrt(x)*cos(x)
## =============================================================================

f2 <- function(x, y, parms)
{
dy  <- y
dy2 <- -1/x * y- (1-1/(4*x^2))*y + 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)

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,(((1+ga)/2 -ks*AA)*y*y-y/y-
(AA/A)*(1-(ga-1)*y^2/2))/(ks*A*y)))
}

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 = 1e-4, func = Prob24, eps = 1e-4)
))

## 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 first-order derivatives
f1st<- function(x, y, S) {
list(c(y,
y,
y,
1/S*(y*y - y*y) ))
}

# 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/S, y/S, y/S, -y/S))
}

# boundary
g2 <- function(i, y, S)  {
if (i == 1) return (y)
if (i == 2) return (y)
if (i == 3) return (y - 1)
if (i == 4) return (y)
}

# 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 higher-order 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 fourth-order derivative
f4th <- function(x, y, S) {
list(1/S * (y*y - y*y))
}

# jacobian of derivative function
df4th <- function(x, y, S)  {
matrix(nrow = 1, ncol = 4, byrow = TRUE, data = c(
-1*y/S, y/S, y/S, -y/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 = 1e-6, order = 4, eps = 1e-6,
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 - 1)/2,
(y*y - x)/mu))
}

bound <- function (i, y, p) {
if (i == 1) y -1    # at x=0.5: y2=1
else y              # 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 - p2)*x + dp1
dx2 <- dp2
dx3 <- x
res <- (x - p1)*(x - exp(t))

list(c(dx1, dx2, dx3, res), res = res)
}

boundfun <- function(i,  x, par, ...) {
if (i == 1) return(x - sin(0))
if (i == 2) return(x - 1)
if (i == 3) return(x - sin(1))
if (i == 4) return((x - sin(1))*(x - 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 = 1e-4, ncomp = 4, leftbc = 2,
dae = list(index = 2,  nalg = 1))

# solved using yini, yend
out1 <- bvpcol (func = bvpdae, x = x, parms = 1e-4,
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)

[Package bvpSolve version 1.4.3 Index]