Fzero for Colebrook Equation

I have looked here: https://uk.mathworks.com/help/matlab/ref/fzero.html

I am having no success in using fzero for

f = -1/sqrt(x)-2.*log10(((e./D)./3.7)+2.51/(Re*sqrt(x))) 

between

[0 0.1]

though.

Can someone please help with this? I want to know the iterations, root, time elapsed, and error.

I would ideally like to use the version that uses a function file, f.m.

[Tony Rankin]

So I have this

function y = f(x)
y = -1/sqrt(x)-2.*log10(((e./D)./3.7)+2.51/(Re*sqrt(x)));
end

which when I use their example works.

But when I call it using

fun = @f; 
p = 0.9*1000; 
mu = 8*0.001; 
e = 0.01; 
D = 4*0.0254;
Q = (2000*42*3.785*10^-3)/(24*60*60);
A = (pi*(D.^2))/4; 
Re = (p*D*Q)/(A*mu); 
x0 = 0.01; 
z = fzero(fun,x0)

I receive an error stating that 'e' is unrecognised.

[Rik]

http://web.archive.org/web/20210316072011/https://www.mathworks.com/matlabcentral/answers/772553-fzero-for-colebrook-equation

Answer 1

[John D'Errico]

What you do not understand is that function have their own workspaces. If you create thse variables outside the function, then call it without passing in these variables, the function will fail, because it cannot "see" those variables. They don't exist in the function workspace. (Nested functions can see variables above them.)

If you insist on using an m-file function, then this should work:

p = 0.9*1000; 
mu = 8*0.001; 
e = 0.01; 
D = 4*0.0254;
Q = (2000*42*3.785*10^-3)/(24*60*60);
A = (pi*(D.^2))/4; 
Re = (p*D*Q)/(A*mu); 
x0 = 0.01; 
fplot(@(x) myfun(x,e,D,Re),[0,10])

I've plotted the function to see if a root does exist. Only now do I want to use fzero.

Also, note the use of fzero where I provide a bracket around the root. When you give it only ONE value, then it must do a search. And negative values for x will cause problems with your function.

options = optimset('Display','iter');
[xval,fval,exitflag] = fzero(@(x) myfun(x,e,D,Re),[eps,10],options)

 
 Func-count    x          f(x)             Procedure
    2              10       2.82898        initial
    3              10       2.82898        interpolation
    4               5       2.69594        bisection
    5             2.5        2.5078        bisection
    6            1.25       2.24175        bisection
    7           0.625       1.86552        bisection
    8          0.3125       1.33353        bisection
    9         0.15625      0.581296        bisection
   10        0.078125     -0.482273        bisection
   11        0.113551      0.136933        interpolation
   12        0.105717     0.0273238        interpolation
   13        0.103868  -0.000326184        interpolation
   14         0.10389   4.29348e-06        interpolation
   15         0.10389   6.66853e-10        interpolation
   16         0.10389  -4.44089e-16        interpolation
   17         0.10389  -4.44089e-16        interpolation
 
Zero found in the interval [2.22045e-16, 10]

xval = 0.1039

fval = -4.4409e-16

exitflag = 1

Save this function as an m-file. Note my use of ./ operator. It seems you know about .* but the ./ operator has escaped your attention.

function y = myfun(x,e,D,Re)
  y = -1./sqrt(x)-2.*log10(((e./D)./3.7)+2.51./(Re*sqrt(x)));
end

As far as being told the time elapsed, you can use tic and toc, I suppose. fzero does not return that. I have no clue what you mean by the "error".

[Tony Rankin]

I have used tic and toc and it is now providing me with the elapsed time for each iteration. When I said error, I was referring to the 'fval' that is computed.

For this I didn't need to paste the parameters into the command window as I had been doing with bisection and Newton-Raphson.

The following is the bisection method for the same:

function [R] = bisection(f,a,b,tol)
% R is an estimation of the root of f using the bisection method
% f is Colebrook equation for turbulent flow
% a is the lower limit of friction factor
% b is the upper limit of friction factor
% tol is the accepted tolerance

if sign(f(a)) == sign(f(b))
    error('Scalars a and b do not bound a root.')
end
m = (a+b)/2;                      % calculate the midpoint of a and b 
if abs(f(m)) < tol
    R = m;                        % end loop, takes m as root of f
elseif sign(f(a)) == sign(f(m))
    R = bisection(f,m,b,tol);     % make a recursive call
elseif sign(f(b)) == sign(f(m))
    R = bisection(f,a,m,tol);     % make a recursive call
end
error = abs(f(m));                % calculate the absolute error
end

When I have been using this function I have had to put the following parameters

error = 1; 
tol = 1e-4; 
a = 0; 
b = 0.1; 
p = 0.9*1000; 
mu = 8*0.001; 
e = 0; 
D = 4*0.0254; 
Q = (2000*42*3.785*10^-3)/(24*60*60);
A = (pi*(D.^2))/4; 
Re = (p*D*Q)/(A*mu); 
f = @(x) -1/sqrt(x)-2.*log10(((e./D)./3.7)+2.51/(Re*sqrt(x)));

is the command window to make it run.

Is there an approach that avoids this like I can do with the fzero function above?

[John D'Errico]

I fail to see the problem. In the call to fzero that I wrote, I passed in the parameters like this:

[xval,fval,exitflag] = fzero(@(x) myfun(x,e,D,Re),[eps,10],options)

You can do identically the same thing, as:

xbis = bisection(@(x) myfun(x,e,D,Re),eps,10,tol)

In either case, I created a function handle on the fly in the call to the optimizer, using that to pass in the parameters as necessary.

Someplace you will need to define those parameters, thus e,D,Re, and that will need to be in the command window or in a script. So wherein lies the confusion?

[Tony Rankin]

I am new to MATLAB. I wanted to put the parameters above the function and then use the one line to call it in the command window. I think the problem is that the file is bisection.m and so the function needs to be the first line in the file. So it doesn't work.

If I change the name and put those parameters above the function, will it work?

[John D'Errico]

I'd suggest you get used to the idea that functions are best kept separate from your scripts. This is far better programming practice. Now you can use a general function, to be used to solve many different problems, not just one. Don't encapsulate your code with your data. That is just poor programming practice.