Get trouble with bisectional method for solving equation

Hi everyone!
I want to solve f(x)=0 by using Bisection’s method. In this way, my function f(x) is any function such as f(x)=A.x^3+Bx^2+Cx+D. So I set an array (x,f(x)) then I use the subroutine having Bisectional method. My expectation is the solution having three roots. If the equation has 2 roots, I set one of the roots is 0.1 for my purpose. My interval is from 0 to 100.
The problem is all the equation I input the roots always appear 0.0 or no roots found.
Here is my code

``````program any_equation
implicit none
double precision A,B,C,D,E,dminroot,dmaxroot
double precision r,dr,r1,r2,r3,ff
integer i,nroots
double precision x(3000),f(3000),roots(10)
! input
print*, "Input  A,B,C,D,E:"
write(*,12) A,B,C,D,E
dr=0.1d0
! creating an array
open(20,file='array_x_fx',status='unknown')
do i=0,1000
r=dr*dfloat(i)
call find_f(A,B,C,D,r,ff)
x(i)=dr*dfloat(i)
f(i)=ff
write(20,*) x(i),f(i)
end do
close (20)
! find roots
call find_roots(roots,nroots,x,f,E)
!----
if (nroots ==0) then
write(*,*) 'no roots found'
stop
end if
!-----
if (nroots.ge.4) stop
!-----
if (nroots.eq.2) then
if(roots(1).lt.roots(2)) then
r2=roots(1)
r3=roots(2)
else
r2=roots(2)
r3=roots(1)
end if
r1=0.1d0
end if
!-----
if (nroots.eq.3) then
dmaxroot=roots(1)
dminroot=roots(1)
do i =1,nroots
if (dmaxroot.lt.roots(i))  dmaxroot = roots(i)
if (dminroot.gt.roots(i))  dminroot = roots(i)
enddo
do i =1,nroots
if((roots(i).lt.dmaxroot).and.(roots(i).gt.dminroot)) then
r2= roots(i)
end if
enddo
r1 =dminroot
r3 =dmaxroot
end if
write(*,131) r1,r2,r3
131     format('Solution are:',/,F7.2,/,F7.2,/,F7.2)

end program
! tinh f(x)
Subroutine find_f(A,B,C,D,r,f)
implicit none
double precision A,B,C,D,f,r
f=A*r**3.0d0+B*r**2.0d0+C*r+D
end subroutine timf

