Problem with Function Handle, I want this to spit out four f values. I'm new to Mathworks too.

2 views (last 30 days)
D = [0.1, 0.2, 0.05, 0.15];
L = [10, 3, 4, 8];
A = 0.25*pi.*D.^2;
e = 0.0015e-3;
rho = 998;
nu = 1.01e-6;
eoverD = e./D;
alpha = 0.1;
Vo = [0.1 0.1 0.1 0.1];
Re = (D.*Vo) / nu;
if Re < 2300
f = 64 / Re;
else
Formula = @(x) 1/sqrt(x)+2*log10(eoverD/3.7 + 2.51/Re/sqrt(x));
f = fzero(Formula,0.01);
end
disp(f);
  2 Comments
Walter Roberson
Walter Roberson on 29 Aug 2013
Do you mean you want one vector of length 4 such that Formula(x) is [0, 0, 0, 0] ? Or do you mean that you want 4 distinct roots, each one a scalar? If you want 4 distinct roots, do they need to be in any particular order, or within any particular interval ?
Jay
Jay on 29 Aug 2013
Hi Walter, What I did was create for IF statements for each Re and Vo value. Presented below. Thanks for your help.
D = [0.1, 0.2, 0.05, 0.15];
L = [10, 3, 4, 8];
A = 0.25*pi.*D.^2;
e = 0.0015e-3;
rho = 998;
nu = 1.01e-6;
eoverD = e./D;
alpha = 0.1;
Vo = [2 2 2 2];
Re = (D.*Vo) / nu;
Re1 = Re(1);
Re2 = Re(2);
Re3 = Re(3);
Re4 = Re(4);
if Re1 < 2300
f1 = 64 / Re1;
else
Formula = @(x) 1/sqrt(x)+2*log10((e/D(1))/3.7 + (2.51/(Re1*sqrt(x))));
f1 = fzero(Formula,0.01);
end
if Re2 < 2300
f2 = 64 / Re2;
else
Formula = @(x) 1/sqrt(x)+2*log10((e/D(2))/3.7 + (2.51/(Re2*sqrt(x))));
f2 = fzero(Formula,0.01);
end
if Re3 < 2300
f3 = 64 / Re3;
else
Formula = @(x) 1/sqrt(x)+2*log10((e/D(3))/3.7 + (2.51/(Re3*sqrt(x))));
f3 = fzero(Formula,0.01);
end
if Re4 < 2300
f4 = 64 / Re4;
else
Formula = @(x) 1/sqrt(x)+2*log10((e/D(4))/3.7 + 2.51/(Re4*sqrt(x)));
f4 = fzero(Formula,0.01);
end

Sign in to comment.

Answers (1)

Walter Roberson
Walter Roberson on 29 Aug 2013
Replace from the first "if" on down to the end with:
Formula = @(x,k) 1/sqrt(x)+2*log10((e/D(k))/3.7 + (2.51/(Re(k)*sqrt(x))));
for K = 1 : 4
if Re(K) < 2300
f(K) = 64 / Re(K);
else
f(K) = fzero( @(x) Formula(x,K), 0.01);
end
end
Note, though, that the formula has a closed form solution: x is
(86248369/4) * log(10)^2 * D(k)^2 / (-9287*LambertW((50/251)*log(10)*Re(k)* exp((500/9287) * log(10) * e * Re(k) / D(k))) * D(k) + 500 * log(10) * e * Re(k))^2
where LambertW here represents the symbolic toolbox lambertw function (which Mathworks does not supply a numeric-only routine for)
You might note a bunch of rational numbers there. I converted the floating point numbers you had into their exact rational equivalents, as algebraic reasoning about floating point numbers and error analysis is frustrating. I need to presume that where you had 2.51 you meant exactly 251/100 -- because if you meant the 2.51 to be the rounded approximation of the real number then the solution could potentially be quite different.
You could vectorize this to do all of them at once, with post-correction for the Re(:) < 2300 situations. I include sym() here to force the lambertw function to be invoked on symbols instead of on double (as there is no lambertw class support for double.)
(86248369/4) .* log(10)^2 .* D.^2 ./ (-9287 .* lambertw( sym((50/251) .* log(10) .* Re .* exp((500/9287) .* log(10) .* e .* Re ./ D) )) .* D + 500 .* log(10) .* e .* Re).^2

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!