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
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.