Solves for n roots of n (nonlinear) equations, created by discretizing ordinary differential equations.
Description
multiroot.1D finds the solution to boundary value problems of ordinary differential equations, which are approximated using the methodoflines approach.
Assumes a banded Jacobian matrix, uses the NewtonRaphson method.
Usage
1 2 3 4 
Arguments
f 
function for which the root is sought; it must return a vector
with as many values as the length of 
start 
vector containing initial guesses for the unknown x;
if 
maxiter 
maximal number of iterations allowed. 
rtol 
relative error tolerance, either a scalar or a vector, one value for each element in the unknown x. 
atol 
absolute error tolerance, either a scalar or a vector, one value for each element in x. 
ctol 
a scalar. If between two iterations, the maximal change in the variable values is less than this amount, then it is assumed that the root is found. 
nspec 
the number of *species* (components) in the model.
If 
dimens 
the number of *boxes* in the model. If NULL, then

verbose 
if 
positive 
if 
names 
the names of the components; used to label the output, which will be written as a matrix. 
parms 
vector or list of parameters used in 
... 
additional arguments passed to function 
Details
multiroot.1D
is similar to steady.1D
, except for the
function specification which is simpler in multiroot.1D
.
It is to be used to solve (simple) boundary value problems of differential equations.
The following differential equation:
0=f(x,y,y',y'')
with boundary conditions
y(x=a) = ya, at the start and y(x=b)=yb at the end of the integration interval [a,b] is approximated
as follows:
1. First the integration interval x is discretized, for instance:
dx < 0.01
x < seq(a,b,by=dx)
where dx
should be small enough to keep numerical errors small.
2. Then the first and secondorder
derivatives are differenced on this numerical
grid. R's diff
function is very efficient in taking numerical
differences, so it is used to approximate the first, and secondorder
derivates as follows.
A firstorder derivative y' can be approximated either as:
 y'=
diff(c(ya,y))/dx
if only the initial condition ya is prescribed,
 y'=
diff(c(y,yb))/dx
if only the final condition, yb is prescribed,
 y'=
0.5*(diff(c(ya,y))/dx+diff(c(y,yb))/dx)
if initial, ya, and final condition, yb are prescribed.
The latter (centered differences) is to be preferred.
A secondorder derivative y” can be approximated by differencing twice.
y”=diff(diff(c(ya,y,yb))/dx)/dx
3. Finally, function multiroot.1D
is used to locate the root.
See the examples
Value
a list containing:
root 
the values of the root. 
f.root 
the value of the function evaluated at the 
iter 
the number of iterations used. 
estim.precis 
the estimated precision for 
Note
multiroot.1D
makes use of function steady.1D
.
It is NOT guaranteed that the method will converge to the root.
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
See Also
stode
, which uses a different function call.
uniroot.all
, to solve for all roots of one (nonlinear) equation
steady
, steady.band
, steady.1D
,
steady.2D
, steady.3D
, steadystate solvers,
which find the roots of ODEs or PDEs. The function call differs from
multiroot
.
jacobian.full
, jacobian.band
, estimates the
Jacobian matrix assuming a full or banded structure.
gradient
, hessian
, estimates the gradient
matrix or the Hessian.
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  ## =======================================================================
## Example 1: simple standard problem
## solve the BVP ODE:
## d2y/dt^2=3py/(p+t^2)^2
## y(t= 0.1)=0.1/sqrt(p+0.01)
## y(t= 0.1)= 0.1/sqrt(p+0.01)
## where p = 1e5
##
## analytical solution y(t) = t/sqrt(p + t^2).
##
## =======================================================================
bvp < function(y) {
dy2 < diff(diff(c(ya, y, yb))/dx)/dx
return(dy2 + 3*p*y/(p+x^2)^2)
}
dx < 0.0001
x < seq(0.1, 0.1, by = dx)
p < 1e5
ya < 0.1/sqrt(p+0.01)
yb < 0.1/sqrt(p+0.01)
print(system.time(
y < multiroot.1D(start = runif(length(x)), f = bvp, nspec = 1)
))
plot(x, y$root, ylab = "y", main = "BVP test problem")
# add analytical solution
curve(x/sqrt(p+x*x), add = TRUE, type = "l", col = "red")
## =======================================================================
## Example 2: bvp test problem 28
## solve:
## xi*y'' + y*y'  y=0
## with boundary conditions:
## y0=1
## y1=3/2
## =======================================================================
prob28 <function(y, xi) {
dy2 < diff(diff(c(ya, y, yb))/dx)/dx # y''
dy < 0.5*(diff(c(ya, y)) +diff(c(y, yb)))/dx # y'  centered differences
xi*dy2 +dy*yy
}
ya < 1
yb < 3/2
dx < 0.001
x < seq(0, 1, by = dx)
N < length(x)
print(system.time(
Y1 < multiroot.1D(f = prob28, start = runif(N),
nspec = 1, xi = 0.1)
))
Y2< multiroot.1D(f = prob28, start = runif(N), nspec = 1, xi = 0.01)
Y3< multiroot.1D(f = prob28, start = runif(N), nspec = 1, xi = 0.001)
plot(x, Y3$root, type = "l", col = "green", main = "bvp test problem 28")
lines(x, Y2$root, col = "red")
lines(x, Y1$root, col = "blue")
