Numerical Methods

The Bisection Method

  Bisection method is a technique to find the roots of algebraic and transcendental equations of the form `f(x)=0` such as:

`xe^x - 1 = 0`.

This technique is based on the Intermediate Value Theorem which states that if `f(x)` is a continuous function in `[a,b]` and if `q` is any number between `f(a)` and `f(b)` ,then, there exists a number `c` in `(a,b)` such that `f(q)=c`. By corollary, if `f(a)` and `f(b)` are of opposite sign, there must be a number `q` in between `a` and `b` such that `f(q)=0` or in other words, `q` is the solution of the equation.


To find the root of `f(x)=0` we have to choose two numbers `x_L` ans `x_H ` such that `f(x_L)` and `f(x_H)` are of opposite sign, for example, let `f(x_L)` is negative and `f(x_H)` is positive. Then, the root must lie between `x_L` and `x_H`, and lets assume that the approximate root is given by `x_"mid" =(x_L+x_H)/2`,ie, the midpoint of the interval. If `f(x_"mid")=0` then, `x_"mid"` is the root and there is no need to proceed further. Otherwise, if `f(x_"mid")>0` the root must lie between `x_L` and `x_"mid"` and we set `x_H = x_"mid"` while `x_L` remains unchanged. Else, if `f(x_"mid") < 0` the root must lie between `x_"mid"` and `x_H` we set `x_L = x_"mid"` and `x_H` remains unchanged. The process is repeated to find the next interval and so on until the width of the interval is within the desired limit or there is negligible difference between the successive midpoints.


Visualising Bisection Method:

Algorithm:

Step 1: Read `x_L, x_H and epsilon` such that `f(x_L)`is negative and `f(x_H)` is positive.

Step 2: Compute `x_"mid"=(x_L+x_H)/2`

Step 3: `previousX= x_"mid"`

Step 4: If `f(x_L)f(x_"mid")<0, x_H=x_"mid"`

Step 5: If `f(x_L)f(x_"mid")>0, x_L=x_"mid"`

Step 6: Compute `x_"mid"=(x_L+x_H)/2`

Step 7: If ` |(previousX-x_"mid")/x_"mid"|>epsilon` go to Step:3

Step 8: OUTPUT :The root is approximately equal to `x_"mid"`

Step 9: Stop

Following are programs in different programming languages to find the root of the equation

`xe^x - 1 = 0`

Program in C
======================================================================
 Progaram to find solutions to equations of the form f(x)=0 using      
 Bisection Method (Bolanzo method)										
 LANGUAGE  : C   (Compiler: GNU GCC)                                  	
=======================================================================

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double func(double x){
    double y;
    y=x*exp(x)-1;
    return y;
}

int main()
{
    double xL,xH,e;     //e is the tolerance
    printf("Enter two real numbers between which the root is to be found\n");
    scanf("%lf %lf",&xL,&xH);
    printf("\nEnter the tolerance\n");
    scanf("%lf",&e);
    double xMid = (xL+xH)/2;
    double prev_xMid;

    do{
        prev_xMid=xMid;
        if(func(xL)*func(xMid)<0){
            xH=xMid;
        }
        if(func(xL)*func(xMid)>0){
            xL=xMid;
        }
        xMid=(xL+xH)/2;
    }while(fabs((prev_xMid-xMid)/xMid)>e);

    printf("\n Root is %lf ",xMid);
    return 0;
}

/*---------------------------  OUTPUT  ----------------------------------/
	Enter two real numbers between which the root is to be found
	0
	1

	Enter the tolerance
	0.000001

	Root is 0.567143
------------------------------------------------------------------------*/

Program in Fortran 95
!======================================================================!
!Progaram to find solutions to equations of the form f(x)=0 using 
!bisection method (Bolanzo method)                                     
!LANGUAGE     :: FORTRAN 95  (compiler: G95 Fortran Compiler)          !
!======================================================================!

