19
$\begingroup$

I calculate a scaled error function defined as

f[x_] := Erfc[x]*Exp[x^2]

but it can not calculate f[30000.]. f[20000.] is not very small (0.0000282). I think Mathematica is supposed to switch to high precision instead of machine precision, but it does not. It says:

General::unfl: Underflow occurred in computation. >>
General::ovfl: Overflow occurred in computation. >>

How can I calculate f for large values of x? Even with N[f[30000], 50], it does not use high precision and fails.

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ Mathematica doesn't switch to arbitrary precision like you seem to believe. If you enter a machine precison number like 30000. all further calculations are done in machine precision. You may want to read some of the tutorials on the bottom of this page. $\endgroup$
    – Sjoerd C. de Vries
    Commented Apr 3, 2012 at 9:38
  • $\begingroup$ The answer mathematica.stackexchange.com/a/161285/534 is exact. It should be the accepted answer. Can you please pick it so it moves to the top? $\endgroup$
    – a06e
    Commented May 22, 2020 at 15:22

5 Answers 5

Reset to default
17
$\begingroup$

If you have an analytic formula for f[x_] := Erfc[x]*Exp[x^2] not using Erfc[x] you could do what you expect. However it is somewhat problematic to do in this form because Erfc[x] < $MinNumber for x == 27300.

$MinNumber

1.887662394852454*10^-323228468

N[Erfc[27280.], 20]

5.680044213569341*10^-323201264

Edit

A very good approximation of your function f[x] for values x > 27280 you can get making use of these bounds ( see e.g. Erfc on Mathworld) :

enter image description here

which hold for x > 0.

Here we find values of the lower and upper bounds with relative errors for various x:

