##########################################
nlsolve <- function(par, fn, method="BFGS", nstarts = 1, ...) {
###########################
#  A solver for finding a zero of a system of non-linear equations
#  Actually minimizes the squared-norm of the set of functions by calling optim()
#  Uses the BFGS algorithm within optim()
#  All the control parameters can be passed as in the call to optim()
#  Uses a random multiple start approach to locate the zeros of the system
#
#  Author:  Ravi Varadhan, Center on Aging and Health, Johns Hopkins University, rvaradhan@jhmi.edu
#
#  June 21, 2007
#
func <- function(x, ...) sum(fn(x, ...)^2)
ans <- optim(par, fn=func, method=method, ...)
err <- sqrt(ans$val/length(par)) 
istart <- 1

while (err > 0.01 & istart < nstarts) {
par <- par * ( 1 + runif(length(par), -1, 1) )  # generating random starting values
ans <- try (optim(par, fn=func, method=method, ...))
if (class(ans) != "try-error") err <- sqrt(ans$val/length(par)) 
istart <- istart + 1
}

if (err > 0.01) cat(" \n Caution:  Zero of the function has not been located!!! \n Increase nstart \n \n")

ans
}

#
##########################################
#
# Some examples illustrating the use of nlsolve()
#
expo1 <- function(p) {
#  From La Cruz and Raydan, Optim Methods and Software 2003, 18 (583-599)
n <- length(p)
f <- rep(NA, n)
f[1] <- exp(p[1] - 1) - 1
f[2:n] <- (2:n) * (exp(p[2:n] - 1) - p[2:n])
f
}

n <- 50
p0 <- runif(n)
nlsolve(par=p0, fn=expo1)

############
expo3 <- function(p) {
#  From La Cruz and Raydan, Optim Methods and Software 2003, 18 (583-599)
n <- length(p)
f <- rep(NA, n)
onm1 <- 1:(n-1) 
f[onm1] <- onm1/10 * (1 - p[onm1]^2 - exp(-p[onm1]^2))
f[n] <- n/10 * (1 - exp(-p[n]^2))
f
}

n <- 15
p0 <- (1:n)/(4*n^2)
nlsolve(par=p0, fn=expo3)

############
func1 <- function(x) {
f <- rep(0,3)
f[1] <- 1/3*cos(x[2]*x[3]) - x[1] + 1/6
f[2] <- 1/9 * sqrt(x[1]^2 + sin(x[3])+ 1.06) - x[2] - 0.1
f[3] <- -1/20*exp(-x[1]*x[2]) - x[3] - (10*pi - 3)/60
f
}

p0 <- runif(3)
nlsolve(par=p0, fn=func1)

############
# This example illustrates the use of multiple random starts to obtain the global minimum of zero (i.e. the roots of the system)
#
froth <- function(p){
# Freudenstein and Roth function (Broyden, Mathematics of Computation 1965, p. 577-593)
f <- rep(NA,length(p))
f[1] <- -13 + p[1] + (p[2]*(5 - p[2]) - 2) * p[2]
f[2] <- -29 + p[1] + (p[2]*(1 + p[2]) - 14) * p[2]
f
}

p0 <- c(3,2)  # this gives the zero of the system
nlsolve(par=p0, fn=froth)

p0 <- c(1,1)  # this gives the local minimum that is not the zero of the system
nlsolve(par=p0, fn=froth)
nlsolve(par=p0, fn=froth, nstart=100)

p0 <- rpois(2,10)
nlsolve(par=p0, fn=froth)
nlsolve(par=p0, fn=froth, nstart=10)