fbpx

Bisection Method: Squeeze that root !

0
Lesson 2: The Bisection Method

Bisection Method: What you need to know

Bisection method and its MATLAB implementation are detailed in this article. In addition, we give a bisection method convergence analysis and explain why the bisection method enjoys a linear type convergence speed.

As in my previous post on the fixed point method, the goal here is to find a root of

    \[f(x) = 0\]

The bisection method works on an interval [a,b] and the initial guess could be the midpoint frac{a+b}{2} = x_0. To see if x_1 is the midpoint of [a,x_0] or [x_0,b], we have to check where the root falls, namely if the root falls in [a,x_0], then we take the midpoint that interval, else we choose [x_0,b], and so on. Stopping criteria could be when vert x_k - x_{k-1} vert < epsilon, where epsilon is a tolerance parameter.

How Bisection Method Works

Well, the only information we need is the starting interval boundaries [a,b] that has to satisfy

    \[f(a)f(b)<0\]

so that we are sure there is a root in [a,b] worth iterating and searching over. The bisection method starts by computing

    \[x_0 = \frac{a+b}{2}\]

which is actually the midpoint of [a,b] followed by a boundary test, i.e. if

    \[f(a)f(x_0) < 0\]

then the new interval would be [a,x_0] else it would be [x_0,b]. Said differently, if

    \[f(a)f(x_0) < 0\]

then replace b by x_0 else replace a by x_0 and repeat this process over and over again until \vert x_k - x_{k-1} \vert < \epsilon where \epsilon is a user pre-defined value.

MATLAB implementation of Bisection method

This is easily achieved on MATLAB. First let’s write a bisection function so that we can re-use it over and over again

function [x,error] = bisection(f,a,b,Niter)

% Input:
%   f - input funtion at [a,b] 
%   a - left limit interval [a,b]
%   b - right limit interval [a,b]
%   Niter - Number of iterations
%% Output:
%   x     - is a vector giving root at each iteration
%   error - is a vector giving  error at each iteration

Now let’s write the bisection main function.

function [x,error] = bisection(f,a,b,Niter)

% Input:
%   f - input funtion at [a,b] 
%   a - left limit interval [a,b]
%   b - right limit interval [a,b]
%   Niter - Number of iterations
%% Output:
%   x     - is a vector giving root at each iteration
%   error - is a vector giving  error at each iteration


if f(a)*f(b)>0 
    disp('You have not respected the Mean Valued Theorem')
else
    % Take midpoint
    x(1) = (a + b)/2;
    %compute error
    error = abs(f(x));
    for n = 2:Niter
   if f(a)*f(x(n-1))<0 
       % use b
       b = x(n-1);
   else
       % else use b
       a = x(n-1);          
   end
   % Take NEW midpoint
    x(n) = (a + b)/2; 
     %compute NEW error
   error(n-1) = abs(f(x(n)));
    end
end

As one can see from the code above, all we did was first test the condition of success of the bisection method, i.e.

    \[f(a)f(b)<0\]

then compute the midpoint

    \[x_k = \frac{a+b}{2}\]

followed by a new interval selection based off

    \[f(x_k)f(a) < 0\]

. Now in the main function, we can initialize some parameters like:

Niter = 20;

Note that in the above, we have assumed the stopping criterion of max number of iterations instead of min tolerance. Next, we can proceed to defining some functions to test, for example:

Progress .. You’re almost done with this article 🙂
%given functions
f1 = @(x)x^3 - 2*x -5; 
f2 = @(x)exp(-x) - x;
f3 = @(x)x*sin(x) -1;
f4 = @(x)x^3 - 3*x^2 +3*x - 1;

Up till now in the main function, we have defined

    \[f_1(x) = x^3 - 2x - 5\]

    \[f_2(x) = e^{-x} - x\]

    \[f_3(x) = x\sin(x) - 1\]

and

    \[f_4(x) = x^3 - 3x^2 + 3x - 1\]

We can now solve each function using

%solve using bisection, newton, secant
[x_1_bisection,error_1_bisection] = bisection(f1,0,5,Niter);
[x_2_bisection,error_2_bisection] = bisection(f2,0,5,Niter);
[x_3_bisection,error_3_bisection] = bisection(f3,0,2,Niter);
[x_4_bisection,error_4_bisection] = bisection(f4,0,5,Niter);

For comparison sake, we will also check whether we converged properly or not. We can either us the error function outputted by our bisection.m function or we can use MATLAB’s fsolve function as

%solving using a library routine
x_1 = fsolve(f1,1);
x_2 = fsolve(f2,1);
x_3 = fsolve(f3,1);
x_4 = fsolve(f4,1);

Right now, we are ready to plot

