Can someone help me fix my bisection function?

The function is just adding up the upper and lower bound and then dividing by 2, it goes into the function but produces an incorrect root with the maximum iterations.

function [Root] = BisectionFunc25(func,Xu,Xl,RelE,MaxIt)
% Input: 
%       Func = the function
%       Xmin and Xmax = the lower and upper bounds
%       RelE = the desired relative absolute error- defaulted to
%       0.0001 percent
%       MaxIt = the maximum number of iterations - defaulted to 50
% Output: 
%       Root = the root
%       Funcroot = the function value at the root
%       AssE = the associated relative percent error
%       Steps = the number of iterations
MaxIt = 50;

% Set format long to show an adequate number of digits in the results. 
format long 
% Check the number of inputs. Need at least the first three, the last two
% default to select values. 
switch nargin 
    case 0 
        error ('Need Func, upper and lower limits; Defaults: desired relative abs error and max number of iterations')
        RelE = 0.0001 
        MaxIt = 50
        disp('Defaults: desired relative abs error and max number of iterations')
    case 1 
        error ('Need upper and lower limits')
        RelE = 0.0001 
        MaxIt = 50
        disp('Defaults: desired relative abs error and max number of iterations')
    case 2 
        RelE = 0.0001; 
        MaxIt = 50;
        Xu = 4;
        Xl = -4;
        fprintf('upper limit is %4d, lower limit is %4d', Xu, Xl)
end 
% Use the Bisection method to find the root to find func that will result in
% a pipe velocity of 5 m/s. 
Mid2 = (Xu + Xl)/2;
NumIt = 0;
% Initialize the number of iterations and set up the for loop for the
% method. 
for i = 1:MaxIt
    NumIt = NumIt + 1
    if func(Xu) * func(Mid2) < 0
        Mid1 = Mid2 
        Mid2 = (Xu + Xl)./2
        E = abs((Mid1 - Mid2)./Mid1) * 100
        if (E <= RelE || NumIt >= MaxIt)
            fprintf('Root found is %10.7f at %4.7f iterations', Mid2,NumIt)
        elseif func(Xu) * func(Mid2) > 0
            Mid1 = Mid2
            Mid2 = (Xu + Xl)./2
            E = abs((Mid2 - Mid1)./Mid2) * 100
            % Check if the new value of x is zero, then check if the error limit or the iteration limit has been reached
        elseif Mid2 == 0 
            E = abs((Mid2 - Mid1)./Mid2) * 100
            fprintf('Root found %10.7f')
            break
        end 
    end 
end 

Answer 1

[John D'Errico]

Good code. Well, at least in terms of the documentation. Excellent in that respect. Readable code is important when you need to debug it. It is also important if you get run down by the crosstown bus, and your colleague needs to use your code. So well done in that respect.

You did not return all of the outputs you listed. But my guess is, you were not done writing it, because you found a problem. The lack of semi-colons on some lines tells me you were not done, and were debugging your code. For that though, LEARN TO USE THE DEBUGGER!

One obvious problem lies at this line:

elseif Mid2 == 0

Surely you wanted to know if func(Mid2) is zero, NOT if Mid2 is zero.

Anyway, I think your problem lies in confusion about the endpoints of the interval. They need to change as necessary. Xl and Xu are not going to be constant. Let me give an example. I'll use bisection to determine the root of

F = @(x) 3*x - 2;

And yes, I'll bet you will guess it happens at x==2/3==0.6666666666... But let me try bisection. I'll set the initial interval of interest as [0,1]. Clearly that brackets a root, since there is a sign change between the endpoints.

Xbounds = [0,1];
Fbounds = F(Xbounds)

Fbounds = 1×2

-2 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Xmid = mean(Xbounds)

Xmid = 0.5000

fmid = F(Xmid)

fmid = -0.5000

Now what is the new interval? We choose the endpoint that will show a sign change from Fmid.

Xbounds = [Xmid,Xbounds(2)]

Xbounds = 1×2

0.5000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

F(Xbounds)

ans = 1×2

-0.5000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Now, we start a new iteration.

Xmid = mean(Xbounds),Fmid = F(Xmid)

Xmid = 0.7500

Fmid = 0.2500

And here, it is the LOWER bound that we will keep. Replace the upper bound with Xmid.

Xbounds = [Xbounds(1),Xmid]

Xbounds = 1×2

0.5000 0.7500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I hope you also see we are getting closer to x = 2/3 as the solution. One more iteration may help.

F(Xbounds),Xmid = mean(Xbounds),Fmid = F(Xmid)

ans = 1×2

-0.5000 0.2500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Xmid = 0.6250

Fmid = -0.1250

Xbounds = [Xmid,Xbounds(2)]

Xbounds = 1×2

0.6250 0.7500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Now, go back, and look at your code. Your code never updates the interval for bisection. It keeps Xu and Xl as unchanged.

Do you need to change your code to look like what I wrote? Well, I think my version is cleaner, by using the variable Xbounds as always the interval of interest, that brackets the root. It means I have fewer variables to keep track of, and that is ALWAYS a good thing. Regardless of what code you do write, you need to look at my example, and see how it tracks the root, updating the lower of upper endpoints of the bracket as necessary. That is the crucial part.