ChemEng stuff followers

The false position method (rtfls) - Numerical recipes in F90

 This page uses has very similar structure to the one for the bisection method rtbis, zbrac, and zbrak pages. So, I recommend the reader to find those pages in this blog in case of  doubt.

A short perspective of the false position method

The False Position method as the bisection is root finding tool that depends on knowing a priori the location of the root. However, this method sems to be better than the bisection one. For many problems, convergence and computational time would not be relevant to choose between the bisection 
and the false position method.

Here, the routine proposed in the Numerical Recipes in F90 book (NR) is being used along with the main program proposed in the Numerical Recipes example book written by the same authors.

The test with $J_0(x)$

For this case the false position method is implemented in F90 by Press, Vetterling and Teukolsky in their NR book as a function. I just added this function into the root_Finding_module.f90 cretaed for this purpose so that it coul be called from the main program.

Another, subroutine is also needed for the rtflsp function. This is maned swap and its purpose is just exchange the two values assigned to two parameters. Since, these two parameters are real I choosed swap_rswap_r was also stored into the nrutil_module.f90 too. I finally, adjust the main program code, in F77 supplied by Press, Vetterling and Teukolsky, for this case to work in F90. The code is as follows:

False_Position_rtflsp.f90
program False_Position_rtflsp
 !Driver for rtflsp 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 False_Position_rtflsp

It runs smoothly and gives the following output




which is in very good agreement with the roots for the zero order Bessel function of the first kind. Convergence is good too and these results are barely the same obtained with the bisection method.

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

=========
Ildebrando

No comments:

Post a Comment

Most popular posts