From 673566bff1c8962b6cf03bda46b1b463e2eea40d Mon Sep 17 00:00:00 2001 From: Michael Reed <18372368+chakravala@users.noreply.github.com> Date: Fri, 13 Dec 2019 22:19:19 -0500 Subject: [PATCH] initialized algebra-of-space tutorial --- README.md | 4 +- docs/make.jl | 1 + docs/src/references.md | 1 + docs/src/tutorials/algebra-of-space.md | 214 +++++++++++++++++++++++++ docs/src/tutorials/quick-start.md | 43 ++++- src/multivectors.jl | 4 +- 6 files changed, 262 insertions(+), 5 deletions(-) create mode 100644 docs/src/tutorials/algebra-of-space.md diff --git a/README.md b/README.md index fbd966b..8d7dd67 100644 --- a/README.md +++ b/README.md @@ -121,7 +121,6 @@ The direct sum operator `⊕` can be used to join spaces (alternatively `+`), an The direct sum of a `VectorBundle` and its dual `V⊕V'` represents the full mother space `V*`. In addition to the direct-sum operation, several other operations are supported, such as `∪,∩,⊆,⊇` for set operations. Due to the design of the `VectorBundle` dispatch, these operations enable code optimizations at compile-time provided by the bit parameters. -**Remark**. Although some of the operations like `∪` and `⊕` are similar and sometimes result in the same values, the `union` and `sum` are entirely different operations in general. Calling manifolds with sets of indices constructs the subspace representations. Given `M(s::Int...)` one can encode `SubManifold{length(s),M,s}` with induced orthogonal space, such that computing unions of submanifolds is done by inspecting the parameter `s`. @@ -145,6 +144,8 @@ Additionally, a universal unit volume element can be specified in terms of `Line The universal interoperability of `LinearAlgebra.UniformScaling` as a pseudoscalar element which takes on the `VectorBundle` form of any other `TensorAlgebra` element is handled globally. This enables the usage of `I` from `LinearAlgebra` as a universal pseudoscalar element. +More information about `AbstractTensors` is available at https://github.com/chakravala/AbstractTensors.jl + # Grassmann elements and geometric algebra Λ(V) The Grassmann `Basis` elements `vₖ` and `wᵏ` are linearly independent vector and covector elements of `V`, while the Leibniz `Operator` elements `∂ₖ` are partial tangent derivations and `ϵᵏ` are dependent functions of the `tangent` manifold. @@ -267,6 +268,7 @@ Full `MultiVector` elements are not representable when `ExtendedAlgebra` is used The sparse representations are a work in progress to be improved with time. ## References +* Michael Reed, [Differential geometric algebra with Leibniz and Grassmann](https://crucialflow.com/grassmann-juliacon-2019.pdf). 2019. * Emil Artin, [Geometric Algebra](https://archive.org/details/geometricalgebra033556mbp). 1957. * John Browne, [Grassmann Algebra, Volume 1: Foundations](https://www.grassmannalgebra.com/). 2011. * C. Doran, D. Hestenes, F. Sommen, and N. Van Acker, [Lie groups as spin groups](http://geocalc.clas.asu.edu/pdf/LGasSG.pdf), J. Math Phys. (1993) diff --git a/docs/make.jl b/docs/make.jl index 4387eeb..8512b71 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -18,6 +18,7 @@ makedocs( "Tutorials" => Any[ "tutorials/install.md", "tutorials/quick-start.md", + "tutorials/algebra-of-space.md", "tutorials/forms.md" ], "References" => "references.md" diff --git a/docs/src/references.md b/docs/src/references.md index f3e4964..82ff0d6 100644 --- a/docs/src/references.md +++ b/docs/src/references.md @@ -7,6 +7,7 @@ [![Gitter](https://badges.gitter.im/Grassmann-jl/community.svg)](https://gitter.im/Grassmann-jl/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge) [![BiVector](https://img.shields.io/badge/bivector.net-Discourse-blueviolet)](https://bivector.net) +* Michael Reed, [Differential geometric algebra with Leibniz and Grassmann](https://crucialflow.com/grassmann-juliacon-2019.pdf). 2019. * Emil Artin, [Geometric Algebra](https://archive.org/details/geometricalgebra033556mbp). 1957. * John Browne, [Grassmann Algebra, Volume 1: Foundations](https://www.grassmannalgebra.com/). 2011. * C. Doran, D. Hestenes, F. Sommen, and N. Van Acker, [Lie groups as spin groups](http://geocalc.clas.asu.edu/pdf/LGasSG.pdf), J. Math Phys. (1993) diff --git a/docs/src/tutorials/algebra-of-space.md b/docs/src/tutorials/algebra-of-space.md new file mode 100644 index 0000000..c256afd --- /dev/null +++ b/docs/src/tutorials/algebra-of-space.md @@ -0,0 +1,214 @@ +# The Algebra of Space (G3) + +This notebook is an adaptation from the [clifford](https://clifford.readthedocs.io/en/latest/TheAlgebraOfSpaceG3.html) python documentation. + +Import `Grassmann` and instantiate a three dimensional geometric algebra + +```Julia +julia> using Grassmann + +julia> basis"3" +(⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃) +``` + +Given a three dimensional GA with the orthonormal basis ``v_i\cdot v_j = \delta_{ij}``, the basis consists of scalars, three vectors, three bivectors, and a trivector. +```math +\{\underbrace{v}_{\mbox{scalar}},\qquad\underbrace{v_{1},v_{2},v_{3}}_{\mbox{vectors}},\qquad\underbrace{v_{12},v_{23},v_{13}}_{\mbox{bivectors}},\qquad\underbrace{v_{123}}_{\mbox{trivector}}\} +``` +The `@basis` macro declares the algebra and assigns the `Basis` elements to local variables. The `Grassmann.Algebra` can also be assigned to `G3` as +```Julia +julia> G3 = Λ(3) +Grassmann.Algebra{⟨+++⟩,8}(v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃) +``` +You may wish to explicitly assign the blades to variables like so, +```Julia +e1 = G3.v1 +e2 = G3.v2 +# etc ... +``` +Or, if you're lazy you can use the macro with different local names +```Julia +julia> @basis ℝ^3 E e +(⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃) + +julia> e3, e123 +(v₃, v₁₂₃) +``` + +## Basics + +The basic products are available + +```Julia +julia> v1*v2 # geometric product +v₁₂ + +julia> v1|v2 # inner product +0v + +julia> v1 ∧ v2 # exterior product +v₁₂ + +julia> v1 ∧ v2 ∧ v3 # even more exterior products +v₁₂₃ +``` + +Multivectors can be defined in terms of the basis blades. For example, you can construct a rotor as a sum of a scalar and a bivector, like so +```Julia +julia> θ = π/4 +0.7853981633974483 + +julia> R = cos(θ) - sin(θ)*v23 +0.7071067811865476 - 0.7071067811865475v₂₃ + +julia> R = ~exp(θ*v23) +0.7071067811865476 - 0.7071067811865475v₂₃ +``` +You can also mix grades without any reason +```Julia +julia> A = 1 + 2v1 + 3v12 + 4v123 +1 + 2v₁ + 3v₁₂ + 4v₁₂₃ +``` +The reversion operator is accomplished with the tilde `~` in front of the `MultiVector` on which it acts +```Julia +julia> ~A +1 + 2v₁ - 3v₁₂ - 4v₁₂₃ +``` +Taking a projection into a specific `grade` of a `MultiVector` is usually written ``\langle A\rangle_n`` and can be done using the soft brackets, like so +```Julia +julia> A(0) +1v + +julia> A(1) +2v₁ + 0v₂ + 0v₃ + +julia> A(2) +3v₁₂ + 0v₁₃ + 0v₂₃ +``` +Using the reversion and grade projection operators, we can define the magnitude of `A` as ``|A|^2 = \langle{\tilde A A\rangle`` +```Julia +julia> ~A*A +30 + 4v₁ + 12v₂ + 24v₃ + +julia> scalar(ans) +30v +``` +This is done in the `abs` and `abs2` operators +```Julia +julia> abs2(A) +30 + 4v₁ + 12v₂ + 24v₃ + +julia> scalar(ans) +30v +``` +The dual of a multivector `A` can be defined as ``\tilde AI``, where `I` is the pseudoscalar for the geometric algebra. In `G3`, the dual of a vector is a bivector: +```Julia +julia> a = 1v1 + 2v2 + 3v3 +1v₁ + 2v₂ + 3v₃ + +julia> ⋆a +3v₁₂ - 2v₁₃ + 1v₂₃ +``` + +## Reflections + +Reflecting a vector ``c`` about a normalized vector ``n`` is pretty simple, ``c\mapsto -ncn`` +```Julia +julia> c = v1+v2+v3 # a vector +1v₁ + 1v₂ + 1v₃ + +julia> n = v1 # the reflector +v₁ + +julia> -n*c*n # reflect a in hyperplane normal to n +0.0 - 1.0v₁ + 1.0v₂ + 1.0v₃ +``` +Because we have the `inv` available, we can equally well reflect in un-normalized vectors using ``a\mapsto n^{-1}an`` +```Julia +julia> a = v1+v2+v3 # the vector +1v₁ + 1v₂ + 1v₃ + +julia> n = 3v1 # the reflector +3v₁ + +julia> inv(n)*a*n +0.0 + 1.0v₁ - 1.0v₂ - 1.0v₃ + +julia> n\a*n +0.0 + 1.0v₁ - 1.0v₂ - 1.0v₃ +``` +Reflections can also be made with respect to the hyperplane normal to the vector, in which case the formula is negated. + +## Rotations + +A vector can be rotated using the formula ``a\mapsto \tilde R aR``, where `R` is a rotor. A rotor can be defined by multiple reflections, ``R = mn`` or by a plane and an angle ``R = e^{\theta B/2}``. +For example, +```Julia +julia> R = exp(π/4*v12) +0.7071067811865476 + 0.7071067811865475v₁₂ + +julia> ~R*v1*R +0.0 + 2.220446049250313e-16v₁ + 1.0v₂ +``` +Maybe we want to define a function which can return rotor of some angle ``\theta`` in the ``v_{12}``-plane, ``R_{12} = e^{\theta v_{12}/2}`` +```Julia +R12(θ) = exp(θ/2*v12) +``` +And use it like this +```Julia +julia> R = R12(π/2) +0.7071067811865476 + 0.7071067811865475v₁₂ + +julia> a = v1+v2+v3 +1v₁ + 1v₂ + 1v₃ + +julia> ~R*a*R +0.0 - 0.9999999999999997v₁ + 1.0v₂ + 1.0v₃ +``` +You might as well make the angle argument a bivector, so that you can control the plane of rotation as well as the angle +```Julia +R_B(B) = exp(B/2) +``` +Then you could do +```Julia +julia> Rxy = R_B(π/4*v12) +0.9238795325112867 + 0.3826834323650898v₁₂ + +julia> Ryz = R_B(π/5*v23) +0.9510565162951535 + 0.3090169943749474v₂₃ +``` +or +```Julia +julia> R_B(π/6*(v23+v12)) +0.9322404424570728 + 0.25585909935689327v₁₂ + 0.25585909935689327v₂₃ +``` +Maybe you want to define a function which returns a *function* that enacts a specified rotation, ``f(B) = a\mapsto e^{B/2}\\ae^{B/2}``. +This just saves you having to write out the sandwich product, which is nice if you are cascading a bunch of rotors, like so +```Julia +R_factory(B) = (R = exp(B/2); a -> ~R*a*R) +Rxy = R_factory(π/3*v12) +Ryz = R_factory(π/3*v23) +Rxz = R_factory(π/3*v13) +``` +Then you can do things like +```Julia +julia> R = R_factory(π/6*(v23+v12)) # this returns a function +#7 (generic function with 1 method) + +julia> R(a) # which acts on a vector +0.0 + 0.5229556000177233v₁ + 0.7381444851051178v₂ + 1.4770443999822769v₃ + 2.7755575615628914e-17v₁₂₃ + +julia> Rxy(Ryz(Rxz(a))) +0.0 + 0.40849364905389035v₁ - 0.6584936490538903v₂ + 1.5490381056766584v₃ +``` +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 some vector ``x``. +```Julia +julia> R_seq(args...) = (R = prod(exp.(args./2)); a -> ~R*a*R) +R_seq (generic function with 1 method) + +julia> R = R_seq(π/2*v23, π/2*v12, v1) +#11 (generic function with 1 method) + +julia> R(v1) +2.220446049250313e-16 + 3.469446951953614e-16v₁ + 0.9999999999999996v₂ - 1.3877787807814457e-17v₃ - 5.551115123125783e-17v₂₃ + 2.7755575615628914e-17v₁₂₃ +``` diff --git a/docs/src/tutorials/quick-start.md b/docs/src/tutorials/quick-start.md index 997b848..18a3641 100644 --- a/docs/src/tutorials/quick-start.md +++ b/docs/src/tutorials/quick-start.md @@ -1,3 +1,42 @@ -# Quick start +# Quick start (G2) -Add some quick start tutorials here. +Import the `Grassmann` package and instantiate a two-dimensional algebra (G2), + +```Julia +julia> using Grassmann + +julia> @basis ℝ^2 +(⟨++⟩, v, v₁, v₂, v₁₂) + +julia> v1*v2 # geometric product +v₁₂ + +julia> v1|v2 # inner product +0v + +julia> v1∧v2 # exterior product +v₁₂ +``` + +## Reflection + +```Julia +julia> a = v1+v2 +1v₁ + 1v₂ + +julia> n = v1 +v₁ + +julia> -n*a/n # reflect a in hyperplane normal to n +0.0 - 1.0v₁ + 1.0v₂ +``` + +## Rotation + +```Julia +julia> R = exp(π/4*v12) +0.7071067811865476 + 0.7071067811865475v₁₂ + +julia> R*v1*~R +0.0 + 2.220446049250313e-16v₁ - 1.0v₂ +``` diff --git a/src/multivectors.jl b/src/multivectors.jl index eb10c92..c790bc2 100644 --- a/src/multivectors.jl +++ b/src/multivectors.jl @@ -227,9 +227,9 @@ setindex!(m::MultiVector{V,T} where V,k::T,i::Int,j::Int) where T = (m[i][j] = k Base.firstindex(m::MultiVector) = 0 Base.lastindex(m::MultiVector{V,T} where T) where V = ndims(V) -(m::MultiVector{V,T})(g::Int,b::Type{B}=Chain) where {T,V,B} = m(Val{g}(),b) +(m::MultiVector{V,T})(g::Int) where {T,V,B} = m(Val{g}()) function (m::MultiVector{V,T})(::Val{g}) where {V,T,g,B} - Chain{V,T,g}(m[g]) + Chain{V,g,T}(m[g]) end function (m::MultiVector{V,T})(g::Int,i::Int) where {V,T,B} Simplex{V,g,Basis{V}(indexbasis(ndims(V),g)[i]),T}(m[g][i])