fzero giving incorrect result

Question

Ted Smith 2023-2-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1919710-fzero-giving-incorrect-result

评论： Ted Smith 2023-2-27

采纳的回答： Daniel Vieira

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Torsten 2023-2-27

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

Stephen23 2023-2-27

在 MATLAB Online 中打开

"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 2023-2-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1919710-fzero-giving-incorrect-result#answer_1180805

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Ted Smith 2023-2-27

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!!

请先登录，再进行评论。