Azzera filtri
Azzera filtri

fzero giving incorrect result

7 visualizzazioni (ultimi 30 giorni)
Ted Smith
Ted Smith il 27 Feb 2023
Commentato: Ted Smith il 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 Commenti
Torsten
Torsten il 27 Feb 2023
Seems MATLAB is better than Excel (see above). :-)
Stephen23
Stephen23 il 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.

Accedi per commentare.

Risposta accettata

Daniel Vieira
Daniel Vieira il 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.
  1 Commento
Ted Smith
Ted Smith il 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!!

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su MATLAB in Help Center e File Exchange

Tag

Prodotti


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by