Your code has several issues, including working with specific polynomials only (presumably cubic spline polynomials), mixing old FORTRAN with modern Fortran, and, most importantly, it’s “vanilla” Bisection which is not recommended (unless it’s a specific college projects, or an introduction to numerical root finding). For any “real-world” project, I strongly recommend Bisection Simplified instead (sometimes called Improved Bisection).
The module below implements Bisection Simplified in Fortran 90 and should work for any function. I wrote it many years ago (in fact this version is strictly Fortran 90 and it’s further simplified to avoid using other modules in order to save space). In most cases, Bisection Simplified works as fast, if not faster, than the popular Newton method, but with the added advantage that it’s still Bisection, meaning it guarantees convergence. Note, however, that Bisection won’t find any root at all if the interval contains an even number of roots.
module Bisection_Simplified
implicit none
private
integer, parameter, public :: single_precision=selected_real_kind(6)
integer, parameter, public :: double_precision=selected_real_kind(14)
integer, parameter, public :: acc=double_precision
public :: Bisect_Root_Simplified
contains
!-------------------------------------------------------------------------------
subroutine Bisect_Root_Simplified(f, x_left_initial, x_right_initial, root,&
estimated_error, accuracy_tolerance,&
max_iterations)
! Returns the estimated root of the equation f(x)=0 that lies in the interval
! [x_left_initial,x_right_initial], and, optionally, the estimated error of the
! root found.
! The intent(in) argument accuracy_tolerance is optional, representing the
! desired accuracy. The default value is accuracy_tolerance=2*epsilon(0.0);
! typically, this is equal to 2.38418579e-7.
! The intent(in) argument max_iterations is also optional, representing the
! maximum number of iterations. The default value is max_iterations=40.
interface
function f(x) result(fx)
import :: acc
real(kind=acc), intent(in) :: x
real(kind=acc) :: fx
end function f
end interface
real(kind=acc), intent(in) :: x_left_initial, x_right_initial
real(kind=acc), intent(out) :: root
real(kind=acc), intent(out), optional :: estimated_error
real(kind=acc), intent(in), optional :: accuracy_tolerance
integer, intent(in), optional :: max_iterations
real(kind=acc) :: eps, half_interval
integer :: iterations_limit, i
real(kind=acc) :: x_previous, x_current, f_current, sign_f_left, sign_f_previous
real(kind=acc), parameter :: zero=2.*epsilon(0._acc)
! Setting up accuracy:
if (present(accuracy_tolerance)) then
if (accuracy_tolerance > 0._acc) then
eps=accuracy_tolerance
else
stop "Bisect_Root_Simplified: Invalid accuracy_tolerance argument."
end if
else
eps=2.*real(epsilon(0.0),acc)
end if
! Setting up maximum number of iterations:
if (present(max_iterations)) then
if (max_iterations >= 2) then
iterations_limit=max_iterations
else
stop "Bisect_Root_Simplified: Invalid max_iterations argument."
end if
else
iterations_limit=40
end if
! Performing interval arguments validity check:
if (x_left_initial >= x_right_initial) stop "Bisect_Root: Invalid initial interval."
! Initializing interval boundary points:
x_previous=x_left_initial
x_current=x_left_initial
sign_f_left=sign(1.0_acc,f(x_left_initial))
sign_f_previous=sign_f_left
half_interval=x_right_initial-x_left_initial
! Performing interval validity check:
if (sign_f_left*sign(1.0_acc,f(x_right_initial)) > 0.0_acc) &
stop "Bisect_Root_Simplified: Initial interval does not contain a root, or contains an even number of roots."
!
do i=1,iterations_limit
! Calculating the mid-point:
half_interval=0.5_acc*half_interval
x_current=x_previous+sign_f_left*sign_f_previous*half_interval
! Checking if desired accuracy tolerance has been achieved:
if (half_interval <= eps) then
root=x_current; if (present(estimated_error)) estimated_error=half_interval
return
end if
f_current=f(x_current)
! Checking the (remote) possibility that f is exactly zero at the middle:
if (abs(f_current) <= zero) then
root=x_current; if (present(estimated_error)) estimated_error=0.0_acc
return
! Preparing next iteration:
else
x_previous=x_current; sign_f_previous=sign(1.0_acc,f_current)
end if
end do
! Normally, the subroutine should return earlier, never reaching this point.
! So return the root, but inform the user that accuracy has not been reached:
root=x_current
if (present(estimated_error)) estimated_error=half_interval
print "(a)","WARNING: Bisect_Root_Simplified: Maximum number of iterations reached,"
print "(a)"," but the accuracy criterion is not statisfied."
end subroutine Bisect_Root_Simplified
!-------------------------------------------------------------------------------
end module Bisection_Simplified
In general, the maximum number of iterations should never be reached (if it does, then most probably something is wrong with the function given). Bisection Simplified should converge after a few iterations, typically less than ten. A simple driver program to test the module could be:
program Bisection_Simplified_driver
use Bisection_Simplified
implicit none
real(kind=acc) :: root
real(kind=acc) :: estimated_error
call Bisect_Root_Simplified(f, -10.0_acc, 0.0_acc, root, estimated_error)
print "(a,f10.7)"," Estimated root:",root
print "(a,es14.7)","Estimated error:",estimated_error
contains
function f(x) result(fx)
! Defines the left-hand side of the equation to be solved.
real(kind=acc), intent(in) :: x
real(kind=acc) :: fx
fx=x+exp(x)
end function
end program Bisection_Simplified_driver
Last but not least, note that any “normal” numerical method will only find one root within a given interval, even if there are many. If you want to find all the roots within a given interval, that’s a completely different task and way more complicated; though doable, it requires special packages designed for that purpose, or a method I mentioned elsewhere. I would not recommend a brute force method such as the one you use in your code.