Back to Algebra of subspaces Cl(3, 0)
# Description: Tutorial code for Cl(3, 0)
if True:
### Initialization
# Import all symbols from the `sympy` and `sympy.galgebra` packages.
from sympy.galgebra import *
from sympy import *
### Defining a basis
# Use the basis e1, e2, e3 for the 1-vectors, and
# choose the symmetric bilinear form so that we
# have Cl_{3, 0}.
e1, e2, e3 = MV.setup('e1 e2 e3', metric = '[1, 1, 1]')
# The e1, e2, e3 on the left are the basis-vectors which
# you use for programming. The e1, e2, e3 on the right
# are the printed names of those basis vectors. They
# need not be the same.
# You can also leave the symmetric bilinear form
# unspecified. This is useful when you want to prove
# that something does not depend on the choice of the
# bilinear form.
print
print 'Algebra of subspaces Cl(3, 0)'
print '-----------------------------'
print
### Grade involution
# Definition:
# Grade-involution is the automorphism of geometric
# algebra determined by negating 1-vectors.
# Unfortunately, the package does not support
# grade-involution. Fortunately, it is easy to
# implement yourself.
def gradeInvolution(v):
return v.even() - v.odd()
a0, a1, a2, a3, a12, a13, a23, a123 = symbols('a0 a1 a2 a3 a12 a13 a23 a123')
v = (a0 +
a1 * e1 + a2 * e2 + a3 * e3 +
a12 * e1 * e2 + a13 * e1 * e3 + a23 * e2 * e3 +
a123 * e1 * e2 * e3)
print 'Grade-involution of'
print v.Fmt(3)
print '='
print gradeInvolution(v).Fmt(3)
print
### Reversion
# Definition:
# Reversion is the anti-automorphism of geometric
# algebra determined by the identity function on
# 1-vectors.
print 'Reversion of'
print v.Fmt(3)
print '='
print v.rev().Fmt(3)
print
### Conjugation
# Definition:
# Conjugation is the anti-automorphism of geometric
# algebra determined by negating 1-vectors.
# The package does not support conjugation either.
# Fortunately it is easy to implement too.
def conjugation(v):
return gradeInvolution(v).rev()
print 'Conjugation of'
print v.Fmt(3)
print '='
print conjugation(v).Fmt(3)
print
### Exterior product of vectors
# Define some 1-vectors.
a = 5 * e1 - 3 * e2
b = 2 * e1 + 9 * e2
# Compute a 2-blade as the exterior product `a ^ b`:
p = a ^ b
# The `p` now represents a weighted, oriented
# 2-dimensional plane spanned by `a` and `b`.
# Theorem:
# a_1 ^ ... ^ a_n = 0 <=> {a_1, ..., a_n} is linear independent.
# We can use this fact to test whether an arbitrary
# vector is on the plane `p` or not.
# Test whether various vectors are on the plane `p`.
print
print 'Plane containment:'
print 'a ^ p =', a ^ p
print 'b ^ p =', b ^ p
print 'e3 ^ p =', e3 ^ p
print
### Reverse of a k-versor
# Definition:
# A k-versor is a geometric product of 1-vectors.
# Define some 1-vectors.
a = 5 * e1 - 3 * e2
b = 2 * e1 + 9 * e2
# Compute a geometric product of vectors to get a 2-versor.
v = a * b
# Compute the reverse of `v`.
print 'Reverse of a 2-versor', v, '=', v.rev()
# Compute the reverse of `v` explicitly.
d = b * a
### Reverse of a k-blade
# Theorem:
# A k-blade is a k-versor.
# Define some 1-vectors.
a = 5 * e1 - 3 * e2
b = 2 * e1 + 9 * e2
# Compute a geometric product of vectors to get a 2-versor.
v = a ^ b
# Compute the reverse of `v`.
print 'Reverse of a 2-blade', v, '=', v.rev()
print
# Compute the reverse of `v` explicitly.
d = b ^ a
### Inverse of a k-versor
# Define some 1-vectors.
a = 5 * e1 - 3 * e2
b = 2 * e1 + 9 * e2
# Compute a 2-versor.
v = a * b
# Theorem:
# A k-versor v is invertible under the geometric product
# if and only if v* v != 0, where * is the reversion.
# Theorem:
# The inverse of an invertible k-versor `v` is given by v* / (v* v),
# where * is the reverse of a multi-vector.
# Unfortunately, the inverse in the package
# is buggy. For example,
# print inv(v)
# crashes.
# Fortunately, the versor inverse is easy to
# implement yourself.
def versorInverse(v):
return v.rev() / (v.rev() * v).scalar()
# Compute the inverse of `v` by the theorem.
print 'k-versor inverse:'
print 'v * versorInverse(v) =', v * versorInverse(v)
print 'versorInverse(v) * v', '=', versorInverse(v) * v
print
### Inverse of a k-blade
# Theorem:
# A k-blade is a k-versor.
# Define some 1-vectors.
a = 5 * e1 - 3 * e2
b = 2 * e1 + 9 * e2
# Compute a 2-blade.
v = a ^ b
# Theorem:
# A k-blade is a k-versor.
# Compute the inverse of `v` by the versor inverse.
print 'k-blade inverse:'
print 'v * versorInverse(v) =', v * versorInverse(v)
print 'versorInverse(v) * v', '=', versorInverse(v) * v
print
### Versor rotation
# Define a rotation 2-blade.
B = 4 * e1 ^ e2
# Define a rotation angle.
alpha = pi / 2
# Construct a versor representing
# rotation around 'p' with angle `alpha`.
q = ((alpha / 2) * (B / B.norm())).exp()
# This is the same as
# q = cos(alpha / 2) + sin(alpha / 2) * (p / p.norm())
# Define a vector to rotate.
v = 3 * e1 + 6 * e2 - 7 * e3
# Rotate the vector.
rv = q * v * inv(q)
# Output the rotated vector.
print 'Rotate', v, 'around', B, '=', rv
# Define a 2-blade to rotate.
A = (2 * e1) ^ (-3 * e2)
# Rotate the 2-blade.
# Note: can't do it with quaternions.
rA = q * A * inv(q)
# Output the rotated 2-blade.
print 'Rotate', A, 'around', B, '=', rA
# You can also rotate a 0-blade (a scalar),
# or a 3-blade (a pseudo-scalar). Both
# remain unaffected by the rotation;
# volume is preserved. Of course, you can
# also rotate multi-vectors by linearity.
### Versor reflection
# Define a vector along which to reflect.
q = e3
# Define a vector to reflect.
v = 3 * e1 + 6 * e2 - 7 * e3
# Reflect along e3.
rv = q * gradeInvolution(v) * inv(q)
# Note that this is a reflection along `q`.
# We can also define reflection through `q`;
# just remove the grade-involution
# Output the reflected vector.
print 'Reflect', v, 'along', q, '=', rv
# Define a 2-blade to reflect.
A = (2 * e1) ^ (-3 * e2 + 5 * e3)
# Reflect A along e3.
rA = q * gradeInvolution(A) * inv(q)
# Output the reflected 2-blade.
print 'Reflect', A, 'along', q, '=', rA
### Orthogonal complement
# Compute a 2-blade `p`, representing a plane.
a = 5 * e1 - 3 * e2
b = 2 * e1 + 9 * e2
p = a^b
# Definition:
# A pseudoscalar is a scalar multiple of a non-zero
# n-vector, where `n` is the dimension of the underlying
# vector space.
# Define a pseudoscalar.
I = 3 * e1 ^ e2 ^ e3
# Theorem:
# An orthogonal complement (a dual) of a k-blade can be
# computed by a contraction with the inverse of a pseudo-scalar.
# Compute the orthogonal complement of the plane `p`;
# Since we are in Cl(3, 0), this is the normal vector of
# the plane.
print 'Complement of', p, 'in', I, 'is', (p < inv(I))
### Contraction and orthogonal projection
# Define a plane `p`.
a = 5 * e1 - 3 * e2
b = 2 * e1 + 9 * e2
p = a^b
# Define a vector outside the plane.
c = -4 * e1 + 7 * e2 - e3
# Compute the contraction of `c` onto `p`.
# This is the orthogonal complement of the
# orthogonal projection of `c` onto `p`,
# on `p`.
d = c < p
print 'Contract ', c, 'onto', p, '=', d
# Compute the orthogonal complement of `d` on `p`.
e = d < inv(p)
print 'Project ', c, 'onto', p, '=', e
# The `e` is the orthogonal projection of `c` onto `p`!
### Contraction properties
a1, a2, a3 = symbols('a1 a2 a3')
b0, b1, b2, b3, b12, b13, b23, b123 = symbols('b0 b1 b2 b3 b12 b13 b23 b123')
c0, c1, c2, c3, c12, c13, c23, c123 = symbols('c0 c1 c2 c3 c12 c13 c23 c123')
a = a1 * e1 + a2 * e2 + a3 * e3
B = (b0 +
b1 * e1 + b2 * e2 + b2 * e3 +
b12 * e1 * e2 + b13 * e1 * e3 + b23 * e2 * e3 +
b123 * e1 * e2 * e3)
C = (c0 +
c1 * e1 + c2 * e2 + c2 * e3 +
c12 * e1 * e2 + c13 * e1 * e3 + c23 * e2 * e3 +
c123 * e1 * e2 * e3)
L = a < (B ^ C)
R = ((a < B) ^ C) + (gradeInvolution(B) ^ (a < C))
print (L - R)