Solves a boundary value problem of a system of ordinary differential
equations using the single shooting method.
This combines the integration routines from package
deSolve
with rootfinding methods from package rootSolve
.
Preferentially bvptwp
or bvpcol
should be used rather than
bvpshoot
, as they give more precise output.
1 2 3 4 5 6 
yini 
either a vector with the initial (state) variable
values for the ODE system, or a function that calculates the
initial condition, or If if If if 
x 
sequence of the independent variable for which output is wanted;
the first value of 
func 
an Rfunction that computes the values of the derivatives in
the ODE system (the model definition) at x. The return value of Note that it is not possible to use 
yend 
either a vector with the final (state) variable values for the
ODE system, a function that calculates the final condition
or if If If if 
parms 
vector or a list with parameters passed to 
order 
the order of each derivative in If 
guess 
guess for the value(s) of the unknown initial conditions; if initial and final conditions are specified by If initial and final conditions are specified by the boundary function

jacfunc 
jacobian (optional)  an Rfunction that evaluates the
jacobian of
If

bound 
boundary function (optional)  only if
if not 
jacbound 
jacobian of the boundary function (optional)  only if

leftbc 
only if 
posbound 
only used if 
ncomp 
only used if the boundaries are specified via the boundary
function 
atol 
absolute error tolerance, either a scalar or a vector, one value for each unknown element  passed to function multiroot  see help of this function. 
rtol 
relative error tolerance, either a scalar or a vector, one value for each unknown element  passed to function multiroot  see help of this function. 
extra 
if too many boundary conditions are given, then it is assumed that an extra parameter has to be estimated.

