bvptwp {bvpSolve} R Documentation

Solves two-point boundary value problems of ordinary differential equations, using a mono-implicit Runge-Kutta formula

Description

Solves a boundary value problem of a system of ordinary differential equations. This is an implementation of the fortran code twpbvpc, based on mono-implicit Runge-Kutta formulae of orders 4, 6 and 8 in a deferred correction framework and that uses conditioning in the mesh selection.

written by J.R. Cash, F. Mazzia and M.H. Wright.

Rather than MIRK, it is also possible to select a lobatto method. This then uses the code 'twpbvplc', written by Cash and Mazzia.

It is possible to solve stiff systems, by using an automatic continuation strategy. This then uses the code 'acdc'.

Usage

bvptwp(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, cond = FALSE, lobatto = FALSE,
allpoints = 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, ...)

Details

This is an implementation of the method twpbvpC, written by Cash, Mazzia and Wright, to solve two-point boundary value problems of ordinary differential equations.

A boundary value problem does not have all initial values of the state variable specified. Rather some conditions are specified at the end of the integration interval. The number of unknown boundary conditions must be equal to the number of equations (or the number of dependent variables y).

The ODEs and boundary conditions are made available through the user-provided routines, func and vectors yini and yend or (optionally) bound.

The corresponding partial derivatives for func and bound are optionally available through the user-provided routines, jacfunc and jacbound. Default is that they are automatically generated by bvptwp, using numerical differences.

The user-requested tolerance is provided through tol. The default is 1e^-6

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.' It may also help to start with a simple version of the model, and use its result as initial guess to solve the more complex problem (continuation strategy, see example 2, and package vignette "bvpSolve").

Models may be defined in compiled C or Fortran code, as well as in an R-function.

This is similar to the initial value problem solvers from package deSolve. See package vignette "bvpSolve" for details about writing compiled code.

The fcontrol argument is a list that can supply any of the following components (conform the definitions in the approx function):

method

specifies the interpolation method to be used. Choices are "linear" or "constant",

rule

an integer describing how interpolation is to take place outside the interval [min(times), max(times)]. If rule is 1 then an error will be triggered and the calculation will stop if times extends the interval of the forcing function data set. If it is 2, the *default*, the value at the closest data extreme is used, a warning will be printed if verbose is TRUE,

Note that the default differs from the approx default

f

For method="constant" a number between 0 and 1 inclusive, indicating a compromise between left- and right-continuous step functions. If y0 and y1 are the values to the left and right of the point then the value is y0*(1-f)+y1*f so that f=0 is right-continuous and f=1 is left-continuous,

ties

Handling of tied times values. Either a function with a single vector argument returning a single number result or the string "ordered".

Note that the default is "ordered", hence the existence of ties will NOT be investigated; in the C code this will mean that -if ties exist, the first value will be used; if the dataset is not ordered, then nonsense will be produced.

Alternative values for ties are mean, min etc

The defaults are:

fcontrol=list(method="linear", rule = 2, f = 0, ties = "ordered")

Note that only ONE specification is allowed, even if there is more than one forcing function data set.

This -advanced- feature is explained in deSolve's package vignette "compiledCode".

Value

A matrix of class bvpSolve, with up to as many rows as elements in x and as many columns as elements in yini or ncomp plus the number of "global" values returned from func, plus an additional column (the first) for the x-value.

Typically, there will be one row for each element in x unless the solver returns with an unrecoverable error. In certain cases, bvptwp will return the solution also in a number of extra points. This will occur if the number of points as in x was not sufficient. In order to not return these extra points, set allpoints equal to FALSE.

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.

Note

When order is not NULL, then it should contain the order of all equations in func. If the order of some equations > 1, then there will be less equations than variables. The number of equations should be equal to the length of order, while the number of variables will be the sum of order.

For instance, if order = c(1,2,3,4), then the first equation will be of order 1, the second one of order 2, ...and the last of order 4.

There will be 1+2+3+4 = 10 variables. For instance, if the derivative function defines (A', B”, C”', D””) respectively, then the variable vector will contain values for A, B, B', C, C', C”, D, D', D”, D”'; in that order. This is also the order in which the initial and end conditions of all variables must be specified.

If neq are the number of equations, and ncomp the number of variables, then the Jacobian of the derivative function as specified in jacfunc must be of dimension (neq, ncomp).

The jacobian of the boundaries, as specified in jacbound should return a vector of length = ncomp

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

Jeff Cash <j.cash@imperial.ac.uk>

Francesca Mazzia <mazzia@dm.uniba.it>

References

