Finding the root of an equation of 1.75th order

Question

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:

Y_crit	U_crit	Q	p	U_max	Y
12327.9	123.2790	57.14286	0.75	198.38	14170

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'"

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.

Answer 1

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}