Функция Планка в Фортране 90

Question

У меня есть следующая формула для функции в Фортране 90, используемая для нахождения спектрального излучения при заданной длине волны и температуре. Я выразил это как:

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

Где я принимаю длину волны в качестве аргумента, который затем преобразуется в nu. Поскольку мои длины волн стремятся к бесконечности, Фортран начинает говорить, что интенсивность бесконечна, возвращая Inf. Однако функция Планка сходится к 0, когда длина волны стремится к бесконечности. Почему это так?

Моя теория такова, что, возможно, поскольку nu = c / lambda, поскольку лямбда стремится к бесконечности, часть функции делится на 0, так что это стремится к бесконечности. Тем не менее, существует также степенной закон -5 о длине волны, который предназначен для общего сохранения функции, сходящейся к нулю, но, похоже, это не так в вычислениях Фортрана. Почему это? Что не так? По сути, он говорит, что ноль раз Inf равно Inf.

Для протокола, я не думал, что отправка блоков моего кода здесь важна - это часть упражнения, в котором эта проблема должна возникнуть Я знаю, что мой код не работает со сбоями.

Answer 1

Я думаю, это потому, что exp (x) - 1 в знаменателе становится точно 0 (с двойной точностью) для очень малого 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

Результат (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

Согласно эта страница , exp (x) - 1 (для маленьких x) банок оценивать более стабильно как 2 * tanh (x / 2) / (1 - tanh (x / 2)). Последнее выражение действительно улучшает результат.

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

Результат:

  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