Bisection function to determine the value of (k) not working.

2 次查看（过去 30 天）
ABDALLA AL KHALEDI
ABDALLA AL KHALEDI 2023-3-18
评论： ABDALLA AL KHALEDI 2023-3-18
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 个评论
ABDALLA AL KHALEDI
ABDALLA AL KHALEDI 2023-3-18
I am trying to call the function in the editor window to find the value of k, all the other variables in the equation are known and the initial bisection limits are 10^(-12) to 1-10^(-12). I tried every thing I knew and found but couldn't get it to run properly.

请先登录，再进行评论。

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

采纳的回答

Torsten
Torsten 2023-3-18
编辑：Torsten 2023-3-18
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by