An improvement for bracketing roots (zbrak) - Numerical recipes in F90

The zbrak routine

This routine supplied in the Numerical Recipes (NR) book is an improvement to the zbrac routine. In the output of zbrac to bracket the roots of the zero order Bessel function of the first kind it can be seen that some roots are bracketed twice which seems to be a waste of effort.

Fortunately, our friends Press, Teukolsky, Vetterling and Flannery proposed some improvements to zbrac through zbrak, which includes more advanced fortrand instructions. In order to check how zbrak works I included it in the previously created Root_Finding_module.f90 so that it could easily used in the new main program.

Also, zbrak requires handling data through arrays another computation is needed to create an array. So the function arth_i was added to nrutil_module.f90 previously created for this porpuse. Notice that zbrak as proposed in the NR book uses the function arth but nrutil has three different options for this one. As you may see, arth_i is the one fitting our needs for integer increments of the arithmetic progression as the array. I just add "_i" to the arth in the zbrak code.

Since for this case the problem function/equation is the zero order Bessel function of the first kind I did use the same bessel_module.f90 as for zbrac.

The main program for zbrak in F90

I took the main program for zbrak in F77 and adapted it to F90. The code is as follows:

bracketing_k.f90

program bracketing_k
!Driver for routine zbrak in F90
    use root_finding_module
    use bessel_module
!
    integer, parameter:: n=100, nbmax=20
!Update here the x-range over which roots shall be bracketed
    real, parameter:: x1=1.0, x2=50.0
    integer:: i, nb
    real, dimension(:), pointer:: xb1, xb2
!
    nb=nbmax
    call zbrak(bessj0,x1,x2,n,xb1,xb2,nb)
    print "(/1x,a/)", 'Brackets for roots of BessJ0:'
    print "(/1x,t17,a,t27,a,t40,a,t50,a/)", &
    'lower', 'upper', 'F(lower)', 'F(upper)'
!
    do i=1,nb
        print "(1x,a,i2,2(4x,2f10.4))", 'Root', &
        i,xb1(i),xb2(i),bessj0(xb1(i)),bessj0(xb2(i))
    end do
!
end program bracketing_k

After running the program the following output is obtained:

You may compare the lower and upper limits bracketing the roots and the corresponding values of J0 at those values of x with crossings shown in the plot below.

Notice that 16 crossings shown in the plot correspond to the 16 bracketed roots in the program output.

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

 =========

Ildebrando