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