Solving a non linear equation using the secant method in Fortran

Question

Sooner than expected, I have yet another problem with my Fortran code. This time I'm trying to find the zeros of a function, which I can normally do, but now I have a function that has more than one input:

Where gamma, mp and k are just numbers; but u and mu are other values that I have to input for external arrays. In particular mu is another non-linear function of T, but I have an array of all mus for the temperatures I'm considering.

REAL*8 FUNCTION f(T, mu, u)
IMPLICIT NONE
REAL*8:: T, mu, u
REAL*8, PARAMETER:: mp = 1.67e-24, gammaa = 1.666666666d0, k_boltz = 1.380649e-16

f = T - ((gammaa - 1) * ((u * mu * mp) / k_boltz))
END FUNCTION f

The way I'm implementing the secant method is:

SUBROUTINE secant(f, x0, n_max, accuracy, mu, u, res)
  IMPLICIT NONE
  REAL*8, EXTERNAL:: f
  REAL*8, INTENT(in):: accuracy
  INTEGER, INTENT(in):: n_max
  REAL*8, INTENT(inout):: x0
  REAL*8, INTENT(out):: res(2)
  REAL*8:: x_next, dx, derivate, mu, u
  INTEGER:: counter
  
  derivate = (f(x0+0.001d0, mu, u) - f(x0, mu, u)) / 0.001d0
 
  counter = 0
  DO
    counter = counter + 1
    IF (counter > n_max) THEN
      PRINT*, "Too many iterations with no result."
      STOP
    END IF
    x_next = x0-f(x0, mu, u)/derivate                 
    dx = ABS(x_next - x0)
    IF (dx<accuracy) THEN
      res(1) = x_next
      res(2) = f(x_next, mu, u)
      STOP
    END IF
    x0 = x_next
  END DO
END SUBROUTINE secante

Which is the same way I implement it for functions of single input functions, only with the extra inputs now hard coded in. I thought this would work fine, but for some reason after calling:

t0 = 1.d0 ! arbitrary starting point
DO i = 1, n
  CALL secante(f, t0, 5, 0.00001d0, mu(i), int_energy(i), temperature(i))
ENDDO

I get a division by zero error in my terminal: "Note: The following floating-point exceptions are signalling: IEEE_DIVIDE_BY_ZERO"

I've tried changing all the values I can change on the fractions in the subroutine, but still no luck.

Any ideas?

Answer 1

FUNCTION f( T, mu, u )
   USE iso_fortran_env
   IMPLICIT NONE
   REAL(real64) f
   REAL(real64) T, mu, u
   REAL(real64), PARAMETER :: mp = 1.67e-24, gammaa = 5.0 / 3.0, k_boltz = 1.380649e-16

   f = T - ( gammaa - 1 ) * u * mu * mp / k_boltz

END FUNCTION f

!=======================================================================

SUBROUTINE secant( func, x0, n_max, accuracy, mu, u, res )
   USE iso_fortran_env
   IMPLICIT NONE
   REAL(real64), EXTERNAL :: func
   REAL(real64), INTENT(inout) :: x0
   INTEGER, INTENT(in) :: n_max
   REAL(real64), INTENT(in) :: accuracy, mu, u
   REAL(real64), INTENT(out) :: res(2)

   REAL(real64) :: x_next, derivative, delta
   INTEGER :: counter
  
   counter = 0
   delta = accuracy
   DO
      counter = counter + 1
      IF ( counter > n_max ) THEN
         PRINT*, "Too many iterations with no result."
         RETURN
      END IF
      derivative = ( func( x0 + delta, mu, u ) - func( x0, mu, u ) ) / delta
      x_next = x0 - func( x0, mu, u ) / derivative
      IF ( ABS( x_next - x0 ) < accuracy ) THEN
         res(1) = x_next
         res(2) = func( x_next, mu, u )
         RETURN
      END IF
      x0 = x_next
   END DO

