Bracketing roots (zbrac) - Numerical Recipes in F90

 

About bracketing roots

This is an important tool to search for possible roots in the real or complex plane. This is, if you know where the root is located you may use this information to narrow the search and/or guide the program by supplying a better guess. It is common practice to do this by plotting the function or equation in question to investigate how many roots and where are these located. In fact, if you have a good software with plotting capabilities it would be a good tool for most cases.

Take for instance the case of the following equation: 

$4x\sin(2x)=0$

can be plotted to look for possible roots between $x=1$ and $x=10$, for example. If we plot it in Maple the following is obtained:

You can see in the plot above the several crossings in the horizontal axis or roots. At a first inspection you would not probably realize the existence of such number of roots. Fortunately, the guys of Numerical Recipes faced this problem decades ago and created a routine capable of bracketing the roots.

The zbrac routine in F90

You may find this routine in Chapter B9 in the NR book. However, you could need a main program to make use of the zbrac routine. Two options appear at this point: you may write your own program or make use of the Numerical Recipes Example Book (NREB). I took the shortest way and take advantage of the main example program for this routine.

Another obstacle that you would need to overcome is that some modules could be needed in order to see what this routine can do. You will need:

1. the nrtype module. You can find it in the appendix C1 of the NR book,
2. the nrutil module. You can find it in the appendix C2 of the NR book,
3. an extra module, that I created to store zbrac, so that the routine could be called into the main program
4. an extra module, that I cretaed to store a fortran function of the function/equation to be studied, so that it could be used in the main program.

If you have the NR book but not the code in electronic format I warned you. You will need to type the whole nrtype module but, for this case, only a part of the nrutil module, which is quite large. From nrutil you will only need the subroutine nrerror.

Modules, in points 3 and 4, listed above are not mentioned either in NR nor in NREB. Here is a short explanation of this:

A module for zbrac

As you can see in the driver program proposed in the NREB for zbrac the subroutine is called and the compiler should find either in the same file for the main program or in a separate file. I preferred the second option to get a more easy-to-review for errors program. th code for this module is as follows

root_finding_module.f90
module root_finding_module
!
    public:: zbrac
!
    contains
!
    subroutine zbrac(func,x1,x2,succes)
    ...
    end subroutine zbrac
!
end module root_finding_module  

In this way, all you have to do is to add a line code to connect this module with the main program as use root_finding_module.

A module for the function/equation to be studied

As you may have noticed in the codes supplied by NR and NREB no function was defined at all. However, for the test program proposed in NREB the Bessel functions of the first kind and order zero are considered since it is known that this functions has several zeros. You can see this graphically as follows:

where several crossings of the horizontal axis can be appreciated.

We know all this because in the NREB a BessJ0 function is used. Then, it should be supplied some how and again I used a module to do this. This is the one i used:

bessel_module.f90
module bessel_module
!
    public:: bessj0
!
    contains
!
    function bessj0(x) result(bessj0_result)
    intrinsic:: bessel_j0
    real, intent(in):: x
    real:: bessj0_result
!
    bessj0_result=bessel_j0(x)
    end function bessj0
!
end module bessel_module

where as you can see I used an intrinsic function of fortran: bessel_j0. In this way, and as before, all you have to do is to add a line code to connect this further module with the main program as use bessel_module.

Summarizing, the beginning of the driver program is just slightly changed as:

Bracketing_c.f90

program Bracketing_c

!Driver program for zbrac routine

    use root_finding_module

    use bessel_module

!

    logical:: succes

    integer:: i

    real:: x1, x2

    write(*,'(1x,t4,a,t29,a/)') 'Bracketing values:', &

    'Function values:'

    write(*,'(1x,t6,a,t16,a,t29,a,t41,a/)') x1, x2, &

    'BessJ0(x1)', 'BessJ0(x2)'

!Change the limit ofr i to expand the search for roots

    do i=1,10

        x1=i

        x2=x1+1.0

        call zbrac(bessj0,x1,x2,succes)

        if (succes) then

            write(*,'(1x,f7.2,f10.2,7x,2f12.6)') x1, x2, &

            bessj0(x1),bessj0(x2)

        end if

    end do

end program Bracketing_c

I changed the name of the main program to look more familiar.

Bracketing the Bessel function J0 between x=1 and x=12 with zbrac

Then, all I did was run Bracketing_c.f90 in Code::Blocks from Fortran Tools. It did what it was expected:

Notice that the change in sign in the columns for BessJ0(x1) and BessJ0(x2) correspond with the previous plot.

Bracketing a random function between x=1 and x=12 with zbrac

I tried another function to test the zbrac routine. It was the one shown at the beginning:

$4x\sin(2x)=0$

and as you could see above, it has several roots. Some changes to the fortran program had to be done.

Since the problem function/equation was different I created another module to stored it. This is it:

my_function_module.f90

module my_function_module

    public:: mfunc0

!

    contains

!

    function mfunc0(x) result(mfunc0_result)

        intrinsic:: sin

        real, intent(in):: x

        real:: mfunc0_result

!

        mfunc0_result=4*x*sin(2*x)

    end function mfunc0

!

end module my_function_module

and the main program had to be changed too as follows:

Bracketing_mine.f90

program Bracketing_mine

!Driver program for zbrac routine

    use root_finding_module

    use my_function_module

!

    logical:: succes

    integer:: i

    real:: x1, x2

    write(*,'(1x,t4,a,t29,a/)') 'Bracketing values:', &

    'Function values:'

    write(*,'(1x,t6,a,t16,a,t29,a,t41,a/)') x1, x2, &

    'f(x1)', 'f(x2)'

    do i=1,10

        x1=i

        x2=x1+1.0

        call zbrac(mfunc0,x1,x2,succes)

        if (succes) then

            write(*,'(1x,f7.2,f10.2,7x,2f12.6)') x1, x2, &

            mfunc0(x1),mfunc0(x2)

        end if

    end do

!

end program Bracketing_mine

After compilation and running of Bracketing_mine.f90 you get

where it is shown that the program correctly gives the two values of $x$ that bracket the roots.

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

=========

Ildebrando