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

Share:

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

1
2
3
4
5
6
7
8
9
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, ...)

Arguments

yini

either a vector with the initial (state) variable values for the ODE system, or NULL.

If yini is a vector, use NA for an initial value which is not specified.

If yini has a names attribute, the names will be available within func and used to label the output matrix.

If yini = NULL, then the boundary conditions must be specified via function bound; if not NULL then yend should also be not NULL.

x

sequence of the independent variable for which output is wanted; the first value of x must be the initial value (at which yini is defined), the final value the end condition (at which yend is defined).

func

either an R-function that computes the values of the derivatives in the ODE system (the model definition) at point x, or a character string giving the name of a compiled function in a dynamically loaded shared library.

If func is an R-function, it must be defined as: func = function(x, y, parms,...). x is the current point of the independent variable in the integration, y is the current estimate of the (state) variables in the ODE system. If the initial values yini has a names attribute, the names will be available inside func. parms is a vector or list of parameters; ... (optional) are any other arguments passed to the function.

The return value of func should be a list, whose first element is a vector containing the derivatives of y with respect to x, and whose next elements are global values that are required at each point in x.

If func is a string, then dllname must give the name of the shared library (without extension) which must be loaded before bvptwp is called. See package vignette "bvpSolve" for more details.

yend

either a vector with the final (state) variable values for the ODE system, or NULL;

if yend is a vector, use NA for a final value which is not specified.

If yend has a names attribute, and yini does not, the names will be available within the functions and used to label the output matrix.

If yend = NULL, then the boundary conditions must be specified via function bound; if not NULL then yini should also be not NULL.

parms

vector or a list with parameters passed to func, jacfunc, bound and jacbound (if present).

If eps is given a value then it should be the **first** element in parms.

epsini

the initial value of the continuation parameter. If NULL and eps is given a value, then epsini takes the default starting value of 0.5. For many singular perturbation type problems, the choice of 0.1 < eps < 1 represents a (fairly) easy problem. The user should attempt to specify an initial problem that is not ‘too’ challenging. epsini must be initialised strictly less than 1 and greater than 0.

eps

the desired value of precision for which the user would like to solve the problem. eps must be less than or equal to epsini. If this is given a value, it must be the first value in parms.

ynames

The names of the variables; used to label the output, and avaliable within the functions.

If ynames is NULL, names can also be passed via yini, yend or yguess.

xguess

Initial grid x, a vector. bvptwp requires the length of xguess to be at least equal to the length of x. If this is not the case, then xguess and yguess will be interpolated to x and a warning printed. If xguess is given, so should yguess be.

Supplying xguess and yguess, based on results from a previous (simpler) BVP-ODE can be used for model continuation, see example 2.

yguess

First guess values of y, corresponding to initial grid xguess; a matrix with number of rows equal to the number of equations, and whose number of columns equals the length of xguess.

if the rows of yguess have a names attribute, the names will be available within the functions and used to label the output matrix.

jacfunc

jacobian (optional) - either an R-function that evaluates the jacobian of func at point x, or a string with the name of a function or subroutine in dllname that computes the Jacobian (see vignette "bvpSolve" for more about this option).

If jacfunc is an R-function, it must be defined as: jacfunc = function(x, y, parms,...). It should return the partial derivatives of func with respect to y, i.e. df(i,j) = dfi/dyj. See last example.

If jacfunc is NULL, then a numerical approximation using differences is used. This is the default.

bound

boundary function (optional) - only if yini and yend are not available. Either an R function that evaluates the i-th boundary element at point x, or a string with the name of a function or subroutine in dllname that computes the boundaries (see vignette "bvpSolve" for more about this option).

If bound is an R-function, it should be defined as: bound = function(i, y, parms, ...). It should return the i-th boundary condition. See last example.

jacbound

jacobian of the boundary function (optional) - only if bound is defined. Either an R function that evaluates the gradient of the i-th boundary element with respect to the state variables, at point x, or a string with the name of a function or subroutine in dllname that computes the boundary jacobian (see vignette "bvpSolve" for more about this option).

