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

In [1]:
from numpy import e,pi
import clifford as cf

layout, blades = cf.Cl(3) # creates a 3-dimensional clifford algebra

Given a three dimensional GA with the orthonormal basis,

\[e_{i}\cdot e_{j}=\delta_{ij}\]

The basis consists of scalars, three vectors, three bivectors, and a trivector.

\[\{\underbrace{\alpha}_{\mbox{scalar}},\qquad\underbrace{e_{1},e_{2},e_{3}}_{\mbox{vectors}},\qquad\underbrace{e_{12},e_{23},e_{13}}_{\mbox{bivectors}},\qquad\underbrace{e_{123}}_{\mbox{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,

In [2]:
blades
Out[2]:
{'e1': (1^e1),
 'e12': (1^e12),
 'e123': (1^e123),
 'e13': (1^e13),
 'e2': (1^e2),
 'e23': (1^e23),
 'e3': (1^e3)}

You may wish to explicitly assign the blades to variables like so,

In [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.

In [4]:
locals().update(blades)

Now, all the blades have been defined in the local namespace

In [5]:
e3, e123
Out[5]:
((1^e3), (1^e123))

Basics

Products

The basic products are available

In [6]:
e1*e2 # geometric product
Out[6]:
(1.0^e12)
In [7]:
e1|e2 # inner product
Out[7]:
0
In [8]:
e1^e2 # outer product
Out[8]:
(1.0^e12)
In [9]:
e1^e2^e3 # even more outer products
Out[9]:
(1.0^e123)

Defects in Precidence

Python’s operator precidence makes the outer product evaluate after addition. This requires the use of parenthesis when using outer products. For example

In [10]:
e1^e2+e2^e3 # fail
Out[10]:
(2.0^e123)
In [11]:
(e1^e2) + (e2^e3) # correct
Out[11]:
(1.0^e12) + (1.0^e23)

Also the inner product of a scalar and a Multivector is 0,

In [12]:
4|e1
Out[12]:
0

So for scalars, use the outer product or geometric product instead

In [13]:
4*e1
Out[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

In [14]:
from scipy import cos, sin

theta = pi/4
R = cos(theta) - sin(theta)*e23
R
Out[14]:
0.70711 - (0.70711^e23)

You can also mix grades without any reason

In [15]:
A = 1 + 2*e1 + 3*e12 + 4*e123
A
Out[15]:
1.0 + (2.0^e1) + (3.0^e12) + (4.0^e123)

Reversion

The reversion operator is accomplished with the tilde ~ in front of the Multivector on which it acts

In [16]:
~A
Out[16]:
1.0 + (2.0^e1) - (3.0^e12) - (4.0^e123)

Grade Projection

Taking a projection onto a specific grade \(n\) of a Multivector is usually written

\[\langle A \rangle _n\]

can be done by using soft brackets, like so

In [17]:
A(0) # get grade-0 elements of R
Out[17]:
1.0
In [18]:
A(1) # get grade-1 elements of R
Out[18]:
(2.0^e1)
In [19]:
A(2)  #  you get it
Out[19]:
(3.0^e12)

Magnitude

Using the reversion and grade projection operators, we can define the magnitude of \(A\)

\[|A|^2 = \langle A\tilde{A}\rangle\]
In [20]:
(A*~A)(0)
Out[20]:
30.0

This is done in the abs() operator

In [21]:
abs(A)**2
Out[21]:
30.0

Inverse

The inverse of a Multivector is defined as \(A^{-1}A=1\)

In [22]:
A.inv()*A
Out[22]:
1.0
In [23]:
A.inv()
Out[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

\[AI^{-1}\]

Where, \(I\) is the psuedoscalar for the GA. In \(G_3\), the dual of a vector is a bivector,

In [24]:
a = 1*e1 + 2*e2 + 3*e3
a.dual()
Out[24]:
-(3.0^e12) + (2.0^e13) - (1.0^e23)

Pretty, Ugly, and Display Precision

You can toggle pretty printing with with pretty() or ugly(). ugly returns an eval-able string.

In [25]:
cf.ugly()
A.inv()
Out[25]:
MultiVector(Layout([1, 1, 1], [(), (1,), (2,), (3,), (1, 2), (1, 3), (2, 3), (1, 2, 3)], firstIdx=1, names=['', 'e1', 'e2', 'e3', 'e12', 'e13', 'e23', 'e123']), value=[0.13414634146341464, 0.12195121951219513, -0.0, -0.14634146341463417, 0.18292682926829273, -2.0816681711721682e-17, 0.097560975609756101, -0.29268292682926833])

You can also change the displayed precision

In [26]:
cf.pretty(precision=2)

A.inv()
Out[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

In [27]:
from IPython.display import Image
Image(url='_static/reflection_on_vector.svg')
Out[27]:

Reflecting a vector \(c\) about a normalized vector \(n\) is pretty simple,

\[c \rightarrow ncn\]
In [28]:
c = e1+e2+e3    # a vector
n = e1          # the reflector
n*c*n          # reflect `a` in hyperplane normal to `n`
Out[28]:
(1.0^e1) - (1.0^e2) - (1.0^e3)

Because we have the inv() available, we can equally well reflect in un-normalized vectors using,

\[a \rightarrow nan^{-1}\]
In [29]:
a = e1+e2+e3    # the vector
n = 3*e1          # the reflector
n*a*n.inv()
Out[29]:
(1.0^e1) - (1.0^e2) - (1.0^e3)

Refelections can also be made with repsect to the a ‘hyperplane normal to the vector \(n\)’, in this case the formula is negated

\[c \rightarrow -ncn^{-1}\]

Rotations

A vector can be rotated using the formula

\[a \rightarrow Ra\tilde{R}\]

Where \(R\) is a rotor. A rotor can be defined by multiple reflections,

\[R=mn\]

or by a plane and an angle,

\[R = e^{-\frac{\theta}{2}\hat{B}}\]

For example

In [30]:
from numpy import pi

R = e**(-pi/4*e12) # enacts rotation by pi/2
R
Out[30]:
0.71 - (0.71^e12)
In [31]:
R*e1*~R    # rotate e1 by pi/2 in the e12-plane
Out[31]:
(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,

\[R_{12} = e^{-\frac{\theta}{2}e_{12}}\]
In [32]:
R12 = lambda theta: e**(-theta/2*e12)
R12(pi/2)
Out[32]:
0.71 - (0.71^e12)

And use it like this

In [33]:
a = e1+e2+e3
R = R12(pi/2)
R*a*~R

Out[33]:
-(1.0^e1) + (1.0^e2) + (1.0^e3)

You might as well make the angle arugment a bivector, so that you can control the plane of rotation as well as the angle

\[R_B = e^{-\frac{B}{2}}\]
In [34]:
R_B = lambda B: e**(-B/2.)

Then you could do

In [35]:
R12 = R_B(pi/4*e12)
R23 = R_B(pi/5*e23)

or

In [36]:
R_B(pi/6*(e23+e12))  # rotor enacting a pi/6-rotation in the e23+e12-plane
Out[36]:
0.93 - (0.26^e12) - (0.26^e23)

Maybe you want to define a function which returns a function that enacts a specified rotation,

\[f(B) \rightarrow \underline{R_B}(a) = R_Ba\tilde{R_B}\]

This just saves you having to write out the sandwhich product, which is nice if you are cascading a bunch of rotors, like so

\[\underline{R_C}( \underline{R_B}( \underline{R_A}(a)))\]
In [37]:
def R_factory( B):
    def dummy_f(a):
        R = e**(-B/2)
        return R*a*~R
    return dummy_f

R = R_factory(pi/6*(e23+e12)) # this returns a function
R(a) # which acts on a vector
Out[37]:
(0.52^e1) + (0.74^e2) + (1.48^e3)

Then you can do things like

In [38]:
R12 = R_factory(pi/3*e12)
R23 = R_factory(pi/3*e23)
R13 = R_factory(pi/3*e13)

R12(R23(R13(a)))
Out[38]:
(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\).

\[f(A,B,C,x) = \underline{R_A} (\underline{R_B} (\underline{R_C}(x)))\]
In [39]:
from functools import reduce

# a sequence of rotations
def R_seq(*args):
    Bs,a = args[:-1],args[-1]
    R_lst =  [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 lambda a: R*a*~R



# rotation sequence by  pi/2-in-e12 THEN pi/2-in-e23
R = R_seq(pi/2*e23, pi/2*e12, e1)

R(e1)
Out[39]:
(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,

In [40]:
layout,blades = cf.Cl(2, names = ['','x','y','i'])

blades
Out[40]:
{'i': (1^i), 'x': (1^x), 'y': (1^y)}
In [41]:
locals().update(blades)
In [42]:
1*x+2*y
Out[42]:
(1^x) + (2^y)
In [43]:
(1+4*i)
Out[43]:
1.0 + (4.0^i)
In [44]: