The bisection method (rtbis) - Numerical Recipes in F90

 

The easiest way to find roots

You may, of course, also find roots graphically but it could be tedious and it could not be connected to larger numerical programs. Fortunately, the bisection method is not difficult to understand and easy to implement. You may take a look to the introductory explanation to this numerical method in the book of Numerical Methods for Engineers by Steven Chapra and Raymond Canale.

Chapra and Canale also provide a step by step calculation to find roots by the bisection method and an algorithm easy to implement in any familiar software. However, as it usually occur in many scientific areas: there is a better and more effective way to do it. The algorthim proposed by Chapra and Canale is a valuable point of start but it can not deal with multiple roots and it is thought to search for only one root.

I found very interesting how Press, Vetterling and Teukolsky in the Numerical Recipes (NR) in F90 book manage to give another advantages to their bisection algorithm by allowing for bracketing roots with the zbrak routine. If you are a beginner I widely recommend you to first read the section devoted to the bisection method writen by Chapra an Canale.

For short, the bisection method takes advange of knowing a priori that there is a root in a given interval of values. The idea is to split that interval in two and look for a function sign change so that one could know in which half the root lies. If you repeat this process several times a very good precision could be achieved. I t is in fact the repetition of a bracketing process.

The bisection function (rtbis)

The fortran function rtbis is provided in the NR book. I choosed to store the function in the module created for this purpose Roof_Finding_module.f90 where zbrak was stored too. Modules nrtype_module.f90 and nrutil_module.f90 are needed as well. 

Next, I adapted the main program in F77 supplied by Press, Vetterling and Teukolsky in the Numerical Recipes Example Book to F90. Again, they use the zero order Bessel function as an example to test the program over a the range $x=1$ to $x=50$. The zbrak routine is then invoked to bracket the roots and the bisection method is used to refine the root under certain accuracy. Notice that the accuracy is not constant during the execution of the program but it changes with the size of the interval.

The code of the main program in F90 is as follows:

Bisection_Method_NR.f90
program Bisection_Method_NR
!Driver for rtbis function
    use root_finding_module
    use bessel_module
!
    integer, parameter:: n=100, nbmax=20
    real, parameter:: x1=1.0, x2=50.0
    integer i, nb
    real:: root, xacc
    real, dimension(:), pointer:: xb1, xb2
!
    nb=nbmax
    call zbrak(bessj0,x1,x2,n,xb1,xb2,nb)
    print "(/1x,a)", 'Roots of Bessel J0:'
    print "(/1x,t19,a,t31,a/)", 'x', 'F(x)'
!Loop for using the bisection function
    do i=1,nb
        xacc=(1.0e-6)*(xb1(i)+xb2(i))/2.0
        root=rtbis(bessj0,xb1(i),xb2(i),xacc)
        print "(1x,a,i2,2x,f12.6,e16.4)", &
        'Root',i,root,bessj0(root)
    end do
!
end program Bisection_Method_NR

The roots of $J_0(x)$ with rtbis

In order to test rtbis the roots of $J_0(x)$ from $x=1$ to $x=50$ are considered. Since the Bessel function $J_0$ is not provided I used one based on the intrinsic function Bessel_J0 in F90 and stored it in the bessel_module.90 as you may have noticed from the code shown above.

It is expected that the program outputs the crossings to the horizontal axis shown next

There 16 roots to be found. Then, after running Bisection_Method_NR.f90 it is obtained:

Finally, checked how good were the zeros in the output by diretly solving $J_0(x)=0$ in Maple. The agreement for every root was satisfactory.

Any questtion? Write in the comments and I shall try to help.

=========

Ildebrando