Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download

This notebook is a tutorial on how to use Clifford Algebras and SageMath. In it, we derive the trigonometric parametrization, curvature, and volume for the sphere.

Project: SageMathExploration

Path: CliffordAlgebraTutorial.ipynb

Views: ⁶⁹
Image: ubuntu2204

Kernel: SageMath 10.2

Clifford Algebra Tutorial

This notebook is a tutorial on how to use Clifford Algebras and SageMath. To run it, you will need either a local sage executable, or to run it on cocalc. I’ve tested it with Sage 10.2. Clifford algebras represent geometric objects such as vectors in an algebraic manner which is amenable to direct computation. SageMath allows me to run these computations symbolically so that we can end up with generic formulas rather than just final numerical answers.

Mathematical Preliminaries

Clifford algebras start with an associative and distributive multiplication operation over vectors and the additional constraint that for all vectors v, $v^2 = v^TMv$ for some fixed symmetric matrix $M$ .

For example, if M is the 3 by 3 identity matrix, then for any basis vector $e_i$ we have $e_i^2 = e_i^TMe_i = 1$ . So then we can rewrite a string such as $e_1e_2e_2e_1\rightarrow e_1e_1 \rightarrow 1$ .

But we have enough information to derive other rewrite rules as well. Consider $M = D$ , an arbitrary diagonal matrix, so that if $i \neq j$ then $e_i^TDe_j = 0$ . Let’s calculate $(e_i + e_j)^2$ . We can simplify this expression two different ways:

(i) $(e_i + e_j)^2 = e_i^2 + e_ie_j + e_je_i + e_j^2$ .

(ii) $(e_i + e_j)^2 = (e_i + e_j)^TD(e_i + e_j) = e_i^TDe_i + e_i^TDe_j + e_j^TDe_ + e_j^TDe_j = e_i^2 + 2e_i^TDe_j + e_j^2 = e_i^2 + e_j^2$

Equating (i) and (ii) tells us that $e_ie_j = -e_ie_j$ . So now we can make calculations like $(e_ie_j)^2 = e_ie_je_ie_j = -e_ie_ie_je_j = -e_i^2e_j^2$ .

Notice that if $e_i^2 = e_j^2 = 1$ , then we have $(e_ie_j)^2=-1$ . In 2D, $e_1e_2$ turns out to be analogous to the complex unit $i$ , and in 3D $e_1e_2$ , $e_2e_3$ , $e_3e_1$ , and $1$ turn out to be analogous to $i$ , $j$ , $k$ , and $1$ of the quaternions. We call $e_ie_j$ bivectors and think of them as oriented areas in the plane spanned by $e_i$ and $e_j$ .

The grade of a Clifford algebra term is the number of basis elements in each reduced term — for instance $1$ is grade 0, $ae_i + be_j$ is grade 1, $e_ie_j$ is grade 2, and so on ( $1+e_i$ is called mixed grade). We can generalize the complex numbers and the quaternions by considering even-grade subalgebras of Clifford algebras with more basis elements. We keep associativity in these higher dimensional algebras by giving up totality of the inverse operation (I think).

Notebook Contents

In this notebook, I will demonstrate the use of a Clifford algebra by deriving:

Rational and irrational parametrizations for the surface of a sphere
Curvature of the sphere
Volume of the sphere

Clifford Algebra Sage Setup

In [1]:

qf = DiagonalQuadraticForm(SR,[1,1,1])
Cl.<e1,e2,e3> = CliffordAlgebra(qf)
Cl

The Clifford algebra of the Quadratic form in 3 variables over Symbolic Ring with coefficients: 
[ 1 0 0 ]
[ * 1 0 ]
[ * * 1 ]

This is some boilerplate to get SageMath to run simplic manipulation of a Clifford algebra term. Using this lift function, trigonometric simplification and differentiation work as expected.

In [2]:

r, theta, phi = var('r theta phi')
def build_ga_term(a,b,c,d,e,f,g,h):
    return a + b*e1 + c*e2 + d*e3 + e*e1*e2 + f*e1*e3 + g*e2*e3 + h*e1*e2*e3
def lift(fxn):
    def fxn_aux(ga_term, *args):
        a,b,c,d,e,f,g,h = [fxn(x, *args) for x in ga_term.dense_coefficient_list()]
        return build_ga_term(a,b,c,d,e,f,g,h)
    return fxn_aux
simp = lift(Expression.simplify_full)
ga_diff = lift(diff)

