How to solve a non-linear equation using nleqslv p

2019-02-16 00:24发布

问题:

I have this equation to solve (e.g. f(x,y) = 0):

library(nleqslv)
target <- function(x)
{
  z = x[1]/(x[1]+x[2])
  y = numeric(2)
  y[1] <- z*exp(-x[2]*(x[2]+z*(1-exp(-x[1]/z))))-0.00680
  y[2] <- z/x[2]*(1-exp(-x[2]))-exp(-x[2])*z/x[1]*(1-exp(-x[1]))-3.43164
  y
}

# Usage
xstart <- c(1,1)
target(xstart)
nleqslv(xstart, target, control=list(ftol=.0001, allowSingular=TRUE),jacobian=TRUE,method="Newton")

using R with nleqslv or another you have others :)

Thanks

回答1:

I have been experimenting with your function. Rewrite the target function to use the a;b constants in your comment as in your second comment as follows:

target <- function(x, a=.00680,b=3.43164)
{
  z <- x[1]/(x[1]+x[2])
  y <- numeric(2)
  y[1] <- z*exp(-x[2]*(x[2]+z*(1-exp(-x[1]/z))))-a
  y[2] <- z/x[2]*(1-exp(-x[2]))-exp(-x[2])*z/x[1]*(1-exp(-x[1]))-b
  y
}

The default values for a and b are what you initially specified. Using them you'll get an ill-conditioned jacobiam.

However if we give some other values to a and b for example like so

nleqslv(xstart, target, control=list(btol=.01),jacobian=TRUE,method="Newton",a=2,b=1)

or

nleqslv(xstart, target, control=list(btol=.01),jacobian=TRUE,method="Newton",a=2,b=2)

then for the first expression the full return value of nleqslv is

$x
[1]  2.4024092 -0.7498464

$fvec
[1] 1.332268e-15 2.220446e-16

$termcd
[1] 1

$message
[1] "Function criterion near zero"

$scalex
[1] 1 1

$nfcnt
[1] 10

$njcnt
[1] 7

$iter
[1] 7

$jac
           [,1]       [,2]
[1,] -0.2930082 -1.2103174
[2,]  0.1801120 -0.6566861

I am inclined to conclude that either your function is incorrect or that you have specified impossible values for a and b. Method Broyden also seems to work nicely.