multiroot.1D: Solves for n roots of n (nonlinear) equations, created by...

In rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations

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 method-of-lines approach.

Assumes a banded Jacobian matrix, uses the Newton-Raphson method.

Usage

multiroot.1D(f, start, maxiter = 100,
             rtol = 1e-6, atol = 1e-8, ctol = 1e-8, 
             nspec = NULL, dimens = NULL, verbose = FALSE, 
             positive = FALSE, names = NULL, parms = NULL, ...)

Arguments

f	function for which the root is sought; it must return a vector with as many values as the length of start. It is called either as f(x, ...) if parms = NULL or as f(x, parms, ...) if parms is not NULL.
start	vector containing initial guesses for the unknown x; if start has a name attribute, the names will be used to label the output vector.
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 NULL, then dimens should be specified.
dimens	the number of *boxes* in the model. If NULL, then nspec should be specified.
verbose	if TRUE: full output to the screen, e.g. will output the steady-state settings.
positive	if TRUE, the unknowns y are forced to be non-negative numbers.
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 f.
...	additional arguments passed to function f.

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,\frac{dy}{dx},\frac{d^2y}{dx^2})

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 second-order 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 second-order derivates as follows.

A first-order 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 second-order 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 root.
iter	the number of iterations used.
estim.precis	the estimated precision for root. It is defined as the mean of the absolute function values (mean(abs(f.root))).

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 <[e-mail]>

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, steady-state 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

## =======================================================================
## 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).
##
## =======================================================================

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  <- 1e-5
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*y-y
}

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")

rootSolve documentation built on Sept. 21, 2023, 5:06 p.m.