Re: Re: Numerical Problem
- To: mathgroup at smc.vnet.net
- Subject: [mg107091] Re: [mg107044] Re: Numerical Problem
- From: Mark McClure <mcmcclur at unca.edu>
- Date: Tue, 2 Feb 2010 03:30:24 -0500 (EST)
- References: <201001291249.HAA29176@smc.vnet.net> <4B6318FC.3020402@wolfram.com>
On Mon, Feb 1, 2010 at 6:13 AM, Richard Fateman <fateman at 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