Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: Numerical Problem

• To: mathgroup at smc.vnet.net
• Subject: [mg107120] Re: [mg107044] Re: Numerical Problem
• From: Richard Fateman <fateman at cs.berkeley.edu>
• Date: Wed, 3 Feb 2010 06:10:43 -0500 (EST)
• References: <201001291249.HAA29176@smc.vnet.net> <4B6318FC.3020402@wolfram.com> <hk17nr$oja$1@smc.vnet.net> <201002011113.GAA22665@smc.vnet.net> <4B686AE7.5030206@wolfram.com>

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

• References:
  • Re: Numerical Problem
    • From: Richard Fateman <fateman@cs.berkeley.edu>