Prev: ListLinePlot Labels
Next: Journals dying?, apparently rather slowly (was
From: schochet123 on 2 Feb 2010 03:25
Actually, the behavior obtained from version 7 is in some sense the
"correct" behavior, and result in version 5.2 is the anomaly. As you
note, the ODE system is stiff, and as Daniel notes, the approximation
scheme is very ill-conditioned. For example, changing the initial data
from {2,-1} to {2,-1}+{0,10^{-16}} yields a problem whose exact
solution (without any round-off) blows up rapidly like the solution
obtained in version 7. The solution you obtained in version 5.2, which
I get in version 6.0, and Daniel gets on some machines and not others
for all versions, happens because the numerical calculation somehow
manages to keep the numerical approximation to be exactly a multiple
of the one eigenvector {-2,1} for which such blow-up does not occur in
the exact solution of the approximation scheme. This perfect accuracy
for the ratio of the two components of the solution is not expected
and should not be relied upon.

Looking into the internals of how Mathematica implements the
approximation scheme numerically is a bad idea, for several reasons:
1) Whatever fix one might find may not work on future versions, other
machines, or other problems. 2) The benefits of optimizations
Mathematica has made in organizing computations may be lost. 3) It is
too much work. In particular, simply calculating with more digits of
accuracy as mentioned by RJF should be expected to merely delay the
problem for a few time steps, not eliminate it.

The correct approach is to use an approximation scheme that is well-
conditioned. Although this could be accomplished by reducing h
drastically, a more efficient method is to switch to to an implicit
scheme or another type of scheme that is well-conditioned for stiff
problems. For the theory involved, look up stiff problems in any text
on numerical solutions of ODEs. Mathematica knows how to do this
automatically. Indeed, even in version 7,

thex[t_]=x[t]/. NDSolve[{x'[t] == {{998, 1998}, {-999, -1999}} . x
[t],
x[0] == {2, -1}}, x, {t, 0, 1}][[1]]

yields an accurate solution.

Steve


