Algebra of subspaces

Back to Geometric algebra in sympy

In this section we study the geometric algebra .

In the other part of the tutorial, we looked at , the even sub-algebra of , and saw that it encodes the quaternions. By using the full algebra , we add the ability to work directly with sub-spaces, rather than just with vectors.

A k-blade, an exterior product of -vectors, is interpreted geometrically as the subspace its -vector exterior-factors span. This subspace can be transformed by -versors, a geometric product of -vectors. Each -vector factor in a -versor corresponds to an orthogonal reflection along that vector. By the Cartan-Dieudonne theorem, the orthogonal transformations are exactly those which can be obtained by composing reflections. Therefore, -versors encode orthogonal transformations as first-class citizens into the algebra. Moreover, -versors are exactly the quaternions.

An analogous discussion could be written about , , and the complex numbers.

# 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)

Files

Tutorial code for Cl(3, 0)