R help - Sep 2010 - Programming: loop versus vector oriented

Dear all,

I am new to R and to it's programming philosophy. The following function
is supposed to work on a vector, but I can't figure out how to do that
without looping through every element of it. Is there a more elegant
way?

Note: I have shortened it, so it is NOT correct from the pipe hydraulics
point of view

# Calculate wall friction factor
# reynolds: Reynolds number
# dk: relative roughness of pipe
lambda_wall <- function (reynolds, dk) {
  z <- 1
  result <- 0

  for (Re in reynolds) {
    if (Re <= 2320) {
      # Laminar flow
      lambda <- 64/Re
    } else if (Re < 65 * dk[z]) {
      # Turbulent flow
      if (Re < 1e+5) {
	lambda <- 0.3164 / sqrt(sqrt(Re))
      } else {
	lambda <- 0.309/(log10(Re/7))^2
      }
    } else {
      # Intermediate area
      lambdanew <- 1 / (2 * log10(3.71 * dk[z]))^2 # Start value
      iter <- 0
    
      repeat {
	lambda <- lambdanew
	lambdanew <- 1 / (2 * log10(2.51/(Re * sqrt(lambda)) + 0.27/dk[z]))^2
	iter <- iter + 1
	if ((abs(lambdanew - lambda) < 0.001) || (iter > 100)) break
      }
    
      lambda = lambdanew
    }
    
    result[z] <- lambda
    z <- z + 1
  }

  result
} # lambda_wall()

Hi

I do not want to go too much deep to internals of your function. What do 
you suppose to get as the result.

If you want to get results of your function for a vector of reynolds and 
dk you can use function outer and probably get rid of for cycle in the 
function.

outer(c(100, 530,2410), c(10, 150,200),lambda_wall)
          [,1]       [,2]       [,3]
[1,] 0.6400000 0.64000000 0.64000000
[2,] 0.1207547 0.12075472 0.12075472
[3,] 0.1081338 0.04515774 0.04515774


If you want to feed vector(s) to your function it would not be so easy as 
you need to check some predefined conditions for laminar or turbulent or 
any other flow.

In this case you cold check function

?switch

together with separate function definitions for each type of flow. But 
three if's could be as good as switch.

Basically your function probably works for vectors but only in certain 
conditions.

Laminar OK
> lambda_wall(10:20)
 [1] 6.400000 5.818182 5.333333 4.923077 4.571429 4.266667 4.000000 
3.764706
 [9] 3.555556 3.368421 3.200000

Turbulent needs dk 
> lambda_wall(2400:2410)
Error: argument "dk" is missing, with no default

But not any dk
> lambda_wall(2400:2410, 10)
Error in if (Re < 65 * dk[z]) { : missing value where TRUE/FALSE needed

only a vector of the same length as reynolds
> lambda_wall(2400:2410, 10:21)
 [1] 0.10815993 0.10325278 0.09912072 0.09558754 0.09252750 0.08984834
 [7] 0.08748069 0.08537137 0.08347884 0.08177018 0.08021892

So before trying to elaborate your function further what shall be inputs 
and what is desired output?

Regards
Petr


r-help-bounces at r-project.org napsal dne 15.09.2010 11:57:28:
> Dear all,
> 
> I am new to R and to it's programming philosophy. The following
function
> is supposed to work on a vector, but I can't figure out how to do that
> without looping through every element of it. Is there a more elegant
> way?
> 
> Note: I have shortened it, so it is NOT correct from the pipe hydraulics
> point of view
> 
> # Calculate wall friction factor
> # reynolds: Reynolds number
> # dk: relative roughness of pipe
> lambda_wall <- function (reynolds, dk) {
>   z <- 1
>   result <- 0
> 
>   for (Re in reynolds) {
>     if (Re <= 2320) {
>       # Laminar flow
>       lambda <- 64/Re
>     } else if (Re < 65 * dk[z]) {
>       # Turbulent flow
>       if (Re < 1e+5) {
>    lambda <- 0.3164 / sqrt(sqrt(Re))
>       } else {
>    lambda <- 0.309/(log10(Re/7))^2
>       }
>     } else {
>       # Intermediate area
>       lambdanew <- 1 / (2 * log10(3.71 * dk[z]))^2 # Start value
>       iter <- 0
> 
>       repeat {
>    lambda <- lambdanew
>    lambdanew <- 1 / (2 * log10(2.51/(Re * sqrt(lambda)) + 0.27/dk[z]))^2
>    iter <- iter + 1
>    if ((abs(lambdanew - lambda) < 0.001) || (iter > 100)) break
>       }
> 
>       lambda = lambdanew
>     }
> 
>     result[z] <- lambda
>     z <- z + 1
>   }
> 
>   result
> } # lambda_wall()
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

Hi Jan

Jan private <jrheinlaender at gmx.de> napsal dne 18.09.2010 12:12:29:
> Hello Petr,
> 
> thank you for your ideas. The split() looks most realistic. 
> 
> What about this idea:
> 
> 1. Define three functions Refun1, Refun2, Refun3 for the three different
> sections of the calculations (same as you suggested)
> 2. lambda = (Re <= 2320) * Refun1(Re)  + ((Re > 2320) && (Re
< 65 * dk))
> * Refun2(Re) etc.
> 
> But my thought is that probably the values of RefunXYZ will be
> calculated for every value of Re, even if the condition (Re <= 2320) is
> FALSE (= 0). So that would give a lot of unnecessary function
> evaluations.
Probably yes but it could be only small overhead and you can check speed 
by

system.time(result1 <- your function1)
system.time(result2 <- your function2)

all.equal(result1, result2)

If you do not perceive differences you can wrap function into some loop 
and let it perform several times inside system.time call

Regards
Petr
> 
> Regards,
>    Jan
> 
>