Accurate computation of exp(-x**2)/erfc(x)

I’m calculating some integrals that include the function

f(x) = \frac{e^{-x^2}}{\mathrm{erfc}(x)}

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:

\mathrm{erfcx(x)} = e^{x^2} \mathrm{erfc} (x)

Using some logarithms, it’s possible to compute f(x) as the expression

f(x) = \exp \left(-x^2 - \underbrace{\ln(\mathrm{erfc}(x))}_{\ln(\mathrm{erfcx}(x)) - x^2}\right)

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:

f(x) = \frac{e^{x^2}}{e^{x^2}} \frac{e^{-x^2}}{\mathrm{erfc}(x)} = \frac{1}{\mathrm{erfcx}(x)}

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
1 Like

You did not state which integral you want to evaluate, but it may be that the well-known asymptotic behavior of the complementary error function may be of use.

The expression √π x exp(x^2 ) erfc(x) approaches 1 as x approaches infinity, and there is an asymptotic series in terms of (1/2x^2) for this expression, see for example https://young.physics.ucsc.edu/116A/erfc.pdf .

2 Likes

Well the integral are slightly more complex. The problem comes from applying the method of weighted residuals to a distribution function W(x) = \sum_{i=0}^N C_i \phi_i(x) expressed on the basis of Laguerre functions \phi_n(x) = e^{-x} L_n(x) where L_n is the n-th Laguerre polynomial. The weighting functions are \psi_j(x) = x^{j-1} e^{-\lambda x}. Next, defining the inner product as

\langle f, g \rangle = \int_0^\infty f(x) \cdot g(x) dx

I am looking for the matrix

\langle \psi_j, x M(x) \phi_n(x)\rangle, \\ j = 1,...,N-1,\; n = 0,...,N

where

M(x) = \frac{2}{\sqrt\pi} \frac{e^{-(x - x_c)^2}}{\mathrm{erfc}(x - x_c)}

I think you can also compute quadrature weights and points exactly for your function \exp(-x^2)\over\mathrm{erfc}(x), you can use sympy or any other arbitrary precision package to get exact double precision points and weights. Then it should work nicely.

Isn’t Method-A less accurate because of the cancellation of large values (-x**2 and logerfc(x)) for large x? According to the page below, erfcx(x) behaves like 1/x for large x, so Method-B may not suffer from such cancellation? (If the built-in routine is accurate enough…)

https://www.intel.com/content/www/us/en/develop/documentation/fortran-compiler-oneapi-dev-guide-and-reference/top/language-reference/a-to-z-reference/e-to-f/erfc-scaled.html

The complementary error function is asymptotic to note graphic a. As such it underflows for note graphicb when using IEEE single precision arithmetic. The exponentially-scaled complementary error function is asymptotic to note graphic c. As such it does not underflow until note graphic d.

1 Like

This would involve finding the recurrence relation of orthogonal polynomials for the custom weight function? Then using the one of the procedures from the OPQ package (or similar tools in SymPy) to find the custom quadrature rule. Sounds doable but complicated.

The computations we are trying to reproduce where originally done in 1969 on an IBM 7044, however the numerical values are not given so we have no means of error-checking. The authors said they verified their Gauss-Laguerre results with a Newton-Cotes rule, the two different integration methods matching to three digits. The accuracy might not have been so crucial at the time, given that the model contains some experimental parameters, and the constraints on computer time (~20 hrs). Still I’d like to be sure we’ve done the integrations as accurately as possible with modern tools.

The authors state they used the 32-node Gauss-Laguerre quadrature rule, with the abscissas taken from the book “Approximate calculation of integrals” by V. I. Krylov, published in 1962 (translation by Arthur Stroud). These are the same values provided above, meaning the domain was [0, 111.75…]

The computations were performed at the Computer Center of IIT Kanpur in India. The paper is given here: On the solution of statistical models of cell populations - ScienceDirect (however it’s missing many details of the calculation procedures).

