Warning
This page was created from a pull request (#389).
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., 1., 7.],
[ -1., 0., 10.],
[ -7., -10., 0.]])
Whats the B? you can read its values straight off the matrix.
[4]:
B
[4]:
-(1^e12) - (7^e13) - (10^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.66607628, 0.73131933, -0.14667795],
[-0.38694602, 0.50691486, 0.77026625],
[ 0.63766383, -0.45629963, 0.62062507]])
The inverse of this is ,
[7]:
rinv = lambda x: ~R*x*R # the inverse rotation
func2Mat(rinv,I=I)
[7]:
array([[ 0.66607628, -0.38694602, 0.63766383],
[ 0.73131933, 0.50691486, -0.45629963],
[-0.14667795, 0.77026625, 0.62062507]])
Orthogonal, Reflection¶
[8]:
a = l.randomIntMV()(1)
n = lambda x: -a*x/a
func2Mat(n,I=I)
[8]:
array([[ 0.07834101, -0.82949309, -0.55299539],
[-0.82949309, 0.25345622, -0.49769585],
[-0.55299539, -0.49769585, 0.66820276]])
[9]:
a
[9]:
-(10^e1) - (9^e2) - (6^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]:
(-1.0, 1.00000000000003)
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([[ 1.57785681, 0.14475449, -0.51576477],
[ 0.14475449, 2.44358288, 0.72478803],
[-0.51576477, 0.72478803, 1.97856032]])
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.223-0.476j, -0.223+0.476j, 0.668+0.j ],
[ 0.639+0.j , 0.639-0.j , 0.427+0.j ],
[-0.204+0.523j, -0.204-0.523j, 0.609+0.j ]])
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.60914^e12) - (0.42725^e13) + (0.66814^e23)
[14]:
B
[14]:
(0.70829^e12) - (0.49679^e13) + (0.77688^e23)
[15]:
vals
[15]:
array([0.39680811+0.91790159j, 0.39680811-0.91790159j,
1. +0.j ])
[16]:
cos(abs(B)), sin(abs(B))
[16]:
(0.3968081091656232, 0.9179015875900874)
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.66608, 0.38695, -0.63766],
[ 0.73132, 0.50691, -0.4563 ],
[-0.14668, 0.77027, 0.62063]])
This is easily found by using the rotation part of the symmetric operator,
[18]:
[R*a*~R for a in I.basis()]
[18]:
[(0.66608^e1) - (0.38695^e2) + (0.63766^e3),
(0.73132^e1) + (0.50691^e2) - (0.4563^e3),
-(0.14668^e1) + (0.77027^e2) + (0.62063^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/envs/389/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 0x7ff4a02bc1c0>
[20]:
func2Mat(f,I =I )
[20]:
array([[ 0.67127791, -0.74120575],
[ 0.74120575, 0.67127791]])