solving a simple (?) system of nonlinear equations

I'm trying to solve a simple system of non-linear equations described in this post.

The system is two equations with two unknowns p and q and a free parameter lambda:

When lambda = 1 the system looks like this:

There is a unique solution and it's in the vicinity of p = 0.3, q = 0.1.

I'm trying to solve it with nleqslv. My objective function is:

library(nleqslv)

fn = function(x, lambda = 1){ 
  # p = x[1]
  # q = x[2]
  pstar = exp(lambda * (1*x[2])) / (exp(lambda * (1*x[2])) + exp(lambda * (1 - x[2])))
  qstar = exp(lambda * (1 - x[1])) / (exp(lambda * ((1 - x[1]))) + exp(lambda * (9*x[1])))
  return(c(pstar,qstar))
}

but the results don't match what the plot:

> xstart = c(0.1, 0.3)
> nleqslv(xstart, fn)$x
[1]  1.994155 -8.921285

My first question is: am I using nleqslv correctly? I thought so after looking at other examples. But now I'm not sure.

My second question: is this a good problem nleqslv? Or am I barking up the wrong tree?

Your function does not reflect properly what you want.

You can see this by evaluating fn(c(0.3,0.1)) as follows.

fn(c(0.3,0.1))
[1] 0.3100255 0.1192029

So the output is very close to the input. You wanted (almost) zero as output.

So you want to solve the system for p and q. What you need to do is to make your function return the difference between the input p and the expression for pstar and the difference between the input q and the expression for qstar.

So rewrite your function as follows

fn <- function(x, lambda = 1){ 
  p <- x[1]
  q <- x[2]
  pstar <- exp(lambda * (1*x[2])) / (exp(lambda * (1*x[2])) + exp(lambda * (1 -    x[2])))
  qstar <- exp(lambda * (1 - x[1])) / (exp(lambda * ((1 - x[1]))) + exp(lambda * (9*x[1])))
  return(c(pstar-p,qstar-q))
}  

and then call nleqslv as follows (PLEASE always show all the code you are using. You left out the library(nleqslv)).

library(nleqslv)
xstart <- c(0.1, 0.3)
nleqslv(xstart, fn)

This will display the full output of the function. Always a good idea to check for succes. Always check $termcd for succes.

$x
[1] 0.3127804 0.1064237

$fvec
[1] 5.070055e-11 6.547240e-09

$termcd
[1] 1

$message
[1] "Function criterion near zero"

$scalex
[1] 1 1

$nfcnt
[1] 7

$njcnt
[1] 1

$iter
[1] 7

The result for $x is more what you expect.

Finally please use <- for assignment. If you don't there will come the day that you will be bitten by R and its magic.

This is nothing wrong in using nleqslv for this problem. You only made a small mistake.