Re: Re: Numerical Problem
- To: mathgroup at smc.vnet.net
- Subject: [mg107115] Re: [mg107044] Re: Numerical Problem
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Wed, 3 Feb 2010 06:09:49 -0500 (EST)
- References: <201001291249.HAA29176@smc.vnet.net> <4B6318FC.3020402@wolfram.com> <hk17nr$oja$1@smc.vnet.net> <201002011113.GAA22665@smc.vnet.net>
Richard Fateman wrote:
> Thanks, Tony, for providing a nice example of how arithmetic matters.
>
> Your 3-line program got pretty plausible results in version 5.0
>
> Clear[epcm,fepcm]
> epcm[h_,y_List,f_List]:=y+0.5 h (f.y+f.(y+h f.y))
> fepcm[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcm[h,#,f]&,y0,10]
> app1=fepcm[.1,0.,1,{2.,-1.},{{998,1998},{-999,-1999}}]
>
> But apparently not in version 7, where DanL suggests that Intel's
> numerical library may have sent it off in the wrong direction.
>
> Let's try using Mathematica's significance arithmetic, and do everything
> "exactly" except for one number, 0.5. Let us use N[1/2,40]. That is,
> 40 decimal digits or so. Observe that 1/2 is exactly representable in
> binary, so that SetPrecision[ [N[1/2,40], 100] is
>
> 0.50000000000000000000000000000000000000000000000000000000000000000000\
> 00000000000000000000000000000000
>
> We change the program then to this, below, using software high-precision
> and therefore significance arithmetic.
>
>
> Clear[epcms,fepcms] (*software floats *)
> epcms[h_,y_List,f_List]:= y+N[1/2,40]*h (f.y+f.(y+h f.y))
> fepcms[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcms[h,#,f]&,y0,10]
> app1=fepcms[1/10,0,1,{2,-1},{{998,1998},{-999,-1999}}]
>
>
> If you look at the trailing few numbers using InputFormat you find that
> for some of them, the numbers that approximate the exactly correct
> answer for about 18 decimal digits are believe wrong after 3 or 4
> decimal digits. And the last two numbers, which print as 0.x10^1 are,
> how can we put this delicately... wrong?
>
> This is using 40 digits -- an extra 22 decimal digits compared to
> machine precision. But machine precision on Mathematica 5.2 apparently
> gets the numbers pretty much right. I am curious to see how version 7
> can be using Intel's library to get different answers.
>
> NOte that if you carry even MORE digits, the problem eventually goes
> away, but for that to work you would have to be looking for this
> problem to appear, and you might not notice it if the results of this
> program were fed in to yet another program.
>
> RJF
Somebody's library has a big effect on it. As for why a library code
would be invoked at all, I think this is how Dot with machine numbers is
handled in the Mathematica kernel (that is, by using e.g. MKL BLAS).
Anyway, I can now emulate the "good" behavior in a way that explains
what is happening. It is, as RJF suggests, the use of extended precision
arithmetic in intermediate computations that makes the difference.
In brief, some platforms support, and some machine double libraries use,
higher precision register arithmetic. Specifically, machine doubles
support mantissas of 53 bits, and extended precision registers use 64,
via 80 bit registers (the other 16 bits are for sign and exponent). But
the results are eventually stored in 64 bit memory, hence the mantissas
get altered (and on't ask what happens when the exponents don't fit
their 11 bit allotments...).
It took me a while to get the emulation correct because I was truncating
the 64 bit results, and zero-padding the emulated "true" doubles (the
things with 53 bits). The latter was appropriate, but the former really
requires rounding, so I was getting junk in the last bits of the
doubles, and the ill conditioning magnified it, and...
Anyway, here it is.
SetAttributes[emulateFloat, Listable];
SetAttributes[emulateExtended, Listable];
emulateFloat[rr_Real] :=
Module[{rf = RealDigits[SetPrecision[rr, $MachinePrecision], 2]},
rf[[1]] = PadRight[rf[[1]], 53]; (*probably no longer necessary*)
Sign[rr]*N[FromDigits[rf, 2], $MachinePrecision]]
emulateExtended[rr_Real] :=
Module[{rf = RealDigits[rr, 2]}, rf[[1]] = PadRight[rf[[1]], 64];
Sign[rr]*N[FromDigits[rf, 2], 64*Log[10., 2]]]
The iteration, with extended precision emulated, involves bumping floats
to 64-bit-mantissa bignums, doing the operation, and forcing
$MachinePrecision digits (converted to bits, that is, 53) to be retained.
epcm[h_, y_List, f_List] := y + 1/2*h*(f.y + f.(y + h*f.y))
app[0] = {2., -1.};
mat = {{998, 1998}, {-999, -1999}};
Do[Print[InputForm[
app[j] =
emulateFloat[epcm[1/10, emulateExtended[app[j - 1]], mat]]]], {j,
1, 16}]
Do[Print[app[j] =
emulateFloat[epcm[1/10, emulateExtended[app[j - 1]], mat]]], {j, 1,
16}]
{1.810000000000000,-0.905000000000000}
{1.638050000000000,-0.8190250000000000}
{1.482435250000000,-0.7412176250000000}
{1.341603901250000,-0.6708019506250000}
{1.214151530631250,-0.6070757653156250}
{1.098807135221281,-0.5494035676106406}
{0.994420457375259,-0.4972102286876297}
{0.8999505139246098,-0.4499752569623049}
{0.8144552151017719,-0.4072276075508859}
{0.7370819696671036,-0.3685409848335518}
{0.6670591825487288,-0.3335295912743644}
{0.6036885602065996,-0.3018442801032998}
{0.5463381469869726,-0.2731690734934863}
{0.4944360230232102,-0.2472180115116051}
{0.4474646008360053,-0.2237323004180026}
{0.4049554637565848,-0.2024777318782924}
One might ask whether this "magic" is really worth the bother. Maybe it
is but I do not know that it is all that valuable for this particular
problem. As Steve Schochet pointed out in a different response, changing
the input by even an ULP will cause the instability to manifest in the
iterations, and that holds true regardless of whether or not extended
precision registers are used. Numeric ill-conditioning trumps clever
arithmetic every time.
Daniel Lichtblau
Wolfram Research
- References:
- Re: Numerical Problem
- From: Richard Fateman <fateman@cs.berkeley.edu>
- Re: Numerical Problem