Yes. I think this task is actually very common, so having an easy to use package (say in SymPy) would be very helpful (I created an issue for this: Gaussian quadrature rule for arbitrary weight function and interval · Issue #22594 · sympy/sympy · GitHub). On the Fortran side, it would be very simple to use the new points and weights.

I don’t see the functionality here: https://github.com/sympy/sympy/tree/master/sympy/integrals, but we should add it.

There is an easy to use function for generating arbitrary precision Gaussian quadrature nodes and weights from an arbitrary weight function in the Julia QuadGK package. See this section of the README.md file and the documentation for the gauss method. It does not require one to first find the recurrence relation of the associated orthogonal polynomials.

1 Like

Indeed. Unfortunately it seems numerical, so it works for smaller intervals, but it’s hard to get it to converge to get high accuracy:

julia> using SpecialFunctions

julia> using QuadGK

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, 0, 20, rtol=1e-9)
([0.38768537755077603, 1.7568858761982542, 3.803903391026576, 6.336780964014735, 9.152644231849653, 12.025290131154804, 14.722945139994678, 17.027550936537043, 18.752707347179335, 19.759208293813685], [1.385352625729292, 6.137392669224711, 16.158630107464138, 30.815101135175087, 47.04694093046928, 60.38252780711042, 66.31501257350119, 61.74935694664822, 46.08036828553858, 21.583335269112947])

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, 0, 30, rtol=1e-9)
ERROR: DomainError with 15.0:
integrand produced NaN in the interval (0, 30)
Stacktrace:
  [1] evalrule(f::var"#3#4", a::Int64, b::Int64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, nrm::typeof(LinearAlgebra.norm))
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/evalrule.jl:41
  [2] #2
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:10 [inlined]
  [3] ntuple
    @ ./ntuple.jl:48 [inlined]
  [4] do_quadgk(f::var"#3#4", s::Tuple{Int64, Int64}, n::Int64, atol::Nothing, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm))
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:10
  [5] (::QuadGK.var"#28#29"{Nothing, Float64, Int64, Int64, typeof(LinearAlgebra.norm)})(f::Function, s::Tuple{Int64, Int64}, #unused#::Function)
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:179
  [6] handle_infinities
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:113 [inlined]
  [7] #quadgk#27
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:177 [inlined]
  [8] _gauss(W::var"#3#4", N::Int64, a::Int64, b::Int64, rtol::Float64, quad::Function)
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/weightedgauss.jl:74
  [9] gauss(W::Function, N::Int64, a::Int64, b::Int64; rtol::Float64, quad::Function)
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/weightedgauss.jl:62
 [10] top-level scope
    @ REPL[4]:1

I think this has to be done symbolically (say in SymPy) to obtain accurate answers.

You can get quite a bit further using extended precision float types:

julia> using SpecialFunctions, QuadGK, Quadmath

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(20), rtol=1e-9)
(Float128[3.87685377550762212241773092779845364e-01, 1.75688587619823879175737923069572830e+00, 3.80390339102656684792616744050199912e+00, 6.33678096401472421945915155325746325e+00, 9.15264423184965727190966650300169598e+00, 1.20252901311547957434572683103166754e+01, 1.47229451399946729266320462093104682e+01, 1.70275509365370377100260435382257308e+01, 1.87527073471793350183947750993046031e+01, 1.97592082938136925783941627417989462e+01], Float128[1.38535262572932963828648215216764593e+00, 6.13739266922469909867619904385367344e+00, 1.61586301074638846605178484544102100e+01, 3.08151011351749197195145122889382423e+01, 4.70469409304694732526895358042952104e+01, 6.03825278071104610553729784281794292e+01, 6.63150125735013793052707456707147675e+01, 6.17493569466472859903879566530466625e+01, 4.60803682855379145652097342145010420e+01, 2.15833352691123531123418668325872807e+01])

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(30), rtol=1e-9)
(Float128[6.34659407224465463142196297858512046e-01, 2.71171768735281275985578726914630361e+00, 5.77172336757198750829571201010836730e+00, 9.55841360601610044207575638130383600e+00, 1.37700767174041187617258255999648973e+01, 1.80676417080122640594540138116366018e+01, 2.21038634233279372174925678952818276e+01, 2.55522054067187704666617804877047030e+01, 2.81336073698157885360344273804487407e+01, 2.96396882894937045382468107761158695e+01], Float128[2.67466117744153237486666477119578476e+00, 1.33106909436555609478674658058197336e+01, 3.60073195706674527496696891058388764e+01, 6.90461357829318821070728103163382012e+01, 1.05537609676403258061350223554796575e+02, 1.35475607326787583431931249161960996e+02, 1.48777295551138076209416609062028994e+02, 1.38521426665211170169908953165660962e+02, 1.03363635159087066450004869773698807e+02, 4.84118206425137524135593371145666744e+01])

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(40), rtol=1e-9)
(Float128[9.00780115683006498306807582292482960e-01, 3.67737225112924308072550023859929769e+00, 7.74713912436374194731955235756185976e+00, 1.27862534586589845322831579891730792e+01, 1.83924496689391570427540879394873073e+01, 2.41136476243153953855474693368597716e+01, 2.94872120320537425920291275042218204e+01, 3.40782336562836041685666433805137078e+01, 3.75150851457707441115157831550554809e+01, 3.95202799518310128320584312641343963e+01], Float128[4.42137875691497530987515277003128497e+00, 2.34095476609499634285705096610784881e+01, 6.38949552209232300492328233141610626e+01, 1.22642636137265081544108285968894329e+02, 1.87438369155112030182459335689592933e+02, 2.40553300756315163130237215575937092e+02, 2.64120951575474784096709157734308366e+02, 2.45880067727169222441715410450082335e+02, 1.83457612796487614807187590772431716e+02, 8.59209671882741275499940977858120291e+01])

