# Numerical underflow for a scaled error function

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, 50], it does not use high precision and fails.

• 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. – Sjoerd C. de Vries Apr 3 '12 at 9:38

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) : 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}]] 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, f[300.], f[30000.], f[-30000.]}

{E^90000 Erfc, 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]}] • But why Mathematica does not switch to arbitrary precision, then? – こじま Apr 3 '12 at 9:15
• @こじま 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. – Sjoerd C. de Vries Apr 3 '12 at 9:24 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}] 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]]}] • This is a good idea I also wanted to work on the same direction after the university hours! Nice job... – PlatoManiac Apr 3 '12 at 10:42 • 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$. – J. M. is away Apr 17 '12 at 14:09 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, 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}, {Exp[z^2] Erfc[z], f[z]}] {0.01128153626532, 0.0112815} N[f, 30] 0.0000188063194411439209981315314042 Plot[f[x], {x, 1, 50}] 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]; 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.

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]$$

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 bound 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] .$$

• Sqrt[2/Pi] Exp[x^2/2] ParabolicCylinderD[-1, Sqrt x] also works. – J. M. is away Sep 27 '18 at 2:59