Round-off: cross- vs. dot-products

160 Views Asked by At

I have three vectors

e0 = numpy.array([1.0, 0.0,  0.0])
e1 = numpy.array([-0.5, 1.0e-4, 0.0])
e2 = numpy.array([0.5, 1.0e-4, 0.0])

(Edges that form a very flat triangle.) I need to compute $\langle e_2\times e_0, e_1\times e_2\rangle$, and

import numpy

e0 = numpy.array([1.0, 0.0,  0.0])
e1 = numpy.array([-0.5, 1.0e-4, 0.0])
e2 = numpy.array([0.5, 1.0e-4, 0.0])

# Cross-product formulation:
e2_cross_e0 = numpy.cross(e2, e0)
e1_cross_e2 = numpy.cross(e1, e2)

for x in e2_cross_e0:
    print('%0.20e' % x)
print('.dot.')
for x in e1_cross_e2:
    print('%0.20e' % x)
print

# perfect!
print('%0.20e' % numpy.dot(e2_cross_e0, e1_cross_e2))
print('------------------')
print

does that to a high precision. Nice!

In my case, I can compute dot-products a lot faster, though; and cross-product can be formulated in terms of a dot-product too (via $\langle a \times b, c \times d\rangle = \langle a, c\rangle \langle b, d\rangle - \langle a, d\rangle \langle b, c\rangle$):

# Dot-product formulation:
e0_dot_e0 = numpy.dot(e0, e0)
e0_dot_e1 = numpy.dot(e0, e1)
e0_dot_e2 = numpy.dot(e0, e2)
e1_dot_e2 = numpy.dot(e1, e2)

print('%0.20e * %0.20e - %0.20e * %0.20e' %
      (e0_dot_e0, e1_dot_e2, e0_dot_e1, e0_dot_e2)
      )
print

alpha = e0_dot_e0 * e1_dot_e2 - e0_dot_e1 * e0_dot_e2
print('%0.20e' % alpha)
print('------------------')

Here, however, the round-off is much larger since we're taking the difference of two almost equal values. Yikes!

Any idea on how to avoid this?