expr = 3*theta + (cos(theta)^2 + sin(theta)^2)*e1
show("(cos(theta)^2 + sin(theta)^2)*e1 + 3*theta -> ", simp(expr))
show("d/d(theta) ((cos(theta)^2 + sin(theta)^2)*e1 + 3*theta): ", ga_diff(expr,theta))

$\displaystyle \verb|(cos(theta)^2|\verb| |\verb|+|\verb| |\verb|sin(theta)^2)*e1|\verb| |\verb|+|\verb| |\verb|3*theta|\verb| |\verb|->| e_{1} + 3 \, \theta$

Since for any vector v we have $v^2 = v^TMv$ , we can invert a vector via $v^{-1} = \frac{v}{v^TMv}$ .

In [3]:

v1,v2,v3 = var('v1 v2 v3')
w1,w2,w3 = var('w1 w2 w3')
v = v1*e1 + v2*e2 + v3*e3
w = w1*e1 + w2*e2 + w3*e3
def to_scalar(t): return t.dense_coefficient_list()[0]
#def norm(t): 
#    expr = sqrt(to_scalar(inner(t, t)))
#    return expr.canonicalize_radical().full_simplify()
def norm(t): # just squaring and adding all coeffs
    return sqrt(sum(map(lambda x: x^2, t.coefficients()))).canonicalize_radical()
def inv(t): return t / to_scalar(t^2)
show("v: ", v); show("inverse(v): ", inv(v)); show("v * inverse(v): ", simp(v*inv(v)))

$\displaystyle \verb|v:| v_{1} e_{1} + v_{2} e_{2} + v_{3} e_{3}$

$\displaystyle \verb|inverse(v):| \left( \frac{v_{1}}{v_{1}^{2} + v_{2}^{2} + v_{3}^{2}} \right) e_{1} + \left( \frac{v_{2}}{v_{1}^{2} + v_{2}^{2} + v_{3}^{2}} \right) e_{2} + \left( \frac{v_{3}}{v_{1}^{2} + v_{2}^{2} + v_{3}^{2}} \right) e_{3}$

The operation $-ve_1v^{-1}$ reflects $e_1$ across the (hyper) plane which is orthogonal to $v$ . Since reflection preserves a vector’s size and $e_1$ is unit length, we get a (rational) parametrization of the unit sphere. Then we show that the sphere’s parametrization is simplified by translating the sphere to touch the origin. We also define a general bivector.

In [4]:

show("Rational param of sphere: ", simp(v*e1*inv(v)))
show("Rational param of core sphere: ", simp(e1 + v*e1*inv(v))) # rational param of core sphere]
x12,x23,x13 = var('x12,x23,x13')
B = x12*e1*e2 + x13*e1*e3 + x23*e2*e3
show("B: ", B)
show("B^2: ", B*B)
show("|B|: ", norm(B))

$\displaystyle \verb|Rational|\verb| |\verb|param|\verb| |\verb|of|\verb| |\verb|sphere:| \left( \frac{v_{1}^{2} - v_{2}^{2} - v_{3}^{2}}{v_{1}^{2} + v_{2}^{2} + v_{3}^{2}} \right) e_{1} + \left( \frac{2 \, v_{1} v_{2}}{v_{1}^{2} + v_{2}^{2} + v_{3}^{2}} \right) e_{2} + \left( \frac{2 \, v_{1} v_{3}}{v_{1}^{2} + v_{2}^{2} + v_{3}^{2}} \right) e_{3}$

$\displaystyle \verb|Rational|\verb| |\verb|param|\verb| |\verb|of|\verb| |\verb|core|\verb| |\verb|sphere:| \left( \frac{2 \, v_{1}^{2}}{v_{1}^{2} + v_{2}^{2} + v_{3}^{2}} \right) e_{1} + \left( \frac{2 \, v_{1} v_{2}}{v_{1}^{2} + v_{2}^{2} + v_{3}^{2}} \right) e_{2} + \left( \frac{2 \, v_{1} v_{3}}{v_{1}^{2} + v_{2}^{2} + v_{3}^{2}} \right) e_{3}$

$\displaystyle \verb|B:| x_{12} e_{1} e_{2} + x_{13} e_{1} e_{3} + x_{23} e_{2} e_{3}$

$\displaystyle \verb|B^2:| -x_{12}^{2} - x_{13}^{2} - x_{23}^{2}$

$\displaystyle \verb"|B|:" \sqrt{x_{12}^{2} + x_{13}^{2} + x_{23}^{2}}$

