fzero giving incorrect result

I am trying to find the root of a rather com plicated equation.

I know that the correct answer is Hin = 2.33 becasue i have used GoalSeek in Excel and the result is correct.

However fzero gives Hin = 2.216.

Any help would be appreciated!

Ted

testfzero()

Hin = 2.1206

ans = 0

ans = -0.0443

function testfzero()
% BASIC DATA
el=0.01;      % m     length of taper
R=0.5;        % m     radius of tip
hmin=1e-6;    % m     min film thickness
e=0.005;       % m     length of tip
pin=0;        % Pa    inlet pressure
pout=0;       % Pa    outlet pessure
slope=-1.4e-4; %       slope of taper
eta=0.05;      % Pa-s  viscosity
U=5;

% temporary values to est routine  htese are normall passed over
Qin=0.65;

% CALCULATE INITIAL CONSTANTS
K=el^2/2/R/hmin;  % constant in film thickness equation in radius
%  H = 1 + K*(X-1)^2
%  H = H1 + slope*X is taper equation
const=6/sqrt(K);  %  used later
h1 = hmin-slope*el;%  inlet at taper
H1 = h1/hmin;
Pin = pin*hmin^2/eta/U/el;
Pout=pout*hmin^2/eta/U/el;
dX=0.01;
n=int64(1.5/dX)+1;      %total number of poiints
m=int64(1/dX)+1;        % number of points in taper
% INITIALIZE VARIABLES

%gamma_out=0.54;         % sets the outlet posistion which will be detremined later by required flow
gamma_out=sqrt(atan(Qin*2-1)); % because dp/dx is zero at outlet so just equal to Qin=Hout/2
gamma_out=double(gamma_out);
LHS=cos(gamma_out)^3*sin(gamma_out)/4+0.75*(gamma_out/2+sin(2*gamma_out)/4);
RHS=(sin(2*gamma_out)/4+gamma_out/2);
Hout=sec(gamma_out)^2;

Resid=@(Hin) -(Pout-(6/sqrt(K))*(sin(2*gamma_out)/4+gamma_out/2-sec(gamma_out)^2*(cos(gamma_out)^3*sin(gamma_out)/4+0.75*(gamma_out/2+sin(2*gamma_out)/4))-6*(1-Hin)/(1-H1)/Hin)*(1-sec(gamma_out)^2*(1+Hin)/2/Hin)-Pin);

Hin=fzero(Resid,H1)  % H1 is first guess
Resid(Hin)
Resid(2.33)
end

[Torsten]

Seems MATLAB is better than Excel (see above). :-)

[Stephen23]

"I know that the correct answer is Hin = 2.33 becasue i have used GoalSeek in Excel and the result is correct. However fzero gives Hin = 2.216."

Lets plot the function you defined:

% BASIC DATA
el=0.01;      % m     length of taper
R=0.5;        % m     radius of tip
hmin=1e-6;    % m     min film thickness
e=0.005;       % m     length of tip
pin=0;        % Pa    inlet pressure
pout=0;       % Pa    outlet pessure
slope=-1.4e-4; %       slope of taper
eta=0.05;      % Pa-s  viscosity
U=5;

% temporary values to est routine  htese are normall passed over
Qin=0.65;

% CALCULATE INITIAL CONSTANTS
K=el^2/2/R/hmin;  % constant in film thickness equation in radius
%  H = 1 + K*(X-1)^2
%  H = H1 + slope*X is taper equation
const=6/sqrt(K);  %  used later
h1 = hmin-slope*el;%  inlet at taper
H1 = h1/hmin;
Pin = pin*hmin^2/eta/U/el;
Pout=pout*hmin^2/eta/U/el;
dX=0.01;
n=int64(1.5/dX)+1;      %total number of poiints
m=int64(1/dX)+1;        % number of points in taper
% INITIALIZE VARIABLES

%gamma_out=0.54;         % sets the outlet posistion which will be detremined later by required flow
gamma_out=sqrt(atan(Qin*2-1)); % because dp/dx is zero at outlet so just equal to Qin=Hout/2
gamma_out=double(gamma_out);
LHS=cos(gamma_out)^3*sin(gamma_out)/4+0.75*(gamma_out/2+sin(2*gamma_out)/4);
RHS=(sin(2*gamma_out)/4+gamma_out/2);
Hout=sec(gamma_out)^2;

Resid=@(Hin) -(Pout-(6/sqrt(K))*(sin(2*gamma_out)/4+gamma_out/2-sec(gamma_out)^2*(cos(gamma_out)^3*sin(gamma_out)/4+0.75*(gamma_out/2+sin(2*gamma_out)/4))-6*(1-Hin)/(1-H1)/Hin)*(1-sec(gamma_out)^2*(1+Hin)/2/Hin)-Pin);

Hin=fzero(Resid,H1)  % H1 is first guess

Hin = 2.1206

Resid(Hin)

ans = 0

fplot(Resid,[1,3])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold on
plot(Hin,0,'r*')

So far I don't see anything to support your title "fzero giving incorrect result", because as far as I can tell, FZERO is returning the correct result for the function you defined. Whether you defined the correct function is another question entirely.

Answer 1

[Daniel Vieira]

I ran this code, got Hin = 2.216, and checked it with Resid(Hin), which gave zero. Testing Resid(2.33) does not give zero. I think you may have you missed some operation or parameter value in your function definition.

[Ted Smith]

Thank you all for helping me. In the end I got the R.H. brackets wrong in the function! Lesson to myself - delete some brackets, then insert them again. Matlab tells you where the corresponding L.H. bracket is!!