Trying to solve for roots

1 view (last 30 days)
David Ziehl on 16 Nov 2021
Commented: David Ziehl on 16 Nov 2021
I am trying to solve for the velocity in this problem. I have attached a photo of my work on paper, but I was trying to set it up in matlab using bisection method, but it calls error when I try to input the large term the includes velocity. I have attached my work on paper and my code. One other thing is my work on paper, in the denominator there should be two parentheses on the right, signifying that the entire denominator is being squared. When I run it, it says my parenthesis are wrong, but I thought I did them correctly. Please help.
% Bisection Method Code #2
max_error=.001;
itr_max=100;
g=9.8
L=80
D=.1
Km=6.5
Re=4826 %coefficient of v in reynolds number
e=.00015
f = @(x) .625-(x^2/(2*g))*(Km+(L/D)*1/(-1.8log((e/(D*3.7))^1.11+6.9/(Re*x)))^2) ; %function
xl=0; %lower bound
xu=100; %upper bound
xr0=0; %initial guess
for i = 1:itr_max
fxl=f(xl);
fxu=f(xu);
xr=(xl+xu)/2;
fxr=f(xr);
err=abs((xr-xr0)/xr)*100;
ferr=abs((fxr-f(xr0))/fxr)*100;
xr0=xr;
fprintf('i=%d \t xl=%f \t xu=%f \t fxl=%f \t fxu=%f \t xr=%f \t fxr=%f \t er=%f\n',i,xl,xu,fxl,fxu,xr,fxr,err)
if err < max_error
break;
end
test1=fxr*fxl;
test2=fxr*fxu;
if test1<0
xu=xr;
elseif test2<0
xl=xr;
end
end

Cris LaPierre on 16 Nov 2021
MATLAB does not do implied multiplicaton. You are getting the error in your first image because you have
f = @(x) .625-(x^2/(2*g))*(Km+(L/D)*1/(-1.8log((e/(D*3.7))^1.11+6.9/(Re*x)))^2); %function
% ^^ implied multiplication
Specify the operation to take place between -1.8 and log to fix this error.
f = @(x) .625-(x^2/(2*g))*(Km+(L/D)*1/(-1.8*log((e/(D*3.7))^1.11+6.9/(Re*x)))^2); %function
% ^
Your code runs without error once that change is made. Just be aware that, in MATLAB, log is the natural log, or ln. If you want log base 10, that is log10.
% Bisection Method Code #2
max_error=.001;
itr_max=100;
g=9.8;
L=80;
D=.1;
Km=6.5;
Re=4826; %coefficient of v in reynolds number
e=.00015;
f = @(x) .625-(x^2/(2*g))*(Km+(L/D)*1/(-1.8*log((e/(D*3.7))^1.11+6.9/(Re*x)))^2); %function
xl=0; %lower bound
xu=100; %upper bound
xr0=0; %initial guess
for i = 1:itr_max
fxl=f(xl);
fxu=f(xu);
xr=(xl+xu)/2;
fxr=f(xr);
err=abs((xr-xr0)/xr)*100;
ferr=abs((fxr-f(xr0))/fxr)*100;
xr0=xr;
fprintf('i=%d \t xl=%f \t xu=%f \t fxl=%f \t fxu=%f \t xr=%f \t fxr=%f \t er=%f\n',i,xl,xu,fxl,fxu,xr,fxr,err)
if err < max_error
break;
end
test1=fxr*fxl;
test2=fxr*fxu;
if test1<0
xu=xr;
elseif test2<0
xl=xr;
end
end
i=1 xl=0.000000 xu=100.000000 fxl=0.625000 fxu=-5023.050705 xr=50.000000 fxr=-1262.750501 er=100.000000 i=2 xl=0.000000 xu=50.000000 fxl=0.625000 fxu=-1262.750501 xr=25.000000 fxr=-318.703577 er=100.000000 i=3 xl=0.000000 xu=25.000000 fxl=0.625000 fxu=-318.703577 xr=12.500000 fxr=-80.759644 er=100.000000 i=4 xl=0.000000 xu=12.500000 fxl=0.625000 fxu=-80.759644 xr=6.250000 fxr=-20.370447 er=100.000000 i=5 xl=0.000000 xu=6.250000 fxl=0.625000 fxu=-20.370447 xr=3.125000 fxr=-4.877890 er=100.000000 i=6 xl=0.000000 xu=3.125000 fxl=0.625000 fxu=-4.877890 xr=1.562500 fxr=-0.845006 er=100.000000 i=7 xl=0.000000 xu=1.562500 fxl=0.625000 fxu=-0.845006 xr=0.781250 fxr=0.223444 er=100.000000 i=8 xl=0.781250 xu=1.562500 fxl=0.223444 fxu=-0.845006 xr=1.171875 fxr=-0.230346 er=33.333333 i=9 xl=0.781250 xu=1.171875 fxl=0.223444 fxu=-0.230346 xr=0.976562 fxr=0.016848 er=20.000000 i=10 xl=0.976562 xu=1.171875 fxl=0.016848 fxu=-0.230346 xr=1.074219 fxr=-0.101702 er=9.090909 i=11 xl=0.976562 xu=1.074219 fxl=0.016848 fxu=-0.101702 xr=1.025391 fxr=-0.041163 er=4.761905 i=12 xl=0.976562 xu=1.025391 fxl=0.016848 fxu=-0.041163 xr=1.000977 fxr=-0.011841 er=2.439024 i=13 xl=0.976562 xu=1.000977 fxl=0.016848 fxu=-0.011841 xr=0.988770 fxr=0.002583 er=1.234568 i=14 xl=0.988770 xu=1.000977 fxl=0.002583 fxu=-0.011841 xr=0.994873 fxr=-0.004609 er=0.613497 i=15 xl=0.988770 xu=0.994873 fxl=0.002583 fxu=-0.004609 xr=0.991821 fxr=-0.001008 er=0.307692 i=16 xl=0.988770 xu=0.991821 fxl=0.002583 fxu=-0.001008 xr=0.990295 fxr=0.000788 er=0.154083 i=17 xl=0.990295 xu=0.991821 fxl=0.000788 fxu=-0.001008 xr=0.991058 fxr=-0.000110 er=0.076982 i=18 xl=0.990295 xu=0.991058 fxl=0.000788 fxu=-0.000110 xr=0.990677 fxr=0.000339 er=0.038506 i=19 xl=0.990677 xu=0.991058 fxl=0.000339 fxu=-0.000110 xr=0.990868 fxr=0.000115 er=0.019249 i=20 xl=0.990868 xu=0.991058 fxl=0.000115 fxu=-0.000110 xr=0.990963 fxr=0.000003 er=0.009624 i=21 xl=0.990963 xu=0.991058 fxl=0.000003 fxu=-0.000110 xr=0.991011 fxr=-0.000053 er=0.004812 i=22 xl=0.990963 xu=0.991011 fxl=0.000003 fxu=-0.000053 xr=0.990987 fxr=-0.000025 er=0.002406 i=23 xl=0.990963 xu=0.990987 fxl=0.000003 fxu=-0.000025 xr=0.990975 fxr=-0.000011 er=0.001203 i=24 xl=0.990963 xu=0.990975 fxl=0.000003 fxu=-0.000011 xr=0.990969 fxr=-0.000004 er=0.000601
David Ziehl on 16 Nov 2021
Thank you very much for the reply. My code works now.
Regards,
David

R2021b

Community Treasure Hunt

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

Start Hunting!