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

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

[Torsten]

How do you call "bisect" ?

[ABDALLA AL KHALEDI]

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.

Answer 1

[Torsten]

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

[ABDALLA AL KHALEDI]

Thanks a lot, It works now