I guess it’ll run out of gas at some point, but perhaps this is far enough for the OP’s application?

Thanks. I tried it but I am getting:

julia> using SpecialFunctions, QuadGK, Quadmath

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(20), rtol=1e-9)
ERROR: MethodError: no method matching eigen!(::LinearAlgebra.SymTridiagonal{Float128, Vector{Float128}})
Closest candidates are:
  eigen!(::LinearAlgebra.SymTridiagonal{var"#s814", V} where {var"#s814"<:Union{Float32, Float64}, V<:AbstractVector{var"#s814"}}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:280
  eigen!(::LinearAlgebra.SymTridiagonal{var"#s814", V} where {var"#s814"<:Union{Float32, Float64}, V<:AbstractVector{var"#s814"}}, ::UnitRange) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:283
  eigen!(::LinearAlgebra.SymTridiagonal{var"#s814", V} where {var"#s814"<:Union{Float32, Float64}, V<:AbstractVector{var"#s814"}}, ::Real, ::Real) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:288
  ...
Stacktrace:
 [1] eigen(A::LinearAlgebra.SymTridiagonal{Float128, Vector{Float128}})
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:281
 [2] _gauss(W::var"#5#6", N::Int64, a::Float128, b::Float128, rtol::Float64, quad::Function)
   @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/weightedgauss.jl:98
 [3] gauss(W::Function, N::Int64, a::Float128, b::Float128; rtol::Float64, quad::Function)
   @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/weightedgauss.jl:62
 [4] top-level scope
   @ REPL[8]:1

I am able to reproduce your error when I start from a clean Julia session. Prior to posting my previous result, I had attempted to use the DoubleFloats package instead of the Quadmath package. The DoubleFloats package worked for an upper limit of 20 for the range of integration but failed when I tried 30 as the upper limit. So I then switched to the higher precision Quadmath package, but in the same Julia session. If I start with a clean Julia session I find that the following works (I’ve added some checks on the results versus standard Float64 for the case with an upper limit of 10):

julia> using SpecialFunctions, QuadGK, DoubleFloats, Quadmath

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, 0, 10, rtol=1e-9)
([0.17002966561186994, 0.8275983996020136, 1.8519610263611523, 3.1268538804877104, 4.544187900979246, 5.989495947115991, 7.346369502028654, 8.5053477390109, 9.372837776505476, 9.878927439535989], [0.5113140688691087, 1.798695873902142, 4.288803672361198, 7.917855524736181, 11.95622529170446, 15.281890838276784, 16.75469046153688, 15.589435499209033, 11.629641986675479, 5.446374524138117])

julia> x2, w2 = gauss(x -> exp(-x^2) / erfc(x), 10, Double64(0), Double64(10), rtol=1e-9)
(Double64[0.17002966561186245, 0.8275983996020075, 1.851961026361147, 3.1268538804877046, 4.544187900979243, 5.9894959471159845, 7.346369502028653, 8.505347739010896, 9.372837776505476, 9.87892743953599], Double64[0.5113140688691216, 1.7986958739021315, 4.288803672361115, 7.91785552473613, 11.956225291704522, 15.28189083827682, 16.754690461536985, 15.589435499208923, 11.629641986675392, 5.446374524138095])