figure
subplot(2,2,1)
plot(x_1_bisection,'r','Linewidth',1)
hold on
plot(x_1*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n')
ylabel('x')
title('Solving f(x) = x^3 - 2x - 5')
legend('Bisection','Library Routine')
grid on
grid minor

subplot(2,2,2)
plot(x_2_bisection,'r','Linewidth',1)
hold on
plot(x_2*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n')
ylabel('x')
title('Solving f(x) = e^{-x} - x')
legend('Bisection','Library Routine')
grid on
grid minor


subplot(2,2,3)
plot(x_3_bisection,'r','Linewidth',1)
hold on
plot(x_3*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n')
ylabel('x')
title('Solving f(x) = xsin(x) - 1')
legend('Bisection','Library Routine')
grid on
grid minor

subplot(2,2,4)
plot(x_4_bisection,'r','Linewidth',1)
hold on
plot(x_4*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n')
ylabel('x')
title('Solving f(x) = x^{3} - 3x^{2} +3x - 1')
legend('Bisection','Library Routine')
grid on
grid minor

Our complete main function looks like this

Niter = 20;

%given functions
f1 = @(x)x^3 - 2*x -5; 
f2 = @(x)exp(-x) - x;
f3 = @(x)x*sin(x) -1;
f4 = @(x)x^3 - 3*x^2 +3*x - 1;

%solve using bisection, newton, secant
[x_1_bisection,error_1_bisection] = bisection(f1,0,5,Niter);
[x_2_bisection,error_2_bisection] = bisection(f2,0,5,Niter);
[x_3_bisection,error_3_bisection] = bisection(f3,0,2,Niter);
[x_4_bisection,error_4_bisection] = bisection(f4,0,5,Niter);

%solving using a library routine
x_1 = fsolve(f1,1);
x_2 = fsolve(f2,1);
x_3 = fsolve(f3,1);
x_4 = fsolve(f4,1);

figure
subplot(2,2,1)
plot(x_1_bisection,'r','Linewidth',1)
hold on
plot(x_1*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n')
ylabel('x')
title('Solving f(x) = x^3 - 2x - 5')
legend('Bisection','Library Routine')
grid on
grid minor

subplot(2,2,2)
plot(x_2_bisection,'r','Linewidth',1)
hold on
plot(x_2*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n')
ylabel('x')
title('Solving f(x) = e^{-x} - x')
legend('Bisection','Library Routine')
grid on
grid minor


subplot(2,2,3)
plot(x_3_bisection,'r','Linewidth',1)
hold on
plot(x_3*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n')
ylabel('x')
title('Solving f(x) = xsin(x) - 1')
legend('Bisection','Library Routine')
grid on
grid minor

subplot(2,2,4)
plot(x_4_bisection,'r','Linewidth',1)
hold on
plot(x_4*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n')
ylabel('x')
title('Solving f(x) = x^{3} - 3x^{2} +3x - 1')
legend('Bisection','Library Routine')
grid on
grid minor
bisection method
bisection method

Bisection Method Caculator

Well, now that the bisection method matlab code is done, all we have to do is run the code on MATLAB or any octave compiler. A small note here would be that this code would have a very similar python equivalent. So with minor tweaks, you can get the code running on python as well.

Convergence Speed

The bisection method has a linear convergence speed. We can also view this on MATLAB, i.e. by plotting convergence error (\vert x_k - x_{k-1} \vert) vs. the number of iterations. Mathematically speaking, this means that

    \[\vert x^* - x_{k}\vert \leq \left(\frac{1}{2}\right)^{k-1}|b-a|\]

This means that at any given iteration k, we are sure that the error of the form | x^* - x_{k} | is always lower bounded by a function which decays exponentially with the iteration number (here I refer to the term \frac{1}{2}\right)^{k-1} multiplied by the interval length. Why do we say linear convergence ? Well, in many convergence analysis studies, the error expressed on a log-scale and hence applying the logs on both sides of the above inequality, we have that

    \[\log | x^* - x_{k} | \leq \log (\frac{1}{2})^{k-1}| b-a |\]

which gives

    \[\log | x^* - x_{k} | \leq \(k-1) \log (\frac{1}{2}) | b-a |\]

hence upper bounded by \mathcal{O}(k)


Summary

In this article, we have detailed the bisection method along with its conditions for proper operation. We also sliced down the MATLAB code so that we explain each and every block appearing in our youtube lecture. Even more, we provide a running example on four functions that are defined above to demonstrate convergence. We also presented a convergence analysis bound that proves that

Follow my other courses such as convex optimization.

Buy me a cup of coffee using the donate link below 😊☕️

[paypal_donation_block email=’myprivateotherstuff@gmail.com’ amount=’15.00′ currency=’USD’ size=’large’ purpose=’Buy me coffee ^^’ mode=’live’]

PS: I’m on twitter. I retweet stuff around algorithms, python, MATLAB and mathematical optimization, mostly convex.

PSS: We are so close to 100K subscribers on YouTube. It would mean so much if you could share the channel and subscribe to help us sustain.