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?