Re: locating overlow/underflow
- To: mathgroup at smc.vnet.net
- Subject: [mg112200] Re: locating overlow/underflow
- From: Bill Rowe <readnews at sbcglobal.net>
- Date: Sat, 4 Sep 2010 04:01:40 -0400 (EDT)
On 9/3/10 at 6:10 AM, nbbienia at cyf-kr.edu.pl (Leslaw Bieniasz) wrote:
>I am not very much familiar with MATHEMATICA, but I have to tabulate
>a certain mathematical expression in a possibly large range of
>independent variable, and with the highest possible number of
>accurate digits. I have a notebook file that does the job (included
>below in text form, together with the results) but during the cell
>evaluation I obtain error messages stating that underflow or
>overflow occurred. My question is: is there any way to rewrite the
>code in such a way so that I can obtain an information where exactly
>(that is for which values of independent variable) the errors
>actually occur? Or, is there any way to rewrite the code or change
>the program settings in such a way so that the errors do not occur
>at all? I am using MATHEMATICA 6.0.
>In[1]:= Table[{N[y],
>SetPrecision[ 2 ((y - 1) Exp[y] Erfc[Sqrt[y]] + (y/3 - 1) 2
>Sqrt[y/Pi] + 1)/(y y), 70], Precision[ 2 ((y - 1) Exp[y]
>Erfc[Sqrt[y]] + (y/3 - 1) 2 Sqrt[y/Pi] + 1)/(y y)]
>},
>{y, { 1/10^19, ..., 5 10^15
>}}]
Where to start?
SetPrecision simply doesn't do what you seem to think it does.
Using the first value for y (1/10^19) I can get 50 digit
precision by doing:
In[18]:= f =
2 ((y - 1) Exp[y] Erfc[Sqrt[y]] + (y/3 - 1) 2 Sqrt[y/Pi] +
1)/(y y);
In[19]:= N[f /. y -> 10^-19 // Simplify, 50]
Out[19]= 0.99999999971454014148222328832683872045949424371145
There is no need to compute this number twice to get the
precision since doing N[expr, 50] will always return a number
with a precision of 50 if all of the values in expr are exact values.
Strictly speaking, I did not need to simplify your expression
before converting it to a number with 50 digit precision. But
experience indicates this is a good idea for a complex
expression with extreme exact values.
While this approach will work for several values of y, it will
not work for all values.
Notice for y = 100
In[44]:= Log[10, Exp[{100, 200, 500}]] // N
Out[44]= {43.4294,86.8589,217.147}
and
In[45]:= Log[10, Erfc[Sqrt[{100, 200, 500}]]] // N
Out[45]= {-44.6802,-88.2591,-218.746}
That is for y values of about 100, the first term
(y - 1) Exp[y] Erfc[Sqrt[y]]
will be about 1/10 the size of the second term
(y/3 - 1) 2 Sqrt[y/Pi] + 1)
But the next decade gives
In[46]:= Log[10, Exp[{1000, 2000, 5000}]] // N
Out[46]= {434.294,868.589,2171.47}
and
In[47]:= Log[10, Erfc[Sqrt[{1000, 2000, 5000}]]] // N
Out[47]= {-436.043,-870.488,-2173.57}
That is the first term decreases about one order of magnitude
while the second one increase by one order of magnitude.
Clearly, the first term can be dropped well before you get to
the maximum value you have for y without any significant error.
The bottom line is if you truly need 70 digits of precision for
this expression with a value of y like 5x10^15 you will need to
have a far more sophisticated approach than what you are using.