Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: LogErfc

  • To: mathgroup at smc.vnet.net
  • Subject: [mg122635] Re: LogErfc
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Fri, 4 Nov 2011 05:59:53 -0500 (EST)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • References: <201111030844.DAA15151@smc.vnet.net>

On 11/03/2011 03:44 AM, Shawn Garbett wrote:
> I need to have a LogErfc function to be able to compute another
> function in a numerically stable manner. I found the function in the
> GSL, I found it in R, and even in another system, there's a copy of code to do
> it floating around. However, I need it in Mathematica and Log[Erfc[x]]
> doesn't cut in due to underflow. Cody solved this in 1969, and I could
> code his solution, but I was hoping that someone had already solved
> this in Mathematica. Anyone?
>

If you are interested mostly in x either positive or not too negative, 
could do Pade approximations at zero and infinity, and use whichever is 
appropriate in given ranges.

padeZero[x_] = PadeApproximant[Log[Erfc[x]], {x,0,{5,5}}];
padeInfinity[x_] = Log[1/x] +
   PadeApproximant[Log[Erfc[x]]-Log[1/x], {x,Infinity,{5,5}}];

logErfcSmallApprox[x_] := Block[{z},padeZero[z]/.z->x]
logErfcLargeApprox[x_] := Block[{z},padeInfinity[z]/.z->x]

Quick check:

InputForm[Table[{Log[Erfc[x]],logErfcSmallApprox[x],logErfcLargeApprox[x]}, 
{x, -3.5, 10.5, 1}]]

Out[238]//InputForm=
{{0.69314680901069, 1.0493788756124784, -14.112164948457105 +
    3.141592653589793*I}, {0.6929436838471712, 0.716213867787085,
   -7.805322341465522 + 3.141592653589793*I},
  {0.6760545027339606, 0.6761720084686005, -3.370687193890008 +
    3.141592653589793*I}, {0.41903914777555973, 0.41903914813762605,
   -0.46255109569808817 + 3.141592653589793*I},
  {-0.7350111298370846, -0.7350111299422797, -0.46255109569808817},
  {-3.384492089551553, -3.3844975258920105, -3.370687193890008},
  {-7.806815272727264, -7.8072778682726955, -7.805322341465522},
  {-14.112437402148174, -14.118997281238098, -14.112164948457105},
  {-22.34976830343594, -22.390749679580356, -22.34969815365446},
  {-32.543008907376965, -32.70519933060154, -32.54298605103614},
  {-44.705670189533116, -45.18860445062361, -44.70566137269984},
  {-58.845967473872314, -60.03403469956622, -58.845963615640876},
  {-74.96923568857265, -77.5233971405431, -74.96923382750941},
  {-93.07912219224849, -98.05396915740401, -93.07912122240553},
  {-113.1782250431496, -122.17524249202071, -113.17822450502089}}

There is clearly an overlap range where both approximations are good to 
three or so places. If this is not enough, could use higher degree Pade 
approximations.

For sufficiently negative x you can approximate by Log[2].

Daniel Lichtblau
Wolfram Research



  • References:
    • LogErfc
      • From: Shawn Garbett <shawn.garbett@gmail.com>
  • Prev by Date: Re: Mathematica 8.0.4 now available
  • Next by Date: Re: Bernoulli Numbers
  • Previous by thread: LogErfc
  • Next by thread: Re: LogErfc