Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

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

• Follow-Ups:
  • Re: Re: Numerical Problem
    • From: Daniel Lichtblau <danl@wolfram.com>