Re: Numerical precision problem
- To: mathgroup at smc.vnet.net
- Subject: [mg43041] Re: Numerical precision problem
- From: "Christos Argyropoulos M.D." <charg at med.upatras.gr>
- Date: Sat, 9 Aug 2003 02:57:41 -0400 (EDT)
- References: <bgv9nf$5ln$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi Gareth !
Sound like you came across a numerical monster (Essex, C., M. Davison, and
C. Schulzky. 2000. Numerical Monsters. SIGSAM Bulletin 134 (Dec.): 16-32. )
i.e. functions
"guaranteed to wreak numerical mayhem across both software packages and
hardware platforms". Your problems boil down to the infamous MachineEpsilon
value and if you have access to the paper you can find out more about the
pecularities of functions that are given as a product of a monotonically
increasing function AND a monotonically decreasing function that depends on
1+MachineEpsilon. Compared to the problems that the 3 authors considered
yours is more interesting,
in the sense that the irregular behaviour only shows up in a limited range
of values.
To see that change {y,0,10} to {y,0,30}in the Mathematica notebook given
below.
There is a Mathematica-way out; cut and paste the following text
corresponding to mathematica notebook and evaluate
Best Regards
Christos Argyropoulos
Notebook[{
Cell[CellGroupData[{
Cell["Original Function : ", "Section"],
Cell[BoxData[
StyleBox[\(\(f[
r_, \ {t1_, \ t2_}, \ {d1_, \
d2_}]\ := \ \((2*E^\((d1*r)\)\ - \
2*E^\((d2*r)\)\ - \ \((d1\ - \ \nd2)\)*E^\((r*t2)\)*
r*\((2\ + \
r*\((d1\ + \ d2\ - \ 2*t2)\))\))\)/\((2*\((d1\ - \
d2)\)*\nE^\((r*t2)\)*
r^2*\((\(-t1\)\ + \ t2)\))\)\)\(\n\)\(\n\)
\),
FormatType->StandardForm]], "Input"]
}, Open ]],
Cell[CellGroupData[{
Cell["\<\
Evaluate Cell to observe wiggly behaviour of the function as we go down to \
zero\
\>", "Section"],
Cell[BoxData[{
\(\(grplot =
Table[Plot[f[x, {2, 4}, {4, 2}], {x, 0, 0.01\ 10^\((\(-y\))\)},
PlotLabel \[Rule] "\<x \[Epsilon] [0,\>" <>
ToString[AccountingForm[0.01\ 10^\((\(-y\))\)]] <> "\<]\>",
DisplayFunction \[Rule] Identity], {y, 0,
9}];\)\), "\[IndentingNewLine]",
\(\(Show[GraphicsArray[Partition[grplot, 3]]];\)\)}], "Input"]
}, Open ]],
Cell[CellGroupData[{
Cell[TextData[{
"We try a number of transformations to arrive at a simpler expression; \na)
\
Evaluate f for a symbolic argument to get :\t \nb) Simplify\t\t\t\t\t\t\t
\t\n\
c) Replace x by y/2\t\t\t\t\t\t\t\nd) Expand \ne) Collect powers of y\nThis
\
transforms the expression f[x,{2,4},{4,2}] = ",
Cell[BoxData[
StyleBox[\(\(\[ExponentialE]\^\(\(-4\)\ x\)\ \((\(-2\)\ \[ExponentialE]\
\^\(2\ x\) + 2\ \[ExponentialE]\^\(4\ x\) -
2\ \[ExponentialE]\^\(4\ x\)\ \((2 -
2\ x)\)\ x)\)\)\/\(8\ x\^2\)\),
FontSize->18,
FontColor->RGBColor[0, 0, 1]]], "Output"],
" \[LongRightArrow] ",
Cell[BoxData[
\(1\/2 + \(1 - \[ExponentialE]\^\(-y\)\)\/y\^2 - 1\/y\)],
FontSize->18,
FontColor->RGBColor[1, 0, 0]],
"\n"
}], "Section",
TextJustification->1],
Cell[BoxData[{
\(f[x, {2, 4}, {4, 2}]\), "\n",
\(Simplify[%]\), "\n",
\(% /. x \[Rule] y/2\), "\n",
\(Expand[%]\), "\n",
\(Collect[%, y]\), "\n",
\(\)}], "Input"]
}, Open ]],
Cell[CellGroupData[{
Cell[TextData[{
"Another transformation y \[LongRightArrow] 1/r (y \[NotEqual] 0) brings \
the original expression to a form that can be analysed according to\n",
StyleBox["Essex, C., M. Davison, and C. Schulzky. 2000. Numerical Monsters.
\
",
FontColor->RGBColor[0, 1, 0],
Background->GrayLevel[0]],
StyleBox["SIGSAM Bulletin",
FontSlant->"Italic",
FontColor->RGBColor[0, 1, 0],
Background->GrayLevel[0]],
StyleBox[" 134 (Dec.): 16-32. \n",
FontColor->RGBColor[0, 1, 0],
Background->GrayLevel[0]],
"i.e. a product of a monotonically increasing function ",
Cell[BoxData[
\(TraditionalForm\`\((r\^2)\)\)]],
"and one ",
Cell[BoxData[
\(\((1 - \[ExponentialE]\^\(\(-1\)/r\))\)\)]],
"that its numerical value as it approaches zero depends on \[Epsilon] , \
the Machine Epsilon value."
}], "Section"],
Cell[BoxData[
\(% /. y \[Rule] 1/r\)], "Input"]
}, Open ]],
Cell[CellGroupData[{
Cell[TextData[{
"Fortunately there is a way out with ",
StyleBox["Mathematica",
FontSlant->"Italic"],
".\na) SetPrecision to a high value (i.e. 100) (in the function definition)
\
\nb) wrap the function with Hold\nc) RelaseHold before plotting the
function\n\
d) Evaluate cells to see the difference !!\n\[GrayCircle] Depending on your
\
application domain you can raise or lower the precision as needed. Note that
\
wiggly behaviour will never go away, but will only be cast away to values \
closer to zero"
}], "Section"],
Cell[BoxData[
\(\(\(fhold[r_, \ {t1_, \ t2_}, \ {d1_, \ d2_}]\ := \ \ Hold[
SetPrecision[\
Unevaluated[\((2*E^\((d1*r)\)\ - \
2*E^\((d2*r)\)\ - \ \((d1\ - \ \nd2)\)*E^\((r*t2)\)*
r*\((2\ + \
r*\((d1\ + \ d2\ - \ 2*t2)\))\))\)/\((2*\((d1\ - \
d2)\)*\nE^\((r*t2)\)*r^2*\((\(-t1\)\ + \ t2)\))\)],
100]]\)\(\n\)
\)\)], "Input"],
Cell[BoxData[{
\(\(grplot2 =
Table[Plot[
ReleaseHold[fhold[x, {2, 4}, {4, 2}]], {x, 0,
0.01\ 10^\((\(-y\))\)},
PlotLabel \[Rule] "\<x \[Epsilon] [0,\>" <>
ToString[AccountingForm[0.01\ 10^\((\(-y\))\)]] <> "\<]\>",
DisplayFunction \[Rule] Identity], {y, 0,
9}];\)\), "\[IndentingNewLine]",
\(\(Show[GraphicsArray[Partition[grplot2, 3]]];\)\)}], "Input"]
}, Open ]]
},
FrontEndVersion->"4.0 for Microsoft Windows",
ScreenRectangle->{{0, 1024}, {0, 693}},
WindowSize->{1016, 666},
WindowMargins->{{0, Automatic}, {Automatic, 0}}
]
----- Original Message -----
From: "Gareth J. Russell" <gjr2008 at columbia.edu>
To: mathgroup at smc.vnet.net
Subject: [mg43041] Numerical precision problem
> Hi All,
>
> I have the following expression (in the form of a function):
>
> f[r_, {t1_, t2_}, {d1_, d2_}] := (2*E^(d1*r) - 2*E^(d2*r) - (d1 -
> d2)*E^(r*t2)*r*(2 + r*(d1 + d2 - 2*t2)))/(2*(d1 - d2)*
> E^(r*t2)*r^2*(-t1 + t2))
>
> r is a rate parameter, with values 0 to infinity. The expression goes to
> 0 in the limit as r goes to 0.
>
> The problem is that as r gets smaller (roughly, smaller than about 0.
> 0001), significant numerical errors appear and then get huge:
>
> In[531]:=
> fDOverlap[1/1000000000., {2, 4}, {2, 4}]
>
> Out[531]=
> 39.8702
>
> The only way to avoid them that I have found is to provide r as an exact
> number, and them simplify the result before getting an approximate
> result
>
> In[525]:=
> fDOverlap[1/1000000000, {2, 4}, {2, 4}]
> Simplify[%]
> % // N
>
> Out[525]=
> (-125000000000000000*(2*E^(1/500000000) - (499999999000000001*E^(1/
> 250000000))/
> 250000000000000000))/E^(1/250000000)
>
> Out[526]=
> 499999999000000001/2 - 250000000000000000/E^(1/500000000)
>
> Out[527]=
> 0.
>
> But of course, I need to be able to use real vales of r.
>
> Any suggestions for how to deal with the precision issue, other than
> arbitrarily setting the value to zero for r less than some constant?
>
> The original expression is as simple as I can make it, but I'm not
> convinced some more terms might not cancel, which could help.
>
> Advice appreciated as always!
>
> Gareth
> Columbia University
>