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