multiroot: Solves for n roots of n (nonlinear) equations.

Description

Given a vector of n variables, and a set of n (nonlinear) equations in these variables,

estimates the root of the equations, i.e. the variable values where all function values = 0.

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

Usage

multiroot(f, start, maxiter = 100, 
          rtol = 1e-6, atol = 1e-8, ctol = 1e-8, 
          useFortran = TRUE, positive = FALSE,
          jacfunc = NULL, jactype = "fullint",
          verbose = FALSE, bandup = 1, banddown = 1, 
          parms = NULL, ...)

Value

a list containing:

root: the location (x-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))).

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.

useFortran

logical, if FALSE, then an R -implementation of the Newton-Raphson method is used - see details.

positive

if TRUE, the unknowns y are forced to be non-negative numbers.

jacfunc

if not NULL, a user-supplied R function that estimates the Jacobian of the system of differential equations dydot(i)/dy(j). In some circumstances, supplying jacfunc can speed up the computations. The R calling sequence for jacfunc is identical to that of f.

If the Jacobian is a full matrix, jacfunc should return a matrix dydot/dy, where the ith row contains the derivative of $dy_i/dt$ with respect to $y_j$, or a vector containing the matrix elements by columns.

If the Jacobian is banded, jacfunc should return a matrix containing only the nonzero bands of the jacobian, (dydot/dy), rotated row-wise.

jactype

the structure of the Jacobian, one of "fullint", "fullusr", "bandusr", "bandint", or "sparse" - either full or banded and estimated internally or by the user, or arbitrary sparse. If the latter, then the solver will call, stodes, else stode

If the Jacobian is arbitrarily "sparse", then it will be calculated by the solver (i.e. it is not possible to also specify jacfunc).

verbose

if TRUE: full output to the screen, e.g. will output the steady-state settings.

bandup

number of non-zero bands above the diagonal, in case the Jacobian is banded.

banddown

number of non-zero bands below the diagonal, in case the jacobian is banded.

parms

vector or list of parameters used in f or jacfunc.

...

additional arguments passed to function f.

Author

Karline Soetaert <karline.soetaert@nioz.nl>

Details

start gives the initial guess for each variable; different initial guesses may return different roots.

The input parameters rtol, and atol determine the error control performed by the solver.

The solver will control the vector e of estimated local errors in f, according to an inequality of the form max-norm of ( e/ewt ) $\leq$ 1, where ewt is a vector of positive error weights. The values of rtol and atol should all be non-negative.

The form of ewt is: $$\mathbf{rtol} \times \mathrm{abs}(\mathbf{f}) + \mathbf{atol}$$

where multiplication of two vectors is element-by-element.

In addition, the solver will stop if between two iterations, the maximal change in the values of x is less than ctol.

There is no checking whether the requested precision exceeds the capabilities of the machine.

Examples

Run this code

## =======================================================================
## example 1  
## 2 simultaneous equations
## =======================================================================

model <- function(x) c(F1 = x[1]^2+ x[2]^2 -1, 
                       F2 = x[1]^2- x[2]^2 +0.5)

(ss <- multiroot(f = model, start = c(1, 1)))

## =======================================================================
## example 2
## 3 equations, two solutions
## =======================================================================

model <- function(x) c(F1 = x[1] + x[2] + x[3]^2 - 12,
                       F2 = x[1]^2 - x[2] + x[3] - 2,
                       F3 = 2 * x[1] - x[2]^2 + x[3] - 1 )

# first solution
(ss <- multiroot(model, c(1, 1, 1), useFortran = FALSE))
(ss <- multiroot(f = model, start = c(1, 1, 1)))

# second solution; use different start values
(ss <- multiroot(model, c(0, 0, 0)))
model(ss$root)

## =======================================================================
## example 2b: same, but with parameters
## 3 equations, two solutions
## =======================================================================

model2 <- function(x, parms) 
      c(F1 = x[1] + x[2] + x[3]^2 - parms[1],
        F2 = x[1]^2 - x[2] + x[3] - parms[2],
        F3 = 2 * x[1] - x[2]^2 + x[3] - parms[3])

# first solution
parms <- c(12, 2, 1)
multiroot(model2, c(1, 1, 1), parms = parms)
multiroot(model2, c(0, 0, 0), parms = parms*2)

## =======================================================================
## example 3: find a matrix
## =======================================================================

f2<-function(x)   {
 X <- matrix(nrow = 5, x)
 X %*% X %*% X -matrix(nrow = 5, data = 1:25, byrow = TRUE)
}
x <- multiroot(f2, start = 1:25 )$root
X <- matrix(nrow = 5, x)

X%*%X%*%X

Run the code above in your browser using DataLab