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?

The best thing you can do here is take a look at what's actually happening inside your code, either with a debugger or with a bunch of print statements. Then you can identify exactly where the divide by zero happens. — veryreverie, Sep 06 '22 at 18:48
What values of T are you looking at, roughly? If you replace your numerical derivative with a central difference formula does it help? But really, without a complete program that shows the issue it may well be difficult to help — Ian Bush, Sep 06 '22 at 19:02
Note also real*8 is not standard Fortran, has never been standard Fortran, might not be supported by a compiler and might not be doing what you want. — Ian Bush, Sep 06 '22 at 19:03
I'm looking at a range of temperatures from 10,000 K to 10,000,000 K. I've tried setting the derivative value to 1 to see if that was misbehaving, but it still ran into that error. Also i'm using real*8 because for some reason that's the standard in my uni labs compilers — Marco Leonardi, Sep 06 '22 at 19:08
*"i'm using real*8 because for some reason that's the standard in my uni labs compilers"* is a very strange expression. It can be standard in your lab's coding conventions - and then I despise them, but will hardly be standard in any compiler. Definitely see https://stackoverflow.com/questions/838310/fortran-90-kind-parameter https://stackoverflow.com/a/10521326/721644 and https://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4?noredirect=1&lq=1 — Vladimir F Героям слава, Sep 06 '22 at 20:12
I don’t know, I only started using fortran like 3 weeks ago and I’m not an experienced programmer to begin with . But I do hear your advice and will drop the real8 from future projects :) — Marco Leonardi, Sep 06 '22 at 21:17
This in the code is not the secant method, but a simplified Newton method. I do not see the re-computation of mu. As you say it is represented by a function table, it would have to be interpolated in every step from that table. So mu should be a function argument that is evaluated in the approximation steps or in `f`. /// There is no guarantee of convergence in such "open" methods like Newton or secant. However, the alternative of bracketing methods requires an initial bracketing interval which might be hard to find reliably. — Lutz Lehmann, Sep 07 '22 at 05:52
REAL*8 has been -and still is- a "de facto standard" in Fortran. It is supported by all decent compilers, it will be supported "forever", and in practice it is always what it is supposed to be: a 64 bits floating point. Yes, one should encourage the use of the now standard kind mechanism in new codes, but keeping using REAL*8 is not such a big issue (and I prefer that over the awful REAL(kind=8) that is a worst practice than REAL*8, although perfectly standard) — PierU, Sep 07 '22 at 08:05
If I am right, the only possible division by zero is when derivate=0. Track this value and see why it cancels. — , Sep 08 '22 at 16:21

score 0 · Answer 1 · answered Sep 22 '22 at 22:13

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

lastchance · Answer 2 · 2022-09-07T09:34:17.063

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

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

What is this code supposed to do? Please comment it somehow. Does it solve some coding error Marco Leonardi made or does it use a better algorithm or what does it do? — Vladimir F Героям слава, Sep 07 '22 at 12:02
The code is to provide a complete compileable version of Marco Leonardi's post. Amongst other things I had to modify his "secant" routine to move the calculation of the derivative into the DO loop. Since the solution is trivial if mu and u are constants, he really needs to tell us what they are (at least, physically). As was pointed out earlier, it's more a Newton-type method (with numerically-estimated derivative) than a secant method (which simply uses a chord). If he is intending to use SI units then his value for the mass of a proton is out by a factor of 1000 - but that's for another day. — lastchance, Sep 07 '22 at 12:24
Some context would really help the answer - otherwise just like Vladimir I was wondering what is this code supposed to do. Also having the kinds of the constants consistent with the kinds of the variables would be good, as would not using external procedures, but sticking them in a module instead, — Ian Bush, Sep 07 '22 at 13:30

score -1 · Answer 3 · answered Sep 07 '22 at 14:46

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

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

I think that Vladimir F Героям слава's comment on your previous answer is applicable here too. Also Ian Bush's. — High Performance Mark, Sep 07 '22 at 15:38
@High Performance Mark: Well, I did put comments in it and I put everything in modules rather than external statements too. So, without further explanation of the physical system from Mr Leonardi ... there's little else I could do. — lastchance, Sep 07 '22 at 15:43
Is this meant to replace your earlier answer? If so, just delete the old one. If not, perhaps explain how they are different? — High Performance Mark, Sep 07 '22 at 15:54
@High Performance Mark - forgive my ignorance of the procedures in Stack Overflow: I do not know whether to delete a previous post and what that would do to other peoples' comments. "How are they different"? Well, as I just said ... I commented the code and put everything in modules. As requested. — lastchance, Sep 07 '22 at 15:57
To be clear I wasn't asking for comments, I asked for context - i.e. a (short) paragraph of text explaining why this answers the question and ideally where you had a choice of ways why you chose the one that you did. Code only answers, even if well commented, can be difficult for new coders to follow and understand the point. — Ian Bush, Sep 07 '22 at 16:22
My apologies @Ian Bush - I was actually responding to Vladimir F's request for "comments", which I misinterpreted as within, not outside, a code. It's my new-to-StackOverflow error - I will try to give more reasoned textual answers in future. — lastchance, Sep 08 '22 at 17:18

Solving a non linear equation using the secant method in Fortran

3 Answers3