If jacbound is an R-function, it should be defined as: jacbound = function(i, y, parms, ...). It should return the gradient of the i-th boundary condition. See last example.

If jacbound is NULL, then a numerical approximation using differences is used. This is the default.

leftbc

only if yini and yend are not available and posbound is not specified: the number of left boundary conditions.

posbound

only used if bound is given: a vector with the position (in the mesh) of the boundary conditions - only the boundary points are allowed. Note that it is simpler to use leftbc.

islin

set to TRUE if the problem is linear - this will speed up the simulation.

nmax

maximal number of subintervals during the calculation.

order

the order of each derivative in func. The default is that all derivatives are 1-st order, in which case order can be set = NULL.

If order is not NULL, the number of equations in func must equal the length of order; the summed values of order must equal the number of variables (ncomp). The jacobian as specified in jacfunc must have number of rows = number of equations and number of columns = number of variables. bound and jacbound remain defined in the number of variables. See example 4 and 4b.

ncomp

used if the model is specified by compiled code, the number of components. See package vignette "bvpSolve".

Also to be used if the boundary conditions are specified by bound, and there is no yguess

atol

error tolerance, a scalar.

cond

if TRUE, uses conditioning in the mesh selection

lobatto

if TRUE, selects a lobatto method.

allpoints

sometimes the solver estimates the solution in a number of extra points, and by default the solutions at these extra points will also be returned.

By setting allpoints equal to FALSE, only output corresponding to the elements in x will be returned.

dllname

a string giving the name of the shared library (without extension) that contains all the compiled function or subroutine definitions referred to in func, jacfunc, bound and jacbound. Note that ALL these subroutines must be defined in the shared library; it is not allowed to merge R-functions with compiled functions.

See package vignette "bvpSolve".

initfunc

if not NULL, the name of the initialisation function (which initialises values of parameters), as provided in ‘dllname’. See package vignette "bvpSolve".

rpar

only when ‘dllname’ is specified: a vector with double precision values passed to the dll-functions whose names are specified by func and jacfunc.

ipar

only when ‘dllname’ is specified: a vector with integer values passed to the dll-functions whose names are specified by func and jacfunc.

nout

only used if dllname is specified and the model is defined in compiled code: the number of output variables calculated in the compiled function func, present in the shared library. Note: it is not automatically checked whether this is indeed the number of output variables calculated in the dll - you have to perform this check in the code. See deSolve's package vignette "compiledCode".

outnames

only used if function is specified in compiled code and nout > 0: the names of output variables calculated in the compiled function. These names will be used to label the output matrix. The length of outnames should be = nout.

forcings

only used if ‘dllname’ is specified: a list with the forcing function data sets, each present as a two-columned matrix, with (time,value); interpolation outside the interval [min(times), max(times)] is done by taking the value at the closest data extreme. This feature is included for consistency with the initial value problem solvers from package deSolve.

See package vignette "compiledCode" from package deSolve.

initforc

if not NULL, the name of the forcing function initialisation function, as provided in ‘dllname’. It MUST be present if forcings has been given a value.

See package vignette "compiledCode" from package deSolve.

fcontrol

A list of control parameters for the forcing functions.

See package vignette "compiledCode" from package deSolve.

verbose

if TRUE then more verbose output will be generated as "warnings".

...

additional arguments passed to the model functions.

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.

See Also

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

  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
## =============================================================================
## 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[2]
  dy2 <- - 3*p*y[1] / (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")

# add analytical solution
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[1] / (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[2]
  dy2 <- -1/x*y[2] - (1-1/(4*x^2))*y[1] + 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")

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

Prob1 <- function(t, y, xi)
   list(c( y[2] , y[1]/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[2], y[3], y[4], R*(y[2]*y[3] - y[1]*y[4]) ))
}

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[4],R*y[3],R*y[2],-R*y[1]))
}

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

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[2]*y[3] - y[1]*y[4]))
}

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

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