Problem with Function Handle, I want this to spit out four f values. I'm new to Mathworks too.
2 views (last 30 days)
Show older comments
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
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 ?
Answers (1)
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
0 Comments
See Also
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!