This page was generated from
docs/tutorials/g3-algebra-of-space.ipynb.
Interactive online version:
.
The Algebra Of Space (G3)¶
In this notebook, we give a more detailed look at how to use clifford
, using the algebra of three dimensional space as a context.
Setup¶
First, we import clifford as cf
, and instantiate a three dimensional geometric algebra using Cl()
(docs).
[1]:
import clifford as cf
layout, blades = cf.Cl(3) # creates a 3-dimensional clifford algebra
Given a three dimensional GA with the orthonormal basis,
The basis consists of scalars, three vectors, three bivectors, and a trivector.
Cl()
creates the algebra and returns a layout
and blades
. The layout
holds information and functions related this instance of G3
, and the blades
is a dictionary which contains the basis blades, indexed by their string representations,
[2]:
blades
[2]:
{'': 1,
'e1': (1^e1),
'e2': (1^e2),
'e3': (1^e3),
'e12': (1^e12),
'e13': (1^e13),
'e23': (1^e23),
'e123': (1^e123)}
You may wish to explicitly assign the blades to variables like so,
[3]:
e1 = blades['e1']
e2 = blades['e2']
# etc ...
Or, if you’re lazy and just working in an interactive session you can use locals()
to update your namespace with all of the blades at once.
[4]:
locals().update(blades)
Now, all the blades have been defined in the local namespace
[5]:
e3, e123
[5]:
((1^e3), (1^e123))
Basics¶
Products¶
The basic products are available
[6]:
e1*e2 # geometric product
[6]:
(1^e12)
[7]:
e1|e2 # inner product
[7]:
0
[8]:
e1^e2 # outer product
[8]:
(1^e12)
[9]:
e1^e2^e3 # even more outer products
[9]:
(1^e123)
Defects in Precedence¶
Python’s operator precedence makes the outer product evaluate after addition. This requires the use of parentheses when using outer products. For example
[10]:
e1^e2 + e2^e3 # fail, evaluates as
[10]:
(2^e123)
[11]:
(e1^e2) + (e2^e3) # correct
[11]:
(1^e12) + (1^e23)
Also the inner product of a scalar and a Multivector is 0,
[12]:
4|e1
[12]:
0
So for scalars, use the outer product or geometric product instead
[13]:
4*e1
[13]:
(4^e1)
Multivectors¶
Multivectors can be defined in terms of the basis blades. For example you can construct a rotor as a sum of a scalar and bivector, like so
[14]:
import math
theta = math.pi/4
R = math.cos(theta) - math.sin(theta)*e23
R
[14]:
0.70711 - (0.70711^e23)
You can also mix grades without any reason
[15]:
A = 1 + 2*e1 + 3*e12 + 4*e123
A
[15]:
1 + (2^e1) + (3^e12) + (4^e123)
Reversion¶
The reversion operator is accomplished with the tilde ~
in front of the Multivector on which it acts
[16]:
~A
[16]:
1 + (2^e1) - (3^e12) - (4^e123)
Grade Projection¶
Taking a projection onto a specific grade \(n\) of a Multivector is usually written
can be done by using soft brackets, like so
[17]:
A(0) # get grade-0 elements of R
[17]:
1
[18]:
A(1) # get grade-1 elements of R
[18]:
(2^e1)
[19]:
A(2) # you get it
[19]:
(3^e12)
Magnitude¶
Using the reversion and grade projection operators, we can define the magnitude of \(A\)
[20]:
(A*~A)(0)
[20]:
30
This is done in the abs()
operator
[21]:
abs(A)**2
[21]:
30.0
Inverse¶
The inverse of a Multivector is defined as \(A^{-1}A=1\)
[22]:
A.inv()*A
[22]:
1.0
[23]:
A.inv()
[23]:
0.13415 + (0.12195^e1) - (0.14634^e3) + (0.18293^e12) + (0.09756^e23) - (0.29268^e123)
Dual¶
The dual of a multivector \(A\) can be defined as
Where, \(I\) is the pseudoscalar for the GA. In \(G_3\), the dual of a vector is a bivector,
[24]:
a = 1*e1 + 2*e2 + 3*e3
a.dual()
[24]:
-(3^e12) + (2^e13) - (1^e23)
Pretty, Ugly, and Display Precision¶
You can toggle pretty printing with with pretty()
or ugly()
. ugly
returns an eval-able string.
[25]:
cf.ugly()
A.inv()
[25]:
MultiVector(Layout([1, 1, 1],
ids=BasisVectorIds.ordered_integers(3),
order=BasisBladeOrder.shortlex(3),
names=['', 'e1', 'e2', 'e3', 'e12', 'e13', 'e23', 'e123']),
[0.13414634146341464, 0.12195121951219512, 0.0, -0.14634146341463414, 0.18292682926829268, 0.0, 0.0975609756097561, -0.2926829268292683])
You can also change the displayed precision
[26]:
cf.pretty(precision=2)
A.inv()
[26]:
0.13 + (0.12^e1) - (0.15^e3) + (0.18^e12) + (0.1^e23) - (0.29^e123)
This does not effect the internal precision used for computations.
Applications¶
Reflections¶
Reflecting a vector \(c\) about a normalized vector \(n\) is pretty simple,
[27]:
c = e1+e2+e3 # a vector
n = e1 # the reflector
n*c*n # reflect `a` in hyperplane normal to `n`
[27]:
(1^e1) - (1^e2) - (1^e3)
Because we have the inv()
available, we can equally well reflect in un-normalized vectors using,
[28]:
a = e1+e2+e3 # the vector
n = 3*e1 # the reflector
n*a*n.inv()
[28]:
(1.0^e1) - (1.0^e2) - (1.0^e3)
Reflections can also be made with respect to the a ‘hyperplane normal to the vector \(n\)’, in this case the formula is negated
Rotations¶
A vector can be rotated using the formula
Where \(R\) is a rotor. A rotor can be defined by multiple reflections,
or by a plane and an angle,
For example
[29]:
R = math.e**(-math.pi/4*e12) # enacts rotation by pi/2
R
[29]:
0.71 - (0.71^e12)
[30]:
R*e1*~R # rotate e1 by pi/2 in the e12-plane
[30]:
(1.0^e2)
Some Ways to use Functions¶
Maybe we want to define a function which can return rotor of some angle \(\theta\) in the \(e_{12}\)-plane,
[31]:
R12 = lambda theta: math.e**(-theta/2*e12)
R12(math.pi/2)
[31]:
0.71 - (0.71^e12)
And use it like this
[32]:
a = e1+e2+e3
R = R12(math.pi/2)
R*a*~R
[32]:
-(1.0^e1) + (1.0^e2) + (1.0^e3)
You might as well make the angle argument a bivector, so that you can control the plane of rotation as well as the angle
[33]:
R_B = lambda B: math.e**(-B/2)
Then you could do
[34]:
R12 = R_B(math.pi/4*e12)
R23 = R_B(math.pi/5*e23)
or
[35]:
R_B(math.pi/6*(e23+e12)) # rotor enacting a pi/6-rotation in the e23+e12-plane
[35]:
0.93 - (0.26^e12) - (0.26^e23)
Maybe you want to define a function which returns a function that enacts a specified rotation,
This just saves you having to write out the sandwich product, which is nice if you are cascading a bunch of rotors, like so
[36]:
def R_factory(B):
def apply_rotation(a):
R = math.e**(-B/2)
return R*a*~R
return apply_rotation
R = R_factory(math.pi/6*(e23+e12)) # this returns a function
R(a) # which acts on a vector
[36]:
(0.52^e1) + (0.74^e2) + (1.48^e3)
Then you can do things like
[37]:
R12 = R_factory(math.pi/3*e12)
R23 = R_factory(math.pi/3*e23)
R13 = R_factory(math.pi/3*e13)
R12(R23(R13(a)))
[37]:
(0.41^e1) - (0.66^e2) + (1.55^e3)
To make cascading a sequence of rotations as concise as possible, we could define a function which takes a list of bivectors \(A,B,C,..\) , and enacts the sequence of rotations which they represent on a some vector \(x\).
[38]:
from functools import reduce
# a sequence of rotations
def R_seq(*args):
*Bs, x = args
R_lst = [math.e**(-B/2) for B in Bs] # create list of Rotors from list of Bivectors
R = reduce(cf.gp, R_lst) # apply the geometric product to list of Rotors
return R*x*~R
# rotation sequence by pi/2-in-e12 THEN pi/2-in-e23
R_seq(math.pi/2*e23, math.pi/2*e12, e1)
[38]:
(1.0^e3)
Changing Basis Names¶
If you want to use different names for your basis as opposed to e’s with numbers, supply the Cl()
with a list of names
. For example for a two dimensional GA,
[39]:
layout, blades = cf.Cl(2, names=['','x','y','i'])
blades
[39]:
{'': 1, 'x': (1^x), 'y': (1^y), 'i': (1^i)}
[40]:
locals().update(blades)
[41]:
1*x + 2*y
[41]:
(1^x) + (2^y)
[42]:
1 + 4*i
[42]:
1 + (4^i)