-1

I have the following formula for a function in Fortran 90 used to find the spectral radiance at a given wavelength and temperature. I've expressed it as:

intensity = (2.0d0*h*nu**(3.0d0))/c**(2.0d0) * 1d0 / (EXP((h*nu)/(k*T))-1d0)

Where I take wavelength as an argument, which is then converted to nu. As my wavelengths tend to infinity, Fortran begins to say that the intensity is infinite, returning Inf. However, the Planck function converges to 0 as the wavelength tends to infinity. Why is this?

My theory is that perhaps, since nu = c/lambda, as lambda tends to infinity, a part of the function is divided by 0, so this tends it to infinity. However, there is also a -5 power law on wavelength as well, which is meant to overall keep the function converging at 0, but this doesn't seem to be the case by Fortran's computation. Why is this? What's going wrong? It's essentially saying zero times Inf is Inf.

For the record, I didn't think submitting blocks of my code was important here -- this is part of an exercise where this problem is meant to arise, so I know my code isn't malfunctioning.

Ian Bush
  • 6,996
  • 1
  • 21
  • 27
genjong
  • 115
  • 1
  • 15
  • I think the Stack Exchange "Mathematics" might be a better fit for this. – Frank Merrow Feb 15 '20 at 23:06
  • I disagree. Mathematically, this function should converge to 0, but not with FORTRAN specifically. – genjong Feb 16 '20 at 00:07
  • The name of the language is Fortran. It has not been spelled as FORTRAN since 1990. Why are you multiplying by 1d0 in the numerator? Why are you using `REAL`-valued exponents. That could be your problem, because `nu**(3.d0)` is likely evaluated as `exp(3.d0*log(nu))`. – evets Feb 16 '20 at 01:02
  • "For the record, I didn't think submitting blocks of my code was important here -- this is part of an exercise where this problem is meant to arise, so I know my code isn't malfunctioning. " So you want us to solve your homework? But still even without a full MWE it is hard to tell why the function diverges – albert Feb 16 '20 at 09:42
  • You might be interested in using [this method to compute `exp(x) - 1`](https://stackoverflow.com/questions/30393928) – kvantour Feb 18 '20 at 14:53

1 Answers1

0

I think it is because exp(x) - 1 in the denominator becomes exactly 0 (in double precision) for very small x (cf. Loss of significance), while the numerator remains nonzero for the same x.

program main
    implicit none
    integer, parameter :: dp = kind(0.0d0)
    real(dp) :: nu, h, kT, ana, lim, numer, denom, expfac

    h = 0.8d0; kT = 1.2d0

    nu = 1.0d0
    do while ( nu > 1.0d-20 )

        expfac = exp((h * nu) / kT)
        numer = 2.0d0 * h * nu**3   ! (nu**3.0d0 also works but nu**3 is recommended)
        denom = expfac - 1.0d0

        !! Intensity (analytical & limiting expressions).
        ana = numer / denom
        lim = 2.0d0 * nu**2 * kT   !! use exp(x) ~ 1 + x

        print "(5(2x,a,es9.2),2x,a,es24.17)", &
                "nu=", nu, "ana=", ana, "lim=", lim, &
                "numer=", numer, "denom=", denom, "exp=", expfac

        nu = nu / 2.0d0
    enddo
end

Result (gfortran-9)

  nu= 1.00E+00  ana= 1.69E+00  lim= 2.40E+00  numer= 1.60E+00  denom= 9.48E-01  exp= 1.94773404105467596E+00
  nu= 5.00E-01  ana= 5.06E-01  lim= 6.00E-01  numer= 2.00E-01  denom= 3.96E-01  exp= 1.39561242508608951E+00
  nu= 2.50E-01  ana= 1.38E-01  lim= 1.50E-01  numer= 2.50E-02  denom= 1.81E-01  exp= 1.18136041286564608E+00
  nu= 1.25E-01  ana= 3.60E-02  lim= 3.75E-02  numer= 3.13E-03  denom= 8.69E-02  exp= 1.08690404952122899E+00
  nu= 6.25E-02  ana= 9.18E-03  lim= 9.37E-03  numer= 3.91E-04  denom= 4.25E-02  exp= 1.04254690518999138E+00
  nu= 3.12E-02  ana= 2.32E-03  lim= 2.34E-03  numer= 4.88E-05  denom= 2.11E-02  exp= 1.02105186214510746E+00
  nu= 1.56E-02  ana= 5.83E-04  lim= 5.86E-04  numer= 6.10E-06  denom= 1.05E-02  exp= 1.01047110901059778E+00
  nu= 7.81E-03  ana= 1.46E-04  lim= 1.46E-04  numer= 7.63E-07  denom= 5.22E-03  exp= 1.00522192027959556E+00
  nu= 3.91E-03  ana= 3.66E-05  lim= 3.66E-05  numer= 9.54E-08  denom= 2.61E-03  exp= 1.00260756045403721E+00
  nu= 1.95E-03  ana= 9.15E-06  lim= 9.16E-06  numer= 1.19E-08  denom= 1.30E-03  exp= 1.00130293141188642E+00
  nu= 9.77E-04  ana= 2.29E-06  lim= 2.29E-06  numer= 1.49E-09  denom= 6.51E-04  exp= 1.00065125364029117E+00
  nu= 4.88E-04  ana= 5.72E-07  lim= 5.72E-07  numer= 1.86E-10  denom= 3.26E-04  exp= 1.00032557382098908E+00
  nu= 2.44E-04  ana= 1.43E-07  lim= 1.43E-07  numer= 2.33E-11  denom= 1.63E-04  exp= 1.00016277366286199E+00
  nu= 1.22E-04  ana= 3.58E-08  lim= 3.58E-08  numer= 2.91E-12  denom= 8.14E-05  exp= 1.00008138351979237E+00
  nu= 6.10E-05  ana= 8.94E-09  lim= 8.94E-09  numer= 3.64E-13  denom= 4.07E-05  exp= 1.00004069093202008E+00
  nu= 3.05E-05  ana= 2.24E-09  lim= 2.24E-09  numer= 4.55E-14  denom= 2.03E-05  exp= 1.00002034525904526E+00
  nu= 1.53E-05  ana= 5.59E-10  lim= 5.59E-10  numer= 5.68E-15  denom= 1.02E-05  exp= 1.00001017257778191E+00
  nu= 7.63E-06  ana= 1.40E-10  lim= 1.40E-10  numer= 7.11E-16  denom= 5.09E-06  exp= 1.00000508627595597E+00
  nu= 3.81E-06  ana= 3.49E-11  lim= 3.49E-11  numer= 8.88E-17  denom= 2.54E-06  exp= 1.00000254313474413E+00
  nu= 1.91E-06  ana= 8.73E-12  lim= 8.73E-12  numer= 1.11E-17  denom= 1.27E-06  exp= 1.00000127156656360E+00
  nu= 9.54E-07  ana= 2.18E-12  lim= 2.18E-12  numer= 1.39E-18  denom= 6.36E-07  exp= 1.00000063578307974E+00
  nu= 4.77E-07  ana= 5.46E-13  lim= 5.46E-13  numer= 1.73E-19  denom= 3.18E-07  exp= 1.00000031789148935E+00
  nu= 2.38E-07  ana= 1.36E-13  lim= 1.36E-13  numer= 2.17E-20  denom= 1.59E-07  exp= 1.00000015894573213E+00
  nu= 1.19E-07  ana= 3.41E-14  lim= 3.41E-14  numer= 2.71E-21  denom= 7.95E-08  exp= 1.00000007947286296E+00
  nu= 5.96E-08  ana= 8.53E-15  lim= 8.53E-15  numer= 3.39E-22  denom= 3.97E-08  exp= 1.00000003973643059E+00
  nu= 2.98E-08  ana= 2.13E-15  lim= 2.13E-15  numer= 4.24E-23  denom= 1.99E-08  exp= 1.00000001986821507E+00
  nu= 1.49E-08  ana= 5.33E-16  lim= 5.33E-16  numer= 5.29E-24  denom= 9.93E-09  exp= 1.00000000993410754E+00
  nu= 7.45E-09  ana= 1.33E-16  lim= 1.33E-16  numer= 6.62E-25  denom= 4.97E-09  exp= 1.00000000496705366E+00
  nu= 3.73E-09  ana= 3.33E-17  lim= 3.33E-17  numer= 8.27E-26  denom= 2.48E-09  exp= 1.00000000248352694E+00
  nu= 1.86E-09  ana= 8.33E-18  lim= 8.33E-18  numer= 1.03E-26  denom= 1.24E-09  exp= 1.00000000124176336E+00
  nu= 9.31E-10  ana= 2.08E-18  lim= 2.08E-18  numer= 1.29E-27  denom= 6.21E-10  exp= 1.00000000062088179E+00
  nu= 4.66E-10  ana= 5.20E-19  lim= 5.20E-19  numer= 1.62E-28  denom= 3.10E-10  exp= 1.00000000031044078E+00
  nu= 2.33E-10  ana= 1.30E-19  lim= 1.30E-19  numer= 2.02E-29  denom= 1.55E-10  exp= 1.00000000015522050E+00
  nu= 1.16E-10  ana= 3.25E-20  lim= 3.25E-20  numer= 2.52E-30  denom= 7.76E-11  exp= 1.00000000007761014E+00
  nu= 5.82E-11  ana= 8.13E-21  lim= 8.13E-21  numer= 3.16E-31  denom= 3.88E-11  exp= 1.00000000003880518E+00
  nu= 2.91E-11  ana= 2.03E-21  lim= 2.03E-21  numer= 3.94E-32  denom= 1.94E-11  exp= 1.00000000001940248E+00
  nu= 1.46E-11  ana= 5.08E-22  lim= 5.08E-22  numer= 4.93E-33  denom= 9.70E-12  exp= 1.00000000000970135E+00
  nu= 7.28E-12  ana= 1.27E-22  lim= 1.27E-22  numer= 6.16E-34  denom= 4.85E-12  exp= 1.00000000000485056E+00
  nu= 3.64E-12  ana= 3.18E-23  lim= 3.18E-23  numer= 7.70E-35  denom= 2.43E-12  exp= 1.00000000000242539E+00
  nu= 1.82E-12  ana= 7.94E-24  lim= 7.94E-24  numer= 9.63E-36  denom= 1.21E-12  exp= 1.00000000000121259E+00
  nu= 9.09E-13  ana= 1.98E-24  lim= 1.99E-24  numer= 1.20E-36  denom= 6.06E-13  exp= 1.00000000000060640E+00
  nu= 4.55E-13  ana= 4.96E-25  lim= 4.96E-25  numer= 1.50E-37  denom= 3.03E-13  exp= 1.00000000000030309E+00
  nu= 2.27E-13  ana= 1.24E-25  lim= 1.24E-25  numer= 1.88E-38  denom= 1.52E-13  exp= 1.00000000000015166E+00
  nu= 1.14E-13  ana= 3.10E-26  lim= 3.10E-26  numer= 2.35E-39  denom= 7.57E-14  exp= 1.00000000000007572E+00
  nu= 5.68E-14  ana= 7.74E-27  lim= 7.75E-27  numer= 2.94E-40  denom= 3.80E-14  exp= 1.00000000000003797E+00
  nu= 2.84E-14  ana= 1.95E-27  lim= 1.94E-27  numer= 3.67E-41  denom= 1.89E-14  exp= 1.00000000000001887E+00
  nu= 1.42E-14  ana= 4.81E-28  lim= 4.85E-28  numer= 4.59E-42  denom= 9.55E-15  exp= 1.00000000000000955E+00
  nu= 7.11E-15  ana= 1.23E-28  lim= 1.21E-28  numer= 5.74E-43  denom= 4.66E-15  exp= 1.00000000000000466E+00
  nu= 3.55E-15  ana= 2.94E-29  lim= 3.03E-29  numer= 7.17E-44  denom= 2.44E-15  exp= 1.00000000000000244E+00
  nu= 1.78E-15  ana= 8.08E-30  lim= 7.57E-30  numer= 8.97E-45  denom= 1.11E-15  exp= 1.00000000000000111E+00
  nu= 8.88E-16  ana= 1.68E-30  lim= 1.89E-30  numer= 1.12E-45  denom= 6.66E-16  exp= 1.00000000000000067E+00
  nu= 4.44E-16  ana= 6.31E-31  lim= 4.73E-31  numer= 1.40E-46  denom= 2.22E-16  exp= 1.00000000000000022E+00
  nu= 2.22E-16  ana= 7.89E-32  lim= 1.18E-31  numer= 1.75E-47  denom= 2.22E-16  exp= 1.00000000000000022E+00
  nu= 1.11E-16  ana= Infinity  lim= 2.96E-32  numer= 2.19E-48  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 5.55E-17  ana= Infinity  lim= 7.40E-33  numer= 2.74E-49  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 2.78E-17  ana= Infinity  lim= 1.85E-33  numer= 3.42E-50  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 1.39E-17  ana= Infinity  lim= 4.62E-34  numer= 4.28E-51  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 6.94E-18  ana= Infinity  lim= 1.16E-34  numer= 5.35E-52  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 3.47E-18  ana= Infinity  lim= 2.89E-35  numer= 6.68E-53  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 1.73E-18  ana= Infinity  lim= 7.22E-36  numer= 8.35E-54  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 8.67E-19  ana= Infinity  lim= 1.81E-36  numer= 1.04E-54  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 4.34E-19  ana= Infinity  lim= 4.51E-37  numer= 1.31E-55  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 2.17E-19  ana= Infinity  lim= 1.13E-37  numer= 1.63E-56  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 1.08E-19  ana= Infinity  lim= 2.82E-38  numer= 2.04E-57  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 5.42E-20  ana= Infinity  lim= 7.05E-39  numer= 2.55E-58  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 2.71E-20  ana= Infinity  lim= 1.76E-39  numer= 3.19E-59  denom= 0.00E+00  exp= 1.00000000000000000E+00
  nu= 1.36E-20  ana= Infinity  lim= 4.41E-40  numer= 3.98E-60  denom= 0.00E+00  exp= 1.00000000000000000E+00

According to this page, exp(x) - 1 (for small x) can be evaluated in a more stable manner as 2 * tanh(x/2) / (1 - tanh(x/2)). The latter expression indeed seems to improve the result.

program main
    implicit none
    integer, parameter :: dp = kind(0.0d0)
    real(dp) :: nu, h, kT, ana, lim, numer, denom, expfac, x

    h = 0.8d0; kT = 1.2d0

    nu = 1.0d0
    do while ( nu > 1.0d-20 )

        x = h * nu / kT
        numer = 2.0d0 * h * nu**3
        denom = (2 * tanh(x / 2)) / (1 - tanh(x / 2))

        ana = numer / denom
        lim = 2.0d0 * nu**2 * kT

        print "(5(2x,a,es9.2))", &
                "nu=", nu, "ana=", ana, "lim=", lim, &
                "numer=", numer, "denom=", denom

        nu = nu / 2.0d0
    enddo
end

Result:

  nu= 1.00E+00  ana= 1.69E+00  lim= 2.40E+00  numer= 1.60E+00  denom= 9.48E-01
  nu= 5.00E-01  ana= 5.06E-01  lim= 6.00E-01  numer= 2.00E-01  denom= 3.96E-01
  nu= 2.50E-01  ana= 1.38E-01  lim= 1.50E-01  numer= 2.50E-02  denom= 1.81E-01
  nu= 1.25E-01  ana= 3.60E-02  lim= 3.75E-02  numer= 3.13E-03  denom= 8.69E-02
  nu= 6.25E-02  ana= 9.18E-03  lim= 9.37E-03  numer= 3.91E-04  denom= 4.25E-02
  nu= 3.12E-02  ana= 2.32E-03  lim= 2.34E-03  numer= 4.88E-05  denom= 2.11E-02
  nu= 1.56E-02  ana= 5.83E-04  lim= 5.86E-04  numer= 6.10E-06  denom= 1.05E-02
  nu= 7.81E-03  ana= 1.46E-04  lim= 1.46E-04  numer= 7.63E-07  denom= 5.22E-03
  nu= 3.91E-03  ana= 3.66E-05  lim= 3.66E-05  numer= 9.54E-08  denom= 2.61E-03
  nu= 1.95E-03  ana= 9.15E-06  lim= 9.16E-06  numer= 1.19E-08  denom= 1.30E-03
  nu= 9.77E-04  ana= 2.29E-06  lim= 2.29E-06  numer= 1.49E-09  denom= 6.51E-04
  nu= 4.88E-04  ana= 5.72E-07  lim= 5.72E-07  numer= 1.86E-10  denom= 3.26E-04
  nu= 2.44E-04  ana= 1.43E-07  lim= 1.43E-07  numer= 2.33E-11  denom= 1.63E-04
  nu= 1.22E-04  ana= 3.58E-08  lim= 3.58E-08  numer= 2.91E-12  denom= 8.14E-05
  nu= 6.10E-05  ana= 8.94E-09  lim= 8.94E-09  numer= 3.64E-13  denom= 4.07E-05
  nu= 3.05E-05  ana= 2.24E-09  lim= 2.24E-09  numer= 4.55E-14  denom= 2.03E-05
  nu= 1.53E-05  ana= 5.59E-10  lim= 5.59E-10  numer= 5.68E-15  denom= 1.02E-05
  nu= 7.63E-06  ana= 1.40E-10  lim= 1.40E-10  numer= 7.11E-16  denom= 5.09E-06
  nu= 3.81E-06  ana= 3.49E-11  lim= 3.49E-11  numer= 8.88E-17  denom= 2.54E-06
  nu= 1.91E-06  ana= 8.73E-12  lim= 8.73E-12  numer= 1.11E-17  denom= 1.27E-06
  nu= 9.54E-07  ana= 2.18E-12  lim= 2.18E-12  numer= 1.39E-18  denom= 6.36E-07
  nu= 4.77E-07  ana= 5.46E-13  lim= 5.46E-13  numer= 1.73E-19  denom= 3.18E-07
  nu= 2.38E-07  ana= 1.36E-13  lim= 1.36E-13  numer= 2.17E-20  denom= 1.59E-07
  nu= 1.19E-07  ana= 3.41E-14  lim= 3.41E-14  numer= 2.71E-21  denom= 7.95E-08
  nu= 5.96E-08  ana= 8.53E-15  lim= 8.53E-15  numer= 3.39E-22  denom= 3.97E-08
  nu= 2.98E-08  ana= 2.13E-15  lim= 2.13E-15  numer= 4.24E-23  denom= 1.99E-08
  nu= 1.49E-08  ana= 5.33E-16  lim= 5.33E-16  numer= 5.29E-24  denom= 9.93E-09
  nu= 7.45E-09  ana= 1.33E-16  lim= 1.33E-16  numer= 6.62E-25  denom= 4.97E-09
  nu= 3.73E-09  ana= 3.33E-17  lim= 3.33E-17  numer= 8.27E-26  denom= 2.48E-09
  nu= 1.86E-09  ana= 8.33E-18  lim= 8.33E-18  numer= 1.03E-26  denom= 1.24E-09
  nu= 9.31E-10  ana= 2.08E-18  lim= 2.08E-18  numer= 1.29E-27  denom= 6.21E-10
  nu= 4.66E-10  ana= 5.20E-19  lim= 5.20E-19  numer= 1.62E-28  denom= 3.10E-10
  nu= 2.33E-10  ana= 1.30E-19  lim= 1.30E-19  numer= 2.02E-29  denom= 1.55E-10
  nu= 1.16E-10  ana= 3.25E-20  lim= 3.25E-20  numer= 2.52E-30  denom= 7.76E-11
  nu= 5.82E-11  ana= 8.13E-21  lim= 8.13E-21  numer= 3.16E-31  denom= 3.88E-11
  nu= 2.91E-11  ana= 2.03E-21  lim= 2.03E-21  numer= 3.94E-32  denom= 1.94E-11
  nu= 1.46E-11  ana= 5.08E-22  lim= 5.08E-22  numer= 4.93E-33  denom= 9.70E-12
  nu= 7.28E-12  ana= 1.27E-22  lim= 1.27E-22  numer= 6.16E-34  denom= 4.85E-12
  nu= 3.64E-12  ana= 3.18E-23  lim= 3.18E-23  numer= 7.70E-35  denom= 2.43E-12
  nu= 1.82E-12  ana= 7.94E-24  lim= 7.94E-24  numer= 9.63E-36  denom= 1.21E-12
  nu= 9.09E-13  ana= 1.99E-24  lim= 1.99E-24  numer= 1.20E-36  denom= 6.06E-13
  nu= 4.55E-13  ana= 4.96E-25  lim= 4.96E-25  numer= 1.50E-37  denom= 3.03E-13
  nu= 2.27E-13  ana= 1.24E-25  lim= 1.24E-25  numer= 1.88E-38  denom= 1.52E-13
  nu= 1.14E-13  ana= 3.10E-26  lim= 3.10E-26  numer= 2.35E-39  denom= 7.58E-14
  nu= 5.68E-14  ana= 7.75E-27  lim= 7.75E-27  numer= 2.94E-40  denom= 3.79E-14
  nu= 2.84E-14  ana= 1.94E-27  lim= 1.94E-27  numer= 3.67E-41  denom= 1.89E-14
  nu= 1.42E-14  ana= 4.85E-28  lim= 4.85E-28  numer= 4.59E-42  denom= 9.47E-15
  nu= 7.11E-15  ana= 1.21E-28  lim= 1.21E-28  numer= 5.74E-43  denom= 4.74E-15
  nu= 3.55E-15  ana= 3.03E-29  lim= 3.03E-29  numer= 7.17E-44  denom= 2.37E-15
  nu= 1.78E-15  ana= 7.57E-30  lim= 7.57E-30  numer= 8.97E-45  denom= 1.18E-15
  nu= 8.88E-16  ana= 1.89E-30  lim= 1.89E-30  numer= 1.12E-45  denom= 5.92E-16
  nu= 4.44E-16  ana= 4.73E-31  lim= 4.73E-31  numer= 1.40E-46  denom= 2.96E-16
  nu= 2.22E-16  ana= 1.18E-31  lim= 1.18E-31  numer= 1.75E-47  denom= 1.48E-16
  nu= 1.11E-16  ana= 2.96E-32  lim= 2.96E-32  numer= 2.19E-48  denom= 7.40E-17
  nu= 5.55E-17  ana= 7.40E-33  lim= 7.40E-33  numer= 2.74E-49  denom= 3.70E-17
  nu= 2.78E-17  ana= 1.85E-33  lim= 1.85E-33  numer= 3.42E-50  denom= 1.85E-17
  nu= 1.39E-17  ana= 4.62E-34  lim= 4.62E-34  numer= 4.28E-51  denom= 9.25E-18
  nu= 6.94E-18  ana= 1.16E-34  lim= 1.16E-34  numer= 5.35E-52  denom= 4.63E-18
  nu= 3.47E-18  ana= 2.89E-35  lim= 2.89E-35  numer= 6.68E-53  denom= 2.31E-18
  nu= 1.73E-18  ana= 7.22E-36  lim= 7.22E-36  numer= 8.35E-54  denom= 1.16E-18
  nu= 8.67E-19  ana= 1.81E-36  lim= 1.81E-36  numer= 1.04E-54  denom= 5.78E-19
  nu= 4.34E-19  ana= 4.51E-37  lim= 4.51E-37  numer= 1.31E-55  denom= 2.89E-19
  nu= 2.17E-19  ana= 1.13E-37  lim= 1.13E-37  numer= 1.63E-56  denom= 1.45E-19
  nu= 1.08E-19  ana= 2.82E-38  lim= 2.82E-38  numer= 2.04E-57  denom= 7.23E-20
  nu= 5.42E-20  ana= 7.05E-39  lim= 7.05E-39  numer= 2.55E-58  denom= 3.61E-20
  nu= 2.71E-20  ana= 1.76E-39  lim= 1.76E-39  numer= 3.19E-59  denom= 1.81E-20
  nu= 1.36E-20  ana= 4.41E-40  lim= 4.41E-40  numer= 3.98E-60  denom= 9.04E-21
roygvib
  • 7,218
  • 2
  • 19
  • 36
  • What about the zero in the numerator? – genjong Feb 16 '20 at 17:33
  • The numerator remains nonzero for small 'nu' (see "numer=" in the output), but the denominator becomes _exactly_ zero for very small 'nu' (because of exp(x) - 1), i.e. I think the denominator loses precision faster than the numerator because the denominator is a difference of values of order 1 (if x is very small). – roygvib Feb 16 '20 at 18:07
  • I think the sentence in the question "It's essentially saying zero times Inf is Inf." is the reason for confusion. We cannot set x to zero exactly in an expression if it approaches 0/0 (which is something like a removable singularity, similar to sin(x)/x). But the limit x -> eps (a small but nonzero number) may suffer from limited accuracy of floating-point numbers (which is the case for the denominator). I think this problem is not specific to Fortran but also occurs in other languages. – roygvib Feb 17 '20 at 02:11