modifying secant method algorithm

Question

my code below uses the secant method to find the root of an analytic function. The analytic function, f must be specified in the function part of my code. The code below works well and has no compilation errors. However, for the problem I want to solve I do not know the analytic function f.

Instead I calculate the function numerically, and its stored as an array. I want now apply my code to find the roots of this function. So how can I modify my code such that the input is not an analytic function, instead just an array which I have already calculated?

My working code is below, I assume I just need to modify the last part where I call the function f, I just am unsure how to go about doing this. Thanks!

program main
  implicit none 
  real :: a = 1.0, b = -1.0
  integer :: m = 8
  interface
  function f(x)
  real, intent(in) :: x
  end function
  end interface

  call secant(f,a,b,m)
  end program main

  subroutine secant(f,a,b,m)
  implicit none
  real, intent(in out) :: a,b
  integer, intent(in) :: m
  real :: fa, fb, temp
  integer :: n
  interface
  function f(x)
  real, intent(in) :: x
  end function f
  end interface

  fa = f(a)
  fb = f(b)
  if (abs(fa) >  abs(fb)) then
     temp = a
     a = b
     b = temp
     temp = fa
     fa = fb
     fb = temp
  end if
  print *,"    n        x(n)         f(x(n))"
  print *," 0 ", a, fa    
  print *," 1 ", b, fb
  do n = 2,m
     if (abs(fa) >  abs(fb)) then
        temp = a
        a = b
        b = temp
        temp = fa
        fa = fb
        fb = temp
     end if
     temp = (b - a)/(fb - fa)
       b = a
     fb = fa
     a = a - fa*temp
     fa = f(a)
     print *,n,a,fa
  end do   
  end subroutine secant

  real function f(x)
  implicit none 
  real, intent(in) :: x
  f = x**5 + x**3 + 3.0 !analytic form of a function, I don't actually have this though, I just have the function stored as an array
  end function f

Answer 1

What I wanted to say in my comments are something as below.

You can modify your secant subroutine to take an object of an abstract class (FAZ) which is guaranteed to have a function f. For example, as following.

solver.f90

!*****************************************************************
MODULE solver
!*****************************************************************
  IMPLICIT NONE
  PRIVATE

  PUBLIC FAZ
  PUBLIC secant

  TYPE, ABSTRACT :: FAZ
   CONTAINS
     PROCEDURE(f),      deferred, pass :: f
  END TYPE FAZ

  ABSTRACT INTERFACE
     FUNCTION f(this, x)
       IMPORT  :: FAZ
       REAL                           :: f
       CLASS(FAZ),         INTENT(IN) :: this
       REAL,               INTENT(IN) :: x
     END FUNCTION f
  END INTERFACE

!=====================================================================
CONTAINS 
!=====================================================================


  subroutine secant(oFAZ,a,b,m)
    CLASS(FAZ)         :: oFAZ
    real, intent(in out) :: a,b
    integer, intent(in) :: m
    real :: fa, fb, temp
    integer :: n

    fa = oFAZ%f(a)
    fb = oFAZ%f(b)
    if (abs(fa) >  abs(fb)) then
       temp = a
       a = b
       b = temp
       temp = fa
       fa = fb
       fb = temp
    end if
    print *,"    n        x(n)         f(x(n))"
    print *," 0 ", a, fa    
    print *," 1 ", b, fb
    do n = 2,m
       if (abs(fa) >  abs(fb)) then
          temp = a
          a = b
          b = temp
          temp = fa
          fa = fb
          fb = temp
       end if
       temp = (b - a)/(fb - fa)
       b = a
       fb = fa
       a = a - fa*temp
       fa = oFAZ%f(a)
       print *,n,a,fa
    end do
  end subroutine secant

END MODULE solver

You can then implement the behavior of the function f in whatever way you like by extending the abstract class FAZ to a concrete class MyFAZ. For example, I wrote it as following.

myfaz.f90

!*******************************************************************
MODULE my_concrete_faz
!*******************************************************************
  USE solver, ONLY : FAZ
  IMPLICIT NONE
  PRIVATE

  PUBLIC MyFAZ
  PUBLIC MyFAZ_constructor

  TYPE, EXTENDS(FAZ) :: MyFAZ
     PRIVATE
     REAL,     DIMENSION(:),   ALLOCATABLE :: xdata, fdata
   CONTAINS
     PROCEDURE :: destructor
     PROCEDURE :: f
  END TYPE MyFAZ