maxiter 
the maximal number of iterations allowed in the root solver. 
positive 
set to 
method 
the integration method used, one of ("lsoda", "lsode", "lsodes", "vode", "euler", "rk4", "ode23" or "ode45"). 
... 
additional arguments passed to the integrator and (possibly) the model functions. 
This is a simple implementation of the shooting method to solve 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 shooting method, is a rootsolving method. There are two strategies:
If initial and end conditions are specified with yini
and yend
then the (unspecified) initial conditions are the unknown values to be solved
for; the function value whose root has to be found are the deviations from
the specified conditions at the end of the integration interval.
Thus, starting with an initial guess of the initial conditions (as
provided in guess
), the ODE model is solved as an initial value
problem, and after termination, the discrepancy of the modeled final
conditions with the known final condition is assessed (the cost function).
The root of this cost function is to be found.
If initial and end conditions are specified with
bound
, then the
unknowns are all initial conditions; the function whose root is to be
found is given by bound
.
Starting from a guess of the initial values, one of the integrators from
package deSolve
(as specified with method
) is used to solve the
resulting initial value problem.
Function multiroot
from package rootSolve
is used to
retrieve the root.
For this method to work, the model should be even determined, i.e. the number of equations should equal the number of unknowns.
bvpshoot
distinguishes two cases:
1. the total number of specified boundary conditions (on both the start and
end of the integration interval) equals the number of boundary value
problem equations (or the number of dependent variables y
).
2. The number of boundary conditions specified exceeds
the number of equations. In this case, extra
parameters have to be solved
for to make the model even determined.
See example nr 4.
A matrix with up to as many rows as elements in x
and as many columns
as the number of state variables in the ODE system plus the number of
"global" values returned in the next elements of the return from func
,
plus an additional column (the first) for the xvalue.
There will be one row for each element in x
unless the solver
returns with an unrecoverable error.
If yini
has a names attribute, it will be used to label the columns
of the output value. If yini
is not named, the solver will try to find
the names in yend
. If the boundaries are specified by bound
then the names from guess
will be used.
The output will have the attribute roots
, which returns the value(s)
of the root(s) solved for (root
), the function value (f.root
),
and the number of iterations (iter
) required to find the root.
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.
Do not specify the jacobian if the maximal order>1.
Karline Soetaert <karline.soetaert@nioz.nl>
bvptwp
for the MIRK method
lsoda
, lsode
,
lsodes
, vode
,
rk
, rkMethod
for details about the integration method
multiroot
, the rootsolving method used
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 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222  ## =============================================================================
## 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 < bvpshoot(yini = init, x = seq(0.1, 0.1, by = 0.001),
func = fun, yend = end, guess = 1)
plot(sol, which = "y", type = "l")
# add analytical solution
curve(x/sqrt(p+x*x), add = TRUE, type = "p")
## =============================================================================
## Example 1b: 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 < bvpshoot(yini = c(y = 1, dy = NA), yend = c(0.5, NA),
x = x, func = f2, guess = 1)
# plot both state variables
plot(sol, type = "l", lwd = 2)
# plot only y and add the analytic solution
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 2  initial condition is a function of the unknown x
## tubular reactor with axial dispersion
## y''=Pe(y'+Ry^n) Pe=1,R=2,n=2
## on the interval [0,1] and with initial conditions:
## y'0=Pe(y(0)1), y'(1)=0
##
## dy=y2
## dy2=Pe(dyRy^n)
## =============================================================================
Reactor<function(x, y, parms)
{
list(c(y[2],
Pe * (y[2]+R*(y[1]^n))))
}
Pe < 1
R < 2
n < 2
x < seq(0, 1, by = 0.01)
# 1. yini is a function here
yini < function (x, parms) c(x, Pe*(x1))
system.time(
sol < bvpshoot(func = Reactor, yend = c(y = NA, dy = 0),
yini = yini, x = x, extra = 1)
)
plot(sol, which = "y", main = "Reactor", type = "l", lwd = 2)
attributes(sol)$roots
# 2. using boundary function rather than yini...
bound < function(i, y, p) {
if (i == 1) return(y[2]  Pe*(y[1]1))
if (i == 2) return(y[2])
}
# need to give number of left boundary conditions + guess of all initial
# conditions (+ names)
system.time(
Sol2< bvpshoot(func = Reactor, x = x, leftbc = 1,
bound = bound, guess = c(y = 1, dy = 1) )
)
attributes(Sol2)$roots
# 3. boundary function with jacobian of boundary function
jacbound < function(i, y, p) {
if (i == 1) return(c(Pe*y[1], 1))
if (i == 2) return(c(0, 1))
}
system.time(
Sol3<bvpshoot(func = Reactor, x = x, leftbc = 1, bound = bound,
jacbound = jacbound, guess = c(y = 1, dy = 1) )
)
attributes(Sol3)$roots
## =============================================================================
## Example 2b  same as 2 but written as higherorder equation
## y''=Pe(y'+Ry^n) Pe=1,R=2,n=2
## on the interval [0,1] and with initial conditions:
## y'0=Pe(y(0)1), y'(1)=0
## =============================================================================
Reactor2<function(x, y, parms)
list (Pe * (y[2]+R*(y[1]^n)))
Pe < 1
R < 2
n < 2
x < seq(0, 1, by = 0.01)
# 1. yini is a function here
yini < function (x, parms) c(x, Pe*(x1))
# need to specify that order = 2
system.time(
sol2 < bvpshoot(func = Reactor2, yend = c(y = NA, dy = 0), order=2,
yini = yini, x = x, extra = 1)
)
max(abs(sol2sol))
## =============================================================================
## Example 3  final condition is a residual function
## The example MUSN from Ascher et al., 1995.
## Numerical Solution of Boundary Value Problems for Ordinary Differential
## Equations, SIAM, Philadelphia, PA, 1995.
##
## The problem is
## u' = 0.5*u*(w  u)/v
## v' = 0.5*(w  u)
## w' = (0.9  1000*(w  y)  0.5*w*(w  u))/z
## z' = 0.5*(w  u)
## y' = 100*(y  w)
## on the interval [0 1] and with boundary conditions:
## u(0) = v(0) = w(0) = 1, z(0) = 10, w(1) = y(1)
## =============================================================================
musn < function(t, Y, pars)
{
with (as.list(Y),
{
du < 0.5 * u * (wu)/v
dv < 0.5 * (wu)
dw < (0.9  1000 * (wy)  0.5 * w * (wu))/z
dz < 0.5 * (wu)
dy < 100 * (yw)
return(list(c(du, dv, dw, dy, dz)))
})
}
#
# Residuals of end conditions
#
res < function (Y, yini, pars) with (as.list(Y), wy)
#
# Initial values; y= NA= not available
#
init < c(u = 1, v = 1, w = 1, y = NA, z = 10)
sol <bvpshoot(y = init, x = seq(0, 1, by = 0.05), func = musn,
yend = res, guess = 1, atol = 1e10, rtol = 0)
pairs(sol, main = "MUSN")
## =============================================================================
## Example 4  solve also for unknown parameter
## Find the 4th eigenvalue of Mathieu's equation:
## y''+(lam10cos2t)y=0 on the interval [0,pi]
## y(0)=1, y'(0)=0 and y'(pi)=0
## The 2nd order problem is rewritten as 2 firstorder problems:
## dy=y2
## dy2= (lam10cos(2t))*y
## =============================================================================
mathieu< function(t, y, lam)
{
list(c(y[2], (lam  10 * cos(2*t)) * y[1]))
}
yini < c(1, 0) # initial condition(y1=1,dy=y2=0)
yend < c(NA, 0) # final condition at pi (y1=NA, dy=0)
# there is one extra parameter to be fitted: "lam"; its initial guess = 15
Sol < bvpshoot(yini = yini, yend = yend, x = seq(0, pi, by = 0.01),
func = mathieu, guess = NULL, extra = 15)
plot(Sol)
attr(Sol, "roots") # root gives the value of "lam" (17.1068)

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.