julia> x3, w3 = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(10), rtol=1e-9)
(Float128[1.70029665611862433659224649784303822e-01, 8.27598399602007522716310792972464948e-01, 1.85196102636114703415280335809190867e+00, 3.12685388048770480613305321518244925e+00, 4.54418790097924340654885464452625230e+00, 5.98949594711598487888835436986946672e+00, 7.34636950202865233803120255124583270e+00, 8.50534773901089557533579797067995003e+00, 9.37283777650547568931649112718171529e+00, 9.87892743953598985096491870518316858e+00], Float128[5.11314068869121514855563410584036841e-01, 1.79869587390213138052062743401453301e+00, 4.28880367236111457651850781821224047e+00, 7.91785552473613082702253498409252333e+00, 1.19562252917045218072471809363012327e+01, 1.52818908382768184272135307269923489e+01, 1.67546904615369832026084526020777775e+01, 1.55894354992089224522041996753377941e+01, 1.16296419866753912481246448811522357e+01, 5.44637452413809532657056757168147170e+00])

julia> all(x .≈ x2 .≈ x3)
true

julia> all(w .≈ w2 .≈ w3)
true

Clearly, it should not be necessary to use the DoubleFloat type first, and the results are therefore suspect, despite the checks above. I will raise this as a question on the Julia discourse group. Meantime, for whatever it’s worth, I’m providing the following result for the OP:

julia> x32_100, w32_100 = gauss(x -> exp(-x^2) / erfc(x), 32, Float128(0), Float128(100), rtol=1e-9);

julia> for (x,w) in zip(x32_100, w32_100)
         println("(", x, ", ", w, ")")
       end
(1.91963135636113087723944127177350547e-01, 5.91881243864281397404682379530372879e-01)
(9.47784705001635673735635624115440130e-01, 2.28103406732697400070025195716587568e+00)
(2.18462250050000638860517247490728286e+00, 6.15270038294108072957879090736266827e+00)
(3.85626859427070854587320930367100973e+00, 1.32715351122281943804061115615437698e+01)
(5.94303187308467266601038365516999644e+00, 2.44376376900740335278372351151039321e+01)
(8.42638901227391390038722405982762608e+00, 4.02246061438534077839506391402296262e+01)
(1.12844428360141741498899057330258462e+01, 6.09839427552017486920159068027314587e+01)
(1.44916867548752248131178061450934793e+01, 8.68305492503164311358005709379575465e+01)
(1.80192958552878944913476484862152326e+01, 1.17631095079496323732167665455715792e+02)
(2.18354546184332115314525436246271114e+01, 1.53001126546825735713393779220942402e+02)
(2.59056840759581005940415884388064161e+01, 1.92311969911848824229307360626933545e+02)
(3.01931752794688898020601038483737474e+01, 2.34707308606169716260678027779734561e+02)
(3.46591343113510252791447051883222276e+01, 2.79128854746171706162643569323787704e+02)
(3.92631402019628602635706496887434664e+01, 3.24350200145202995678699215690886389e+02)
(4.39635148406255416897901927267930996e+01, 3.69017632568504248917732540166495400e+02)
(4.87177025562547582957593255134124877e+01, 4.11696339913465915607232700011629950e+02)
(5.34826585903904232532942371599805010e+01, 4.50921293491815542882676339208858114e+02)
(5.82152359116053165875036252841666312e+01, 4.85248084612860171671828307608346226e+02)
(6.28725812916780241616703248427819922e+01, 5.13306732970190961236736744924939985e+02)
(6.74125209660953831834501594401474889e+01, 5.33851694739406618351103770667634646e+02)
(7.17939436152740935644911222492100765e+01, 5.45809077392199778016709139664051254e+02)
(7.59771728543853796379063808248813219e+01, 5.48318448338041974701335519348192086e+02)
(7.99243267019881804634226842838910505e+01, 5.40767807257881291911081577349390172e+02)
(8.35996607583697538931447220700354058e+01, 5.22820421115705296761046480758239117e+02)
(8.69698920177915542133292792344160127e+01, 4.94432508255697671603529107031386430e+02)
(9.00045004009175853681883379968963752e+01, 4.55861098417666542493676776663789586e+02)
(9.26760052938985683299002166380371211e+01, 4.07661727779874996996623675559100551e+02)
(9.49602147232202962186959757488623784e+01, 3.50675997908408904323976953687314454e+02)
(9.68364453107542932707453209942753953e+01, 2.86009489026088228757267541583554843e+02)
(9.82877127420278085522421230085364442e+01, 2.15001274170254064551428620628945615e+02)
(9.93009023803250404838548836044167475e+01, 1.39190361543864742545054280921176578e+02)
(9.98670662667736393395987132357698739e+01, 6.03633386886598401030275140720523847e+01)

