bvpcol: Solves multipoint boundary value problems of ordinary...
In bvpSolve: Solvers for Boundary Value Problems of Differential Equations
Description
Usage
Arguments
Details
Value
Note
Author(s)
References
See Also
Examples
Description
Solves Boundary Value Problems For Ordinary Differential Equations (ODE)
or semi-explicit Differential-Algebraic Equations (DAE) with index at most 2.
It is possible to solve stiff ODE systems, by using an automatic continuation strategy
This is an implementation of the fortran codes colsys.f, colnew.f and coldae.f written by
respectively U. Ascher, J. christiansen and R.D. Russell (colsys),
U. Ascher and G. Bader (colnew) and U. Ascher and C. Spiteri.
The continuation strategy is an implementation of the fortran code colmod
written by J.R. Cash, M.H. Wright and F. Mazzia.
Usage
1
2
3
4
5
6
7
8
9
bvpcol (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, colp = NULL, bspline = FALSE,
fullOut = 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, dae = NULL, ...)
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 values of the equations. In case where the equations
are first-order, this will be the derivatives of y with respect to
x. After this can come global values that are required at
each point in x.
If the problem is a DAE, then the algebraic equations should be the last.
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.
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 of
bvptwp.
yguess
First guess values of y, corresponding to initial grid
xguess; a matrix with number of rows equal to the number of
variables, 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.
It is also allowed to pass the output of a previous run for continuation.
This will use the information that is stored in the attributes
istate and rstate. It will only work when for the previous
run, fullOut was set equal to TRUE (the default).
In this case, xguess need not be provided.
See example 3b.
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 examples.
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 - its values should be sorted - and
it should be within the range of x;
(posbound corresponds to fortran input "Zeta" in the colnew/colsys
FORTRAN codes. )
See last example.
Note that two-point boundary value problems can also be specified via
leftbc (which is simpler).
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.
For higher-order derivatives, specifying the order can improve computational
efficiency, but this interface is more complex.
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 3 and 3b.
ncomp
used if the model is specified by compiled code, the number of
components (or equations). 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.
colp
number of collocation points per subinterval.
bspline
if FALSE, then code colnew is used the default,
if TRUE, then fortran code colsys is used. Code colnew
incorporates a new basis representation, while colsys uses b-splines.
fullOut
if set to TRUE, then the collocation output required
e.g. for continuation will be returned in attributes rwork and
iwork. Use attributes(out)\$rwork, attributes(out)\$rwork
to see their contents
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" or deSolve's package vignette
"compiledCode".
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
when TRUE, writes messages (warnings) to the screen.
dae
if the problem is a DAE, should be a list containing the
index of the problem and the number of algebraic equations nalg.
See example 5
...
additional arguments passed to the model functions.
Details
If eps does not have a value and dae = NULL,
then the method is based on an implementation
of the Collocation methods called "colnew" and
"colsys" to solve multi-point boundary value problems of ordinary
differential equations.
The ODEs and boundary conditions are made available through the
user-provided routines, func and vectors yini and yend
or (optionally) bound. bvpcol can also solve multipoint
boundary value problems (see one but last example).
The corresponding partial derivatives are optionally available through the
user-provided routines, jacfunc and jacbound. Default is that
they are automatically generated by R, using numerical differences.
The user-requested tolerance is provided through atol.
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.'
If eps does have a value, then the method is based on an implementation
of the Collocation methods called "colmod".
The type of problems which this is designed to solve typically
involve a small positive parameter 0 < eps << 1.
As eps becomes progressively smaller, the problem normally becomes
increasingly difficult to approximate numerically (for example, due
to the appearance of narrow transition layers in the profile of the
analytic solution).
The idea of continuation is to solve a chain of problems in which the
parameter eps decreases monotonically towards some desired value.
That is, a sequence of problems is attempted to be solved:
epsini > eps1 > eps2 > eps3 > ..... > eps > 0
where epsini is a user provided starting value and eps is a
user desired final value for the parameter.
If dae is not NULL, then it is assumed that a DAE has to be solved.
In that case, dae should contain give the index of the DAE and the number of algebraic
equations (nalg).
(this part comes from the comments in the code coldae).
With respect to the dae, it should be noted that the code does not explicitly
check the index of the problem, so if the index is > 2 then the code will
not work well.
The number of boundary conditions required is independent
of the index. it is the user's responsibility to ensure
that these conditions are consistent with the constraints.
The conditions at the left end point must include
a subset equivalent to specifying the index-2
constraints there.
For an index-2 problem in hessenberg form, the projected
collocation method of Ascher and Petzold [2] is used.
Value
A matrix of class bvpSolve, with up to as many rows as elements in
x and as many columns
as elements in yini plus the number of "global" values returned
in the second element 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 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.
The output will also have attributes istate and rstate
which contain the collocation output required e.g. for continuation of a
problem, unless fullOutput is FALSE
Note
colnew.f (Bader and Ascher, 1987), is a modification of the code
colsys.f (Ascher, Christiansen and Russell, 1981), which incorporates
a new basis representation replacing b-splines, and improvements for
the linear and nonlinear algebraic equation solvers. To toggle on/off
colsys, set bspline = TRUE/FALSE
colmod is a revised version of the package colnew by
Bader and Ascher (1987), which in turn is a modification of the package
colsys by Ascher, Christiansen and Russell (1981). Colmod has been
adapted to allow an automatic continuation strategy to
be used (Cash et al., 1995).
The mesh selection algorithm used in colmod differs from that used in colnew
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
References
U. Ascher, J. Christiansen and R. D. Russell, (1981)
collocation software for boundary-value odes,
acm trans. math software 7, 209-222.
G. Bader and U. Ascher, (1987)
a new basis implementation for a mixed order
boundary value ode solver,
siam j. scient. stat. comput. 8, 487-483.
U. Ascher, J. Christiansen and R.D. Russell, (1979)
a collocation solver for mixed order
systems of boundary value problems,
math. comp. 33, 659-679.
U. Ascher, J. Christiansen and R.D. Russell, (1979)
colsys - a collocation code for boundary value problems,
lecture notes comp.sc. 76, springer verlag,
B. Childs et. al. (eds.), 164-185.
J. R. Cash, G. Moore and R. W. Wright, (1995)
an automatic continuation strategy for the solution of
singularly perturbed linear two-point boundary value problems,
j. comp. phys. 122, 266-279.
U. Ascher and R. Spiteri, 1994.
collocation software for boundary value differential-algebraic equations,
siam j. scient. stat. comput. 15, 938-952.
U. Ascher and L. Petzold, 1991.
projected implicit runge-kutta methods for differential-
algebraic equations,
siam j. num. anal. 28 (1991), 1097-1120.
See Also
bvpshoot for the shooting method
bvptwp for a MIRK formula
diagnostics.bvpSolve, for a description of diagnostic messages
approx.bvpSolve, for approximating solution in new values
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
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
## =============================================================================
## 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(-0.1 / sqrt(p+0.01), NA)
end <- c( 0.1 / sqrt(p+0.01), NA)
# Solve bvp
sol <- bvpcol(yini = init, yend = end,
x = seq(-0.1, 0.1, by = 0.001), func = fun)
plot(sol, which = 1)
# add analytical solution
curve(x/sqrt(p+x*x), add = TRUE, type = "p")
diagnostics(sol)
zoom <- approx(sol, xout = seq(-0.005, 0.005, by = 0.0001))
plot(zoom, which = 1, main = "zoom in on [-0.0005,0.0005]")
## =============================================================================
## Example 1b:
## Same problem, now solved as a second-order equation
## and with different value of "p".
## =============================================================================
fun2 <- function(t, y, pars)
{ dy <- - 3 * p * y[1] / (p+t*t)^2
list(dy)
}
p <- 1e-4
sol2 <- bvpcol(yini = init, yend = end, order = 2,
x = seq(-0.1, 0.1, by = 0.001), func = fun2)
# plot both runs at once:
plot(sol, sol2, which = 1)
## =============================================================================
## Example 1c: 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 <- bvpcol(yini = c(1, NA), yend = c(-0.5, NA), bspline = TRUE,
x = x, func = f2)
plot(sol, which = 1)
# 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 2. Uses continuation
## Test problem 24
## =============================================================================
Prob24<- function(t, y, ks) { #eps is called ks here
A <- 1+t*t
AA <- 2*t
ga <- 1.4
list(c(y[2],(((1+ga)/2 -ks*AA)*y[1]*y[2]-y[2]/y[1]-
(AA/A)*(1-(ga-1)*y[1]^2/2))/(ks*A*y[1])))
}
ini <- c(0.9129, NA)
end <- c(0.375, NA)
xguess <- c(0, 1)
yguess <- matrix(nrow = 2, ncol = 2, 0.9 )
# bvpcol works with eps NOT too small, and good initial condition ...
sol <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
xguess = xguess, yguess = yguess,
parms = 0.1, func = Prob24, verbose = FALSE)
# when continuation is used: does not need a good initial condition
sol2 <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 0.05, func = Prob24,
eps = 0.05)
#zoom <- approx(sol2, xout = seq(0.01, 0.02, by = 0.0001))
#plot(zoom, which = 1, main = "zoom in on [0.01, 0.02]")
sol3 <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 0.01, func = Prob24 , eps = 0.01)
sol4 <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 0.001, func = Prob24, eps = 0.001)
# This takes a long time
## Not run:
print(system.time(
sol5 <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 1e-4, func = Prob24, eps = 1e-4)
))
## End(Not run)
plot(sol, sol2, sol3, sol4, which = 1, main = "test problem 24",
lwd = 2)
legend("topright", col = 1:4, lty = 1:4, lwd = 2,
legend = c("0.1", "0.05", "0.01", "0.001"), title = "eps")
## =============================================================================
## Example 3 - 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)
## =============================================================================
# derivative function: 4 first-order derivatives
f1st<- function(x, y, S) {
list(c(y[2],
y[3],
y[4],
1/S*(y[2]*y[3] - y[1]*y[4]) ))
}
# jacobian of derivative function
df1st <- function(x, y, S) {
matrix(nrow = 4, ncol = 4, byrow = TRUE, data = c(
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1,
-1*y[4]/S, y[3]/S, y[2]/S, -y[1]/S))
}
# boundary
g2 <- function(i, y, S) {
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])
}
# jacobian of boundary
dg2 <- function(i, y, S) {
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))
}
# we use posbound to specify the position of boundary conditions
# we can also use leftbc = 2 rather than posbound = c(0,0,1,1)
S <- 1/100
sol <- bvpcol(x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f1st, parms = S, eps = S,
bound = g2, jacfunc = df1st, jacbound = dg2)
plot(sol)
## =============================================================================
## Example 3b - 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, S) {
list(1/S * (y[2]*y[3] - y[1]*y[4]))
}
# jacobian of derivative function
df4th <- function(x, y, S) {
matrix(nrow = 1, ncol = 4, byrow = TRUE, data = c(
-1*y[4]/S, y[3]/S, y[2]/S, -y[1]/S))
}
# boundary function - same as previous example
# jacobian of boundary - same as previous
# order = 4 specifies the equation to be 4th order
# solve with bspline false
S <- 1/100
sol <- bvpcol (x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f4th, order = 4,
parms = S, eps = S, bound = g2, jacfunc = df4th,
jacbound = dg2 )
plot(sol)
# Use (manual) continuation to find solution of a more difficult example
# Previous solution collocation from sol passed ("guess = sol")
sol2 <- bvpcol(x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f4th,
parms = 1e-6, order = 4, eps = 1e-6,
bound = g2, jacfunc = df4th, jacbound = dg2 )
# plot both at same time
plot(sol, sol2, lwd = 2)
legend("bottomright", leg = c(100, 10000), title = "R = ",
col = 1:2, lty = 1:2, lwd = 2)
## =============================================================================
## Example 4 - a multipoint bvp
## dy1 = (y2 - 1)/2
## dy2 = (y1*y2 - x)/mu
## over interval [0,1]
## y1(1) = 0; y2(0.5) = 1
## =============================================================================
multip <- function (x, y, p) {
list(c((y[2] - 1)/2,
(y[1]*y[2] - x)/mu))
}
bound <- function (i, y, p) {
if (i == 1) y[2] -1 # at x=0.5: y2=1
else y[1] # at x= 1: y1=0
}
mu <- 0.1
sol <- bvpcol(func = multip, bound = bound,
x = seq(0, 1, 0.01), posbound = c(0.5, 1))
plot(sol)
# check boundary value
sol[sol[,1] == 0.5,]
## =============================================================================
## Example 5 - a bvp DAE
## =============================================================================
bvpdae <- function(t, x, ks, ...) {
p1 <- p2 <- sin(t)
dp1 <- dp2 <- cos(t)
dx1 <- (ks + x[2] - p2)*x[4] + dp1
dx2 <- dp2
dx3 <- x[4]
res <- (x[1] - p1)*(x[4] - exp(t))
list(c(dx1, dx2, dx3, res), res = res)
}
boundfun <- function(i, x, par, ...) {
if (i == 1) return(x[1] - sin(0))
if (i == 2) return(x[3] - 1)
if (i == 3) return(x[2] - sin(1))
if (i == 4) return((x[1] - sin(1))*(x[4] - exp(1))) # Not used here..
}
x <- seq(0, 1, by = 0.01)
mass <- diag(nrow = 4) ; mass[4, 4] <- 0
# solved using boundfun
out <- bvpcol (func = bvpdae, bound = boundfun, x = x,
parms = 1e-4, ncomp = 4, leftbc = 2,
dae = list(index = 2, nalg = 1))
# solved using yini, yend
out1 <- bvpcol (func = bvpdae, x = x, parms = 1e-4,
yini = c(sin(0), NA, 1, NA),
yend = c(NA, sin(1), NA, NA),
dae = list(index = 2, nalg = 1))
# the analytic solution
ana <- cbind(x, "1" = sin(x), "2" = sin(x), "3" = 1, "4" = 0, res = 0)
plot(out, out1, obs = ana)
bvpSolve documentation built on Sept. 23, 2021, 3 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Note Author(s) References See Also Examples
Description
Solves Boundary Value Problems For Ordinary Differential Equations (ODE) or semi-explicit Differential-Algebraic Equations (DAE) with index at most 2.
It is possible to solve stiff ODE systems, by using an automatic continuation strategy
This is an implementation of the fortran codes colsys.f, colnew.f and coldae.f written by respectively U. Ascher, J. christiansen and R.D. Russell (colsys), U. Ascher and G. Bader (colnew) and U. Ascher and C. Spiteri.
The continuation strategy is an implementation of the fortran code colmod written by J.R. Cash, M.H. Wright and F. Mazzia.
Usage
1 2 3 4 5 6 7 8 9 | bvpcol (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, colp = NULL, bspline = FALSE,
fullOut = 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, dae = NULL, ...)
|
Arguments
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 R-function that computes the values of the
derivatives in the ODE system (the model definition) at point If The return value of If the problem is a DAE, then the algebraic equations should be the last. 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 It is also allowed to pass the output of a previous run for continuation.
This will use the information that is stored in the attributes
See example 3b. |
jacfunc |
jacobian (optional) - either an R-function 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 For higher-order derivatives, specifying the order can improve computational efficiency, but this interface is more complex. If |
ncomp |
used if the model is specified by compiled code, the number of
components (or equations). See package vignette Also to be used if the boundary conditions are specified by |
atol |
error tolerance, a scalar. |
colp |
number of collocation points per subinterval. |
bspline |
if |
fullOut |
if set to |
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 dll-functions whose names are
specified by |
ipar |
only when ‘dllname’ is specified: a vector with
integer values passed to the dll-functions 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 two-columned 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 |
when |
dae |
if the problem is a DAE, should be a list containing the
See example 5 |
... |
additional arguments passed to the model functions. |
Details
If eps does not have a value and dae = NULL,
then the method is based on an implementation
of the Collocation methods called "colnew" and
"colsys" to solve multi-point boundary value problems of ordinary
differential equations.
The ODEs and boundary conditions are made available through the
user-provided routines, func and vectors yini and yend
or (optionally) bound. bvpcol can also solve multipoint
boundary value problems (see one but last example).
The corresponding partial derivatives are optionally available through the
user-provided routines, jacfunc and jacbound. Default is that
they are automatically generated by R, using numerical differences.
The user-requested tolerance is provided through atol.
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.'
If eps does have a value, then the method is based on an implementation
of the Collocation methods called "colmod".
The type of problems which this is designed to solve typically
involve a small positive parameter 0 < eps << 1.
As eps becomes progressively smaller, the problem normally becomes
increasingly difficult to approximate numerically (for example, due
to the appearance of narrow transition layers in the profile of the
analytic solution).
The idea of continuation is to solve a chain of problems in which the parameter eps decreases monotonically towards some desired value. That is, a sequence of problems is attempted to be solved:
epsini > eps1 > eps2 > eps3 > ..... > eps > 0
where epsini is a user provided starting value and eps is a
user desired final value for the parameter.
If dae is not NULL, then it is assumed that a DAE has to be solved.
In that case, dae should contain give the index of the DAE and the number of algebraic
equations (nalg).
(this part comes from the comments in the code coldae). With respect to the dae, it should be noted that the code does not explicitly check the index of the problem, so if the index is > 2 then the code will not work well. The number of boundary conditions required is independent of the index. it is the user's responsibility to ensure that these conditions are consistent with the constraints. The conditions at the left end point must include a subset equivalent to specifying the index-2 constraints there. For an index-2 problem in hessenberg form, the projected collocation method of Ascher and Petzold [2] is used.
Value
A matrix of class bvpSolve, with up to as many rows as elements in
x and as many columns
as elements in yini plus the number of "global" values returned
in the second element 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 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.
The output will also have attributes istate and rstate
which contain the collocation output required e.g. for continuation of a
problem, unless fullOutput is FALSE
Note
colnew.f (Bader and Ascher, 1987), is a modification of the code
colsys.f (Ascher, Christiansen and Russell, 1981), which incorporates
a new basis representation replacing b-splines, and improvements for
the linear and nonlinear algebraic equation solvers. To toggle on/off
colsys, set bspline = TRUE/FALSE
colmod is a revised version of the package colnew by Bader and Ascher (1987), which in turn is a modification of the package colsys by Ascher, Christiansen and Russell (1981). Colmod has been adapted to allow an automatic continuation strategy to be used (Cash et al., 1995).
The mesh selection algorithm used in colmod differs from that used in colnew
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
References
U. Ascher, J. Christiansen and R. D. Russell, (1981) collocation software for boundary-value odes, acm trans. math software 7, 209-222.
G. Bader and U. Ascher, (1987) a new basis implementation for a mixed order boundary value ode solver, siam j. scient. stat. comput. 8, 487-483.
U. Ascher, J. Christiansen and R.D. Russell, (1979) a collocation solver for mixed order systems of boundary value problems, math. comp. 33, 659-679.
U. Ascher, J. Christiansen and R.D. Russell, (1979) colsys - a collocation code for boundary value problems, lecture notes comp.sc. 76, springer verlag, B. Childs et. al. (eds.), 164-185.
J. R. Cash, G. Moore and R. W. Wright, (1995) an automatic continuation strategy for the solution of singularly perturbed linear two-point boundary value problems, j. comp. phys. 122, 266-279.
U. Ascher and R. Spiteri, 1994. collocation software for boundary value differential-algebraic equations, siam j. scient. stat. comput. 15, 938-952.
U. Ascher and L. Petzold, 1991. projected implicit runge-kutta methods for differential- algebraic equations, siam j. num. anal. 28 (1991), 1097-1120.
See Also
bvpshoot for the shooting method
bvptwp for a MIRK formula
diagnostics.bvpSolve, for a description of diagnostic messages
approx.bvpSolve, for approximating solution in new values
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 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 | ## =============================================================================
## 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(-0.1 / sqrt(p+0.01), NA)
end <- c( 0.1 / sqrt(p+0.01), NA)
# Solve bvp
sol <- bvpcol(yini = init, yend = end,
x = seq(-0.1, 0.1, by = 0.001), func = fun)
plot(sol, which = 1)
# add analytical solution
curve(x/sqrt(p+x*x), add = TRUE, type = "p")
diagnostics(sol)
zoom <- approx(sol, xout = seq(-0.005, 0.005, by = 0.0001))
plot(zoom, which = 1, main = "zoom in on [-0.0005,0.0005]")
## =============================================================================
## Example 1b:
## Same problem, now solved as a second-order equation
## and with different value of "p".
## =============================================================================
fun2 <- function(t, y, pars)
{ dy <- - 3 * p * y[1] / (p+t*t)^2
list(dy)
}
p <- 1e-4
sol2 <- bvpcol(yini = init, yend = end, order = 2,
x = seq(-0.1, 0.1, by = 0.001), func = fun2)
# plot both runs at once:
plot(sol, sol2, which = 1)
## =============================================================================
## Example 1c: 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 <- bvpcol(yini = c(1, NA), yend = c(-0.5, NA), bspline = TRUE,
x = x, func = f2)
plot(sol, which = 1)
# 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 2. Uses continuation
## Test problem 24
## =============================================================================
Prob24<- function(t, y, ks) { #eps is called ks here
A <- 1+t*t
AA <- 2*t
ga <- 1.4
list(c(y[2],(((1+ga)/2 -ks*AA)*y[1]*y[2]-y[2]/y[1]-
(AA/A)*(1-(ga-1)*y[1]^2/2))/(ks*A*y[1])))
}
ini <- c(0.9129, NA)
end <- c(0.375, NA)
xguess <- c(0, 1)
yguess <- matrix(nrow = 2, ncol = 2, 0.9 )
# bvpcol works with eps NOT too small, and good initial condition ...
sol <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
xguess = xguess, yguess = yguess,
parms = 0.1, func = Prob24, verbose = FALSE)
# when continuation is used: does not need a good initial condition
sol2 <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 0.05, func = Prob24,
eps = 0.05)
#zoom <- approx(sol2, xout = seq(0.01, 0.02, by = 0.0001))
#plot(zoom, which = 1, main = "zoom in on [0.01, 0.02]")
sol3 <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 0.01, func = Prob24 , eps = 0.01)
sol4 <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 0.001, func = Prob24, eps = 0.001)
# This takes a long time
## Not run:
print(system.time(
sol5 <- bvpcol(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
parms = 1e-4, func = Prob24, eps = 1e-4)
))
## End(Not run)
plot(sol, sol2, sol3, sol4, which = 1, main = "test problem 24",
lwd = 2)
legend("topright", col = 1:4, lty = 1:4, lwd = 2,
legend = c("0.1", "0.05", "0.01", "0.001"), title = "eps")
## =============================================================================
## Example 3 - 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)
## =============================================================================
# derivative function: 4 first-order derivatives
f1st<- function(x, y, S) {
list(c(y[2],
y[3],
y[4],
1/S*(y[2]*y[3] - y[1]*y[4]) ))
}
# jacobian of derivative function
df1st <- function(x, y, S) {
matrix(nrow = 4, ncol = 4, byrow = TRUE, data = c(
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1,
-1*y[4]/S, y[3]/S, y[2]/S, -y[1]/S))
}
# boundary
g2 <- function(i, y, S) {
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])
}
# jacobian of boundary
dg2 <- function(i, y, S) {
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))
}
# we use posbound to specify the position of boundary conditions
# we can also use leftbc = 2 rather than posbound = c(0,0,1,1)
S <- 1/100
sol <- bvpcol(x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f1st, parms = S, eps = S,
bound = g2, jacfunc = df1st, jacbound = dg2)
plot(sol)
## =============================================================================
## Example 3b - 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, S) {
list(1/S * (y[2]*y[3] - y[1]*y[4]))
}
# jacobian of derivative function
df4th <- function(x, y, S) {
matrix(nrow = 1, ncol = 4, byrow = TRUE, data = c(
-1*y[4]/S, y[3]/S, y[2]/S, -y[1]/S))
}
# boundary function - same as previous example
# jacobian of boundary - same as previous
# order = 4 specifies the equation to be 4th order
# solve with bspline false
S <- 1/100
sol <- bvpcol (x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f4th, order = 4,
parms = S, eps = S, bound = g2, jacfunc = df4th,
jacbound = dg2 )
plot(sol)
# Use (manual) continuation to find solution of a more difficult example
# Previous solution collocation from sol passed ("guess = sol")
sol2 <- bvpcol(x = seq(0, 1, by = 0.01),
ynames = c("y", "dy", "d2y", "d3y"),
posbound = c(0, 0, 1, 1), func = f4th,
parms = 1e-6, order = 4, eps = 1e-6,
bound = g2, jacfunc = df4th, jacbound = dg2 )
# plot both at same time
plot(sol, sol2, lwd = 2)
legend("bottomright", leg = c(100, 10000), title = "R = ",
col = 1:2, lty = 1:2, lwd = 2)
## =============================================================================
## Example 4 - a multipoint bvp
## dy1 = (y2 - 1)/2
## dy2 = (y1*y2 - x)/mu
## over interval [0,1]
## y1(1) = 0; y2(0.5) = 1
## =============================================================================
multip <- function (x, y, p) {
list(c((y[2] - 1)/2,
(y[1]*y[2] - x)/mu))
}
bound <- function (i, y, p) {
if (i == 1) y[2] -1 # at x=0.5: y2=1
else y[1] # at x= 1: y1=0
}
mu <- 0.1
sol <- bvpcol(func = multip, bound = bound,
x = seq(0, 1, 0.01), posbound = c(0.5, 1))
plot(sol)
# check boundary value
sol[sol[,1] == 0.5,]
## =============================================================================
## Example 5 - a bvp DAE
## =============================================================================
bvpdae <- function(t, x, ks, ...) {
p1 <- p2 <- sin(t)
dp1 <- dp2 <- cos(t)
dx1 <- (ks + x[2] - p2)*x[4] + dp1
dx2 <- dp2
dx3 <- x[4]
res <- (x[1] - p1)*(x[4] - exp(t))
list(c(dx1, dx2, dx3, res), res = res)
}
boundfun <- function(i, x, par, ...) {
if (i == 1) return(x[1] - sin(0))
if (i == 2) return(x[3] - 1)
if (i == 3) return(x[2] - sin(1))
if (i == 4) return((x[1] - sin(1))*(x[4] - exp(1))) # Not used here..
}
x <- seq(0, 1, by = 0.01)
mass <- diag(nrow = 4) ; mass[4, 4] <- 0
# solved using boundfun
out <- bvpcol (func = bvpdae, bound = boundfun, x = x,
parms = 1e-4, ncomp = 4, leftbc = 2,
dae = list(index = 2, nalg = 1))
# solved using yini, yend
out1 <- bvpcol (func = bvpdae, x = x, parms = 1e-4,
yini = c(sin(0), NA, 1, NA),
yend = c(NA, sin(1), NA, NA),
dae = list(index = 2, nalg = 1))
# the analytic solution
ana <- cbind(x, "1" = sin(x), "2" = sin(x), "3" = 1, "4" = 0, res = 0)
plot(out, out1, obs = ana)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.