solving a simple (?) system of nonlinear equations

Question

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?

score 2 · Accepted Answer · answered Jan 17 '21 at 08:44

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.

ah i see. why must `fn` return the **difference** between p and pstar? — invictus, Jan 17 '21 at 16:29
You want to solve a function `f(p,q)=0`. You have not defined the function `f`. However you have defined `pstar = expression(p,q)` and `qstar = expression(p,q)` where `p` and `q` are the values `nleqslv` thinks should be a solution of `pstar=0` and `qstar=0`. Your function returns `c(pstar,qstar)` and `nleqslv` tries to make that zero. You want the left `p` and `q` to be equal to the input `p` and `q`. That is why your function should return the difference. — Bhas, Jan 17 '21 at 19:21

Sandipan Dey · Answer 2 · 2023-08-12T00:05:53.937

We can solve the system of nonlinear equations using Newton's method too, where each iteration step is shown below:

Here x_n = [p_n, q_n], the solution obtained in the n-th iteration, with the function F and its Jacobian J_F defined properly, as done in the code below:

# the system of nonlinear equations
# p(1 + exp(λ(1-2q)) - 1 = 0    (1)
# q(1 + exp(λ(10p-1)) - 1 = 0   (2)
# with F = [p(1 + exp(λ(1-2q)) - 1, q(1 + exp(λ(10p-1)) - 1]
# with the Jacobian matrix J_F = [[1+exp(λ(1-2q)), -2pλ.exp(λ(1-2q))],
#                                [10qλ.exp(λ(10p-1)), 1+exp(λ(10p-1))]]

f <- function(x) {
  p <- x[1]
  q <- x[2]
  return (c(p + p*exp(lambda*(1-2*q))-1, q + q*exp(lambda*(10*p-1))-1))
}

jacobian_f <- function(x) {
  p <- x[1]
  q <- x[2]
  return (matrix(c(1+exp(lambda*(1-2*q)), -2*p*lambda*exp(lambda*(1-2*q)), 10*q*lambda*exp(lambda*(10*p-1)), 1+exp(lambda*(10*p-1))), nrow=2, byrow=T))
}

Now, the function Newton() implements the Newton's method given the initial value x0, the function F and the Jacobian JF, as shown below:

Newton <- function(x0, F, JF, niter=10) {
  x <- x0  # Set Newton initial solution to x0
  for (k in range(niter)) {
    x <- x - solve(JF(x), F(x))
  }
  return (x)
}

Finally, let's solve the system of equations:

lambda <- 1
Newton(c(0.1,0.3), f, jacobian_f)
# [1] 0.3157453 0.1084092

The following animation shows how the Newton's method converges to the solution (red points represent the updated solution for the corresponding iteration):

The convergence of the Newton's method is obtained only in 5 iterations:

iter         p         q
1    0 0.1000000 0.3000000
2    1 0.3757359 0.0863961
3    2 0.3157453 0.1084092
4    3 0.3128046 0.1064944
5    4 0.3127804 0.1064237
6    5 0.3127804 0.1064237

solving a simple (?) system of nonlinear equations

2 Answers2