Re: Numerical precision problem
- To: mathgroup at smc.vnet.net
- Subject: [mg43059] Re: [mg43016] Numerical precision problem
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sun, 10 Aug 2003 01:46:48 -0400 (EDT)
- References: <200308080426.AAA05607@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
"Gareth J. Russell" wrote:
>
> 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
Presumably fDOverlap is the same as f.
You might want to "force" high precision arithmetic. A quick and dirty
way is as below.
fhighprecision[r_, {t1_, t2_}, {d1_, d2_}, prec_] := With[
{bigr=Rationalize[r,0], bigt1=Rationalize[t1,0],
bigt2=Rationalize[t2,0], bigd1=Rationalize[d1,0],
bigd2=Rationalize[d2,0]},
N[f[bigr, {bigt1,bigt2}, {bigd1,bigd2}], prec]
]
In[14]:= InputForm[fhighprecision[1/1000000000, {2, 4}, {2, 4}, 25]]
Out[14]//InputForm= 3.33333333166666666733333333311111111168`25.*^-10
If you know a priori certain ranges of the inputs wherein the machine
arithmetic result is trustworthy, you can use that information to bypass
the higher precision approach, as machine arithmetic tends to be faster.
This is of interest if you need to do this sort of computation a million
times, say.
Daniel Lichtblau
Wolfram Research
- References:
- Numerical precision problem
- From: "Gareth J. Russell" <gjr2008@columbia.edu>
- Numerical precision problem