Edit: The usage of Quadmath types in a call to gauss following a call using DoubleFloats types is valid, as explained in this post.

2 Likes

@ivanpribec would you mind trying the points and weights that @PeterSimon provided?

I tested the Julia code and it works for me.

Why are the weights decaying so slowly? I would expect the weights for x=100 to go to zero much faster.

Why are the weights decaying so slowly? I would expect the weights for x=100 to go to zero much faster.

Looking at this a bit more, the weight function exp(-x^2) / erfc(x) grows with increasing x, which is probably why the weights aren’t decaying. I don’t really know much about Gaussian integration theory but I would be surprised if such a weight function is even valid. I like the OP’s original thought of using a Gauss-Laguerre scheme that uses the decaying exponential as the weight function. Some of the other ideas presented in this thread can allow him to evaluate the remaining part of the integrand at large x values.

1 Like

Yeah, right after I posted my comment, I noticed that exp(-x^2) / erfc(x) ~ O(x) for x->oo, so that explains it. I would choose a weight function that decays by incorporating into this some part of the function being integrated, such as exp(-x) or something. Then the weights should look good.

@PeterSimon, @certik, thanks for the demonstration of the custom weight calculation. This is one of the big advantages of Julia, that you can leverage custom types across the whole library ecosystem.

As you’ve both noticed, the function f(x) (or M(x)) grows linearly as x approaches infinity. It’s supposed to model the probability of cell division. Above the critical mass x_c the function slopes upward.

image

However the model cell distribution behaves like exp(-x) for large x values (large cells are very rare, because they divide), so the weight in the integral actually looks something like this (depending on the trial function and weighting functions used in the method of weighted residuals):

image

My current plan was to first compare the results of Gauss-Laguerre quadrature with some adaptive methods available in QUADPACK, and MATLAB, to determine the accuracy. If not sufficient then I know custom quadrature rules are an option.

Thanks for the function using the continued fraction. It appears that in the region x < 20, the product exp(-x**2) * erfc(x) is also slightly more accurate than erf_scaled.

I’ve tested the routine and also found the absolute errors of order 10**(-14), in close agreement with the MPFR results. However if I evaluate at the quadrature nodes, instead of the whole numbers (5,10,15,…) the accuracy appears to drop back to 10**(-12) on par with Method A and B. I guess this is just a property of floating-point math?

Is the effective weight function then something like exp(-x^2) / erfc(x) * exp(-x)? That will behave nicely. If so, then you can use the well behaved points and weights for this weight function, and then when you integrate your actual function f(x) simply integrate f(x)*exp(x), because we include exp(-x) in the new weight function, so that things cancel.

Yes, precisely.

So after shuffling the terms around, we get something like:

\int_0^\infty \underbrace{e^{-(\lambda + 1)x} \frac{e^{-(x-x_c)^2}}{\mathrm{erfc}(x-x_c)}}_{\omega(x)\text{ - effective weight}} x^{i + k - 1} dx

where the i + k - 1 is a non-negative integer constant. There’s an additional \lambda factor in the exponent, which provides some control over the conditioning of the moment system.

The factor x^k originates from the terms of the Laguerre polynomial expansion. In the original work they decided that since the Laguerre polynomials are oscillating functions it’s better to integrate the monomial parts individually and then recombine them by multiplying with the coefficients of the closed-form Laguerre polynomial,

L_n(x) = \sum_{k=0}^n {n \choose k}\frac{(-1)^k}{k!} x^k.
1 Like