RE: Big Loss of Precision

[Montgomery-Smith, Stephen]

(Sorry for top posting, but my email client encourages it.)

In the example you give at the bottom of this email, the numbers created happen 
to be exactly representable in double precision floating point.  (That is, 
assuming the result 0 was what you also had.)

Try using A = 10*rand(1,3) instead, and you get (probably) similar errors as in 
your first example.  Well, I did.


________________________________________
From: Help-octave address@hidden on behalf of Thomas D. Dean address@hidden
Sent: Saturday, July 16, 2016 4:11 PM
To: address@hidden
Subject: Re: Big Loss of Precision

On 07/16/2016 12:05 PM, Montgomery-Smith, Stephen wrote:
> R is of the same order of magnitude as 1e4.  So dot(R,cross(R,V)) is of the 
> same order of magnitude as 1e8.  1e8 times octave eps is about 1e-8.
> ________________________________________
> From: Help-octave address@hidden on behalf of Thomas D. Dean address@hidden
> Sent: Saturday, July 16, 2016 1:35 PM
> To: address@hidden
> Subject: Big Loss of Precision
>
> The cross product of two vectors is orthogonal to both.  The dot product
> of two orthogonal vectors is zero.
>
> This is not an octave problem, just an ill-formed situation.
>
> octave:1405> printf(" %.20f\n",V);
>    3.99120302628332668249
>    -0.91842329903575536942
>    1.29246846883179156151
>
> octave:1406> printf(" %.20f\n",R);
>    -1518.35794621729837672319
>    -4017.17627231122196462820
>    4816.75939627815478161210
>
> octave:1407> dot(cross(R,V),R)
> ans =    1.4901e-08
> octave:1408> dot(cross(R,V),V)
> ans =    3.6380e-12
>
> cross has 6 multiplications and 3 add/subtract
> dot has 3 multiplications and 3 add/subtract
>
> Octave eps is 2.2204e-16
>
> machine eps is 1.11022e-16 (quadmath 1.925930e-34).
>
> This case has a very serious loss of precision!
>
> Tom Dean

A similar problem with similar magnitudes.  Very different outcome.
A = [2 4 8]; norm(A)
B = A.^4; norm(B)
C=cross(A,B)
dot(C,A)
dot(C,B)

tomdean

_______________________________________________
Help-octave mailing list
address@hidden
https://lists.gnu.org/mailman/listinfo/help-octave