END SUBROUTINE secant

!=======================================================================

PROGRAM main
   USE iso_fortran_env
   IMPLICIT NONE
   EXTERNAL secant
   REAL(real64), EXTERNAL :: f
   REAL(real64) t0, mu, u, accuracy
   REAL(real64) res(2)
   INTEGER nmx
   INTEGER i
   CHARACTER(len=*), PARAMETER :: fmt = "( A, '(', i0, ') = ', es11.4 )"


   t0 = 1.0
   mu = 1.0                               ! **** These should be functions,
   u = 1.0e12                             ! **** not constants
   accuracy = 1.0e-6
   nmx = 100

   CALL secant( f, t0, nmx, accuracy, mu, u, res )

   WRITE( *, fmt ) ( "res", i, res(i), i = 1, size( res ) )

END PROGRAM main

!=======================================================================

Answer 2

MODULE Properties
   USE iso_fortran_env
   IMPLICIT NONE

CONTAINS

   !----------------------------------------------

   REAL(real64) FUNCTION u( T )                            ! Returns the value of u(T) ... whatever that is (internal energy?)
      REAL(real64), INTENT(in) :: T
   
      u = 1.0e12                                           ! Put the definition of u as a function of T here
   
   END FUNCTION u

   !----------------------------------------------

   REAL(real64) FUNCTION mu( T )                           ! Returns the value of mu(T) ... whatever that is
      REAL(real64), INTENT(in) :: T
   
      mu = 1.0                                             ! Put the definition of mu as a function of T here
   
   END FUNCTION mu

   !----------------------------------------------

   REAL(real64) FUNCTION f( T )                            ! The function whose root is to be found
      REAL(real64), INTENT(in) :: T                        ! Thermodynamic temperature (presumably in Kelvin)
      REAL(real64), PARAMETER :: mp = 1.67e-24             ! Mass of a proton (funny units?)
      REAL(real64), PARAMETER :: gammaa = 5.0/3.0          ! Ratio of specific heats
      REAL(real64), PARAMETER :: k_boltz = 1.380649e-16    ! Boltzmann constant (also, funny units)

      f = T - ( gammaa - 1 ) * u(T) * mu(T) * mp / k_boltz
   
   END FUNCTION f
   
   !----------------------------------------------

END MODULE Properties

!=======================================================================

MODULE Solver                                              ! Routines for finding roots (only one at present)
   USE iso_fortran_env
   IMPLICIT NONE

CONTAINS
   SUBROUTINE secant( func, x, n_max, accuracy, res )      ! Attempts to solve func(x)=0 by local linear approximation
      REAL(real64), EXTERNAL :: func                       ! Function whose root is to be found
      REAL(real64), INTENT(inout) :: x                     ! Start (and end) value of x
      INTEGER, INTENT(in) :: n_max                         ! Maximum number of iterations allowed
      REAL(real64), INTENT(in) :: accuracy                 ! Tolerance for change in one iteration
      REAL(real64), INTENT(out) :: res(2)                  ! (Not really necessary) x and func(x) at end
   
      REAL(real64) :: x_next, derivative, delta
      INTEGER :: counter
     
      res = 0
      counter = 0
      delta = accuracy
      DO
         counter = counter + 1
         IF ( counter > n_max ) THEN
            PRINT*, "Too many iterations with no result."
            RETURN
         END IF
         derivative = ( func( x + delta ) - func( x ) ) / delta      ! Estimate of local slope
         x_next = x - func( x ) / derivative                         ! Linear update
         IF ( ABS( x_next - x ) < accuracy ) THEN                    ! Test for convergence
            res(1) = x_next
            res(2) = func( x_next )
            RETURN
         END IF
         x = x_next                                                  ! Update current position
      END DO
   
   END SUBROUTINE secant

END MODULE Solver

!=======================================================================

