# The Bisection (Bozano) Equation Solver

### Problem Statement

Given a continuous equation f(x)=0 and two values a and b (a < b), if f(a)*f(b) < 0 (i.e., f(a) and f(b) have opposite signs), it can be proved that there exists a root of f(x)=0 between a and b. More precisely, there exists a c, a <= c <= b, such that f(c)=0 holds.

This result provides us with a method for solving equations. If we take the midpoint of a and b, c=(a+b)/2, and computes its function value f(c), we have the following cases:

1. If f(c) is very small (i.e., smaller than a tolerance value), then c can be considered as a root of f(x)=0 and we are done.
2. Otherwise, since f(a) and f(b) have opposite signs, the sign of f(c) is either identical to that of f(a) or that of f(b).
• If the sign of f(c) is different from that of f(a), then since f(a) and f(c) have opposite signs, f(x)=0 has a root in the range of a and c.
• If the sign of f(c) is different from that of f(b), then since f(b) and f(c) have opposite signs, f(x)=0 has a root in the range of c and b.
3. Therefore, if the original interval [a,b] is replaced with [a,c] in the first case or replaced with [c,b] in the second case, the length of the interval is reduced by half. If continue this process, after a number of steps, the length of the interval could be very small. However, since this small interval still contains a root of f(x)=0, we actually find an approximation of that root!
Write a program that contains two functions: (1) Funct() - the function f(x) and (2) Solve() - the equation solver based on the above theory. Then, reads in a and b and uses function Solve() to find a root in the range of a and b. Note that before calling Solve, your program should check if f(a)*f(b)<0 holds.

Your function Funct() is given below:

This function has a root near -0.89.

Since this process keeps dividing the intervals into two equal halves, it is usually referred to as the bisection method. It is also known as Bozano's method.

### Solution

```! --------------------------------------------------------------------
!    This program solves equations with the Bisection Method.  Given
! a function f(x) = 0.  The bisection method starts with two values,
! a and b such that f(a) and f(b) have opposite signs.  That is,
! f(a)*f(b) < 0.  Then, it is guaranteed that f(x)=0 has a root in
! the range of a and b.  This program reads in a and b (Left and Right
! in this program) and find the root in [a,b].
!    In the following, function f() is REAL FUNCTION Funct() and
! solve() is the function for solving the equation.
! --------------------------------------------------------------------

PROGRAM  Bisection
IMPLICIT  NONE

REAL, PARAMETER :: Tolerance = 0.00001
REAL            :: Left,  fLeft
REAL            :: Right, fRight
REAL            :: Root

WRITE(*,*)  'This program can solves equation F(x) = 0'
WRITE(*,*)  'Please enter two values Left and Right such that '
WRITE(*,*)  'F(Left) and F(Right) have opposite signs.'
WRITE(*,*)
WRITE(*,*)  'Left and Right please --> '
READ(*,*)   Left, Right              ! read in Left and Right

fLeft  = Funct(Left)                 ! compute their function values
fRight = Funct(Right)
WRITE(*,*)
WRITE(*,*)  'Left = ', Left, '    f(Left) = ', fLeft
WRITE(*,*)  'Right = ', Right, '   f(Right) = ', fRight
WRITE(*,*)
IF (fLeft*fRight > 0.0)  THEN
WRITE(*,*)  '*** ERROR: f(Left)*f(Right) must be negative ***'
ELSE
Root = Solve(Left, Right, Tolerance)
WRITE(*,*)  'A root is ', Root
END IF

CONTAINS

! --------------------------------------------------------------------
! REAL FUNCTION  Funct()
!    This is for function f(x).  It takes a REAL formal argument and
! returns the value of f() at x.  The following is sample function
! with a root in the range of -10.0 and 0.0.  You can change the
! expression with your own function.
! --------------------------------------------------------------------

REAL FUNCTION  Funct(x)
IMPLICIT  NONE
REAL, INTENT(IN) :: x
REAL, PARAMETER  :: PI = 3.1415926
REAL, PARAMETER  :: a  = 0.8475

Funct = SQRT(PI/2.0)*EXP(a*x) + x/(a*a + x*x)

END FUNCTION  Funct

! --------------------------------------------------------------------
! REAL FUNCTION  Solve()
!    This function takes Left - the left end, Right - the right end,
! and Tolerance - a tolerance value such that f(Left)*f(Right) < 0
! and find a root in the range of Left and Right.
!    This function works as follows.  Because of INTENT(IN), this
! function cannot change the values of Left and Right and therefore
! the values of Left and Right are saved to a and b.
!    Then, the middle point c=(a+b)/2 and its function value f(c)
! is computed.  If f(a)*f(c) < 0, then a root is in [a,c]; otherwise,
! a root is in [c,b].  In the former case, replacing b and f(b) with
! c and f(c), we still maintain that a root in [a,b].  In the latter,
! replacing a and f(a) with c and f(c) will keep a root in [a,b].
! This process will continue until |f(c)| is less than Tolerance and
! hence c can be considered as a root.
! --------------------------------------------------------------------

REAL FUNCTION  Solve(Left, Right, Tolerance)
IMPLICIT  NONE
REAL, INTENT(IN) :: Left, Right, Tolerance
REAL             :: a, Fa, b, Fb, c, Fc

a = Left                          ! save Left and Right
b = Right

Fa = Funct(a)                     ! compute the function values
Fb = Funct(b)
IF (ABS(Fa) < Tolerance) THEN     ! if f(a) is already small
Solve = a                      ! then a is a root
ELSE IF (ABS(Fb) < Tolerance) THEN     ! is f(b) is small
Solve = b                      ! then b is a root
ELSE                              ! otherwise,
DO                             ! iterate ....
c  = (a + b)/2.0            !   compute the middle point
Fc = Funct(c)               !   and its function value
IF (ABS(Fc) < Tolerance) THEN    ! is it very small?
Solve = c                ! yes, c is a root
EXIT
ELSE IF (Fa*Fc < 0.0) THEN  ! do f(a)*f(c) < 0 ?
b  = c                   ! replace b with c
Fb = Fc                  ! and f(b) with f(c)
ELSE                        ! then f(c)*f(b) < 0 holds
a  = c                   ! replace a with c
Fa = Fc                  ! and f(a) with f(c)
END IF
END DO                         ! go back and do it again
END IF
END FUNCTION  Solve

END PROGRAM  Bisection
```
Click here to download this program.