! ================================================================
CONTAINS
! ================================================================

  ! ****************************************************************
  FUNCTION MyFAZ_constructor(xdata_arg, fdata_arg) RESULT(oMyFAZ)
  ! ****************************************************************
    TYPE(MyFAZ)                         :: oMyFAZ
    REAL,     DIMENSION(:), INTENT(IN)  :: xdata_arg, fdata_arg
    INTEGER :: ndata, jj

    ndata = size(xdata_arg)
    if (size(fdata_arg) /= ndata) then
       stop 'MyFAZ_constructor: array size mismatch .. ndata'
    end if
    do jj=1,ndata-1
       if (xdata_arg(jj)>xdata_arg(jj+1)) then
          stop 'MyFAZ_constructor: expecting a sorted xdata. I am lazy.'
       end if
    end do
    allocate(oMyFAZ%xdata(ndata))
    allocate(oMyFAZ%fdata(ndata))
    oMyFAZ%xdata = xdata_arg
    oMyFAZ%fdata = fdata_arg

  END FUNCTION MyFAZ_constructor


  ! ****************************************************************
  SUBROUTINE destructor(this)
  ! ****************************************************************
    CLASS(MyFAZ),           INTENT(INOUT) :: this

    deallocate(this%xdata)
    deallocate(this%fdata)

  END SUBROUTINE destructor


  ! ****************************************************************
  FUNCTION f(this, x)
  ! ****************************************************************
    ! evaluates the function.  
    ! Linear interpolation is used here, but this will not make sense
    ! in actual application. Everything is written in a very inefficient way. 
    REAL                         :: f
    CLASS(MyFAZ),     INTENT(IN) :: this
    REAL,             INTENT(IN) :: x
    !
    INTEGER  :: jj
    REAL     :: rr

    do jj=1, size(this%xdata)-1
       if (this%xdata(jj)<=x .and. x<=this%xdata(jj+1)) then
          exit
       end if
    end do

    rr = (this%fdata(jj+1) - this%fdata(jj))/(this%xdata(jj+1) - this%xdata(jj))
    f  = rr*(x - this%xdata(jj)) + this%fdata(jj)

  END FUNCTION f

END MODULE my_concrete_faz

I used the linear interpolation, just for demonstration. Actually, if f(x) = r x + s, then you know the solution without using the secant method.
You will have your own appropriate method to evaluate f(x) between data points.

You can use the above two modules as following.

main.f90

PROGRAM demo
  USE solver, ONLY : secant
  USE my_concrete_faz, ONLY : MyFAZ, MyFAZ_constructor
  IMPLICIT NONE

  REAL,  DIMENSION(:), ALLOCATABLE :: xdata, fdata
  INTEGER                          :: ndata
  INTEGER                          :: niter_max
  REAL                             :: xa, xb
  TYPE(MyFAZ)                      :: oMyFAZ

  niter_max = 10
  xa = -2.0
  xb =  3.0

  ! prepare data
  ndata = 4
  allocate(xdata(ndata))
  allocate(fdata(ndata))

  xdata(1) = -3.0
  xdata(2) = -1.1
  xdata(3) =  1.2
  xdata(4) =  3.8

  fdata(1) = -1.5
  fdata(2) = -0.9
  fdata(3) =  0.1
  fdata(4) =  0.8

  ! prepare the function
  oMyFAZ = MyFAZ_constructor(xdata, fdata) 
  deallocate(xdata)
  deallocate(fdata)

  ! solve
  call secant(oMyFAZ,xa,xb,niter_max)


  write(*,*) '**************'
  write(*,*) 'normal end'
  write(*,*) '**************'

END PROGRAM demo

I compiled, built, and got output as following.

$ ifort -c solver.f90
$ ifort -c myfaz.f90
$ ifort -c main.f90
$ ifort -o demo *.o 
$ ./demo 
     n        x(n)         f(x(n))
  0    3.000000      0.5846154    
  1   -2.000000      -1.184211    
           2   1.347448      0.1396975    
           3  0.8285716     -6.1490655E-02
           4  0.9871597      7.4606538E-03
           5  0.9700001      0.0000000E+00
           6  0.9700001      0.0000000E+00
           7            NaN            NaN
           8            NaN            NaN
           9            NaN            NaN
          10            NaN            NaN
 **************
 normal end
 **************
$ 

The NaNs are there because your secant subroutine reached to the solution before the maximum iteration, but had no way to exit in the middle of the loop.

Here is a plot of the data.