       Re: Numerical Problem

• To: mathgroup at smc.vnet.net
• Subject: [mg107044] Re: Numerical Problem
• From: Richard Fateman <fateman at cs.berkeley.edu>
• Date: Mon, 1 Feb 2010 06:13:18 -0500 (EST)
• References: <201001291249.HAA29176@smc.vnet.net> <4B6318FC.3020402@wolfram.com> <hk17nr\$oja\$1@smc.vnet.net>

```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

```

• Prev by Date: Re: Re: Re: position of sequence of
• Next by Date: Re: position of sequence of numbers in list
• Previous by thread: Re: Re: Re: position of sequence of
• Next by thread: Re: Re: Numerical Problem