0

I have the following nonlinear equation system:

r1*r2*(0.25*0.6061-0.5)+r1+(0.25*r2) = 0.6061
r1*r2*(0.25*0.6429-0.5)+r1+(0.25*r2) = 0.6429

I need to find the r1 and r2. Their solution should be equal to: 

r1 = 0.5754
r2 = 0.6201

which resulted from Fortran.
I don't know how can I find the solutions using nleqslv package.

I will appreciate it if anybody can help me.

Thanks.

Karolis Koncevičius
  • 9,417
  • 9
  • 56
  • 89
Amir Aliakbari
  • 3
  • 3
  • Are you sure of those solutions? I am getting values close to `c(r1, r2) == c(1, 4)`. – Rui Barradas Apr 29 '18 at 17:40
  • The solution obtained using NAG (C05PBF) which is Fortran program. By the way, I couldn't understand how you performed the Jacobian – Amir Aliakbari Apr 29 '18 at 19:02
  • The solutions you gave are not correct.. Try solving the second equation with the solutions you have.. They do not work – Onyambu Apr 29 '18 at 19:22

1 Answers1

2

The way to use nleqslv would be to define a function returning a 2-dim vector, just like c(r1, r2).
That function takes a 2-dim argument x. I have replaced r1 and r2 by x[1] and x[2], respectively.
I have also replaced your 1st equation constant (0.25*0.6061-0.5) by its value, -0.348475. And in the 2nd equation, (0.25*0.6429-0.5) by -0.339275.
As for the initial values, I have first tried c(0, 0) but it didn't work. If you run the code below with initial values c(0.5, 0.5) you will get the same solution, within floating point precision.

fn <- function(x){
    c(-0.348475*x[1]*x[2] + x[1] + 0.25*x[2] - 0.6061,
      -0.339275*x[1]*x[2] + x[1] + 0.25*x[2] - 0.6429)
}

nleqslv(c(1, 1), fn)
#$`x`
#[1] 1.000000 3.999999
#
#$fvec
#[1] 4.773959e-14 4.607426e-14
#
#$termcd
#[1] 1
#
#$message
#[1] "Function criterion near zero"
#
#$scalex
#[1] 1 1
#
#$nfcnt
#[1] 2
#
#$njcnt
#[1] 1
#
#$iter
#[1] 2

Note that $fvec is near zero, and that $termcd is 1. This means that the algorithm converged. To have just the solution you can run

sol <- nleqslv(c(1, 1), fn)
sol$x
#[1] 1.000000 3.999999

If you can compute the jacobian, and in this case it is even very easy, the algorithm is more accurate. 

jacfn <- function(x) {
    n <- length(x)
    Df <- matrix(numeric(n*n), n, n)
    Df[1,1] <- -0.348475*x[2] + 1
    Df[1,2] <- -0.348475*x[1] + 0.25
    Df[2,1] <- -0.339275*x[2] + 1
    Df[2,2] <- -0.339275*x[1] + 0.25

    Df
}

sol2 <- nleqslv(c(1, 1), fn, jac = jacfn)
sol2$x
#[1] 1 4

The solution is the same.

Rui Barradas
  • 70,273
  • 8
  • 34
  • 66
  • The solution obtained using NAG (C05PBF) which is Fortran program. By the way, I couldn't understand how you performed the Jacobian. – Amir Aliakbari Apr 29 '18 at 18:25
  • @AmirAliakbari The Jacobian is simple, create a matrix of dimension equal to `length(x)`, then fill it with the partial derivatives of the objective function. As for your solution, try it in the equation. `fn(c(r1, r2))` outputs `[1] -1.281055e-05 -3.353020e-02`, far from zero. – Rui Barradas Apr 29 '18 at 21:17