We would like to have a notion of inner and outer product, but so far we only have a notion of inner product on vectors specifically (via $v^TMw$ ). We can more general inner and outer products on arbitrary vectors via $\langle v,w \rangle = \frac{1}{2}(vw + wv)$ and $v\wedge w = \frac{1}{2}(vw - wv)$ , respectively, but it turns out that this is not the cleanest generalization (in terms of equational laws / how much it acts like our expectations of an inner and outer product). Instead, Dorst suggests using the outer product and left contraction operations as defined below:

In [6]:

def restrict(term, grade):
    a,b,c,d,e,f,g,h = term.dense_coefficient_list()
    if grade == 0:
        return build_ga_term(a,0,0,0,0,0,0,0)
    elif grade == 1:
        return build_ga_term(0,b,c,d,0,0,0,0)
    elif grade == 2:
        return build_ga_term(0,0,0,0,e,f,g,0)
    elif grade == 3:
        return build_ga_term(0,0,0,0,0,0,0,h)
    return build_ga_term(0,0,0,0,0,0,0,0)

def outer_product(C,D):
    out = 0*e1
    for r in range(4):
        for s in range(4):
            out += restrict(restrict(C,r)*restrict(D,s), r+s)
    return out

def left_contraction(C,D):
    out = 0*e1
    for r in range(4):
        for s in range(4):
            out += restrict(restrict(C,r)*restrict(D,s), s-r)
    return out

def right_contraction(C,D):
    out = 0*e1
    for r in range(4):
        for s in range(4):
            out += restrict(restrict(C,r)*restrict(D,s), r-s)
    return out

def scalar_product(C,D):
    out = 0*e1
    for r in range(4):
        for s in range(4):
            out += restrict(restrict(C,r)*restrict(D,s), 0)
    return out

show("v_|v: ", left_contraction(v,v))
show("B_|B: ", left_contraction(B,B))
show("B^2: ", B^2)
show("v cross v: ", outer_product(v,v))
show("v cross w: ", outer_product(v,w))

$\displaystyle \verb"v_|v:" v_{1}^{2} + v_{2}^{2} + v_{3}^{2}$

$\displaystyle \verb"B_|B:" -x_{12}^{2} - x_{13}^{2} - x_{23}^{2}$

$\displaystyle \verb|B^2:| -x_{12}^{2} - x_{13}^{2} - x_{23}^{2}$

$\displaystyle \verb|v|\verb| |\verb|cross|\verb| |\verb|w:| \left( -v_{2} w_{1} + v_{1} w_{2} \right) e_{1} e_{2} + \left( -v_{3} w_{1} + v_{1} w_{3} \right) e_{1} e_{3} + \left( -v_{3} w_{2} + v_{2} w_{3} \right) e_{2} e_{3}$

We mentioned earlier that $-ve_1v^{-1}$ will reflect $e_1$ across the hyperplane orthogonal to $v$ . We can then express rotations via an even number of reflections, such as $vwe_1w^{-1}v^{-1}$ . This product will expand into a bivector $B$ , which in general can be written as $B = x_{12}e_1e_2 + x_{23}e_2e_3 + x_{31}e_3e_1$ . So we can express rotations via $Be_1B^{-1}$ .

However, consider $B=e_1e_2$ . Then $Be_1B^{-1} = e_1e_2e_1e_2e_1 = -e_1^2e_2^2e_1 = -e_1$ . Each application of $e_1e_2$ or $e_2e_1$ rotated the target vector by a quarter turn, Making a collective turn of $\frac{1}{2}$ . But what if we wanted to do this same rotation incrementally — as in by $\theta$ degrees in the $e_1e_2$ plane?

We can use that $(e_1e_2)^2 = -1$ and taylor expansion to derive an analogue for Euler’s formula but with $e_ie_j$ instead of $i$ .

Below we check that $e^{\theta e_1e_2} = cos(\theta) + e_1e_2sin(\theta)$ up to degree 7.

In [7]:

def ga_exp(t, p):
    if p == 0: return 1
    if p % 2 == 1: return t * ga_exp(t, p-1)
    if p % 2 == 0: return ga_exp(t, p/2)^2 
    assert False
    
show("e^(theta*e1*e2): ")
show(sum([ga_exp(theta*e1*e2,i)/factorial(i) for i in range(8)]))
show("cos(theta) + e1*e2*sin(theta): ")
show(taylor(cos(theta),theta,0,7)+e1*e2*taylor(sin(theta),theta,0,8))

$\displaystyle \verb|e^(theta*e1*e2):|$

$\displaystyle \left( -\frac{1}{5040} \, \theta^{7} + \frac{1}{120} \, \theta^{5} - \frac{1}{6} \, \theta^{3} + \theta \right) e_{1} e_{2} - \frac{1}{720} \, \theta^{6} + \frac{1}{24} \, \theta^{4} - \frac{1}{2} \, \theta^{2} + 1$

