# fzero giving incorrect result

5 views (last 30 days)
Ted Smith on 27 Feb 2023
Commented: Ted Smith on 27 Feb 2023
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
##### 2 CommentsShow 1 older commentHide 1 older comment
Stephen23 on 27 Feb 2023
"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.

Daniel Vieira on 27 Feb 2023
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 on 27 Feb 2023
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!!

### Categories

Find more on Programming in Help Center and File Exchange

R2020b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!