bonifaz - 3 months ago 29

R Question

I want to solve a system of non-linear equations

`x1 = f(x1, x2)`

x2 = g(x1, x2)

subject to

`x1 >= 0`

`x2 >= 0`

Using nleqslv, I set up the return vector of my function to be optimized as follows:

`y[1] = x[1] - f(x[1], x[2])`

y[2] = x[2] - g(x[1], x[2])

y[3] = -x[1]

y[4] = -x[2]

where the last two should reflect the non-negativy constraints.

Calling

"Jacobian is singular (1/condition=0.0e+000) (see allowSingular option)"

and calling with allowSingular = T yields:

"x-values within tolerance 'xtol'"

which makes sense, since y[4] doesn't change as x[4] is changed (by construction, doesn't respond to it at all).

How can I do this correctly?

Answer

Try this

```
f <- function(z) {
x <- z^2 # this will force x to be >= 0; other tranformations are possible
y <- numeric(2)
y[1] <- x[1] - f(x[1], x[2])
y[2] <- x[2] - g(x[1], x[2])
y
}
```

and then use `nleqslv`

to solve function `f`

with appropriate starting values for the transformed variables.
So if your starting values are now in `xstart`

use

```
zstart <- sqrt(xstart)
```

for the transformed problem.
And to get the result of solving in the correct units just use `^2`

on the result given by `nleqslv`

.