Using a similar derivation, we can extend Euler’s formula to an arbitrary bivector $B$ .

Namely, we check that $e^{B\theta} = cos(|B|\theta) + \frac{B}{|B|}sin(|B|\theta)$ .

In [8]:

# Generalize further -- e^(B*theta) = cos(|B|theta) + (B/|B|)*sin(|B|theta) 
show("e^(B*theta)")
show(simp(sum([ga_exp(B*theta,i)/factorial(i) for i in range(5)])))
show("cos(|B|theta) + (B/|B|)*sin(|B|theta):")
show(simp(taylor(cos(norm(B)*theta), theta,0,4)+(B/norm(B))*(taylor(sin(norm(B)*theta),theta,0,3))))

$\displaystyle \verb|e^(B*theta)|$

$\displaystyle \left( -\frac{1}{6} \, \theta^{3} x_{12}^{3} - \frac{1}{6} \, \theta^{3} x_{12} x_{13}^{2} - \frac{1}{6} \, \theta^{3} x_{12} x_{23}^{2} + \theta x_{12} \right) e_{1} e_{2} + \left( -\frac{1}{6} \, \theta^{3} x_{13}^{3} - \frac{1}{6} \, \theta^{3} x_{13} x_{23}^{2} - \frac{1}{6} \, {\left(\theta^{3} x_{12}^{2} - 6 \, \theta\right)} x_{13} \right) e_{1} e_{3} + \left( -\frac{1}{6} \, \theta^{3} x_{23}^{3} - \frac{1}{6} \, {\left(\theta^{3} x_{12}^{2} + \theta^{3} x_{13}^{2} - 6 \, \theta\right)} x_{23} \right) e_{2} e_{3} + \frac{1}{24} \, \theta^{4} x_{12}^{4} + \frac{1}{24} \, \theta^{4} x_{13}^{4} + \frac{1}{24} \, \theta^{4} x_{23}^{4} - \frac{1}{2} \, \theta^{2} x_{12}^{2} + \frac{1}{12} \, {\left(\theta^{4} x_{12}^{2} - 6 \, \theta^{2}\right)} x_{13}^{2} + \frac{1}{12} \, {\left(\theta^{4} x_{12}^{2} + \theta^{4} x_{13}^{2} - 6 \, \theta^{2}\right)} x_{23}^{2} + 1$

$\displaystyle \verb"cos(|B|theta)"\verb" "\verb"+"\verb" "\verb"(B/|B|)*sin(|B|theta):"$

I. Parameterization of the Sphere

Now that we have an analogue for Euler’s formula in the Clifford algebra, we can automatically compute rotations.

First we calculate $e^{e_1e_2\frac{\theta}{2}}e_1e^{e_2e_1\frac{\theta}{2}} = cos(\theta)e_1 + sin(\theta)e_2$ .

Then we calculate a trigonometric parametrization for the radius $r$ sphere by rotating $re_1$ by $\theta$ in the $e_1e_2$ plane, followed by a rotation by $\phi$ in the $e_2e_3$ plane. Notice that this derivation will proceed identically in higher dimensions by continuing to apply rotations in the $e_ie_{i+1}$ planes.

In [13]:

def ga_eulers(B, angle): 
    return cos(norm(B)*angle)+(B/norm(B))*sin(norm(B)*angle)
def ga_eulers_unit_B(B, angle): 
    return simp(cos(angle)+B*sin(angle))
def rotate(B, angle, v): 
    expr = simp(ga_eulers(B, -angle/2)*v*ga_eulers(B, angle/2))
    return lift(lambda x: x.reduce_trig())(expr)

show("e1 rotated by theta in e1*e2 plane: ", simp(rotate(e1*e2, theta, e1)))
s = simp(rotate(e2*e3, phi, rotate(e1*e2, theta, r*e1)))
show("Sphere parametrization: ", s)

$\displaystyle \verb|e1|\verb| |\verb|rotated|\verb| |\verb|by|\verb| |\verb|theta|\verb| |\verb|in|\verb| |\verb|e1*e2|\verb| |\verb|plane:| \cos\left(\theta\right) e_{1} + \sin\left(\theta\right) e_{2}$

$\displaystyle \verb|Sphere|\verb| |\verb|parametrization:| r \cos\left(\theta\right) e_{1} + r \cos\left(\phi\right) \sin\left(\theta\right) e_{2} + r \sin\left(\phi\right) \sin\left(\theta\right) e_{3}$

II. Gaussian Curvature

