Get trouble with bisectional method for solving equation

Help
1

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:"
        read *, A,B,C,D,E
        write(*,12) A,B,C,D,E
12      FORMAT('Your equation: ',F7.2,'X^3+',F7.2,'X^2+',F7.2,'X+',F7.2,'=',F7.2)
        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)
! Comments: 
! 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
2

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 :slight_smile:)

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
3

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
4

I really appreciate your advice. Thank you for your help.
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.

5

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.

6

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
7

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
8

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
9

Of course he can do that. There are many options when it comes to ready-made libraries, and sometimes this is the way to go. But since neither root finding or simple cubic splines are complicated topics, I opted to recommend a home-made custom program. For me at least, going that way whenever I could was the best decision I ever made. Writing your own program to solve a problem is not just about programming. You need to read books and papers about the numerical methods involved, then roll up sleeves and implement the methods in Fortran. It is a very educative process. Way slower, sure, you will make a lot of mistakes that you will need to correct, no doubt about it; but it will prove very useful in the long run. You learn a lot, both about the numerical methods you use and about programming. You gain experience for your next, more complicated project, and that’s an invaluable asset.

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
10

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.

11

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.

12

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.