Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I calculate 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.

share|improve this question
2  
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
add comment

4 Answers

up vote 13 down vote accepted

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

share|improve this answer
    
But why Mathematica does not switch to arbitrary precision, then? –  こじま Apr 3 '12 at 9:15
3  
@こじま 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
add comment

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

share|improve this answer
    
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. Apr 17 '12 at 14:09
add comment

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)

share|improve this answer
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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