bvpshoot: Solver for two-point boundary value problems of ordinary...
In bvpSolve: Solvers for Boundary Value Problems of Differential Equations
Description
Usage
Arguments
Details
Value
Note
Author(s)
See Also
Examples
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
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
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)
Example output
Loading required package: deSolve
Attaching package: 'bvpSolve'
The following object is masked from 'package:stats':
approx
user system elapsed
0.017 0.000 0.018
root f.root iter
dy 0.6367841 -8.426954e-09 6
user system elapsed
0.026 0.000 0.026
root f.root iter
y 0.6367841 0.000000e+00 8
dy -0.3632159 -3.053113e-16 8
user system elapsed
0.02 0.00 0.02
root f.root iter
y 0.6367841 0.000000e+00 8
dy -0.3632159 -3.053113e-16 8
user system elapsed
0.023 0.000 0.023
[1] 0
root f.root iter
2 17.10683 -5.205281e-13 6
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) See Also Examples
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 |
Arguments
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 R-function 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 R-function 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. |
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
yiniandyendthen 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 bybound.
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
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)
|
Example output
Loading required package: deSolve
Attaching package: 'bvpSolve'
The following object is masked from 'package:stats':
approx
user system elapsed
0.017 0.000 0.018
root f.root iter
dy 0.6367841 -8.426954e-09 6
user system elapsed
0.026 0.000 0.026
root f.root iter
y 0.6367841 0.000000e+00 8
dy -0.3632159 -3.053113e-16 8
user system elapsed
0.02 0.00 0.02
root f.root iter
y 0.6367841 0.000000e+00 8
dy -0.3632159 -3.053113e-16 8
user system elapsed
0.023 0.000 0.023
[1] 0
root f.root iter
2 17.10683 -5.205281e-13 6
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.