Filter löschen
Filter löschen

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

1 Ansicht (letzte 30 Tage)
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
Torsten
Torsten am 18 Mär. 2023
How do you call "bisect" ?
ABDALLA AL KHALEDI
ABDALLA AL KHALEDI am 18 Mär. 2023
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.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Torsten
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)

Kategorien

Mehr zu Creating and Concatenating Matrices finden Sie in Help Center und File Exchange

Produkte


Version

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by