Bisection Method Using MATLAB

Bisection method is based on Intermediate Value Theorem which states as,
 Let f (x) be a continuous function on the interval [a, b]. If d  [f (a), f (b)], then there is a c  [a, b] such that f (c) = d.

We are using a sub-form  form of this theorem, In our case we are given with two values at which function changes its value form positive to negative and we have to find a value where function is zero or in other words root of the function. We then use a binary search type algorithm that slowly decreases this interval and finally gets solution with some error. 

So if we are given with a interval [xl , xu] as lower and upper limit. Value of function with xl is negative and with  xu is positive then there must be at least one point in between these two values which is  root of the function  or in other words function becomes zero. Function should be real and continuous.

Algorithm:

1. We check weather function changes its value from positive to negative within given interval, If it changes we proceed other wise we can't find solution.
2. We find mid point of the interval, xm. 
3. If function at  xm is negative, it becomes our new xl, if its positive its our new xu and if it happens to be zero, we find our root (It rarely happens). 
4. We find error by seeing the difference between last and new xm. If error is smaller than provided one, we stops other wise step 2 and 3 repeats.
5. If we also limit iterations by 1000. Means our program can execute more than 1000 iterations.

MATLAB Code:

%Author: Muhammad Awais
%UET Lahore, Electrical Engineering
%fb/awais12506
function [iteration,xm]=Bisection(xl,xu,f,e)

%Find Root of a equation by method of BiSection
%Input: Lower Limit, Upper Limit, Function, Error Tolerance
%Please Insert f as f=@(x)x.^2+9*x+3
%Output: Root of the equation with given precison
%Exception: Give error if Equation is not convergent or roots are dont give

%postive and negative values of the function

iteration=0;     %To see how many iterations, equation took to solve
xm=xl;
error=1;
if  (f(xl)*f(xu)>0)
    disp('Interval have some error')

else
     while ( abs(error) > abs(e) ||iteration<=1000 )
        iteration=iteration+1;
        xmOld=xm;
        xm=(xl+xu)/2;   %Mid point
        xmNew=xm;
        if f(xl)*f(xm)<0   %if f(xm) is neg, equate xu with xm
            xu=xm;
        else if f(xl)*f(xm)>0 %if f(xm) is pos, equate xl with xm
            xl=xm;
            else            %it means xm is our root
            break
            end
        error=(xmNew-xmOld)/xmNew;   %Error
        end

    end

end
end

Output:

References:

http://www.sosmath.com/calculus/limcon/limcon07/limcon07.html