Bisection function to determine the value of (k) not working.
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
ABDALLA AL KHALEDI
am 18 Mär. 2023
Kommentiert: ABDALLA AL KHALEDI
am 18 Mär. 2023
Here is the bisection function
function [k,e] = bisect(f,a,b,n)
% function [k e] = mybisect(f,a,b,n)
% Does n iterations of the bisection method for a function f
% Inputs: f -- a function
% a,b -- left and right edges of the interval
% n -- the number of bisections to do.
% Outputs: x -- the estimated solution of f(x) = 0
% e -- an upper bound on the error
% evaluate at the ends and make sure there is a sign change
c = f(a);
d = f(b);
if c*d > 0
error('Function has same sign at both endpoints.')
end
disp(' k y')
for i = 1:n
% find the middle and evaluate there
k = (a + b)/2;
y = f(k);
disp([ k y])
if y == 0.0 % solved the equation exactly
a = k;
b = k;
break % jumps out of the for loop
end
% decide which half to keep, so that the signs at the ends differ
if c*y < 0
b=k;
else
a=k;
end
end
% set the best estimate for x and the error bound
k = (a + b)/2;
e = (b-a)/2;
end
And this is the function I am trying to apply it to (everything is known except for k, and the equation is equal to zero). For example, d = 0.5, H = 0.05, and L = 11.07.
(L/d)-((4*ellipticK(k))/(3*H/d)^0.5)*...
(1+(H/d)*((5/8)-(3/2)*(ellipticE(k)/ellipticK(k)))...
+((H/d)^2)*((-21/128)+(1/16)*(ellipticE(k)/ellipticK(k))+(3/8)*(ellipticE(k)/ellipticK(k))^2)...
+((H/d)^3)*((20127/179200)-(409/6400)*(ellipticE(k)/ellipticK(k))+(7/64)*(ellipticE(k)/ellipticK(k))^2+(1/16)*(ellipticE(k)/ellipticK(k))^3)...
+((H/d)^4)*((-1575087/28672000)+(1086367/1792000)*(ellipticE(k)/ellipticK(k))...
-(2679/25600)*(ellipticE(k)/ellipticK(k))^2+(13/128)*(ellipticE(k)/ellipticK(k))^3+(3/128)*(ellipticE(k)/ellipticK(k))^4))
2 Kommentare
Akzeptierte Antwort
Torsten
am 18 Mär. 2023
Bearbeitet: Torsten
am 18 Mär. 2023
Please fill in the missing values and try again.
If you don't succeed, it usually helps to plot f in the range 1e-12:1-1e-12 for k before calling a root finder.
L = ...
H = ...
d = ...
f = @(k) (L/d)-((4*ellipticK(k))/(3*H/d)^0.5)*...
(1+(H/d)*((5/8)-(3/2)*(ellipticE(k)/ellipticK(k)))...
+((H/d)^2)*((-21/128)+(1/16)*(ellipticE(k)/ellipticK(k))+(3/8)*(ellipticE(k)/ellipticK(k))^2)...
+((H/d)^3)*((20127/179200)-(409/6400)*(ellipticE(k)/ellipticK(k))+(7/64)*(ellipticE(k)/ellipticK(k))^2+(1/16)*(ellipticE(k)/ellipticK(k))^3)...
+((H/d)^4)*((-1575087/28672000)+(1086367/1792000)*(ellipticE(k)/ellipticK(k))...
-(2679/25600)*(ellipticE(k)/ellipticK(k))^2+(13/128)*(ellipticE(k)/ellipticK(k))^3+(3/128)*(ellipticE(k)/ellipticK(k))^4));
a = ...
b = ...
n = ...
[k e] = bisect(f,a,b,n)
function [k,e] = bisect(f,a,b,n)
% function [k e] = mybisect(f,a,b,n)
% Does n iterations of the bisection method for a function f
% Inputs: f -- a function
% a,b -- left and right edges of the interval
% n -- the number of bisections to do.
% Outputs: x -- the estimated solution of f(x) = 0
% e -- an upper bound on the error
% evaluate at the ends and make sure there is a sign change
c = f(a);
d = f(b);
if c*d > 0
error('Function has same sign at both endpoints.')
end
disp(' k y')
for i = 1:n
% find the middle and evaluate there
k = (a + b)/2;
y = f(k);
disp([ k y])
if y == 0.0 % solved the equation exactly
a = k;
b = k;
break % jumps out of the for loop
end
% decide which half to keep, so that the signs at the ends differ
if c*y < 0
b=k;
else
a=k;
end
end
% set the best estimate for x and the error bound
k = (a + b)/2;
e = (b-a)/2;
end
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Creating and Concatenating Matrices finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!