Hi Xin,
You can use the nlsolve() function that I have written to solve a nonlinear
system of equations. It converts a root finding problem into a minimization
problem, and uses optim() to find the minimizer.
A well-known problem with this approach to root-finding is that a local
minimum of the squared residual function is not necessarily a root of the
nonlinear system. A nice feature of nlsolve() is that it will try "really
hard" (by increasing the "nstarts" argument) to find the solution to the
system by identifying multiple local minima of the squared residual
function, if they exist.
Here is a solution to your problem (I rearranged your equations a bit).
func <- function(x, A, B, pz) {
f <- rep(NA,length(x))
f[1] <- A/(1-pz^10) - x[2]*(1-x[1])/(x[1]*(1-x[1]^x[2]))
f[2] <- A + B/A - 1/x[1] - x[2]*(1-x[1])/x[1]
f
}
> p0 <- c(0.2,10)
> set.seed(123)
> nlsolve(par=p0, fn=func, A=327.727, B=9517.336, pz=0.114^10, nstart=1000)
$par
[1] 0.03443476 11.68745035
$value
[1] 7.01806e-05
$counts
function gradient
148 15
$convergence
[1] 0
$message
NULL
Best,
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
----------------------------------------------------------------------------
--------
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Xin
Sent: Wednesday, December 19, 2007 5:12 PM
To: R-Help
Subject: [R] can optimize solve paired euqations?
I used the command below, but R gives me the error message--syntax error.
can anyone see the mistakes I made?
optimize(function(x,y)
+ ((327.727-(1-0.114^10)*y*(1-x)/x/(1-x^y))+(9517.336-327.727
*(1+(1-x)*(1+y)/x-327.727)))^2
+ interval=c(0,1))
At the same time, I use nlm() but R gives me the code
$code
[1] 3
function(vals) {
x <- vals[1]
y <- vals[2]
sum(c((199.458913633542-(1-0.114^10)*y*(1-x)/x/(1-x^y)),(5234.11964684527-19
9.458913633542*(1+(1-x)*(1+y)/x-199.458913633542)))^2)
}
nlm(f, c(2,2))
Any one know how to deal with this case?
Thanks
Xin
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org 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.
##########################################
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()
#
# Author: Ravi Varadhan, Center on Aging and Health, Johns Hopkins
University, [EMAIL PROTECTED]
#
# June 21, 2007
#
func <- function(x, ...) sum(fn(x, ...)^2)
ans <- try(optim(par, fn=func, method=method, ...), silent=T)
if (class(ans) != "try-error") {
err <- sqrt(ans$val/length(par))
istart <- 1
} else {
err <- Inf
istart <- 0
}
p0 <- par
ntry <- 1
while (err > 0.01 & istart < nstarts & ntry < nstarts*10) {
par <- p0 * ( 1 + runif(length(par), -1, 1) ) # generating random starting
values
ans <- try(optim(par, fn=func, method=method, ...), silent=T)
ntry <- ntry + 1
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)
########################
# An example of Xin
#
func2 <- function(x,A, B, pz) {
f <- rep(NA,length(x))
f[1] <- A/(1-pz^10) - x[2]*(1-x[1])/(x[1]*(1-x[1]^x[2]))
f[2] <- A + B/A - 1/x[1] - x[2]*(1-x[1])/x[1]
f
}
p0 <- c(0.2,10)
set.seed(123)
nlsolve(par=p0, fn=func2, A=327.727, B=9517.336, pz=0.114^10, nstart=1000)
______________________________________________
R-help@r-project.org 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.
