CODEWITHSUNDEEP
requation-solving
Amir Aliakbari
Amir Aliakbari

Reputation: 3

nleqslv, R, nonlinear equation system

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.

Upvotes: 0

Views: 2878

Answers (1)

Rui Barradas
Rui Barradas

Reputation: 76651

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.

Upvotes: 2

Related Questions