PROGRAM bisection
    IMPLICIT NONE
    REAL :: xL,xH,tolerance,d,xMid
    INTEGER :: iterations = 0
    WRITE (*,*) 'Enter the two numbers between which root is  to be found'
    READ (*,*) xL,xH
    WRITE (*,*) 'The root should be correct upto how many decimal places?'
    READ (*,*) tolerance

    IF ((f(xL)*f(xH))== 0) THEN
        IF ((f(xL)==0).AND.(f(xH)==0)) THEN
            WRITE (*,*) 'The desired roots are',xL,' and ',xH,' themselves'
        ELSE IF (f(xL)==0) THEN
            WRITE (*,*) 'The desired root is',xL
        ELSE IF (f(xH)==0) THEN
            WRITE (*,*) 'The desired root is',xH
        END IF
    ELSE IF ((f(xL)*f(xH))>0) THEN
        WRITE (*,*) 'The root does not lie between ',xL,' and ',xH
    ELSE IF ((f(xL)*f(xH))<0) THEN
        DO
            xMid=(xL+xH)/2
            IF (f(xMid)==0) EXIT
            d=ABS((xMid-xL)/xMid)
            IF (d<((0.1)**tolerance)) EXIT             !accuracy
            IF ((f(xMid)*f(xL))<0) THEN            !root lies between xMid and xL
                xH=xL
                xL=xMid
                !PRINT 10,xL,xH,d,'g'
                !10 FORMAT (F12.7,3X,F12.7,3X,F12.7,3X,A2)
            ELSE IF ((f(xMid)*f(xH))<0) THEN       !root lies between xMid and xH
                xL=xMid
                xH=xH
            END IF
            iterations = iterations + 1
        END DO
        WRITE (*,*) 'The root is',xMid
        WRITE (*,*) 'Total iterations =',iterations
   END IF

   CONTAINS
   
		FUNCTION f(y)
			REAL :: f
			REAL,INTENT (IN)::y
			f=y*EXP(y)-1
		END FUNCTION f

END PROGRAM

!-------------------------OUTPUT-----------------------------!
! Enter the two numbers between which root is  to be found   !
! 0.1														 !
! 0.9														 !
! The root should be correct upto how many decimal places?   !
! 6							   							     !
! The root is 0.56714296  									 !
! Total iterations = 20										 !
!------------------------------------------------------------!
Program in Scilab
======================================================================
 Progaram to find solutions to equations of the form f(x)=0 using      
 Bisection Method (Bolanzo method)										
 LANGUAGE  : SCILAB                                  	
=======================================================================

iterations = 0
xL = input("Enter number at which f(x) is negative: ")
xH = input("Enter number at which f(x) is positive: ")
tolerance = input("The root should be correct upto how many decimal places?: ")
if ((f(xL)*f(xH))== 0) then
    if ((f(xL)==0)&(f(xH)==0)) then
        printf("The desired roots are %f and %f themselves",xL,xH)
    elseif (f(xL)==0) then
        printf ("The desired root is %f",xL)
    elseif (f(xH)==0) then
        printf ("The desired root is %f",xH)
    end
elseif ((f(xL)*f(xH))>0) then
    printf ("The root does not lie between %f and %f",xL,xH) 
elseif ((f(xL)*f(xH))<0) then
    while (iterations>=0)
        xMid=(xL+xH)/2
        if (f(xMid)==0) break;end
        d=abs((xMid-xL)/xMid)
        if (d<((0.1)**tolerance)) break;end            //accuracy
        if ((f(xMid)*f(xL))<0) then            //root lies between xMid and xL
            xH = xMid
        elseif ((f(xMid)*f(xL))>0) then       //root lies between xMid and xH
            xL=xMid
        end
        iterations = iterations+1
    end
    printf ("The root is %f",xMid) 
    printf ("\nTotal iterations =%d",iterations)
end

function y = f(x)
    y=x*exp(x)-1
endfunction

======================= OUTPUT =============================
Enter number at which f(x) is negative: 0.1

Enter number at which f(x) is positive: 0.9

The root should be correct upto how many decimal places?: 6

The root is 0.567144
Total iterations =20
===========================================================