Broyden's Method — broyden • pracma

Broyden's method for the numerical solution of nonlinear systems of n equations in n variables.

broyden(Ffun, x0, J0 = NULL, ...,
        maxiter = 100, tol = .Machine$double.eps^(1/2))

Arguments

Ffun: n functions of n variables.
x0: Numeric vector of length n.
J0: Jacobian of the function at x0.
...: additional parameters passed to the function.
maxiter: Maximum number of iterations.
tol: Tolerance, relative accuracy.

Details

F as a function must return a vector of length n, and accept an n-dim. vector or column vector as input. F must not be univariate, that is n must be greater than 1.

Broyden's method computes the Jacobian and its inverse only at the first iteration, and does a rank-one update thereafter, applying the so-called Sherman-Morrison formula that computes the inverse of the sum of an invertible matrix A and the dyadic product, uv', of a column vector u and a row vector v'.

Value

List with components: zero the best root found so far, fnorm the square root of sum of squares of the values of f, and niter the number of iterations needed.

References

Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.

Note

Applied to a system of n linear equations it will stop in 2n steps

See also

newtonsys, fsolve

Examples

##  Example from Quarteroni & Saleri
F1 <- function(x) c(x[1]^2 + x[2]^2 - 1, sin(pi*x[1]/2) + x[2]^3)
broyden(F1, x0 = c(1, 1))
#> $zero
#> [1]  0.4760958 -0.8793934
#> 
#> $fnorm
#> [1] 9.092626e-09
#> 
#> $niter
#> [1] 13
#> 
# zero: 0.4760958 -0.8793934; fnorm: 9.092626e-09; niter: 13

F <- function(x) {
    x1 <- x[1]; x2 <- x[2]; x3 <- x[3]
    as.matrix(c(x1^2 + x2^2 + x3^2 - 1,
                x1^2 + x3^2 - 0.25,
                x1^2 + x2^2 - 4*x3), ncol = 1)
}
x0 <- as.matrix(c(1, 1, 1))
broyden(F, x0)
#> $zero
#> [1] 0.4407629 0.8660254 0.2360680
#> 
#> $fnorm
#> [1] 1.34325e-08
#> 
#> $niter
#> [1] 8
#> 
# zero: 0.4407629 0.8660254 0.2360680; fnorm: 1.34325e-08; niter: 8

##  Find the roots of the complex function sin(z)^2 + sqrt(z) - log(z)
F2 <- function(x) {
    z  <- x[1] + x[2]*1i
    fz <- sin(z)^2 + sqrt(z) - log(z)
    c(Re(fz), Im(fz))
}
broyden(F2, c(1, 1))
#> $zero
#> [1] 0.2555197 0.8948303
#> 
#> $fnorm
#> [1] 7.284373e-10
#> 
#> $niter
#> [1] 13
#> 
# zero   0.2555197 0.8948303 , i.e.  z0 = 0.2555 + 0.8948i
# fnorm  7.284374e-10
# niter  13

##  Two more problematic examples
F3 <- function(x)
        c(2*x[1] - x[2] - exp(-x[1]), -x[1] + 2*x[2] - exp(-x[2]))
broyden(F3, c(0, 0))
#> $zero
#> [1] 0.5671433 0.5671433
#> 
#> $fnorm
#> [1] 4.677305e-12
#> 
#> $niter
#> [1] 5
#> 
# $zero   0.5671433 0.5671433   # x = exp(-x)

F4 <- function(x)   # Dennis Schnabel
        c(x[1]^2 + x[2]^2 - 2, exp(x[1] - 1) + x[2]^3 - 2)
broyden(F4, c(2.0, 0.5), maxiter = 100)
#> $zero
#> [1] 1 1
#> 
#> $fnorm
#> [1] 1.188144e-10
#> 
#> $niter
#> [1] 43
#>