Hi All,
Here is a modification of nlsolve() that I had submitted before. The new
version of nlsolve() has the option to generate multiple random starting
values in order to improve the chances of locating the zeros of the
nonlinear system. The last test example in the attached file, the
Freudenstein-Roth function, shows the usefulness of the multiple random
starts.
Hope this is useful,
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
##########################################
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, [EMAIL PROTECTED]
#
# 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)
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
