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

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,

blades

{'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,

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.

locals().update(blades)

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

e3, e123

((1^e3), (1^e123))

Basics

Products

The basic products are available

e1*e2 # geometric product

(1.0^e12)

e1|e2 # inner product

0

e1^e2 # outer product

(1.0^e12)

e1^e2^e3 # even more outer products

(1.0^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

e1^e2+e2^e3 # fail

(2.0^e123)

(e1^e2) + (e2^e3) # correct

(1.0^e12) + (1.0^e23)

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

4|e1

0

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

4*e1

(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

from scipy import cos, sin

theta = pi/4
R = cos(theta) - sin(theta)*e23
R

0.70711 - (0.70711^e23)

You can also mix grades without any reason

A = 1 + 2*e1 + 3*e12 + 4*e123
A

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

~A

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

A(0) # get grade-0 elements of R

1.0

A(1) # get grade-1 elements of R

(2.0^e1)

A(2)  #  you get it

(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\]

(A*~A)(0)

30.0

This is done in the abs() operator

abs(A)**2

30.0

Inverse

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

A.inv()*A

1.0

A.inv()

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 pseudoscalar for the GA. In \(G_3\), the dual of a vector is a bivector,

a = 1*e1 + 2*e2 + 3*e3
a.dual()

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

cf.ugly()
A.inv()

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.1829268292682927, -7.80625564189563e-18, 0.09756097560975611, -0.29268292682926833])

You can also change the displayed precision

cf.pretty(precision=2)

A.inv()

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

from IPython.display import Image
Image(url='_static/reflection_on_vector.svg')

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

\[c \rightarrow ncn\]

c = e1+e2+e3    # a vector
n = e1          # the reflector
n*c*n          # reflect `a` in hyperplane normal to `n`

(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}\]

a = e1+e2+e3    # the vector
n = 3*e1          # the reflector
n*a*n.inv()

(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

\[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

from numpy import pi

R = e**(-pi/4*e12) # enacts rotation by pi/2
R

0.71 - (0.71^e12)

R*e1*~R    # rotate e1 by pi/2 in the e12-plane

(0.0^e1) + (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}}\]

R12 = lambda theta: e**(-theta/2*e12)
R12(pi/2)

0.71 - (0.71^e12)

And use it like this

a = e1+e2+e3
R = R12(pi/2)
R*a*~R

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

\[R_B = e^{-\frac{B}{2}}\]

R_B = lambda B: e**(-B/2.)

Then you could do

R12 = R_B(pi/4*e12)
R23 = R_B(pi/5*e23)

or

R_B(pi/6*(e23+e12))  # rotor enacting a pi/6-rotation in the e23+e12-plane

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 sandwich product, which is nice if you are cascading a bunch of rotors, like so

\[\underline{R_C}( \underline{R_B}( \underline{R_A}(a)))\]

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

(0.52^e1) + (0.74^e2) + (1.48^e3)

Then you can do things like

R12 = R_factory(pi/3*e12)
R23 = R_factory(pi/3*e23)
R13 = R_factory(pi/3*e13)

R12(R23(R13(a)))

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

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)

(0.0^e1) + (0.0^e2) + (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,

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

blades

{'x': (1^x), 'y': (1^y), 'i': (1^i)}

locals().update(blades)

1*x+2*y

(1^x) + (2^y)

(1+4*i)

1.0 + (4.0^i)