T = Table[ 
          N[#, 15]& @ {2 /(Sqrt[Pi] (x + Sqrt[2 + x^2])), 
                       2 /(Sqrt[Pi] ( x + Sqrt[x^2 + 4/Pi])), 
                       1 - ( x + Sqrt[x^2 + 4/Pi])/(x + Sqrt[2 + x^2]),
          {x, 30000 Table[10^k, {k, 0, 5}]}];

Grid[ Array[ InputField[ Dynamic[T[[#1, #2]]], FieldSize -> 13] &, {6, 3}]]

enter image description here

Therefore we propose this definition of the function f (namely the arithetic mean of its bounds for x > 27280 ) :

f[x_]/; x >= 0 := Piecewise[ { { Erfc[x]*Exp[x^2],                      x < 27280 },

                               { 1 /( Sqrt[Pi] ( x + Sqrt[2 + x^2])) 
                               + 1 /( Sqrt[Pi] ( x + Sqrt[x^2 + 4/Pi])), x >= 27280}}
                           ]
f[x_] /; x < 0 := 2 - f[-x]

I.e. we use the original definition of the function f for 0 < x < 27280, the approximation for x > 27280 and for x < 0 we use the known symmetry of the Erfc function, which is relevant when we'd like to calculate f[x] for x < - 27280. Now we can safely use this new definition for a much larger domain :

{f[300], f[300.], f[30000.], f[-30000.]}

{E^90000 Erfc[300], 0.0018806214974, 0.0000188063, 1.99998}

and now we can make plots of f around of the gluing point ( x = 27280.)

GraphicsRow[{ Plot[ f[x], {x, 2000, 55000}, 
                      Epilog -> {PointSize[0.02], Red, Point[{27280., f[27280.]}]}, 
                      PlotStyle -> Thick, AxesOrigin -> {0, 0}], 
              Plot[ f[x], {x, 27270, 27290}, 
                      Epilog -> {PointSize[0.02], Red, Point[{27280., f[27280.]}]}, 
                      PlotStyle -> Thick]}]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ But why Mathematica does not switch to arbitrary precision, then? $\endgroup$
    – こじま
    Commented Apr 3, 2012 at 9:15
  • 4
    $\begingroup$ @こじま From the docs: $MinNumber gives the minimum positive arbitrary-precision number that can be represented on a particular computer system. By the way: could you perhaps choose a user name that one can easily type using a standard keyboard? It's pretty difficult to @-refer to you in this way. $\endgroup$
    – Sjoerd C. de Vries
    Commented Apr 3, 2012 at 9:24
15
$\begingroup$

I've had to work with that kind of function (relying on cancellation of large terms) before, and the most practical workaround I could figure out to be able to evaluate the function numerically is to use its power expansion near the point of trouble (here, $+\infty$).

So, get a good look at the series expansion and find out how it works (or derive it on paper):

Series[f[x], {x, ∞, 50}]

So, that means that you can use

$$f(x) \approx \sum_{n=0}^N \left(-\frac{1}{2}\right)^n \frac{(2n-1)!!}{x^{2n+1}}$$

for $x$ larger than some value $X$. All that remains is to find a couple of values $(X,N)$ suitable to your needs. Because we have a series with signs alternating and terms decrease, the error is bounded by the $(n+1)$th term. Assuming we want to choose a value of $X$ in the range [20,1000], we plot the relative error after the $N$th term as a function of $x$ in this range:

LogLogPlot[Table[(2 n - 1)!!/x^(2 n + 1)/f[x], {n, 5, 10}], {x, 20, 1000}]

enter image description here

So, say we want to have a relative accuracy of $10^{-30}$ (which is better than machine precision), we can for example take $X=100$ and $N=10$. This gives us the following definition for your f function:

f[x_] := Piecewise[{{Erfc[x]*Exp[x^2], x < 100}}, 
   1/Sqrt[π]*Sum[(-1/2)^n*(2 n - 1)!!/x^(2 n + 1), {n, 0, 10}]];

LogLogPlot[{f[x], Erfc[x]*Exp[x^2]}, {x, 10, 10^6}, 
 PlotStyle -> {Red, Directive[Blue, Thick, Opacity[0.4]]}]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ This is a good idea I also wanted to work on the same direction after the university hours! Nice job... $\endgroup$
    – PlatoManiac
    Commented Apr 3, 2012 at 10:42
  • $\begingroup$ There is also the possibility of constructing a Padé approximant from your series. For instance, PadeApproximant[Erfc[x] Exp[x^2], {x, Infinity, {4, 4}}] yields an approximant with absolute relative error $< 10^{-13}$ for $x > 50$. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Apr 17, 2012 at 14:09
14
$\begingroup$

For numerical evaluation, there is the rapidly-converging continued fraction (due to Jones and Thron):

$$\exp(x^2)\mathrm{erfc}(x)=\frac{2x}{\sqrt \pi}\cfrac{1}{2x^2+1-\cfrac{1\cdot2}{2x^2+5-\cfrac{3\cdot4}{2x^2+9-\cdots}}},\qquad x > 0$$

One can use the built-in function ContinuedFractionK[] with a suitable cut-off:

With[{x = N[30000], n = 10}, -2 x /(1 + 2 x^2 + 
      ContinuedFractionK[2 j (1 - 2 j), 1 + 4 j + 2 x^2, {j, 1, n}])/
   Sqrt[Pi]]
0.0000188063

or, even better, use the Lentz-Thompson-Barnett algorithm for evaluating this continued fraction, avoiding unneeded evaluation effort:

f[z_?InexactNumberQ] := Module[{c, d, h, k, u, v, y},
   y = v = 2 z^2 + 1;
   c = y; d = 0; k = 1;
   While[True,
    u = k (k + 1); v += 4;
    c = v - u/c; d = 1/(v - u d);
    h = c*d; y *= h;
    If[Abs[h - 1] <= 10^-Precision[z], Break[]];
    k += 2];
   2 z/y/Sqrt[Pi]] /; Re[z] > 0

With[{z = N[50]}, {Exp[z^2] Erfc[z], f[z]}]
{0.01128153626532, 0.0112815}

N[f[30000], 30]
0.0000188063194411439209981315314042

Plot[f[x], {x, 1, 50}]

plot of exp(x^2)erfc(x)

Improve this answer
$\endgroup$
10
$\begingroup$

In the interest of showing that there's more than one way to skin a cat, I present a method suitable for large positive arguments, due to Chiarella and Reichel. The method uses the approximation

$$\exp(z^2)\mathrm{erfc}(z)\approx\frac{hz}{\pi}\left(z^{-2}+2\sum_{k \geq 1}\frac{\exp(-h^2 k^2)}{z^2+h^2 k^2}\right)$$

where $h$ is a suitably chosen parameter, based on the precision needed.

f[z_?InexactNumberQ] := 
 Module[{prec = Precision[z], y = z^2, e, h, j, k, s, t},
   h = Pi/Sqrt[(Round[prec] + 1) Log[10]]; e = h^2;
   s = 0; j = k = 1;
   While[True,
    t = Exp[-e k]/(y + e k);
    s += t;
    If[Abs[t] <= Abs[s] 10^-prec, Break[]];
    j += 2; k += j];
   h z (1/y + 2 s)/Pi] /; TrueQ[Quiet[z > 0]]

Again, this works best for large positive $z$, which seems to be the arguments of interest for the OP anyway. If evaluation for small $z$ is needed, a correction term has to be added to the Chiarella-Reichel approximation; see their paper for details.

Improve this answer
$\endgroup$
9
$\begingroup$

You need to avoid the underflow or overflow inside the formula $f[x\_]$. And that can be done simply by:

$$ g[x\_] := \frac{2}{\sqrt{\pi}}\ \text{HermiteH}[-1, x].$$

Use this $g[x]$ in place of the original $f[x]$, e.g. $ g[30000.]\approx 0.0000188063$.

To see why, look the general form of Hermite Polynomial:

$$ H_n(x) = (-1)^n e^{x^2} \frac{\text{d}^n}{\text{d}x^n} e^{-x^2}.$$

When $n=-1$, the derivative becomes an integral, and the integral bounds happen to be the same as Erfc in Mathematica. Therefore,

$$ f[x] = \text{Erfc}[x] * \text{Exp}[x^2] = \frac{2}{\sqrt{\pi}} \int_{x}^\infty e^{-x^2} \text{d}x \ e^{x^2} = g[x] .$$

i.e. we transformed $f[x]$ to a single function $\text{HermiteH}[-1, x]$, so Mathematica can see the whole purpose and achieve the said accuracy.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Sqrt[2/Pi] Exp[x^2/2] ParabolicCylinderD[-1, Sqrt[2] x] also works. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Sep 27, 2018 at 2:59
  • $\begingroup$ @J.M. That still has the exponential factor that can overflow. $\endgroup$
    – a06e
    Commented May 22, 2020 at 15:15
  • $\begingroup$ This is exact. It should be the accepted answer $\endgroup$
    – a06e
    Commented May 22, 2020 at 15:16
  • $\begingroup$ @becko, yes, that no longer works since 11.3; and somewhat similar to the Dawson integral question you asked some time back. It's all the same root problem. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 22, 2020 at 15:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.