This page was generated from docs/tutorials/MatrixRepresentationsOfGeometricFunctions.ipynb. Interactive online version: .
Matrix Representations of Geometric Functions¶
More info can be found in Chapter 5 of “New foundations for classical mechanics”, by David Hestenes.
Creating a geometric function
Apply it to an orthonormal frame
Convert the resultant frame into a matrix
The matrix is defined as the inner product of each basis element of original and transformed frame.
(or vice-versa with the i,j, you get the point). Since we are going to do this repeatedly, define a func2Mat()
[1]:
import numpy as np
def func2Mat(f,I):
'''
Convert a function acting on a vector into a matrix, given
the space defined by psuedoscalar I
'''
A = I.basis()
B = [f(a) for a in A]
M = [float(b | a) for a in A for b in B]
return np.array(M).reshape(len(B), len(B))
Start with initializing a euclidean N-dimensional algebra and assign our pseudoscalar to \(I\), pretty standard.
[2]:
from clifford import Cl
from math import *
l,b = Cl(3) # returns (layout,blades). you can change dimesion here
I = l.pseudoScalar
Anti-symmetric¶
This is so easy,
[3]:
B = l.randomIntMV()(2) # we use randIntMV because its easier to read
f = lambda x:x|B
func2Mat(f,I=I)
[3]:
array([[ 0., -3., 9.],
[ 3., 0., -6.],
[-9., 6., 0.]])
Whats the B? you can read its values straight off the matrix.
[4]:
B
[4]:
(3^e12) - (9^e13) + (6^e23)
Diagonal ( Directional Scaling)¶
A bit awkward this one, but its made by projection onto each basis vector, then scaling the component by some amount.
[5]:
ls = range(1,len(I.basis())+1) # some dilation values (eigenvalues)
A = I.basis()
d = lambda x: sum([(x|a)/a*l for a,l in zip(A,ls)])
func2Mat(d,I=I)
[5]:
array([[1., 0., 0.],
[0., 2., 0.],
[0., 0., 3.]])
Orthgonal, Rotation¶
where
[6]:
B = l.randomMV()(2)
R = e**(B/2)
r = lambda x: R*x*~R
func2Mat(r,I=I)
[6]:
array([[-0.02877844, 0.83278728, 0.5528446 ],
[-0.2119728 , -0.54557927, 0.81080873],
[ 0.97685175, -0.09385421, 0.19222917]])
The inverse of this is ,
[7]:
rinv = lambda x: ~R*x*R # the inverse rotation
func2Mat(rinv,I=I)
[7]:
array([[-0.02877844, -0.2119728 , 0.97685175],
[ 0.83278728, -0.54557927, -0.09385421],
[ 0.5528446 , 0.81080873, 0.19222917]])
Orthogonal, Reflection¶
[8]:
a = l.randomIntMV()(1)
n = lambda x: -a*x/a
func2Mat(n,I=I)
[8]:
array([[ 0.97530864, -0.09876543, -0.19753086],
[-0.09876543, 0.60493827, -0.79012346],
[-0.19753086, -0.79012346, -0.58024691]])
[9]:
a
[9]:
-(1^e1) - (4^e2) - (8^e3)
Notice the determinant for reflection is -1, and for rotation is +1.
[10]:
from numpy.linalg import det
det(func2Mat(n,I=I)), det(func2Mat(r,I=I))
[10]:
(-0.9999999999999999, 1.0000000000411782)
Symmetric¶
This can be built up from the functions we just defined ie Rotation*Dilation/Rotation
which if you write it out, looks kind of dumb
So, the antisymmetric matrix is interpreted as a set dilations about a some orthogonal frame rotated from the basis (what basis,eh? exactly what basis!).
More generally we could include reflection in the \(R\) too.
[11]:
g = lambda x: r(d(rinv(x)))
func2Mat(g,I=I)
[11]:
array([[2.30480895, 0.44215098, 0.13438513],
[0.44215098, 2.61247833, 0.36292709],
[0.13438513, 0.36292709, 1.08271272]])
Eigen stuffs¶
By definition the eigen-stuff is the invariants of the transformation, sometimes this is a vector, and other times it is a plane.
Rotation¶
The eigen blades of a rotation are really the axis and plane of rotation.
[12]:
from numpy.linalg import eig
vals, vecs = eig(func2Mat(r,I=I))
np.round(vecs,3)
[12]:
array([[ 0.626+0.j , -0.136-0.535j, -0.136+0.535j],
[ 0.293+0.j , 0.676+0.j , 0.676-0.j ],
[ 0.723+0.j , -0.157+0.463j, -0.157-0.463j]])
If you checkout the real column, and compare this to the bivector which generated this rotation (aka the generator), after its been normalized
[13]:
B/(abs(B))
[13]:
(0.72272^e12) - (0.29331^e13) + (0.62581^e23)
[14]:
B
[14]:
(1.68666^e12) - (0.68452^e13) + (1.46049^e23)
[15]:
vals
[15]:
array([ 1. +0.j , -0.69106427+0.72279332j,
-0.69106427-0.72279332j])
[16]:
cos(abs(B)), sin(abs(B))
[16]:
(-0.6910642689158123, 0.7227933150132574)
Symmetric¶
For the symmetric matrix, the invariant thing the orthonormal frame along which the dilations take place
[17]:
vals, vecs = eig(func2Mat(g,I=I))
np.round(vecs,5).T
[17]:
array([[ 0.55284, 0.81081, 0.19223],
[ 0.83279, -0.54558, -0.09385],
[-0.02878, -0.21197, 0.97685]])
This is easily found by using the rotation part of the symmetric operator,
[18]:
[R*a*~R for a in I.basis()]
[18]:
[-(0.02878^e1) - (0.21197^e2) + (0.97685^e3),
(0.83279^e1) - (0.54558^e2) - (0.09385^e3),
(0.55284^e1) + (0.81081^e2) + (0.19223^e3)]
Primitive Visualization in 2D¶
[19]:
from pylab import linspace, plot,axis,legend
def plot_ps(ps,**kw):
x = [p[e1] for p in ps]
y = [p[e2] for p in ps]
plot(x,y, marker='o',ls='',**kw)
l,b = Cl(2)
locals().update(b)
I = l.pseudoScalar
## define function of interest
B = l.randomMV()(2)
R = e**(B/2)
f = lambda x: R*x*~R
## loop though cartesian grid and apply f,
ps,qs=[],[]
for x in linspace(-1,1,11):
for y in linspace(-1,1,11):
p = x*e1 + y*e2
q = f(p)
ps.append(p)
qs.append(q)
plot_ps(ps,label='before')
plot_ps(qs,label='after')
axis('equal')
legend()
/home/docs/checkouts/readthedocs.org/user_builds/clifford-alex/envs/latest/lib/python3.8/site-packages/traitlets/traitlets.py:3030: FutureWarning: --rc={'figure.dpi': 96, 'savefig.transparent': True} for dict-traits is deprecated in traitlets 5.0. You can pass --rc <key=value> ... multiple times to add items to a dict.
warn(
[19]:
<matplotlib.legend.Legend at 0x7f6a014c8a30>
[20]:
func2Mat(f,I =I )
[20]:
array([[ 0.82037719, 0.57182276],
[-0.57182276, 0.82037719]])