CODEWITHSUNDEEP
rloops
Gerry
Gerry

Reputation: 117

Looping uniroot() on two arguments with sapply()

This function that calculates the pressure loss in a pipe given its diameter, flowrate, and length. 

hazwil2 <- function(diam, flow, leng){
  psi2=((1/(2.31*100))*1050*((flow/140)^1.852)*leng*diam^-4.87)
  return(psi2)
}

I'm looking for the smallest diameter that will keep the pressure loss within 2 psi. The range of diameters is between 2 and 12 inches. Using uniroot() and two arbitrary values for flow and length:

intercept   = 2L
uniroot(
  function(x) hazwil2(x, 100, 400) - intercept ,
  interval = c(2, 12)
)$root

Now I'm trying to loop uniroot on this dataset of values of flow and length

data <- read.csv(
  text = 
  "leng,flow
  100, 100
  200, 100
  300, 100
  100, 150
  200, 150
  300, 150
  100, 200
  200, 200
  300, 200",
  stringsAsFactors = FALSE
)

Upvotes: 3

Views: 488

Answers (2)

Parfait
Parfait

Reputation: 107687

Consider generalizing your unitroot call for f and l params, and pass it in mapply (multiple apply) to iterate elementwise down the equal length vectors (i.e., columns, flow and length of dataframe):

find_root <- function(f, l) {
  intercept <- 2L
  uniroot(
    function(x) hazwil2(x, f, l) - intercept ,
    interval = c(2, 12)
  )$root
}

data$root <- mapply(find_root, data$flow, data$leng)

data
#   leng flow     root
# 1  100  100 2.681145
# 2  200  100 3.091243
# 3  300  100 3.359606
# 4  100  150 3.128141
# 5  200  150 3.606602
# 6  300  150 3.919727
# 7  100  200 3.489780
# 8  200  200 4.023564
# 9  300  200 4.372916

For a larger set which can potentially yield root issues, consider wrapping function call in tryCatch to return NA:

data <- expand.grid(leng = seq(100, 1000, by=100), flow = seq(10, 200, by=10))
data$root <- mapply(function(l,f) tryCatch(find_root(l,f), error=function(e) NA), 
                    data$flow, data$leng)

Output

head(data, 20)

#    leng flow     root
# 1   100   10       NA
# 2   200   10       NA
# 3   300   10       NA
# 4   400   10       NA
# 5   500   10       NA
# 6   600   10       NA
# 7   700   10       NA
# 8   800   10       NA
# 9   900   10       NA
# 10 1000   10       NA
# 11  100   20       NA
# 12  200   20       NA
# 13  300   20       NA
# 14  400   20       NA
# 15  500   20 2.023187
# 16  600   20 2.100367
# 17  700   20 2.167915
# 18  800   20 2.228178
# 19  900   20 2.282705
# 20 1000   20 2.332653

Upvotes: 2

wibeasley
wibeasley

Reputation: 5287

This based on @Parfait's answer, and includes some graphs that came to mind during the previous question. The first graph sets the x & y coordinates to flow and leng and displays the min diameter required to drop below 2 PSI.

data$diam_min        <- purrr::map2_dbl(.x=data$flow, .y=data$leng, ~find_root(f=.x, l=.y))
data$diam_min_pretty <- sprintf("diam min:\n%0.2f", data$diam_min)

library(ggplot2);  library(magrittr)
ggplot(data, aes(x=flow, y=leng, label=diam_min_pretty)) +
  geom_text()

13-points

The second graph make each value of leng a contour line, and the min diameter replaces it on the y axis. That open dot represents the solution to Optimize function in R with constraints.

tidyr::crossing(leng=c(100, 200, 300, 400), flow=100:200) %>% 
  dplyr::mutate(
    diam_min = purrr::map2_dbl(.x=flow, .y=leng, ~find_root(.x, .y)),
    leng     = factor(leng)
  ) %>% 
  ggplot(aes(x=flow, y=diam_min, color=leng)) +
  geom_line() +
  annotate("point", x=100, y=find_root(f=100, l=400), size=4, shape=1)

level-curves

I'm still not sure how you're going to choose the most desirable value (out of this set that achieves <2 psi), but maybe these visual relationships will help.

Upvotes: 1

Related Questions