Search code examples
rsolvernon-linear

Solve non-linear equation using nleqslv package in R


I am new to R and I am trying to solve a non-linear equation. The equation is a function of 4 variables: Q, D, E2+z, and y1.new. I have values for Q, D, and E2+z stored in individual vectors with length=3364. For each of these cases, I want to find y1.new that satisfies the equation below

0 =  `y1.new` + (`Q`/((`D`^2)/8)*(2*acos(1-2*`y1.new`/`D`) - sin(2*acos(1-2*`y1.new`/`D`))))^2/(2*9.81) - `E2+z`.

I set a function as

target.y1.new <- function(y1.new, Q=Q,D=D,`E2+z`=`E2+z`,g=9.81)
{
  y <- numeric(1)
  y[1] <- `y1.new` + (`Q`/((`D`^2)/8)*(2*acos(1-2*`y1.new`/`D`) - sin(2*acos(1-2*`y1.new`/`D`))))^2/(2*`g`) - `E2+z`
}

my initial guesses are stored in a vector named y1, also with 3364 values

I tried to use the function

nleqslv(`y1`,`target.y1.new`,control=list(btol=0.01),jacobian=TRUE)

but it results in an error

Error in fn(par, ...) :
promise already under evaluation: recursive default argument reference or earlier problems?

Can anyone advise what I am doing wrong here?

Thanks in advance,
Benny


Solution

  • Rui is right. You are using backticks where they are not needed.

    Simplifying and correction your function target.y1.new results in this

    target.y1.new <- function(y1.new,g=9.81)
    {
      y <- numeric(length(y1.new))
      y <- y1.new + (Q/((D^2)/8)*(2*acos(1-2*y1.new/D) - sin(2*acos(1-2*y1.new/D))))^2/(2*g) - E2pz
      y
    }
    

    You say you have vectors of length 3364. Your function works with vectors. This means that your function should return a vector of length 3364 which is why the definition of the return value y should have the same length as the input vectors.

    With some trial values we can try this

    Q <- c(1,2)
    D <- c(1,2)
    E2pz <- c(1,2)
    y1 <- rep(.1,2) 
    

    You can now do the following

    nleqslv(y1,target.y1.new,control=list(btol=0.01),jacobian=TRUE)  
    

    This will give you a solution. But it is not the most efficient. You system of equations is diagonal in the sense that the input y of any equation does not occur in any other equation.

    nleqslv has an option for a specifying a band (including diagonal) jacobian (see the manual). As follows

    nleqslv(y1,target.y1.new,control=list(btol=0.01,dsub=0,dsuper=0),
            jacobian=TRUE,method="Newton")
    

    In this case it is better to use the Newton method since the default Broyden method will do too much work. And it does not work with a diagonal jacobian.