My bisection method is freezing when I put smaller error margin.

When I put error margin smaller than 0.1, it keeps going but never stops.

clc; clear all; close all;

D= 4*0.0254;     %inner diameter (inch to meter)
radius=D/2;
Q= (2000 * 42 * 3.785e-3) / 86400 ;     %Volumetric flow rate (bbl to gal to m^3) / (day to second)
area= pi*((radius)^2);   %cross sectional area 
p= 0.9 * 1000;    %density (g/cm^3 to kg/m^3)
Mu= 8 * 0.001    %viscosity (cp to pa*s)
e_D=[0,0.002,0.004,0.006,0.008];

Re=linspace(5000,100000,5);  %Reynolds number


z=@(x)(-2.0).*log10(((e_D)./3.7)+(2.51./(Re.*sqrt(x))))-(1./sqrt(x)); %friction factor
%bisection method
error=1e-8;
a=0;
b=1;
c=(a+b)./2;
while abs(z(c))>error
    if z(a).*z(c)<0
        b=c;
    else
        a=c;
    end
    c=(a+b)./2;
end

Answer 1

[Chunru]

Your function z does not return a scaler output for a scalar input x since e_D and Re are vectors. for bisection method to work, z(x) is a scalar function.

D= 4*0.0254;     %inner diameter (inch to meter)
radius=D/2;
Q= (2000 * 42 * 3.785e-3) / 86400 ;     %Volumetric flow rate (bbl to gal to m^3) / (day to second)
area= pi*((radius)^2);   %cross sectional area 
p= 0.9 * 1000;    %density (g/cm^3 to kg/m^3)
Mu= 8 * 0.001    %viscosity (cp to pa*s)

Mu = 0.0080

e_D=[0,0.002,0.004,0.006,0.008];
e_D = 0;
Re=linspace(5000,100000,5);  %Reynolds number
Re = 5000;
z=@(x)(-2.0).*log10(((e_D)./3.7)+(2.51./(Re.*sqrt(x))))-(1./sqrt(x)); %friction factor
z(1)

ans = 5.5986

%bisection method
error=1e-8;
a=0;
b=1;
c=(a+b)./2;
i= 0;
while abs(z(c))>error
    if z(a).*z(c)<0
        b=c;
    else
        a=c;
    end
    
    fprintf("a=%f  b=%f  c=%f   z(c)=%f\n", a, b, c, z(c));
    
    %if i>50, break; end
    c=(a+b)./2;
    i = i+1;
end

a=0.000000  b=0.500000  c=0.500000   z(c)=4.883349
a=0.000000  b=0.250000  c=0.250000   z(c)=3.996533
a=0.000000  b=0.125000  c=0.125000   z(c)=2.867075
a=0.000000  b=0.062500  c=0.062500   z(c)=1.394473
a=0.031250  b=0.062500  c=0.031250   z(c)=-0.563412
a=0.031250  b=0.046875  c=0.046875   z(c)=0.650732
a=0.031250  b=0.039062  c=0.039062   z(c)=0.130708
a=0.035156  b=0.039062  c=0.035156   z(c)=-0.188738
a=0.037109  b=0.039062  c=0.037109   z(c)=-0.023009
a=0.037109  b=0.038086  c=0.038086   z(c)=0.055256
a=0.037109  b=0.037598  c=0.037598   z(c)=0.016486
a=0.037354  b=0.037598  c=0.037354   z(c)=-0.003169
a=0.037354  b=0.037476  c=0.037476   z(c)=0.006681
a=0.037354  b=0.037415  c=0.037415   z(c)=0.001762
a=0.037384  b=0.037415  c=0.037384   z(c)=-0.000702
a=0.037384  b=0.037399  c=0.037399   z(c)=0.000530
a=0.037392  b=0.037399  c=0.037392   z(c)=-0.000086
a=0.037392  b=0.037395  c=0.037395   z(c)=0.000222
a=0.037392  b=0.037394  c=0.037394   z(c)=0.000068
a=0.037393  b=0.037394  c=0.037393   z(c)=-0.000009
a=0.037393  b=0.037393  c=0.037393   z(c)=0.000030
a=0.037393  b=0.037393  c=0.037393   z(c)=0.000010
a=0.037393  b=0.037393  c=0.037393   z(c)=0.000001
a=0.037393  b=0.037393  c=0.037393   z(c)=-0.000004
a=0.037393  b=0.037393  c=0.037393   z(c)=-0.000002
a=0.037393  b=0.037393  c=0.037393   z(c)=-0.000001
a=0.037393  b=0.037393  c=0.037393   z(c)=0.000000
a=0.037393  b=0.037393  c=0.037393   z(c)=-0.000000
a=0.037393  b=0.037393  c=0.037393   z(c)=-0.000000
a=0.037393  b=0.037393  c=0.037393   z(c)=-0.000000

fplot(z, [-1 1]); grid on