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

```

• Prev by Date: Is this a bug? FromAdjacencyMatrix[]
• Next by Date: Re: Re: What does & mean?
• Previous by thread: Re: Re: Numerical Problem
• Next by thread: Re: Re: Numerical Problem