Prev: ListLinePlot Labels
Next: Journals dying?, apparently rather slowly (was
From: Tony Harker on 29 Jan 2010 07:49
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 pair
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.99442,-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.899951,
-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.00000000000
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: Tony Harker on 30 Jan 2010 07:12


]-> -----Original Message-----
]-> From: Daniel Lichtblau [mailto:danl(a)wolfram.com]
]-> Sent: 29 January 2010 17:21
]-> To: Tony Harker
]-> Cc: mathgroup(a)smc.vnet.net
]-> Subject: Re: Numerical Problem
]->
]-> Tony Harker 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 pair
]-> > 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.7412
]-> 17},{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
]-> .99442,-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.7412
]-> 18},{1.3416,-0.6
]-> >
]-> 70802},{1.21415,-0.607076},{1.09881,-0.549404},{0.99442,-0.4
]-> 9721},{0.899951,
]-> > -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.0000000
0000000000,{2.00
]-> >
]-> 000000000000000,-1.00000000000000000},{{998.000000000000000,
]-> 1998.00000000000
]-> > 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.*1
0^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
]->
]-> It's not the Mathematica version that matters, it's the
]-> processor or
]-> perhaps MKL library. On one machine I get the "bad" results in both
]-> versions, and the "good" results in both versions on
]-> another amchine.

Ah, that's interesting. I was getting different results from the different
versions on the same machine. This suggests to me that the change is in the
MKL library. But, in that all that's being done is multiplication and
addition, this seems like a rather fundamental change.

]->
]-> Here is a pared down example. First the bad result.
]->
]-> epcm[h_, y_List, f_List] := y + N[1/2* h (f.y + f.(y + h f.y))]
]->
]-> app[0] = {2, -1};
]->
]-> Do[Print[InputForm[
]-> app[j] = epcm[.1, app[j - 1], {{998, 1998}, {-999,
]-> -1999}}]]], {j,
]-> 1, 6}]
]->
]-> {1.8100000000000023, -0.9049999999999955}
]-> {1.638049999944292, -0.8190249999442927}
]-> {1.4824349769833363, -0.7412173519833346}
]-> {1.3402658465670585, -0.6694638959420446}
]-> {-5.3436544706007, 5.9507302359163425}
]-> {-32138.70840490274, 32139.257808470356}
]->
]-> Here is the more desired result:
]->
]-> In[1]:= epcm[h_, y_List, f_List] := y + N[1/2* h (f.y +
]-> f.(y + h f.y))]
]->
]-> In[2]:= app[0] = {2, -1};
]->
]-> In[3]:= Do[Print[InputForm[
]-> app[j] = epcm[.1, app[j - 1], {{998, 1998}, {-999,
]-> -1999}}]]], {j,
]-> 1, 6}]
]-> {1.81, -0.905}
]-> {1.63805, -0.819025}
]-> {1.48243525, -0.741217625}
]-> {1.34160390125, -0.670801950625}
]-> {1.21415153063125, -0.607075765315625}
]-> {1.0988071352212814, -0.5494035676106407}
]->
]-> In[4]:= $Version
]-> Out[4]= 8.0 for Linux x86 (32-bit) (January 28, 2010)
]->
]-> This is clearly a very ill-conditioned iteration (small changes to
]-> input, on either machine, lead to relatively larger changes
]-> in output).
]-> There may also be bad library code somewhere. To see it is ill
]-> conditioned (other than observing what is already shown),
]-> one can take
]-> derivatives with respect to the y inputs (the ones that are
]-> changing in
]-> the iteration). For simplicity we preset fh to 1/10 and f
]-> to {-999, -1999}
]->
]-> In[27]:= symbepcm[y : {y1_, y2_}] =
]-> y + 1/20*{-999, -1999}.y + {-999, -1999}.(y + 1/10*{-999,
]-> -1999}.y)
]-> Out[27]= y +
]-> 1/20 {-999, -1999}.y + {-999, -1999}.(y + 1/10 {-999, -1999}.y)
]->
]-> In[33]:= D[symbepcm[{y1, y2}], y1] // N
]-> Out[33]= {298452., 298451.}
]->
]-> In[34]:= D[symbepcm[{y1, y2}], y2] // N
]-> Out[34]= {597201., 597202.}
]->
]-> These indicate that a change in input can be magnified
]-> several orders of
]-> magnitude.

Indeed, this is to be expected, as this represents using a na=EFve numerical
scheme to solve a stiff pair of first-order differential equations in which
the required solution goes like exp(-t) and the parasitic solution is
exp(-1000 t).

]->
]-> One thing that seems to help stabilize the iteration is to
]-> use h=1/10
]-> rather than 0.1 in calling epcm[]. I doubt this is really a
]-> "fix" since,
]-> as noted, this iteration is ill conditioned.
]->

I suppose this is picking up where in my integer-only version working in
integers is having the greatest effect.

Tony


]-> Daniel Lichtblau
]-> Wolfram Research
]->
]->
]->


From: Richard Fateman on 1 Feb 2010 06:13
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

From: Tony Harker on 2 Feb 2010 03:23


For interest, I've tried the same iterations in Excel and in
Fortran (double precision), on the same Intel cpu as that on which I was
running the Mathematica cases. They both give results which are consistent
with Mathematica version 7. It looks as if Mathematica version 5 might have
had some additional subtlety which took it beyond IEEE standard floating
point arithmetic, but that this has been stripped out of version 7.

You're right, of course, that the loss of significance in the iterates is
to be expected, and in this sense the version 7 results are more in keeping
with what we expect from limited precision floating point calculations than
are the version 5 results.

Tony

]-> -----Original Message-----
]-> From: Mark McClure [mailto:mcmcclur(a)unca.edu]
]-> Sent: 01 February 2010 12:49
]-> To: mathgroup(a)smc.vnet.net; fateman(a)cs.berkeley.edu
]-> Cc: Daniel Lichtblau; a.harker(a)ucl.ac.uk
]-> Subject: Re: Re: Numerical Problem
]->
]-> 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: Richard Fateman on 2 Feb 2010 03:24
Mark McClure wrote:
> 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.
>
I think it is not hardware, but that version 5 somehow used higher
precision. It may be that
it used extended-double, supported partly in hardware.
The results on version 7 are the same as you get from executing the
specified calculations
in double precision using other programming languages.


....
>
> 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.
Yes
> I would also expect the iterates to
> lose significance, which is exactly what we see.
Yes, and eventually Mathematica produces nothing but noise.

The exact answer, written out to 40 decimal digits is
([7.3708196966710360351173828125e-1],[-3.68540984833551801755869140625e-1])

In fact, carrying 40 digits of precision is enough to compute this result as
([7.370783.....e-1],[-3.68537...e-1])]

> Seems to be a good
> example of significance arithmetic doing exactly what we want.
>

Mathematica's version of arithmetic thinks that it has lost 5 or 6 extra
decimal digits. (In fact claims to have
lost all of them.) So in one sense if the user gets to look at the
output and notices 0.10^-1 or some such thing,
and knows how to interpret that, yes, that would be good; significance
arithmetic has given a warning, albeit
overly cautious. (Other ways of estimating or bounding the accumulated
error in the iteration would require more thought.)

But if the user does not look at the output but merely has a program
compare this result to
0, then the result would be, why yes! your answer IS zero! or possibly
why yes! your answer is 1! Maybe even both!
This is a hazard, and I wonder if this is exactly what people want.

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.



RJF


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