Solver for two-point boundary value problems of ordinary differential equations, using the single shooting method

Share:

Description

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 root-finding methods from package rootSolve.

Preferentially bvptwp or bvpcol should be used rather than bvpshoot, as they give more precise output.

Usage

1
2
3
4
5
6
bvpshoot(yini = NULL, x, func, yend = NULL, parms = NULL, 
         order = NULL, guess = NULL,
         jacfunc = NULL, bound = NULL, jacbound = NULL, 
         leftbc = NULL, posbound = NULL, ncomp = NULL, 
         atol = 1e-8, rtol = 1e-8, extra = NULL, 
         maxiter = 100, positive = FALSE, method = "lsoda",...)

Arguments

yini

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

If yini is a function, it should be defined as: yini <- function(y, parms,...); where y are the initial values, and parms the parameters.

if yini is a vector then use NA for an initial value which is not available.

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

if yini = NULL then bound should be specified; 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

an R-function that computes the values of the derivatives in the ODE system (the model definition) at x. func 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 or yend 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.

Note that it is not possible to use bvpshoot with functions defined in compiled code. Use bvptwp instead.

yend

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

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

If yend is a function, it should be defined as: yend <- function (y, yini, parms, ...); where y are the final values, yini the initial values and parms the parameters.

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 bound should be specified; 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).

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 examples.

guess

guess for the value(s) of the unknown initial conditions;

if initial and final conditions are specified by yini and yend, then guess should contain one value for each NA in yini. The length of guess should thus equal the number of unknown initial conditions (=NAs in yini). If guess is not provided, a value = 0 is assumed for each NA in yini and a warning printed.

If initial and final conditions are specified by the boundary function bound, then guess should contain the initial guess for all initial conditions, i.e. its length should equal the number of state variables in the ODE system; if in this case guess has a names attribute, the names will be available within the functions and used to label the output matrix. If guess is not provided, then ncomp should specify the total number of variables, a value = 0 will be assumed for the initial conditions and a warning printed.

jacfunc

jacobian (optional) - an R-function that evaluates the jacobian of func at point x.

jacfunc 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.

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

jacfunc is passed to the initial value problem solver.

bound

boundary function (optional) - only if yini and yend are not available. An R function that evaluates the i-th boundary element at point x.

bound should be defined as: bound = function(i, y, parms, ...). It should return the i-th boundary condition.

if not NULL, bound defines the root to be solved by the root solving algorithm.

jacbound

jacobian of the boundary function (optional) - only if bound is defined. An R function that evaluates the gradient of the i-th boundary element with respect to the state variables, at point x.

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

jacbound is passed to the root solver.

leftbc

only if yini and yend are not available: 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 points that are in x are allowed. Note that, if the boundary conditions are at the ends of the integration interval, it is simpler to use leftbc.

ncomp

only used if the boundaries are specified via the boundary function bound and guess is not specified. The number of components.

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.

extra should contain the initial guess of this extra parameter.

maxiter

the maximal number of iterations allowed in the root solver.

positive

set to TRUE if dependent variables (y) have to be positive numbers.

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.

Details

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 root-solving method. There are two strategies:

yini and yend specified

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.

bound specified

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.

Value

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 x-value.

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.

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.

Do not specify the jacobian if the maximal order>1.

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>

See Also

bvptwp for the MIRK method

lsoda, lsode, lsodes, vode,

rk, rkMethod for details about the integration method

multiroot, the root-solving method used

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
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 = 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  <- 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 + (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  <- 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(dy-Ry^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*(x-1))

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 higher-order 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*(x-1))

# 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(sol2-sol))

## =============================================================================
## 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 * (w-u)/v
   dv <- -0.5 * (w-u)
   dw <- (0.9 - 1000 * (w-y) - 0.5 * w * (w-u))/z
   dz <- 0.5 * (w-u)
   dy <- -100 * (y-w)
   return(list(c(du, dv, dw, dy, dz)))
  })
}

#--------------------------------
# Residuals of end conditions
#--------------------------------
res  <- function (Y, yini, pars)  with (as.list(Y), w-y)

#--------------------------------
# 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 = 1e-10, rtol = 0)
pairs(sol, main = "MUSN")

## =============================================================================
## Example 4 - solve also for unknown parameter
## Find the 4th eigenvalue of Mathieu's equation:
## y''+(lam-10cos2t)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 first-order problems:
## dy=y2
## dy2= -(lam-10cos(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)