PROGRAM main
   USE iso_fortran_env
   USE Properties, ONLY: f
   USE Solver, ONLY: secant
   IMPLICIT NONE
   REAL(real64) t0, accuracy
   REAL(real64) res(2)
   INTEGER nmx
   INTEGER i
   CHARACTER(len=*), PARAMETER :: fmt = "( A, '(', i0, ') = ', es11.4 )"


   t0 = 273.0                                              ! Temperature in Kelvin - starting guess
   accuracy = 1.0e-6                                       ! Tolerance (in T) for convergence
   nmx = 100                                               ! Maximum allowed iterations

   CALL secant( f, t0, nmx, accuracy, res )

   WRITE( *, fmt ) ( "res", i, res(i), i = 1, size( res ) )

END PROGRAM main

!=======================================================================

Answer 3

The divide by zero could be related to your implementation of the secant method. Perhaps you can try with one of the solvers from the roots-fortran library (e.g. brent or TOMS748). An example is included below.

As a side observation, is it okay that you restart the solution for each temperature(i) at the previously found root t0? This could be another source of error. Also for the resulting scalar item temperature(i), it looks like it is getting associated with the output array res(2) of length 2, which is probably not what you want? (It might have gone unnoticed because it only overwrites memory beyond the current loop iteration.)

---

Proposed implementation with the roots-fortran package:

! nle_f.f90

!> Nonlinear equation for f(T)
module nle_f

implicit none
private

public :: wp
public :: solve_f, f

integer, parameter :: wp = kind(1.0d0)

contains

real(wp) function f(T, mu, u)
real(wp), intent(in) :: T, mu, u
real(wp), parameter :: mp = 1.67e-24_wp, &
                       gammaa = 1.666666666_wp, &
                       k_boltz = 1.380649e-16_wp

f = T - ((gammaa - 1) * ((u * mu * mp) / k_boltz))
end function f

subroutine solve_f(T,fT,TA,TB,nmax,atol,mu,u,iflag)

  use roots_module, only: root_scalar
    !> library of root solvers, available at:
    !> https://github.com/jacobwilliams/roots-fortran
    !> make sure to build with compatible precision kind

  real(wp), intent(out) :: T, fT
    !> root and function value results
  real(wp), intent(in) :: TA, TB
    !> root search interval, must satisfy f(TA)*f(TB) < 0
  integer, intent(in) :: nmax
    !> maximum number of iterations
  real(wp), intent(in) :: atol
    !> absolute tolerance for zero interval
  real(wp), intent(in) :: mu, u
    !> problem parameters
  integer, intent(out), optional :: iflag
    !> integer status flag, use for troubleshooting
  
  character(len=*), parameter :: method = 'toms748'

  call root_scalar('toms748',hf,TA,TB,T,fT,iflag, &
    atol = atol, &  ! you can specify other options too
    maxiter = nmax)

contains

  !> adapter function used to match procedure interface expected 
  !> by root_scalar
  !>
  !> Nota bene: mu and u accesible through host association, i.e
  !>            hf (the internal subprogram), has access to mu and u
  !>            in solve_f (the host)
  !>
  real(wp) function hf(T) 
    real(wp), intent(in) :: T
    hf = f(T,mu,u)
  end function

end subroutine

end module

Driver code:

! nlsolve.f90

program nlsolve

use nle_f
implicit none

! ... declarations &
!     initializations ...

! Specify search interval
ta = 10.0e3_wp
tb = 10.0e6_wp

do i = 1, n
   call solve_f(t0,ft0,ta,tb, &
      nmax=5,atol=0.00001_wp,mu=mu(i),u=int_energy(i), &
      iflag=iflag)

    if (iflag /= 0) then
      ! ... something went wrong, react appropriately ...
    end if
    temperature(i) = t0  ! store root, could be inlined in the call above

    ! optionally, adjust interval [ta,tb] based on knowledge of t0
    ! e.g. assuming the next root might be close ...

end do

! ... continue 

end program