I have found the yield of a crop (Y) as a function of its nitrogen offtake (U) i.e., Y(U).

The rest of the values for this particular crop are:

I want to solve for U.

I tried solving this using a binary search algorithm, using the uniroot() and polyroot(), all to no avail :(

I tried defining it as

fn <- function(U) 
{
 Y - Y_crit - Q * (U-U_Crit)   ((Q/(p 1)) * ((U - U_crit)/(U_max - U_crit))^(p 1) * (U_max - U_crit)
}
        
U <- polyroot(fn)
print(U)

but it says: "Error in polyroot(fn) : unimplemented type 'closure' in 'polyroot'"

Please help. TIA!

EDIT: I had first presented the value of Y as 14170 (=Y_max) but then confusing it with data for another crop, changed it to 11000. I have now changed it back.

CodePudding user response：

I have rewritten the function in order to make it more clear what it is computing, there is an error in the question's version (wrong parenthesis).

• Pass the arguments values as one object and in the function coerce it to list. This will make argument passing easier and less error prone;
• repeating terms are pre-computed and reused.
• I have plotted the function with values starting with U = 123.79, the value in the data.frame, until a visual inspection found an interval where the root is.

fn <- function(U, args) {
  with(as.list(args), {
    term1 <- U - U_crit
    term2 <- U_max - U_crit
    lhs <- Y_crit   Q*term1 - Q/(p 1) * (term1/term2)^(p 1) * term2
    rhs <- Y
    return(lhs - rhs)
  })
}

U <- uniroot(fn, c(123.279, 350), args = args)
U
#> $root
#> [1] 308.6662
#> 
#> $f.root
#> [1] 0.0004746999
#> 
#> $iter
#> [1] 7
#> 
#> $init.it
#> [1] NA
#> 
#> $estim.prec
#> [1] 6.103516e-05

curve(fn(x, args), 123.3, 350, lwd = 2)
abline(h = 0)
points(U$root, U$f.root, col = "red", pch = 19)

^{Created on 2022-12-22 with reprex v2.0.2}

Edit

According to its documentation, package optimx

Provides a replacement and extension of the optim() function to call to several function minimization codes in R in a single statement.

But it only minimises the objective function, so write a wrapper around it, gn below.

``` r
library(optimx)

gn <- function(x0, args) {
  with(as.list(x0), {
    args$Y <- Y
    -fn(U, args)
  })
}

x0 <- c(U = 124, Y = 10000)
optimx(par = x0, gn, 
       method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B"), 
       args = args)
#>                        U             Y          value fevals gevals niter
#> Nelder-Mead 1.887090e 19 -7.002469e 34  -6.310914e 34    501     NA    NA
#> BFGS        1.917764e 02  8.128266e 03  -6.026305e 03    100    100    NA
#> CG          1.983800e 02  9.853717e 03  -4.315391e 03    201    101    NA
#> L-BFGS-B              NA            NA  8.988466e 307     NA     NA    NA
#>             convcode kkt1  kkt2 xtime
#> Nelder-Mead        1 TRUE FALSE  0.00
#> BFGS               1 TRUE FALSE  0.06
#> CG                 1 TRUE FALSE  0.02
#> L-BFGS-B        9999   NA    NA  0.01

optimx(par = x0, gn, method = c("BFGS", "CG"), args = args)
#>             U        Y     value fevals gevals niter convcode kkt1  kkt2 xtime
#> BFGS 191.7764 8128.266 -6026.305    100    100    NA        1 TRUE FALSE  0.04
#> CG   198.3800 9853.717 -4315.391    201    101    NA        1 TRUE FALSE  0.02

^{Created on 2022-12-23 with reprex v2.0.2}

The first run with 4 methods gives similar results for methods BFGS and CG. The second run keeps only these two methods.
The function's values are the symmetric of the values in column value.

Data

Here is the posted arguments data set as copy&paste able code.

args <- "Y_crit     U_crit  Q   p   U_max   Y
12327.9     123.2790    57.14286    0.75    198.38  11000"
args <- read.table(textConnection(args), header = TRUE)

^{Created on 2022-12-22 with reprex v2.0.2}