J.R. Cash and M.H. Wright, A deferred correction method for nonlinear two-point boundary value problems: implementation and numerical evaluation, SIAM J. Sci. Stat. Comput., 12 (1991) 971 989.

Cash, J. R. and F. Mazzia, A new mesh selection algorithm, based on conditioning, for two-point boundary value codes. J. Comput. Appl. Math. 184 (2005), no. 2, 362–381.

Cash, J. R. and F. Mazzia, in press. Hybrid Mesh Selection Algorithms Based on Conditioning for Two-Point Boundary Value Problems, Journal of Numerical Analysis, Industrial and Applied Mathematics.

bvpshoot for the shooting method

bvpcol for the collocation method

diagnostics.bvpSolve, for a description of diagnostic messages

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(y = -0.1 / sqrt(p+0.01), dy=NA)
end  <- c(     0.1 / sqrt(p+0.01), NA)

# Solve bvp
sol  <- as.data.frame(bvptwp(yini = init, x = seq(-0.1, 0.1, by = 0.001),
func = fun, yend = end))
plot(sol\$x, sol\$y, type = "l")

curve(x/sqrt(p+x*x), add = TRUE, type = "p")

## =============================================================================
## Example 1b:
## Same problem, now solved as a second-order equation.
## =============================================================================

fun2 <- function(t, y, pars)  {
dy <- - 3 * p * y / (p+t*t)^2
list(dy)
}
sol2  <- bvptwp(yini = init, yend = end, order = 2,
x = seq(-0.1, 0.1, by = 0.001), func = fun2)

## =============================================================================
## Example 2: 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  <- bvptwp(yini = c(y = 1, dy = NA),
yend = c(-0.5, NA), x = x, func = f2)
plot(sol, which = "y")

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 3  - solved with automatic continuation
## d2y/dx2 = y/xi
## =============================================================================

Prob1 <- function(t, y, xi)
list(c( y , y/xi ))

x <-  seq(0, 1, by = 0.01)
xi <- 0.1
twp <- bvptwp(yini = c(1, NA), yend = c(0, NA), x = x,
islin = TRUE, func = Prob1, parms = xi, eps = xi)

xi <-0.00001
twp2 <- bvptwp(yini = c(1, NA), yend = c(0, NA), x = x,
islin = TRUE, func = Prob1, parms = xi, eps = xi)

plot(twp, twp2, which = 1, main = "test problem 1")

# exact solution
curve(exp(-x/sqrt(xi))-exp((x-2)/sqrt(xi))/(1-exp(-2/sqrt(xi))),
0, 1, add = TRUE, type = "p")

curve(exp(-x/sqrt(0.1))-exp((x-2)/sqrt(0.1))/(1-exp(-2/sqrt(0.1))),
0, 1, add = TRUE, type = "p")

## =============================================================================
## Example 4  - 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)
## =============================================================================

f2<- function(x, y, parms, R) {
list(c(y, y, y, R*(y*y - y*y) ))
}

df2 <- function(x, y, parms, R) {
matrix(nrow = 4, ncol = 4, byrow = TRUE, data = c(
0,        1,     0,     0,
0,        0,     1,     0,
0,        0,     0,     1,
-1*R*y,R*y,R*y,-R*y))
}

g2 <- function(i, y, parms, R) {
if (i == 1) return(y)
if (i == 2) return(y)
if (i == 3) return(y-1)
if (i == 4) return(y)
}

dg2 <- function(i, y, parms, R) {
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))
}

init <- c(1, NA)
R    <- 100
sol  <- bvptwp(x = seq(0, 1, by = 0.01), leftbc = 2,
func = f2, R = R, ncomp = 4,
bound = g2, jacfunc = df2, jacbound = dg2)
plot(sol[,1:2])  # columns do not have names

mf <- par ("mfrow")
sol  <- bvptwp(x = seq(0, 1, by = 0.01), leftbc = 2,
func = f2, ynames = c("y", "dy", "d2y", "d3y"), R=R,
bound = g2, jacfunc = df2, jacbound = dg2)
plot(sol)        # here they do
par(mfrow = mf)

## =============================================================================
## Example 4b - 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, parms, R) {
list(R * (y*y - y*y))
}

# jacobian of derivative function
df4th <- function(x, y, parms, R)  {
df <- matrix(nrow = 1, ncol = 4, byrow = TRUE, data = c(
-1*R*y, R*y, R*y, -R*y))
}

# boundary function - same as previous example

# jacobian of boundary - same as previous

# order = 4 specifies the equation to be 4th order
sol2 <- bvptwp(x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f4th, R = R, order = 4,
bound = g2, jacfunc = df4th, jacbound = dg2)

max(abs(sol2-sol))

[Package bvpSolve version 1.4.3 Index]