### Program Input and Output

The following is the output from the above program with correct input:
```This program solves equation F(x) = 0
Please enter two values Left and Right such that
F(Left) and F(Right) have opposite signs.

Left and Right please -->
-10.0  0.0

Left = -10.    f(Left) = -9.902540594E-2
Right = 0.E+0   f(Right) = 1.25331414

A root is -0.89050293
```
The following output shows that the function values of the input do not have opposite signs and hence program stops.
```This program solves equation F(x) = 0
Please enter two values Left and Right such that
F(Left) and F(Right) have opposite signs.

Left and Right please -->
-10.0  -1.0

Left = -10.    f(Left) = -9.902540594E-2
Right = -1.   f(Right) = -4.495930672E-2

*** ERROR: f(Left)*f(Right) must be negative ***
```

### Discussion

• Function Funct(x) returns the function value at x.
• The heart of this program is function Solve(). It takes three arguments Left, Right and Tolerance and finds a root in the range of Left and Right.
• Since arguments Left and Right are declared with INTENT(IN), their values cannot be changed. As a result, their values are copied to variables a and b.
• The function values of a and b are stored in Fa and Fb.
• If the absolute value of Fa (resp., Fb) is very small, one can consider a (resp., b) as a root of the given equation.
• Otherwise, the midpoint c of a and b and its function value Fc are computed. If Fc is very small, c is considered as a root.
• Then if Fa and Fc have opposite signs, replacing b and Fb with c and Fc, respectively.
• If Fa and Fc have the same sign, then Fc and Fb must have the opposite sign. In this case, a and Fa are replaced with c and Fc.
• Either way, the original interval [a,b] is replaced with a new one with half length. This process continues until the absolute value of the function value at c, Fc, is smaller than the given tolerance value. In this case, c is considered a root of the given equation.