Solves a boundary value problem of a system of ordinary differential equations. This is an implementation of the fortran code twpbvpc, based on monoimplicit RungeKutta 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'.
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 = 1e8, 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, ...)

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 
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 
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 If 
ncomp 
used if the model is specified by compiled code, the number of
components. See package vignette Also to be used if the boundary conditions are specified by 
atol 
error tolerance, a scalar. 
cond 
if 
lobatto 
if 
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 
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 
if 
... 
additional arguments passed to the model functions. 
This is an implementation of the method twpbvpC, written by Cash, Mazzia and Wright, to solve twopoint 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
userprovided routines, func
and vectors yini
and yend
or (optionally) bound
.
The corresponding partial derivatives for func
and bound
are optionally available through the
userprovided routines, jacfunc
and jacbound
. Default is that
they are automatically generated by bvptwp
, using numerical differences.
The userrequested 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 Rfunction.
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):
specifies the interpolation method to be used. Choices are "linear" or "constant",
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
For method="constant"
a number between 0
and
1
inclusive, indicating a compromise between left and
rightcontinuous step functions. If y0
and y1
are the
values to the left and right of the point then the value is
y0*(1f)+y1*f
so that f=0
is rightcontinuous and
f=1
is leftcontinuous,
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".
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.
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
Karline Soetaert <karline.soetaert@nioz.nl>
Jeff Cash <j.cash@imperial.ac.uk>
Francesca Mazzia <mazzia@dm.uniba.it>
J.R. Cash and M.H. Wright, A deferred correction method for nonlinear twopoint 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 twopoint 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 TwoPoint 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.
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 = 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(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 secondorder 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 + (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 < 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((x2)/sqrt(xi))/(1exp(2/sqrt(xi))),
0, 1, add = TRUE, type = "p")
curve(exp(x/sqrt(0.1))exp((x2)/sqrt(0.1))/(1exp(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 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, 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(sol2sol))

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.