Subroutine find_roots(Roots,nroots,x,f,E)
!=====================================================================
! Multiple roots of equation f(x)=0 on [x1,x2] interval
! Method: Brute force with one of closed domain methods
! Close domain methods: bisectional or false position
! Alex G. January 2010
!---------------------------------------------------------------------
! input ...
! f   - function - evaluates f(x) for any x in [x1,x2]
! x1  - left endpoint of initial interval
! x2  - right endpoint of initial interval
! eps - desired uncertainity of the root as |b-a|<eps
! K   - number of subintervals for [x1,x2]
! output ...
! Root(n) - roots of the equation f(x)=0 on [x1,x2]
! nroots  - number of roots (nroots<=n)
! The program divide [x1,x2] into n subintervals
! Max number of iterations for every subinterval - 200
!======================================================================
implicit none
integer, parameter:: iter=200
integer nroots,i,j
double precision x1,x2,dx,e
double precision xmin,xmax,fmin,fmax,deriv3
double precision fa,fae,fb,fbe,fc,fce,froot,froote
double precision a,b,c,root,eps
double precision roots(10), y2(3000), x(3000), f(3000)
! initializayion
eps = 1.0e-6
x1=0.0d0
x2=100.0d0
dx = (x2-x1)/60. !the number of interval
nroots = 0
xmin=x(1)
xmax=x(1000)
fmin= deriv3(xmin,x,f,1000,1)
fmax= deriv3(xmax,x,f,1000,1)
call SPLINE(X,F,1000,fmin,fmax,y2)
! loop over subintervals
do j=1,60
a = x1  + real(j-1)*dx
b = a + dx
CALL SPLINT(x,f,y2,1000,a,fa)
fae=fa-e
CALL SPLINT(x,f,y2,1000,b,fb)
fbe=fb-e
!======================================================================
! check the closed domain condition f(a)*f(b)<0
if(fae*fbe>0) cycle
! Iterative refining the solution
do i=1,iter
c=(b+a)/2.0
! find c
CALL SPLINT(x,f,y2,1000,c,fc)
fce=fc-e
if(fae*fce.le.0.0) then
b=c
else
a=c
end if
!   condition(s) to stop iterations
if(abs(b-a)<= eps)  EXIT
end do
! check if it is a root or singularity
root = (b+a)/2.0d0
CALL SPLINT(x,f,y2,1000,root,froot)
froote=froot-e
if (abs(froote) < 1.0) then
nroots = nroots+1
Roots(nroots)=root
end if
end do
end subroutine find_roots
c------------------------------------------------------spline
SUBROUTINE SPLINE(x,y,n,yp1,ypn,y2)
INTEGER n,NMAX
REAL*8 yp1,ypn,x(n),y(n),y2(n)
PARAMETER (NMAX=20000)
INTEGER i,k
REAL*8 p,qn,sig,un,u(NMAX)
if (yp1.gt.0.99D40) then
y2(1)=0.0D0
u(1)=0.0D0
else
y2(1)=-0.50D0
u(1)=(3.0D0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
endif
do 11 i=2,n-1
sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
p=sig*y2(i-1)+2.0D0
y2(i)=(sig-1.0D0)/p
u(i)=(6.0D0*((y(i+1)-y(i))/(x(i+1)-x(i))-
& (y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
11    continue
if (ypn.gt.0.99D40) then
qn=0.0D0
un=0.0D0
else
qn=0.5D0
un=(3.0D0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
endif
y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.0D0)
do 12 k=n-1,1,-1
y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
return
END
c------------------------------------------------------splint
SUBROUTINE SPLINT(xa,ya,y2a,n,x,y)
INTEGER n
REAL*8 x,y,xa(n),y2a(n),ya(n)
INTEGER k,khi,klo
REAL*8 a,b,h
klo=1
khi=n
1     if (khi-klo.gt.1) then
k=(khi+klo)/2
if(xa(k).gt.x)then
khi=k
else
klo=k
endif
goto 1
endif
h=xa(khi)-xa(klo)
if(h.eq.0.00D0) stop
a=(xa(khi)-x)/h
b=(x-xa(klo))/h
y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*
&(h**2)/6.0D0
return
END
c------------------------------------------------------deriv3
function deriv3(xx, xi, yi, ni, m)
!====================================================================
! Evaluate first- or second-order derivatives
! using three-point Lagrange interpolation
! written by: Alex Godunov (October 2009)
!--------------------------------------------------------------------
! input ...
! xx    - the abscissa at which the interpolation is to be evaluated
! xi()  - the arrays of data abscissas
! yi()  - the arrays of data ordinates
! ni - size of the arrays xi() and yi()
! m  - order of a derivative (1 or 2)
! output ...
! deriv3  - interpolated value
!============================================================================*/
implicit none
integer, parameter :: n=3
double precision deriv3, xx
integer ni, m
double precision xi(ni), yi(ni)
double precision x(n), f(n)
integer i, j, k, ix
! exit if too high-order derivative was needed,
if (m > 2) then
deriv3 = 0.0
return
end if
! if x is ouside the xi(1)-xi(ni) interval set deriv3=0.0
if (xx < xi(1) .or. xx > xi(ni)) then
deriv3 = 0.0
return
end if
! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
i = 1
j = ni
do while (j > i+1)
k = (i+j)/2
if (xx < xi(k)) then
j = k
else
i = k
end if
end do
! shift i that will correspond to n-th order of interpolation
! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
i = i + 1 - n/2
! check boundaries: if i is ouside of the range [1, ... n] -> shift i
if (i < 1) i=1
if (i + n > ni) i=ni-n+1
! just wanted to use index i
ix = i
! initialization of f(n) and x(n)
do i=1,n
f(i) = yi(ix+i-1)
x(i) = xi(ix+i-1)
end do
! calculate the first-order derivative using Lagrange interpolation
if (m == 1) then
deriv3 =(2.0*xx - (x(2)+x(3)))*f(1)/((x(1)-x(2))*(x(1)-x(3)))
deriv3 = deriv3 + (2.0*xx - (x(1)+x(3)))*f(2)/((x(2)-
1  x(1))*(x(2)-x(3)))
deriv3 = deriv3 + (2.0*xx - (x(1)+x(2)))*f(3)/((x(3)-
2  x(1))*(x(3)-x(2)))
! calculate the second-order derivative using Lagrange interpolation
else
deriv3 =          2.0*f(1)/((x(1)-x(2))*(x(1)-x(3)))
deriv3 = deriv3 + 2.0*f(2)/((x(2)-x(1))*(x(2)-x(3)))
deriv3 = deriv3 + 2.0*f(3)/((x(3)-x(1))*(x(3)-x(2)))
end if
end function deriv3

!-------------------------------------------------------
Function gauss8(xx,fx,dx,a,b)
!==========================================================
! Integration of f(x) on [a,b]
! Method: Gauss 8 points
! written by: Alex Godunov (October 2009)
!----------------------------------------------------------
! IN:
! f   - Function to integrate (supplied by a user)
! a	  - Lower limit of integration
! b	  - Upper limit of integration
! OUT:
! gauss8 - Result of integration
!==========================================================
implicit none
integer, parameter :: n=4
double precision gauss8, a, b, mi,fmi,mti,fmti
double precision ti(n), ci(n)
data ti/0.1834346424,0.5255324099,0.7966664774,0.9602898564/
data ci/0.3626837833,0.3137066458,0.2223810344,0.1012285362/
double precision r, m, c
integer i,dx
double precision xmin,xmax,fmin,fmax,deriv3
double precision xx(1000),fx(1000),y2(1000)
xmin=xx(1)
xmax=xx(dx)
fmin=deriv3(xmin,xx,fx,dx,1)
fmax=deriv3(xmax,xx,fx,dx,1)
call SPLINE(xx,fx,dx,fmin,fmax,y2)
r = 0.0;
m = (b-a)/2.0;
c = (b+a)/2.0;
do i = 1,n
mi= m*(-1.0)*ti(i) + c
call SPLINT(xx,fx,y2,dx,mi,fmi)
mti=m*ti(i) + c
call SPLINT(xx,fx,y2,dx,mti,fmti)
r = r + ci(i)*(fmi + fmti)
end do
gauss8 = r*m
return
end function gauss8

``````

I do not have a solution for you, at least not right away. But just a small thing: in the function evalaution you use real exponents. There is no need for that:

`````` f=A*r**3+B*r**2+C*r+D
``````

works fine too and possibly a bit faster (though in this context it will not be measurable )

The code you supply is pretty long. Have you seen where the roots probably lie from the table you print?

You mention a bisection method, but I see that the function is first tabulated and then splines are used for interpolation. That is not the usual way to go about it.

One way to understand what is going on, is to print some intermediate results or to use a debugger. The code is a trifle lengthy to see what is going wrong by straightforward inspection.

1 Like

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.

4 Likes

At first, I think the way to solve any function f(x)=0 can interpret to this way:

1. I create a tabular including x and f(x)
2. I infer these values from the table to a function by using SPLINE and SPLINE.
3. Then I use some models of solving my equation to finish my work.

So, it doesn’t seem to work. Thank you for pointing out my faults.

Thank you for giving me a lot of details. I learned many things from your advice. I read the Fortran recipe then I came up with this idea to use SPLINE and SPLINT, and forgot the condition to use them. And many thanks for your source code to find the solution by using the bisectional method. Finally, I want to find some roots in a given interval, not all the roots. So I need to study more about your code.
Many thanks! I appreciate your help in this issue.

There is nothing wrong with your idea. You can do that if you don’t have an analytic expression for the function, but only its value at specific points `[xi,f(xi)], i=0,1,2,...,n` - typically data from an experiment, or data from another numerical method. And even then, you do some kind of interpolation, whichever is appropriate for the function at hand (not necessarily splines). Cubic splines work well for that purpose provided the function “behaves”, at least in the interval you are interested (continuous and differentiable function with continuous first and second derivatives).

In case you can use cubic splines, keep in mind the so called “natural” cubic splines are far from being natural by any means, unless you know in advance that the second derivative is zero at boundary points - which is rarely the case. Why natural cubic splines are so popular or even named like that, I have no idea, but I know they should not be used except in special cases (which are usually staged examples).
If you don’t have any additional information, you should use “not-a-knot” cubic spline interpolation instead (see, e.g., this for a short description and an algorithm concerning not-a-knot boundary constraints). John Burkardt’s spline.f90 can be used, but since it covers much more than just cubic splines I would go for a simpler code just for the task, preferably a custom code (you learn much more about a method when you have to write its Fortran implementation). If nothing else, I could upload my module somewhere, which I wrote years ago and uses not-a-knot by default (but you can optionally use other constraints as well).

1 Like

You don’t have to roll your own root finder. Use a library:

Rarely is bisection what you want. Try one of the other faster methods.

4 Likes

And for splines there’s a library for that too:

Both of these use the Fortran Package Manager so it’s trivial to add them to your code.

4 Likes

Look at @haidangwy’s strategy, for example: calculate tabulated data, apply spine interpolation, then use the resulting interpolating function to find its roots. It is actually a plausible strategy (provided the required conditions are there). He probably ended up to this on his own, and I am guessing he wouldn’t if he used a ready-made solution for the whole problem.
In my humble opinion, being spoiled with ready-made solutions while you are at the learning stage is a big mistake. It makes things easier today, but makes things harder tomorrow since you get little to no experience when you use “black boxes” all the time. Interpolation and root finding are very good starting points to “get your hands dirty” and concoct a custom solution. Later on, using existing packages for more complicated interpolation will certainly help, but in the meantime trying a custom, home-made solution makes you knowing what you are doing. I have seen many people using packages such as Matlab, Scilab, Octave all the time, and I realized they use numerical methods without really knowing what they are doing, because they didn’t have to - they are just addicted to ready-made libraries, using them as black boxes for just about anything.

2 Likes

Can you give us some references for your bisection simplified method? How is it different from the standard bisection method?
Also, one strategy you can adopt to find all roots within a given interval is deflation. It is usually employed with the muller method. The routine dzanly from the IMSL library employs it. It is described in Conte & de Boor’s book, Elementary Numerical Analysis, sec. 3.7.

Thank your for your sharing, I’m going to study the way you work with the function and your idea to solve it. I really appreciate your support.

There are many ways to simplify Bisection. The example code above uses (and does the necessary safekeeping) with signs instead of real numbers, and uses subinterval length instead of both its boundaries. But there are other ways to simplify Bisection (see, e.g., this paper). I use the term “simplified” because that’s how it was called the first time I’ve read about it, years ago; I have seen it as “improved” or “enhanced” Bisection too. As far I am aware of, re-inventing Bisection started with Bisection is optimal, a paper with a rather surprising (at first glance) title.

I’m sure there are others too, no doubt about it. Implementations like zanly (dzanly), zbren, zreal, etc are sometimes mentioned in the bibliography (see, e,g, R. L. Burden and J.D. Faires, Numerical Analysis - p. 98 in my edition). However, I never had full access to IMSL and I can’t say much about those, let alone recommend them. I also have to admit I never liked Muller’s method because it requires three initial guesses of the root, although I am aware there are variations of the method promising they can deal with poor initial guesses.