I’m calculating some integrals that include the function
The full integral will be performed using Gauss-Laguerre quadrature, with the following set of abscissa’s:
real(dp), parameter :: x(32) = [&
0.0444893658332670184_dp, &
0.234526109519618537_dp, &
0.576884629301886426_dp, &
1.07244875381781763_dp, &
1.72240877644464544_dp, &
2.52833670642579488_dp, &
3.49221327302199449_dp, &
4.61645676974976739_dp, &
5.90395850417424395_dp, &
7.35812673318624111_dp, &
8.98294092421259610_dp, &
10.7830186325399721_dp, &
12.7636979867427251_dp, &
14.9311397555225573_dp, &
17.2924543367153148_dp, &
19.8558609403360547_dp, &
22.6308890131967745_dp, &
25.6286360224592478_dp, &
28.8621018163234747_dp, &
32.3466291539647370_dp, &
36.1004948057519738_dp, &
40.1457197715394415_dp, &
44.5092079957549380_dp, &
49.2243949873086392_dp, &
54.3337213333969073_dp, &
59.8925091621340182_dp, &
65.9753772879350528_dp, &
72.6876280906627086_dp, &
80.1874469779135231_dp, &
88.7353404178923987_dp, &
98.8295428682839726_dp, &
111.751398097937695_dp]
Unfortunately, for larger x I start getting NaN’s. This is due to an arithmetic underflow in the calculation of erfc, which can be caught at compile time if the array of abscissas is given the parameter attribute.
Later I discovered the Fortran intrinsic erfc_scaled sometimes abbreviated as \mathrm{erfcx} which is equal to:
Using some logarithms, it’s possible to compute f(x) as the expression
I’ll call this method A; in terms of Fortran code this becomes:
elemental function logerfc(x) result(y)
real(dp), intent(in) :: x
real(dp) :: y
y = log(erfc_scaled(x)) - x**2
end function
! Method A
elemental function f(x) result(y)
real(dp), intent(in) :: x
real(dp) :: y
y = exp(-x**2 - logerfc(x))
end function
Another alternative is to compute:
which I’ll call method B, in Fortran code
! Method B
elemental function f(x) result(y)
real(dp), intent(in) :: x
real(dp) :: y
y = 1.0_dp/erfc_scaled(x)
end function
I’ve tested this and found that A and B agree very well at smaller values, however at the larger quadrature nodes, the differences are of order 1.0E-11. (One caveat, actually I’m computing f(x - x_c), where x_c \approx 7)
Which computation method should I prefer, A or B?
Should I perform a comparison of A and B with a quadruple precision calculation of the original function and take the method with the lower error?
Edit: I am using gfortran v10.3, presumably linked with glibc v2.27.
Edit 2: Full code is provided in the details box below
Details
module lag32
implicit none
private
public :: x
integer, parameter :: dp = kind(1.0d0)
real(dp), parameter :: x(32) = [&
0.0444893658332670184_dp, &
0.234526109519618537_dp, &
0.576884629301886426_dp, &
1.07244875381781763_dp, &
1.72240877644464544_dp, &
2.52833670642579488_dp, &
3.49221327302199449_dp, &
4.61645676974976739_dp, &
5.90395850417424395_dp, &
7.35812673318624111_dp, &
8.98294092421259610_dp, &
10.7830186325399721_dp, &
12.7636979867427251_dp, &
14.9311397555225573_dp, &
17.2924543367153148_dp, &
19.8558609403360547_dp, &
22.6308890131967745_dp, &
25.6286360224592478_dp, &
28.8621018163234747_dp, &
32.3466291539647370_dp, &
36.1004948057519738_dp, &
40.1457197715394415_dp, &
44.5092079957549380_dp, &
49.2243949873086392_dp, &
54.3337213333969073_dp, &
59.8925091621340182_dp, &
65.9753772879350528_dp, &
72.6876280906627086_dp, &
80.1874469779135231_dp, &
88.7353404178923987_dp, &
98.8295428682839726_dp, &
111.751398097937695_dp]
end module
program main
use lag32, only: x
implicit none
integer, parameter :: dp = kind(1.0d0)
real(dp), parameter :: mc = 3.0e-12_dp
real(dp), parameter :: epsilon = 3*sqrt(2.0_dp)*1.0E-13_dp
real(dp), parameter :: xc = mc/epsilon
real(dp), allocatable :: ya(:), yb(:)
integer :: i
! Original
! y = exp(-(x - xc)**2)/erfc(x - xc)
! Method A
ya = exp(-(x - xc)**2 - logerfc(x - xc))
! Method B
yb = 1.0_dp/erfc_scaled(x - xc)
do i = 1, size(x)
print *, x(i), ya(i) ,ya(i) - yb(i)
end do
contains
elemental function logerfc(x) result(y)
real(dp), intent(in) :: x
real(dp) :: y
y = log(erfc_scaled(x)) - x**2
end function
end program
. As such it underflows for 
when using IEEE single precision arithmetic. The exponentially-scaled complementary error function is asymptotic to 
. As such it does not underflow until 
.
