Re: //N bug, but WHY?
- To: mathgroup at smc.vnet.net
- Subject: [mg58655] Re: //N bug, but WHY?
- From: Jean-Marc Gulliet <jeanmarc.gulliet at gmail.com>
- Date: Tue, 12 Jul 2005 05:21:38 -0400 (EDT)
- Organization: The Open University, Milton Keynes, England
- References: <data3n$mec$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
symbio wrote:
> Evaluating (using //N) two exact same expressions, gives wrong answer unless
> fullsimplify is used first, I spent 2 days on a problem thinking my answer
> was wrong, but turned out Mathematica 5 was giving me buggy answers, I
> debugged it to this point, but WHY in the world is this happening? Please
> help!!!
>
> cut and paste below to see the test case:
>
> In[243]:=
> \!\(\(Cosh[\(43\ \[Pi]\)\/\@2] + \((1 - Cosh[43\ \@2\ \[Pi]])\)\ Csch[
> 43\ \@2\ \[Pi]]\ Sinh[\(43\ \[Pi]\)\/\@2] // FullSimplify\) //
> N\[IndentingNewLine]
> Cosh[\(43\ \[Pi]\)\/\@2] + \((1 - Cosh[43\ \@2\ \[Pi]])\)\ Csch[
> 43\ \@2\ \[Pi]]\ Sinh[\(43\ \[Pi]\)\/\@2] // N\)
> Out[243]=
> \!\(6.551787517854307`*^-42\)
> Out[244]=
> \!\(\(-1.9342813113834067`*^25\)\)
No bug and nothing magic here; just the misunderstanding of what
machine-precision arithmetic is and its pitfalls in some circumstances.
In[1]:=
expr = Cosh[(43*Pi)/Sqrt[2]] + (1 - Cosh[43*Sqrt[2]*Pi])*
Csch[43*Sqrt[2]*Pi]*Sinh[(43*Pi)/Sqrt[2]]
Out[1]=
Cosh[(43*Pi)/Sqrt[2]] + (1 - Cosh[43*Sqrt[2]*Pi])*Csch[43*Sqrt[2]*Pi]*
Sinh[(43*Pi)/Sqrt[2]]
In[2]:=
exprfs = FullSimplify[expr]
Out[2]=
Sech[(43*Pi)/Sqrt[2]]
Both expressions do not seem to be "two exact same expressions." Let try
to expand the first one to see what happens:
In[3]:=
TrigExpand[expr]
Out[3]=
1/(Cosh[Pi/Sqrt[2]]^43 + 903*Cosh[Pi/Sqrt[2]]^41*Sinh[Pi/Sqrt[2]]^2 +
123410*Cosh[Pi/Sqrt[2]]^39*Sinh[Pi/Sqrt[2]]^4 +
6096454*Cosh[Pi/Sqrt[2]]^37*Sinh[Pi/Sqrt[2]]^6 +
145008513*Cosh[Pi/Sqrt[2]]^35*Sinh[Pi/Sqrt[2]]^8 +
1917334783*Cosh[Pi/Sqrt[2]]^33*Sinh[Pi/Sqrt[2]]^10 +
15338678264*Cosh[Pi/Sqrt[2]]^31*Sinh[Pi/Sqrt[2]]^12 +
78378960360*Cosh[Pi/Sqrt[2]]^29*Sinh[Pi/Sqrt[2]]^14 +
265182149218*Cosh[Pi/Sqrt[2]]^27*Sinh[Pi/Sqrt[2]]^16 +
608359048206*Cosh[Pi/Sqrt[2]]^25*Sinh[Pi/Sqrt[2]]^18 +
960566918220*Cosh[Pi/Sqrt[2]]^23*Sinh[Pi/Sqrt[2]]^20 +
1052049481860*Cosh[Pi/Sqrt[2]]^21*Sinh[Pi/Sqrt[2]]^22 +
800472431850*Cosh[Pi/Sqrt[2]]^19*Sinh[Pi/Sqrt[2]]^24 +
421171648758*Cosh[Pi/Sqrt[2]]^17*Sinh[Pi/Sqrt[2]]^26 +
151532656696*Cosh[Pi/Sqrt[2]]^15*Sinh[Pi/Sqrt[2]]^28 +
36576848168*Cosh[Pi/Sqrt[2]]^13*Sinh[Pi/Sqrt[2]]^30 +
5752004349*Cosh[Pi/Sqrt[2]]^11*Sinh[Pi/Sqrt[2]]^32 +
563921995*Cosh[Pi/Sqrt[2]]^9*Sinh[Pi/Sqrt[2]]^34 +
32224114*Cosh[Pi/Sqrt[2]]^7*Sinh[Pi/Sqrt[2]]^36 +
962598*Cosh[Pi/Sqrt[2]]^5*Sinh[Pi/Sqrt[2]]^38 +
12341*Cosh[Pi/Sqrt[2]]^3*Sinh[Pi/Sqrt[2]]^40 +
43*Cosh[Pi/Sqrt[2]]*Sinh[Pi/Sqrt[2]]^42)
Wow! Things are getting worse :-) except if the second expression can be
expanded in a similar way. Let's check
In[4]:=
TrigExpand[expr] == TrigExpand[exprfs]
Out[4]=
True
Good. Now let's check that the expressions by themselves are equal:
In[5]:=
expr == exprfs
From In[5]:=
\!\(\*
RowBox[{\(N::"meprec"\),
":", "\<\"Internal precision
limit $MaxExtraPrecision = \\!\\(49.99999999999999`\\) reached
while \
evaluating \\!\\(\\(\\(Cosh[\\(\\(\\(43\\\\ Ï?\\)\\/\\@2\\)\\)]\\)\\) - \
\\(\\(Sech[\\(\\(\\(43\\\\ Ï?\\)\\/\\@2\\)\\)]\\)\\) + \\(\\(\\(\\((1 - \
\\(\\(Cosh[\\(\\(43\\\\ \\@2\\\\ Ï?\\)\\)]\\)\\))\\)\\)\\\\ \
\\(\\(Csch[\\(\\(43\\\\ \\@2\\\\ Ï?\\)\\)]\\)\\)\\\\ \
\\(\\(Sinh[\\(\\(\\(43\\\\ Ï?\\)\\/\\@2\\)\\)]\\)\\)\\)\\)\\). \
\\!\\(\\*ButtonBox[\\\"Moreâ?¦\\\", ButtonStyle->\\\"RefGuideLinkText\\\", \
ButtonFrame->None, ButtonData:>\\\"General::meprec\\\"]\\)\"\>"}]\)
Out[5]=
Cosh[(43*Pi)/Sqrt[2]] + (1 - Cosh[43*Sqrt[2]*Pi])*Csch[43*Sqrt[2]*Pi]*
Sinh[(43*Pi)/Sqrt[2]] == Sech[(43*Pi)/Sqrt[2]]
Okay. Mathematica is not able to tell us whether both expressions are
equal; however, the warning message gives us a crucial clue. Thus we now
try to use 50 digits, say.
In[6]:=
N[expr, 50]
From In[6]:=
\!\(\*
RowBox[{\(N::"meprec"\),
":", "\<\"Internal precision
limit $MaxExtraPrecision = \\!\\(49.99999999999999`\\) reached
while \
evaluating \\!\\(\\(\\(Cosh[\\(\\(\\(43\\\\ Ï?\\)\\/\\@2\\)\\)]\\)\\) +
\\(\\(\
\\(\\((1 - \\(\\(Cosh[\\(\\(43\\\\ \\@2\\\\ Ï?\\)\\)]\\)\\))\\)\\)\\\\ \
\\(\\(Csch[\\(\\(43\\\\ \\@2\\\\ Ï?\\)\\)]\\)\\)\\\\ \
\\(\\(Sinh[\\(\\(\\(43\\\\ Ï?\\)\\/\\@2\\)\\)]\\)\\)\\)\\)\\). \
\\!\\(\\*ButtonBox[\\\"Moreâ?¦\\\", ButtonStyle->\\\"RefGuideLinkText\\\", \
ButtonFrame->None, ButtonData:>\\\"General::meprec\\\"]\\)\"\>"}]\)
Out[6]=
6.551787517854344014050058`14.31640268976547*^-42
We still get a warning message but things are getting definitely better.
Having a look at the specified documentation, we try
In[7]:=
Block[{$MaxExtraPrecision = 1000}, N[expr, 50]]
Out[7]=
6.55178751785434401405005796823578794510303347370361442804220253`50.*^\
-42
In[8]:=
N[exprfs, 50]
Out[8]=
6.55178751785434401405005796823578794510303347370361442804220253`50.*^\
-42
and we can see that both values are equal (up to 50 digits in our
example). We have just done is using arbitrary-precision arithmetic
(safe but "slow") rather than machine-size precision arithmetic (fast
but sometimes "surprising", hardware dependent).
Regards,
/J.M.