solving a non linear equation

9 次查看（过去 30 天）
Robert
Robert 2014-9-6
评论： Walter Roberson 2022-11-18
I am trying to solve this equation (Colebrook). I have created a problem to solve it directly. However, that is not the method I can use for this exercise. This is program I solved it initially.
function [f]= moody(ed,Re)
if Re<0
error(sprintf('Reynolds number = %f cannot be negative',Re));
elseif Re<2000
f = 64/Re; return % laminar flow
end
if ed>0.05
warning(sprintf('epsilon/diameter ratio = %f is not on Moody chart',ed));
end
findf = inline('1.0/sqrt(f) + 2.0*log10( ed/3.7 + 2.51/( Re*sqrt(f)) )','f','ed','Re');
fi = 1/(1.8*log10(6.9/Re + (ed/3.7)^1.11))^2; % initial guess at f
I need to use this bisection method to solve the Colebrook equation. which this is the bisection method.
function bisection(f,a,b,tolerance,maxiter)
%bisection method for finding root of f(x)=0
%calculate function at the initial interval and mid-point
a0=a; b0=b;%initial bracket
iter=0;
error=abs(b0-a0);
f_a0=feval(f,a0);
f_b0=feval(f,b0);
if f_a0*f_b0>0
disp ('method does not work, change the initial bracket ')
return
else
%repeat the calculation until error<tolerance
while error>tolerance
f_a=feval(f,a);
f_b=feval(f,b);
c=(a+b)*0.5;
f_c=feval(f,c);
if f_c*f_a<0
b=c;
elseif f_b*f_c<0
a=c;
end
error=abs(b-a);
iter=iter+1;
if iter>maxiter
disp('no solution within maxiter. increase maxiter')
break
end
end
end
%
if iter<=maxiter
disp(['iter= ' num2str(iter)])
disp(['root is ' num2str(c)])
end
fplot(f,[a0,b0])
xlabel('x'),ylabel('f(x)');grid on
dfTol = 5e-6;
f = fzero(findf,fi,optimset('TolX',dfTol,'Display','off'),ed,Re);
Here is where I am at now.
function [f]=friction_factor(f)
Re=10e4;
ed=0.03;
f=(1.0/((sqrt(f)+2.0*log10((ed/3.7)+(2.51/(Re*sqrt(f)))))))^2;
end
I think there is something small that I am not seeing. I just went blank after writing the much. I would appreciate any help to get to the next steps.
thank you
  1 个评论
Roger Stafford
Roger Stafford 2014-9-7
It seems a shame that you should be messing around with iterative methods when you can obtain an explicit solution using the 'lambertw' function of matlab.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Brand
Brand 2022-11-18
编辑：Walter Roberson 2022-11-18
% Bisection Method in MATLAB
a=input('Enter function with right hand side zero:','s');
f=inline(a);
xl=input('Enter the first value of guess interval:') ;
xu=input('Enter the end value of guess interval:');
tol=input('Enter the allowed error:');
if f(xu)*f(xl)<0
else
fprintf('The guess is incorrect! Enter new guesses\n');
xl=input('Enter the first value of guess interval:\n') ;
xu=input('Enter the end value of guess interval:\n');
end
fori=2:1000
xr=(xu+xl)/2;
if f(xu)*f(xr)<0
xl=xr;
else
xu=xr;
end
if f(xl)*f(xr)<0
xu=xr;
else
xl=xr;
end
xnew(1)=0;
xnew(i)=xr;
if abs((xnew(i)-xnew(i-1))/xnew(i))<tol,break,end
end
str = ['The required root of the equation is: ', num2str(xr), '']
Brand
Brand 2022-11-18
编辑：Walter Roberson 2022-11-18
% Bisection Method in MATLAB
a=input('Enter function with right hand side zero:','s');
f=inline(a);
xl=input('Enter the first value of guess interval:') ;
xu=input('Enter the end value of guess interval:');
tol=input('Enter the allowed error:');
if f(xu)*f(xl)<0
else
fprintf('The guess is incorrect! Enter new guesses\n');
xl=input('Enter the first value of guess interval:\n') ;
xu=input('Enter the end value of guess interval:\n');
end
fori=2:1000
xr=(xu+xl)/2;
if f(xu)*f(xr)<0
xl=xr;
else
xu=xr;
end
if f(xl)*f(xr)<0
xu=xr;
else
xl=xr;
end
xnew(1)=0;
xnew(i)=xr;
if abs((xnew(i)-xnew(i-1))/xnew(i))<tol,break,end
end
str = ['The required root of the equation is: ', num2str(xr), '']

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by