One-Log Call Iterative Solution of the Colebrook Equation for Flow Friction Based on Pade Polynomials

The 80 year-old empirical Colebrook function, widely used as an informal standard for hydraulic resistance, relates implicitly the unknown flow friction factor, with the known Reynolds number and the known relative roughness of a pipe inner surface. It is based on logarithmic law in the form that captures the unknown flow friction factor in a way that it cannot be extracted analytically. As an alternative to the explicit approximations or to the iterative procedures that require at least a few evaluations of computationally expensive logarithmic function or non-integer powers, this paper offers an accurate and computationally cheap iterative algorithm based on Pade polynomials with only one log-call in total for the whole procedure (expensive log-calls are substituted with Pade polynomials in each iteration with the exception of the first). The proposed modification is computationally less demanding compared with the standard approaches of engineering practice, but does not influence the accuracy or the number of iterations required to reach the final balanced solution.


Introduction
The empirical Colebrook equation (Colebrook 1939) implicitly relates the unknown flow friction factor with the known Reynolds number and the know relative roughness of inner pipe surface, ; , where is functional symbol, Equation (1). (1) In Equation (1) is Reynolds number; , is relative roughness of inner pipe surface; , and is Darcy flow friction factor; (all three quantities are dimensionless).
The Colebrook equation is based on experiments performed by Colebrook and White in 1937 with the flow of air through a set of artificially roughened pipes (Colebrook and White 1937). The accuracy of this eighty year old equation is disputed many times but it is still accepted in engineering practice as an informal standard for hydraulic resistance. Therefore to repeat results and for comparisons, it is required to solve the Colebrook equation accurately. Numerous evaluations of flow friction factor such as in the case of complex networks of pipes pose extensive burden for computers, so not only an accurate but also a simplified solution is required. Calculation of complex water or gas distribution networks (Brkić 2009, Brkić 2011ab, Praks et al. 2015, Praks et al. 2017 which requires few evaluations of logarithmic function for each pipe, presents a significant and extensive burden which available computer resources hardly can easily manage (Clamond 2009, Giustolisi et al. 2011, Danish et al. 2011, Winning and Coole 2013, Vatankhah 2018.
The Colebrook equation is based on logarithmic law where the unknown flow friction factor is given implicitly, i.e., it appears on both sides of Equation (1) in form , from which it cannot be extracted analytically; an exception is through the Lambert -function (Keady 1998, Sonnad and Goudar 2004, Brkić 2011cd, Brkić 2012ab, Biberg 2017, Brkić 2017a. The common way to solve it is to guess an initial value for friction factor and then to try to balance it using the iterative algorithm (Brkić 2017b) which needs to be terminated after the certain number of iterations when the final balanced value is reached. As an alternative to the iterative procedure, numerous approximate formulas are available (Gregory and Fogarasi 1985, Zigrang and Sylvester 1985, Brkić 2011e, Brkić 2012c, Brkić and Ćojbašić 2017, Pimenta et al. 2018. Usually, more complex approximations are more accurate, but also more computationally expensive because they contain at least a few logarithmic expressions and/or terms with non-integer powers which require use of demanding algorithms (namely an evaluation of fractional exponential and natural logarithm) to be evaluated in processor units of computers and to be stored in registers (Clamond 2009, Giustolisi et al. 2011, Danish et al. 2011, Winning and Coole 2013, Vatankhah 2018, Sonnad and Goudar 2004.
The presented scheme for solving Colebrook equation requires only one single call of the logarithmic function in respect to the whole iterative procedure. It is equally accurate as a standard iterative approach and does not require additional iterations to reach the same accuracy. Instead of computationally expensive logarithmic function, its Padé polynomial equivalent (Baker and Graves-Morris 1996) is used in all iterations, exception the first. The Padé approximant is the approximation of a function by a rational algebraic fraction where both the numerator and the denominator are polynomials. Because these rational functions only use the elementary arithmetic operations, they can be evaluated numerically very easily. In the computer environment, they required less basic floating-point operations compared with the logarithmic function (Kropa 1978, Rising 2007, Pineiro et al. 2004, Al-Mohy 2012.
The presented simplified iterative method can be profitable for future computing software in terms of having a high level of accuracy and speed with a decreased computational burden.

Evaluation of Logarithmic Function through Padé Polynomials
Basic floating-point operations such as addition and multiplication are carried out directly in the Central Processor Unit (CPU) while logarithmic functions, exponents or square roots require expensive operations based on more complex algorithms (Kropa 1978, Rising 2007, Pineiro et al. 2004, Al-Mohy 2012. In addition to logarithms and non-integer powers, Biberg (2017) adds also division in the group of more costly functions for evaluation; addition, subtraction and multiplication are low-cost operations according to Biberg. Evaluated with various compilers and executed on various platforms, integer addition, subtraction, or multiplication requires less than 1 floating-point operation, float addition, subtraction, or multiplication about 1, float division 2-6, integer division 4-10, square root 5-20, while functions , , , as well and functions 10-40 floating-point operations. Winning and Coole (2015) report average time for 100 million operation in seconds and relative effort, respectively as follows: addition 23.40s and 1, subtraction 27.50s and 1.18, division 31.70s and 1.35, multiplication 36.20s and 1.55, squared 51.10s and 2.18, square root 53.70s and 2.29, cubed 55.58s and 2.38, natural log 63.00s and 2.69, cubed root 63.40s and 2.71, fractional exponential 77.60s and 3.32, and log to 10-base 78.80s and 3.37.
Regarding the Colebrook equation, in order to simplify the iterative procedure which is in common use in engineering practice, the logarithmic function is replaced with its relevant Padé polynomial equivalent in all iterations, with exception to the first. Padé polynomials can accurately approximate logarithmic function only in a limited domain. For example, if it is known that , calculation of for example can be evaluated using the fact that where is close to 1. Logarithmic function can be replaced by its Padé polynomial equivalent very accurately in a limited domain, and therefore instead of , already calculated . Evaluation of 10-base logarithmic function in many computing languages goes through natural logarithm where and where is constant, and therefore the Padé polynomials that approximate accurately for are shown; Equations (2-7). Padé polynomials of different orders can be used for approximation of here all accurate for arguments close to 1; . As the expansion point is a root of the accuracy of the Padé approximant decreases. Setting the OrderMode option in Matlab padé command to relative compensates for the decreased accuracy. Thus, here, the Pade approximant of order uses the form , where α and β are coefficients (the coefficients of the polynomials need not be rational numbers). For example, Padé polynomial of order (2, 3) is with polynomial of order 2 in numerator and of order 3 in denominator; Equation (6). Of course, low order formulas are simpler, but they have larger errors than high order formulas and vice versa. As can be seen from Figure 1, even very simple form of Padé polynomials (1,2) and (2,1) are of high accuracy in respect of domain of interest for solving the Colebrook equation which is . Horner algorithm transforms polynomials into a computationally efficient form and therefore, Horner nested polynomial representations of the Padé polynomials of different orders for where are shown here; Equations (2-7). Higher order of Padé approximants are more accurate, but more complex.
Relative error introduced by them; Equations (2-5) compared with is shown in Figure 1 and for Equation (6) in Table 1. The relative error of Padé approximants (2,2) for of is negligible for . Thus, relative error of the used Padé approximants (2,3) of in the proposed iterative procedure is even more negligible and therefore it is not presented in Figure  1, but is available in Table 1.  To illustrate the complexity of computing in modern computers it should be noted that even such a relatively simple equation such as Colebrook's can make a numerical problem in computer registers due to overflow error. Its transformed version in term of the Lambert -function can give such large numbers for some pairs of the Reynolds number Re and the relative roughness of inner pipe surface which are from the practical domain of applicability in engineering practice which cannot be stored in 32-or 64-bit registers of modern computers (Sonnad andGoudar 2004, Brkić 2012a).

Initial Starting Point for the Proposed Iterative Method
In case of the Colebrook equation, practical experience shows that trying to get a good initial starting point has limited value until it is chosen in the domain of applicability of the equation which is . Every initial starting point chosen from the domain of applicability of the Colebrook equation will lead to the final accurate solution surely, with the only difference that in some cases more additional iterations would be needed. Usually, with the initial guess that is close to the exact solution, the iterative procedure converges to it in five or fewer iterations. To date, cases which lead to divergence, fluctuation, or convergence to a  (Brkić 2017). Here the goal is to avoid use of logarithmic functions and therefore, this starting point is not suitable.
2. An efficient procedure for finding a sufficiently good initial starting point is proposed by Yun (2008) in the integral form; Equation (8): In Equation (8), , represents the Colebrook equation, is the lower while is the upper limit from which an initial starting point should be chosen; = 3.68 and = 12.47 because the domain of applicability of the Colebrook equation that is between 3.68 and 12.47 in respect to , is signum function: if , , and , while is hyperbolic tangent which is defined through the exponential function with non-integer power the use of which is as computationally expensive as the use of the logarithmic function and which therefore cannot be recommended. 3. Every explicit approximation of the Colebrook equation (Gregory and Fogarasi 1985, Zigrang and Sylvester 1985, Brkić 2011, Brkić and Ćojbašić 2017, Pimenta et al. 2018); , where is the functional symbol, can be used to choose an initial starting point . On the other hand, almost all available approximations contain logarithmic or/and terms with non-integer powers, which makes them unsuitable for use in the developed approach. On the other hand, having previous experience with training Artificial Neural Networks (ANN) to simulate the Colebrook equation (Özger and Yıldırım 2009, Brkić and Ćojbašić 2016, Bardestani et al. 2017, i.e. to use ability of artificial intelligence to simulate the Colebrook equation not knowing its logarithmic nature but only knowing raw input and corresponding output datasets , a computationally cheap explicit approximation of the Colebrook equation is developed through genetic programming (Giustolisi and Savić 2006, Ćojbašić and Brkić 2013, Brkić 2014, Brkić and Ćojbašić 2016. The developed approximation is computationally efficient because of its polynomial structure; Equation (9): (9) Eureqa [computer software] by Nutonian, Inc., Boston, MA. (Schmidt andLipson 2009, Dubčáková 2011) is used to generate Equation (9). The Eureqa-polynomial approximation; Equation (9) has up to 40% relative error, but it is very cheap and sufficiently accurate to serve as an initial starting point . 4. Extensive tests over the domain of applicability of the Colebrook equation shows that one fixed value also can be used as the initial starting point for the iterative procedure in all cases. Results indicate that the proposed Padé approach works in all cases, as the argument z of ln(z) is always close to one. When Equation (9)  To avoid using a computationally expensive logarithmic function in the initial stage of the iterative procedure, the recommendation is to start calculation with fixed-value starting point or to use a polynomial expression; Equation (9). Power-law formulas from Russian practice which do not contain non-integer powers also can be used (Альтшуль 1982).

Proposed Iterative Method
The Colebrook equation is usually solved iteratively using the Newton-Raphson method (Ypma 1995, Abbasbandy 2003 or even more using a simplified Newton-Raphson method known as the fixed-point method (Brkić 2017b). Recently, hybrid three-point methods have been proposed (Brkić and Praks 201x).
Here is presented an adjusted very accurate, fast and computationally cheap version of the Newton-Raphson method suitable for the Colebrook equation in which the logarithmic function is replaced after the first iteration with the Padé approximant in polynomial form (Baker and Graves-Morris 1996).
Knowing that the Colebrook equation is based on logarithmic law White 1937, Colebrook 1939), the achievement with this simplified approach is more significant. Numerical examples are shown in Section 5 of this paper.

Iteration 0:
In order to avoid use of computationally expensive logarithmic functions or functions with non-integer powers, a required initial starting point should be chosen using recommendations from Section 3 of this paper; points 3 or 4.

Iteration 1:
Having provided an initial starting point , new value can be calculated using Equation In Equation (11), which will be used also in the next iteration (in an additional variant of the proposed method is used in all subsequent iterations), while in Equation (10), the first derivative of in respect to ; is from Equation (12): In Equation (12), is with constant value, and therefore only from Equation (11) requires evaluation of the logarithmic function.
In many programming languages evaluation of logarithmic function of any base is processed by natural logarithm (Vatankhah 2018). Change of 10-base logarithm from the Colebrook equation to e-based natural logarithm where and where is implemented as .

Iteration 2:
New value should be calculated using Equation (13): In Equation (13), is not calculated by , where , but using Padé polynomial replacement for logarithmic function which is accurate for and using the already calculated value of from the previous iteration; Equation (14): In Equation (14), , , and . In the first iteration, is already known; Equation (11). The Padé polynomial used in Equation (14) is of order (2,3) which means that the polynomial in the numerator is of the order of 2 while in the denominator of order 3. The Padé polynomials are also known as Padé approximants and here the maximal relative error of the polynomial expression term in Equation (14) in domain ; is minor as shown in Table 1. Value of z for the procedure shown in practice is and therefore the error of the used Padé approximant can be neglected in the case shown.
The first derivative does not contain any logarithmic functions and should be evaluated using Equation (12), where should be replaced with the new value or knowing that the value of the derivative does not change significantly between two iterations, can be reused in all subsequent iterative cycles. Even knowing that the value of the first derivate in the procedure shown is always near 1; for rough calculations it can be assumed that which gives the fixed-point method as a special case of the Newton-Raphson scheme.

Iteration 3:
New value is again evaluated in the same way using Equation (15): In Equation (15) with the reference to the preceding iteration (here to the second iteration); Equation (17).
The Padé polynomial is a very accurate approximation of logarithmic function, so knowing that is evaluated directly through the logarithmic function, while , , , etc. is based on its Padé polynomial equivalent, it is obvious that the sequence , , , etc. is slightly more accurate compared with the sequence , , , etc. which accumulates error introduced with Padé approximations. On the other hand, the error is minimized when the argument is closer to 1 which is case for the second sequence , , , etc. In both cases, the introduced error can be neglected.

Iteration i:
All indexes in respect the third iteration should be updated as with exemption of index 0 in Equation (16). The calculation is finished when . The algorithm for the proposed improved procedure is given in Figure 2. Only a one-off evaluation of the logarithmic function is needed in the proposed algorithm from Figure 2, which is clearly marked in red; . On the other hand, calculated in iteration 1 is reused in all next steps and it is marked in green in Figure 2.
The proposed procedure can be simplified assuming that which gives the simple fixed-point procedure (Brkić 2017b) instead of the Newton-Raphson.

Numerical examples
Here

Conclusions
An efficient algorithm for the iterative calculation of the Colebrook equation by both an accurate and computationally efficient Padé approximation is presented in this paper. It requires only one evaluation of the logarithmic function in respect to the whole iterative procedure and more specifically only in the first iteration, while the common procedures from current engineering practice require at least one evaluation of logarithmic function for every single iteration. The logarithmic function in the proposed procedure is replaced in all iterations (except the first), with simple, accurate and efficient Padé polynomials (Baker and Graves-Morris 1996). In this way the same accuracy is reached through the proposed less demanding procedure, after the same number of iterations as in the standard algorithm which uses -call in each iterative step. This is a good achievement knowing that the nature of the Colebrook equation is logarithmic. For their evaluation in the Central Processor Unit (CPU) of computers, Padé polynomials require a lower number of floating-point operations to be executed compared with the logarithmic function (Clamond 2009, Giustolisi et al. 2011, Danish et al. 2011, Winning and Coole 2013, Vatankhah 2018, Sonnad and Goudar 2004, Brkić 2012a, Winning and Coole 2015. The here presented iterative approach only introduces a computationally cheaper alternative to the standard iterative procedure. It does not reduce the number of required iterations to reach the final desired accuracy nor provide more accurate results. The proposed method reduces the burden for the Central Processing Unit (CPU) as less floating-point operations need to be executed. In that way, the presented approach also increases speed of computation. On the other hand, many explicit non-iterative approximations to the Colebrook equation are available in literature (Gregory and Fogarasi 1985, Zigrang and Sylvester 1985, Brkić 2011e, Pimenta et al. 2018) which initially appear simple for computation, but are not. They are widely used, but although some of them are very accurate, they contain relatively complex internal iterative steps and also a number of computationally demanding functions. For example, the widely used Haaland approximation introduces relative error up to 1.5% (Haaland 1983, Wood andHaaland 1983), but with the cost of evaluation of one logarithmic expression and one non-integer power. Also, the approximation by Romeo et al. (2002) reaches extremely high accuracy with the relative error of up to 0.14%, but with a cost of evaluation of even three logarithmic expressions and two non-integer powers. Regarding alternative iterative procedures, Clamond (2009) provides an accurate iterative approach using function, but this algorithm requires at least two -calls; one for initialization and one in the first iteration, which is more expensive compared with the here presented approach.
The procedure proposed in this paper can significantly reduce the computational burden for evaluating complex distribution networks with various applications (water, gas) (Brkić 2009, Brkić 2011a, Praks et al. 2015, Praks et al. 2017. For example, a probabilistic approach using time dependent modeling of distribution or transmission networks requires many millions of evaluations of Colebrook's equation, which means that it is not a computationally cheap task at all. For such kinds of computations is always good to have a cheaper but still accurate approach in order to speed up the process.