# Solves for n roots of n (nonlinear) equations. 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. ```r 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, ...) ``` ## 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 -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 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 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`. ## Returns 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`))). ## Author(s) Karline Soetaert ## Examples ```r ## ======================================================================= ## 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 ``` ## 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 ) $<=$ 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}\bold{rtol} * abs(\bold{x}) + \bold{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. ## 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. ## Note The Fortran implementation of the Newton-Raphson method function (the default) is generally faster than the implementation. The implementation has been included for didactic purposes. `multiroot` makes use of function `stode`. Technically, it is just a wrapper around function `stode`. If the sparsity structure of the Jacobian is known, it may be more efficiently to call `stode, stodes, steady, steady.1D, steady.2D, steady.3D`. It is NOT guaranteed that the method will converge to the root.

multiroot function