On Jan 29, 2:49 pm, "Tony Harker" <a.har...(a)ucl.ac.uk> wrote:
> I expect I'm missing something in the documentation, but I can't see why
> this simple example of Euler's predictor-corrector applied to a stiff pai=
r
> of equations should behave differently in Version 7 and Version 5.
> =========================
==========================
=======================
> $Version
> 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}}]
>
> 7.0 for Microsoft Windows (32-bit) (February 18, 2009)
>
> {{2.,-1.},{1.81,-0.905},{1.63805,-0.819025},{1.48243,-0.741217},{1.34027,=
-0.
> 669464},{-5.34365,5.95073},{-32138.7,32139.3},{-1.57517*10^8,1.57517*10^8=
},{
> -7.71992*10^11,7.71992*10^11},{-3.78353*10^15,3.78353*10^15},{-1.85431*10=
^19
> ,1.85431*10^19}}
> =========================
==========================
==========================
==
> $Version
> 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}}]
>
> 5.2 for Microsoft Windows (June 20, 2005)
>
> {{2.,-1.},{1.81,-0.905},{1.63805,-0.819025},{1.48244,-0.741218},{
> 1.3416,-0.670802},{1.21415,-0.607076},{1.09881,-0.549404},{0.9944=
2,-0.\
> 49721},{0.899951,-0.449975},{0.814455,-0.407228},{0.737082,-0.368541}}
> =========================
==========================
==========================
==
> It's clear that the version 5 result is more acceptable, as we can do
> everything in integers and convert to reals at the end:
>
> $Version
> epcm[h_,y_List,f_List]:=y+ h (f.y+f.(y+h f.y))/2
> fepcm[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcm[h,#,f]&,y0,10]
> N[app1=fepcm[1/10,0,1,{2,-1},{{998,1998},{-999,-1999}}]]
>
> 7.0 for Microsoft Windows (32-bit) (February 18, 2009)
>
> {{2.,-1.},{1.81,-0.905},{1.63805,-0.819025},{1.48244,-0.741218},{1.3416,-=
0.6
> 70802},{1.21415,-0.607076},{1.09881,-0.549404},{0.99442,-0.49721},{0.8999=
51,
> -0.449975},{0.814455,-0.407228},{0.737082,-0.368541}}
> =========================
==========================
==========================
==
>
> What I don't understand is why this discrepancy between Versions 5 and 7
> occurs, and whether I can get version 7 to behave in the same way as 5.
>
> For something even more exciting, there's
>
> =========================
==========================
==========================
==
>
> $Version
> epcm[h_,y_List,f_List]:=y+0.500000000000000000 h (f.y+f.(y+h f.y))
> fepcm[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcm[h,#,f]&,y0,10]
> app1=fepcm[.100000000000000000,0.00000000000000000,1.00000000000000000,=
{2.00
> 000000000000000,-1.00000000000000000},{{998.000000000000000,1998.00000000=
000
> 0000},{-999.00000000000000000,-1999.00000000000000000}}]
>
> 7.0 for Microsoft Windows (32-bit) (February 18, 2009)
>
> {{2.0000000000000000,-1.0000000000000000},{1.810000000000,-0.905000000000=
},{
> 1.6380500,-0.8190250},{1.482,-0.741},{0.*10^1,0.*10^1},{0.*10^6,0.*10^6},=
{0.
> *10^10,0.*10^10},{0.*10^15,0.*10^15},{0.*10^20,0.*10^20},{0.*10^24,0.*10^=
24}
> ,{0.*10^29,0.*10^29}}
> =========================
==========================
==========================
==
>
> Tony Harker


From: Mark McClure on 2 Feb 2010 03:30
On Mon, Feb 1, 2010 at 6:13 AM, Richard Fateman <fateman(a)cs.berkeley.edu> wrote:
> Thanks, Tony, for providing a nice example of how arithmetic matters.
> Your 3-line program got pretty plausible results in version 5.0
> ...
> But apparently not in version 7, where DanL suggests that Intel's
> numerical library may have sent it off in the wrong direction.

I cannot reproduce these results on my MacBook Pro, i.e. I get (on V7
and V5.2) the same erroneous results that Tony gets with V7. I
suppose this supports Dan's suggestion that this is a hardware issue.

> 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]. ...

The last few terms in the sequence I get look something like so:
In[1]:= epcms[h_, y_List, f_List] := y +
N[1/2, 40]*h (f.y + f.(y + h f.y));
In[2]:= fepcms[h_, t0_, tmax_, y0_List, f_List] :=
NestList[epcms[h, #, f] &, y0, 10];
In[3]:= app1 = fepcms[1/10, 0, 1, {2, -1},
{{998, 1998}, {-999, -1999}}];
In[4]:= Column[app1[[-4 ;;]]]
Out[4]=
{0.994420457375, -0.497210228688}
{0.89995051, -0.44997526}
{0.814, -0.407}
{0.10^2 , 0.10^2}

I haven't examined the equations closely but, as I understand it, they
arise via Euler's method applied to a stiff system. Thus, I would
expect numerical instability. I would also expect the iterates to
lose significance, which is exactly what we see. Seems to be a good
example of significance arithmetic doing exactly what we want.

Mark McClure

From: Andrzej Kozlowski on 3 Feb 2010 06:07

On 2 Feb 2010, at 09:25, Richard Fateman wrote:

> And again, just doing 40 decimal digits (maybe 135 bit floating-point
> fraction) you can get a pretty good answer,
> even in Mathematica, except that Mathematica will tell you it can't.
> Just insert the occasional SetPrecision.

But there is a standard and easy way to do this in Mathematica, that you have been told about many times, and which does not require inserting "SetPrecision". And 20 digits (rather than 40) is enough.


epcms[h_, y_List, f_List] := y + N[1/2, 20]*h (f.y + f.(y + h f.y))
fepcms[h_, t0_, tmax_, y0_List, f_List] :=
NestList[epcms[h, #, f] &, y0, 10]

app1 = Block[{$MinPrecision = 20, $MaxPrecision = 20},
fepcms[1/10, 0, 1, {2, -1}, {{998, 1998}, {-999, -1999}}]]

Andrzej Kozlowski

From: Daniel Lichtblau on 3 Feb 2010 06:09
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





From: Richard Fateman on 3 Feb 2010 06:10
Daniel Lichtblau wrote:
>
> 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).

I agree with DanL. I was able to get Exactly the same results as
Mathematica 5.2 in Tony Harker's example by using the same program but
doing this.

Using software floats in two different precisions, 53 bits and 64 bits.
I got the same numbers as Tony; It looks to me like DanL's differed
slightly.

According to DanL's prescription,
All arithmetic is done with 53 bit fractions except within "Dot", when
the arithmetic is done with 64 bit fractions, and any results from Dot
are rounded (to nearest even) down to 53 bits.

Now what is interesting is that the system I used needed only 3 lines to
do it, and did it just right.

in Mathematica pseudo-code it would look something like this ..

ThePrecision=53;
Define "~" to have the same syntax precedence as "."
(a_~b_) := Block[{ ThePrecision=64}, a . b)

We also had to rewrite one expression in the program to be this ..

y+ h/2*(f~y+f~(y+(h*f)~y)) (* use ~ instead of "." *)

Now this doesn't detract from DanL's main point that this computation is
ill-conditioned and one should not expect
high-quality results from it.

But it does suggest that doing careful arithmetic in Mathematica
requires a fair amount of work and expertise
to UNDO the features of significance arithmetic.

DanL's program is about 11 lines of pretty non-obvious code.

In case you want to see the 3 line program, email me.
RJF




First  |  Prev  |  Next  |  Last
Pages: 1 2 3
Prev: ListLinePlot Labels
Next: Journals dying?, apparently rather slowly (was