Next we calculate Gaussian Curvature by following the formulas for the first and second fundamental forms of a surface. Notice that where the instructions call for an inner product, we instead use the left contraction that we defined earlier. Also we replace the cross product with the negative pseudoscalar ( $-e_1e_2…e_n$ ) times our generalized outer product.

In [10]:

def calculate_curvature(X,u,v):
    X_u = ga_diff(X, u)
    X_v = ga_diff(X, v)
    cross_Xu_Xv = simp(-e1*e2*e3*outer_product(X_u, X_v))
    n = simp(cross_Xu_Xv/norm(cross_Xu_Xv))
    show("Normal Vector: ", n)
    X_uu, X_uv, X_vv = ga_diff(X_u, u), ga_diff(X_u, v), ga_diff(X_v, v)
    L = simp(left_contraction(n, X_uu))
    M = simp(left_contraction(n, X_uv))
    N = simp(left_contraction(n, X_vv))
    detII = simp(L*N - M^2)
    show("Det 2nd Fundamental Form: ", detII)
    E = simp(left_contraction(X_u, X_u))
    F = simp(left_contraction(X_u, X_v))
    G = simp(left_contraction(X_v, X_v))
    detI = simp(E*G - F^2)
    show("Det 1st Fundamental Form: ", detI)
    curvature = to_scalar(detII) / to_scalar(detI)
    return curvature

show("Curvature: ", calculate_curvature(s,theta,phi))

$\displaystyle \verb|Normal|\verb| |\verb|Vector:| \cos\left(\theta\right) e_{1} + \cos\left(\phi\right) \sin\left(\theta\right) e_{2} + \sin\left(\phi\right) \sin\left(\theta\right) e_{3}$

$\displaystyle \verb|Det|\verb| |\verb|2nd|\verb| |\verb|Fundamental|\verb| |\verb|Form:| r^{2} \sin\left(\theta\right)^{2}$

$\displaystyle \verb|Det|\verb| |\verb|1st|\verb| |\verb|Fundamental|\verb| |\verb|Form:| r^{4} \sin\left(\theta\right)^{2}$

$\displaystyle \verb|Curvature:| \frac{1}{r^{2}}$

III. Volume of the Sphere

We want to calculate the volume of the sphere by integrating in spherical coordinates. In order to do so, we can integrate with an approriate extra factor for the relative shape of the infinitesimal volume elements we are integrating over. That extra factor is the determinant of the Jacobian of the surface parameterization. We can express this determinant using the outer product of the surface’s gradient vectors, and like earlier we convert it into a vector by multiplying by the negative pseudoscalar.

In [11]:

def calc_det_jacobian(S, v1, v2, v3):
    S_1 = ga_diff(S,v1)
    S_2 = ga_diff(S,v2)
    S_3 = ga_diff(S,v3)
    expr1 = outer_product(S_1, S_2)
    expr2 = outer_product(expr1, S_3)
    return simp(-e1*e2*e3*expr2)

detJ = to_scalar(calc_det_jacobian(s, r, theta, phi))
show("Det Jacobian of Coord System: ", detJ)

$\displaystyle \verb|Det|\verb| |\verb|Jacobian|\verb| |\verb|of|\verb| |\verb|Coord|\verb| |\verb|System:| r^{2} \sin\left(\theta\right)$

Now we can calculate the volume of the sphere by integrating!

In [12]:

def calc_volume(detJ, v1, v2, v3):
    expr1 = integral(detJ,v=r)
    expr2 = integral(expr1, v=theta, a=0, b=pi)
    expr3 = integral(expr2, v=phi, a=-pi, b=pi)
    return expr3

show("Volume of Sphere: ", calc_volume(detJ, r, theta, phi))

$\displaystyle \verb|Volume|\verb| |\verb|of|\verb| |\verb|Sphere:| \frac{4}{3} \, \pi r^{3}$

Future Work

In the future I plan to extend this notebook to obtain:

Formulas for surface of ellipse and torus
Volumes of these shapes by deriving the determinant of the Jacobian of the surface parameterization
Curvature for these shapes

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

Clifford Algebra Tutorial

Mathematical Preliminaries

Notebook Contents

Clifford Algebra Sage Setup

I. Parameterization of the Sphere

II. Gaussian Curvature

III. Volume of the Sphere

Future Work

Product

Resources

Company

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more, all in one place. Commercial Alternative to JupyterHub.

Clifford Algebra Tutorial

Mathematical Preliminaries

Notebook Contents

Clifford Algebra Sage Setup

I. Parameterization of the Sphere

II. Gaussian Curvature

III. Volume of the Sphere

Future Work

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.