Symbolic and Numeric Geometric Algebra

Terje Vold

Continuum Technology, Inc.

Introduction

Algebraic Systems for Geometric Information

There are many algebraic systems that are used to represent geometric information.  Some of the better known ones are  

Hamilton's quaternions (1843)
Grassman's exterior algebra (1844)
Gibbs' "vector algebra" (c 1884)
Clifford's "geometric algebra" (c 1885)

Gibbs was a towering figure in physics and chemistry (especially statistical mechanics!), and his vector algebra has become the best known.

Clifford's promising career was cut short by an early death at age 34, but his geometric algebra contains all three of the others, is simple to use and understand, and is practical.

Geometric Algebra in a Nutshell

Elements of the algebra

The fundamental elements are scalars (real or complex numbers) and vectors. 
Secondary or derived quantities are exterior products.

For example, for 3D space with basis vectors  (Subscript[e, 1]), (Subscript[e, 2]), (Subscript[e, 3]), 
we also have bivectors  (Subscript[e, 1])⋀(Subscript[e, 2]), (Subscript[e, 2])⋀(Subscript[e, 3]), (Subscript[e, 3])⋀(Subscript[e, 1])
and the trivector  I = (Subscript[e, 1])⋀(Subscript[e, 2])⋀(Subscript[e, 3]).  
Elements of the algebra are generically called multivectors.  
The number of vector factors is called the grade of the multivector.

Elements of the algebra form a linear space of multivectors over the real (or complex) numbers:  
	if    a  and  b    are multivectors, 
	so is  (Subscript[c, 1])a + (Subscript[c, 2])b
	for any real (or complex) numbers (Subscript[c, 1])and (Subscript[c, 2]).

The geometric product

The product of a geometric algebra - the "geometric product" - is not generally commutative,

		ab ≠ ba 

but is generally associative, 

		a(bc) = (ab)c .

Many results from real analysis can be used if modified to retain orders of factors.

The geometric product between a scalar and anything else is the normal product.

The geometric product between a vector and a vector is defined to be the sum of the conventional inner product and the conventional exterior or outer product: 

		ab = a·b + a⋀b .

From the definition of the geometric product between vectors 
and the associativity of the geometric product for all elements,
the geometric product between any two elements is uniquely determined.

Allowing sums of dissimilar object is similar to allowing sums of real and imaginary numbers in complex algebra.  In fact, various subalgebras of geometric algebras are isomorphic to the conventional complex numbers.

Inner and Outer Products

The geometric product is taken to be the fundamental product,
and the inner and outer products written in terms of the geometric product. 

In particular, for any vectors  a  and  b,  
		 a·b = (1/2)(ab + ba)
and		a⋀b = (1/2)(ab - ba)

For multivectors a and b with one or both grades greater than 1, a more general definition applies.

Gibbs' Vector Algebra in a Nutshell

The basic elements of Gibbs' vector algebra are vectors only; they form a linear space of vectors over the reals.   

There are two distinct products between vectors a and b:  the inner product a·b and the cross product a×b.  

One product is commutative, one is anticommutative. 

Neither product is associative, e.g., a×(b×c) ≠ (a×b)×c.   The dot product is not even defined for more than two operands, e.g.,  a·b·c.

Why Use Geometric Algebra?

Conventional Gibbs' vector algebra is contained within geometric algebra: 
it has the same dot product  a·b, 
and the cross product is related to the outer product by  a × b = -I a⋀b

But associativity makes geometric algebra more powerful, systematic, and compact.

Examples of additional capabilities from elementary geometry:   
a) For any unit vector  (OverHat[a])  and vector (OverVector[v]),  (OverHat[a]) (OverVector[v]) (OverHat[a])  is the result of reflecting (OverVector[v]) in the plane perpendicular to  a.   
b) For any unit vectors (OverHat[a]) and (OverHat[b]) and vector (OverVector[v]),  (OverHat[b]) (OverHat[a]) (OverVector[v]) (OverHat[a]) (OverHat[b])  is the result of rotating (OverVector[v]) in the plane containing a and b by an angle equal to twice the angle from a to b.    Matrix elements of the rotation for a component representation are (Subscript[R, ij])=(Subscript[(OverHat[e]), i])·((OverHat[b]) (OverHat[a]) (Subscript[(OverHat[e]), j]) (OverHat[a]) (OverHat[b])).

Other geometries - such as projective and conformal - can be represented with geometric algebra with useful results.

Geometric differentiation and integration are powerful generalizations of the conventional operations in real analysis.

Some of Nature's most basic physical laws - in particular, Maxwell's equations - are exceptionally simple in geometric algebra.

Suggested reference:  Geometric Algebra for Physicists, by Doran and Lasenby, Cambridge University Press, 2003.

Maxwell's Equation

Maxwell in 4D Spacetime

The geometric derivative in 4D spacetime with orthonormal basis vectors {(Subscript[(OverHat[γ]), t]), (Subscript[(OverHat[γ]), x]), (Subscript[(OverHat[γ]), y]), (Subscript[(OverHat[γ]), z])}, for a function F[r] of the vector variable  r = t(Subscript[(OverHat[γ]), t]) + x(Subscript[(OverHat[γ]), x]) + y(Subscript[(OverHat[γ]), y]) + z(Subscript[(OverHat[γ]), z]) ,  is a straightforward generalization of the derivative of a scalar function of a scalar variable x:

□F[r] = (((OverHat[γ])^t)(∂_t) + ((OverHat[γ])^x)(∂_x) + ((OverHat[γ])^y)(∂_y) + ((OverHat[γ])^z)(∂_z))F[r]

The simplest possible differential equation for an unknown mutlivector field F and a specified multivector field S in spacetime is

	□F[r] = S[r]

For a bivector field F and a vector field S, this is Maxwell's equation in spacetime.

A Spacetime Split

We can project into parts parallel to and perpendicular to any chosen unit vector u by multipling by 1 = uu.  For example, multiplying on the left using  u = (Subscript[(OverHat[γ]), t])  we have 

	□F[r] =((Subscript[(OverHat[γ]), t])(Subscript[(OverHat[γ]), t]))□F[r] = (Subscript[(OverHat[γ]), t])((Subscript[(OverHat[γ]), t])□F[r]) = (Subscript[(OverHat[γ]), t])((∂_t) + (Subscript[(OverHat[γ]), x])⋀(Subscript[(OverHat[γ]), t])(∂_x) + (Subscript[(OverHat[γ]), y])⋀(Subscript[(OverHat[γ]), t])(∂_y) + (Subscript[(OverHat[γ]), z])⋀(Subscript[(OverHat[γ]), t])(∂_z))F[r]	
and
	S[r] =((Subscript[(OverHat[γ]), t])(Subscript[(OverHat[γ]), t]))S[r] = (Subscript[(OverHat[γ]), t])((Subscript[(OverHat[γ]), t])S[r]) = (Subscript[(OverHat[γ]), t])((Subscript[S, t]) + (Subscript[(OverHat[γ]), x])⋀(Subscript[(OverHat[γ]), t])(Subscript[S, x]) + (Subscript[(OverHat[γ]), y])⋀(Subscript[(OverHat[γ]), t])(Subscript[S, y]) + (Subscript[(OverHat[γ]), z])⋀(Subscript[(OverHat[γ]), t])(Subscript[S, z]))

It turns out that the even subalgebra of 4D spacetime is isomorphic to the algebra of 3D space, with the identification 
	(Subscript[(OverHat[e]), x]) = (Subscript[(OverHat[γ]), x])⋀(Subscript[(OverHat[γ]), t]),   (Subscript[(OverHat[e]), y]) = (Subscript[(OverHat[γ]), y])⋀(Subscript[(OverHat[γ]), t]),  (Subscript[(OverHat[e]), z]) = (Subscript[(OverHat[γ]), z])⋀(Subscript[(OverHat[γ]), t]).

Then
	□F[r] = (Subscript[(OverHat[γ]), t])((∂_t) + (Subscript[(OverHat[e]), x])(∂_x) + (Subscript[(OverHat[e]), y])(∂_y) + (Subscript[(OverHat[e]), z])(∂_z))F[r]		
and
	S[r] = ((Subscript[(OverHat[γ]), t])(Subscript[(OverHat[γ]), t]))S[r] = (Subscript[(OverHat[γ]), t])((Subscript[(OverHat[γ]), t])S[r]) = (Subscript[(OverHat[γ]), t])((Subscript[(OverHat[γ]), t])·S[r]+(Subscript[(OverHat[γ]), t])⋀S[r])= (Subscript[(OverHat[γ]), t])(ρ - (OverVector[J]))
	
where 	 (OverVector[J])≡S⋀(Subscript[γ, t])=-(Subscript[γ, t])⋀S=(Subscript[(OverHat[e]), x])(Subscript[J, x])+(Subscript[(OverHat[e]), y])(Subscript[J, y])+(Subscript[(OverHat[e]), z])(Subscript[J, z])
and	 ρ≡S·(Subscript[γ, t])=(Subscript[γ, t])·S

Maxwell in 3D Space

We define the geometric derivative (OverVector[∇]) in 3D space in the usual way,

(OverVector[∇])F[r] = ((Subscript[(OverHat[e]), x])(∂_x) + (Subscript[(OverHat[e]), y])(∂_y) + (Subscript[(OverHat[e]), z])(∂_z))F[r]

With these results, we can write Maxwell's equation in 4D spacetime, times (Subscript[γ, t]), as  

	((∂_t)+(OverVector[∇]))F[r] = ρ[r] - (OverVector[J])[r]	

where
	F[r] =(OverVector[E])[r]+I(OverVector[B])[r] .

Comparison of Maxwell's Equations in Homogeneous Linear Polarizable Media

We compare Maxwell's equation(s) in the important case of homogeneous linear polarizable media, 
with v equal to the speed of light in the medium and ϵ equal to the permittivity, and using SI units.

Conventional Cartesian components and real differential calculus:

	 (∂_x)(Subscript[E, x])+(∂_y)(Subscript[E, y])+(∂_z)(Subscript[E, z]) = (1/ϵ)ρ
	 (∂_y)(Subscript[B, z])-(∂_z)(Subscript[B, y])  = (1/(v^2))((1/ϵ) (Subscript[J, x])+ (∂_t)(Subscript[E, x]))
	 (∂_z)(Subscript[B, x])-(∂_x)(Subscript[B, z])  = (1/(v^2))((1/ϵ) (Subscript[J, y])+ (∂_t)(Subscript[E, y]))
	 (∂_x)(Subscript[B, y])-(∂_y)(Subscript[B, x])  = (1/(v^2))((1/ϵ) (Subscript[J, z])+ (∂_t)(Subscript[E, z]))
	 (∂_y)(Subscript[E, z])-(∂_z)(Subscript[E, y])  - (∂_t)(Subscript[B, x]) = 0
	 (∂_z)(Subscript[E, x])-(∂_x)(Subscript[E, z])  - (∂_t)(Subscript[B, y]) = 0
	 (∂_x)(Subscript[E, y])-(∂_y)(Subscript[E, x])  - (∂_t)(Subscript[B, z]) = 0
	 (∂_x)(Subscript[B, x])+(∂_y)(Subscript[B, y])+(∂_z)(Subscript[B, z]) = 0

Gibbs' vector algebra and vector calculus:  

	 (OverVector[∇])·(OverVector[E]) = (1/ϵ)ρ
	 (OverVector[∇])×(OverVector[B])  =(1/(v^2))((1/ϵ) (OverVector[J])+ (∂_t)(OverVector[E]))
	 (OverVector[∇])×(OverVector[E]) + (∂_t)(OverVector[B]) = 0
	 (OverVector[∇])·(OverVector[B]) = 0

where  (OverVector[E])=(Subscript[E, x])(Subscript[(OverHat[e]), x])+(Subscript[E, y])(Subscript[(OverHat[e]), y])+(Subscript[E, z])(Subscript[(OverHat[e]), z]),  (OverVector[B])=(Subscript[B, x])(Subscript[(OverHat[e]), x])+(Subscript[B, y])(Subscript[(OverHat[e]), y])+(Subscript[B, z])(Subscript[(OverHat[e]), z]),  and  (OverVector[J])=(Subscript[J, x])(Subscript[(OverHat[e]), x])+(Subscript[J, y])(Subscript[(OverHat[e]), y])+(Subscript[J, z])(Subscript[(OverHat[e]), z]).

Geometric algebra and calculus:
 
		 ((OverVector[∇])+(1/v)(∂_t))F= S  
		 
where	 F≡(OverVector[E])+I v (OverVector[B])  and  S ≡(1/ϵ)((Subscript[ρ, f])-(1/v)(Subscript[(OverVector[J]), f])) .

Green's Function for Maxwell in 3D Space

Green's theorem, applied to Maxwell's equation, relates the field at any point r' in or on the surface of volume V to an integral of the field over points r on the surface of the volume:

F[r']=(∫_(∂V)) (Subscript[g, k])[r-r'] da (OverHat[n]) F[r]  -  (∫_V) (Subscript[g, k])[r-r'] dv S[r]

where
		(Subscript[g, k])[r]=((OverVector[∇])-i k)(Subscript[ϕ, k])[(OverVector[r])]
with
		(Subscript[ϕ, k])[r]=-(1/(4π r))(e^(ikr))

or equivalently
		(Subscript[g, k])[r]=(1/(4π (r^2)))(e^(ikr))((OverHat[r])+ikr(1-(OverHat[r]))) .

For simplicity, consider the case of no source (S=0) in volume V. 

Then Green's theorem for Maxwell's equation is written using geometric algebra as 

	F[r']=(∫_(∂V)) (Subscript[g, k])[r-r'] da (OverHat[n]) F[r]		

where		(Subscript[g, k])[r]=(1/(4π (r^2)))(e^(ikr))((OverHat[r])+ikr(1-(OverHat[r]))).

Green's theorem for Maxwell's equation may be written using vector algebra as the two equations 

	(OverVector[E])[x]=(∫_(∂V)) (ik((OverVector[n])×(OverVector[B]))G[x-x']+((OverVector[n])'×(OverVector[E]))×(OverVector[∇])'G[x-x']+((OverVector[n])'·(OverVector[E]))(OverVector[∇])'G[x-x'])da' 
and
	(OverVector[B])[x]=(∫_(∂V)) (-ik((OverVector[n])×(OverVector[E]))G[x-x']+((OverVector[n])'×(OverVector[B]))×(OverVector[∇])'G[x-x']+((OverVector[n])'·(OverVector[B]))(OverVector[∇])'G[x-x'])da'

where		G = (1/(4π r))(e^(ikr)) .

The expression using geometric algebra embodies or uses known patterns or regularities, making it easier to discover other patterns that are more likely to be masked by the irregularities of vector notation.

Computational Geometric Algebra

Geometric Product

We will use SmallCircle to represent the geometric product, as in

a∘b //FullForm

Attributes[SmallCircle]

We set the Attributes to be the same as for NonCommutativeMultiply:

a**b //FullForm

Attributes[NonCommutativeMultiply]

SetAttributes[SmallCircle,{Flat,OneIdentity}]

Attributes[SmallCircle]

Basis Elements

We define a set of orthonormal basis vectors of 3D space as

ex::usage="ex is one of a set of orthonormal basis vectors for 3D space.";
ey::usage="ey is one of a set of orthonormal basis vectors for 3D space.";
ez::usage="ez is one of a set of orthonormal basis vectors for 3D space.";

spaceBasisVectors::usage="spaceBasisVectors equals the List {ex,ey,ez} of orthonormal basis vectors for 3D space.";

spaceBasisVectors={ex,ey,ez};

A basis for the entire geometric algebra of 3D space is

spaceBasis::usage="spaceBasis is a List of basis elements of the geometric algebra of 3D space.";

spaceBasis={1}~Join~(If[Length[#]>1,SmallCircle@@#,#[[1]]]& /@Union[Flatten[Table[Union[{a1,a2,a3}],{a1,spaceBasisVectors},{a2,spaceBasisVectors},{a3,spaceBasisVectors}],2]])

spaceBasis={1,ex,ey,ez,ey∘ez,ez∘ex,ex∘ey,ex∘ey∘ez};

Geometric Products between Basis Elements

Rules for the geometric product between basis elements

We tell the geometric product to place factors in a standard order:

Clear[SmallCircle]

SmallCircle[a___,b_,c_,d___]:=-SmallCircle[a,c,b,d]/;MemberQ[spaceBasisVectors,b]∧MemberQ[spaceBasisVectors,c]∧Order[c,b]==1

ex∘ey∘ex

Now we use the metric of space to define products between a basis vector and itself:

a_∘a_:=1/;MemberQ[spaceBasisVectors,a]

ex∘ey∘ex

We want the geometric product to default to the usual Times if at least one of two factors is a number:

SmallCircle[a_,b_]:=a b /; NumberQ[a]∨NumberQ[b]

ex∘ey∘ex

But it doesn't yet work if one factor has Head Times:

ex∘(7ey)∘ey

We want scalar (real number) coefficients of vectors to be commuted out to the left of all vector factors:

SmallCircle[a___,Times[b_,c_],d___]:=b SmallCircle[a,c,d]/;NumberQ[b]

ex∘(7ey)∘ey

Table of geometric products between basis elements

Now a Table of geometric products between basis elements can be calculated:

gpTable=Table[a∘b,{a,spaceBasis},{b,spaceBasis}];

TableForm[ gpTable,TableHeadings->{spaceBasis,spaceBasis}]

The Geometric Product between List Representations of Multivectors

Representation of multivectors by Lists of coefficients of basis vectors

We will represent a multivector m as a List of coefficients of basis vectors.  For example,

multivectorA=As+Avx ex+Avy ey+Avz ez+Abx ey∘ez+Aby ez∘ex+Abz ex∘ey+At ex∘ey∘ez

It is generally convenient to group together basis elements, or the corresponding coefficients, with the same grade:

groupGrades[{s_,x_,y_,z_,a_,b_,c_,t_}]:={s,{x,y,z},{a,b,c},t}

This grouping explicitly shows the four grades,  {scalar, vector, bivector, trivector}.

Then, for example, the List of basis elements in the standard order,

spaceBasis

may be displayed as

spaceBasis//groupGrades

coefficientList[mv_]:=groupGrades[{1,1,1,1,1,-1,1,1}*Flatten[CoefficientArrays[Flatten[multivectorA],Drop[{1,1,1,1,1,-1,1,1}*spaceBasis,1]][[{1,2}]]]]

coefficientList[multivectorA]

Examples:  
A position r in space is a vector.  All multivector components are zero except for it grade-1 (vector) components.  It is represented as a List of multivector basis coefficients by

r={0,{Xx,Yy,Zz},{0,0,0},0};

The value of an electromagnetic field at some particular point and time is a sum of vector (electric) and bivector (magnetic) components, with scalar or trivector components always equal to zero:

F={0,{Ex,Ey,Ez},{Bx,By,Bz},0};

The geometric product for Flat basis coefficient Lists

The geometric product is linear.  So we can represent it with an array with 3 indices, that we will call gp. 
Then if a and b are List-representations of multivectors, the matrix product

	a.gp.b 
	
is the List-representation of the geometric product ab of a and b.

We can generate gp from gpTable by applying an appropriate List of rules for mapping from the basis element symbols to the associated Lists:

symbolsToLists=
{1->{1,0,0,0,0,0,0,0},
-1->{-1,0,0,0,0,0,0,0},
ex∘ey∘ez->{0,0,0,0,0,0,0,1},
-ex∘ey∘ez->{0,0,0,0,0,0,0,-1},
ey∘ez->{0,0,0,0,1,0,0,0},
ez∘ex->{0,0,0,0,0,1,0,0},
ex∘ey->{0,0,0,0,0,0,1,0},
-ey∘ez->{0,0,0,0,-1,0,0,0},
-ez∘ex->{0,0,0,0,0,-1,0,0},
-ex∘ey->{0,0,0,0,0,0,-1,0},
ex->{0,1,0,0,0,0,0,0},
ey->{0,0,1,0,0,0,0,0},
ez->{0,0,0,1,0,0,0,0},
-ex->{0,-1,0,0,0,0,0,0},
-ey->{0,0,-1,0,0,0,0,0},
-ez->{0,0,0,-1,0,0,0,0}
};

Now we can define the "tensor" gp:

gp=Transpose[Replace[gpTable,symbolsToLists,{2}],{1,3,2}]

Spot check:

{0,0,0,0,1,0,0,0}.gp.{0,0,0,0,1,0,0,0}

The geometric product for grouped basis coefficient Lists

Now we can define the geometric product to act on basis coefficient Lists that represent multivectors:

SmallCircle[List[a__],List[b__]]:=groupGrades[Flatten[List[a]].gp.Flatten[List[b]]]

For example, for two vectors a and b,

a={0,{ax,ay,az},{0,0,0},0};
b={0,{bx,by,bz},{0,0,0},0};

a//MatrixForm

b//MatrixForm

the geometric product is

a∘b//MatrixForm

The geometric product between two general multivectors is

{sa,{xa,ya,za},{aa,ba,ca},ta}∘{sb,{xb,yb,zb},{ab,bb,cb},tb}//MatrixForm

Inner and Outer Products

Definitions

The inner product will be represented by CenterDot,

a·b//FullForm

and the outer product by Wedge,

a⋀b//FullForm

The outer product needs Attributes set:

SetAttributes[Wedge, {Flat, OneIdentity}]

We need to select grades of a general multivector; we do this with the help of

gradeMasks={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}

Now we can define the inner and outer products:

CenterDot[a_,b_]:=Sum[((a*gradeMasks[[i]])∘(b*gradeMasks[[j]]))*gradeMasks[[Abs[i-j]+1]],{i,4},{j,4}]

Wedge[a_, b_] := Sum[((a*gradeMasks[[i]])∘(b*gradeMasks[[g - i + 1]]))*gradeMasks[[g]], {g, 4}, {i, g}]

Examples

For example, for any two vectors

a//MatrixForm
b//MatrixForm

the inner product between them is

a·b//MatrixForm

and the outer product is

a⋀b//MatrixForm

As another example, the inner product  between a vector and a bivector is

{0,{ax,ay,az},{0,0,0},0}·{0,{0,0,0},{Bx,By,Bz},0}

The cross product

It is convenient to define the unit trivector III :

III = {0, {0, 0, 0}, {0, 0, 0}, 1};

Now the cross product  a×b  may be written

-III∘(a⋀b)  //MatrixForm

Greens function for Maxwell

Greens function in traditional notation

The Greens function for Maxwell's equation is

g[r]=(1/(4π (r^3)))(e^(i k r)) ((OverVector[r]) + i k r (r-(OverVector[r])))

To express this in Mathematica using a basis coefficient representation, we first write a few definitions :

Some functions generalized to act on multivectors

The "Reverse" in geometric algebra refers to reversing the order of all vector factors.  This is analogous to complex conjugation.  In 3D space, the result of this is that the bivector components change sign.  So we define

gaReverse[{s_,{x_,y_,z_},{a_,b_,c_},t_}]:={s,{x,y,z},{-a,-b,-c},-t}

For example, for the vector

a

the geometric product of a with its geometric algebraic Reverse is

a∘gaReverse[a]

The Magnitude is a generalization of the magnitude of a complex number:

Clear[gaMag]

gaMag[a_]:=((((a∘gaReverse[a])[[1]]))^(1/2))*{1,{0,0,0},{0,0,0},0}

gaMag[{s,{0,0,0},{bx,by,bz},0}]

gaPower[a_,0]:={1,{0,0,0},{0,0,0},0};
gaPower[a_,n_]:=a∘gaPower[a,n-1]/;n>0;
gaPower[a_,n_]:=gaPower[gaReverse[a],-n](((gaMag[a][[1]]))^(2n))/;n<0;

Look at a few examples:

a

gaPower[a,-1]

gaPower[a,-2]

gaPower[a,-3]

gaMag[a]

gaPower[gaMag[a],-3]

The Exponential function is uniquely defined in geometric algebra by its power series expression; however, there is no general finite expression for all possible multivector arguments.  But for the simple cases of a scalar argument, such as any real or complex number c time the magnitude gaMag of a multivector m, it is equal to the conventional definition :

gaExp[{a_,{0,0,0},{0,0,0},0}]:={Exp[a],{0,0,0},{0,0,0},0}

gaExp[c_ gaMag[m_]]:={Exp[c gaMag[m][[1]]],{0,0,0},{0,0,0},0}

gaExp[{s_,{0,0,0},bv_,0}]:=Exp[s]{Cos[Sqrt[bv.bv]],{0,0,0},(bv/Sqrt[bv.bv])Sin[Sqrt[bv.bv]],0}/;bv.bv≠0

gaExp[{s_,{0,0,0},bv_,0}]:=Exp[s]{1,{0,0,0},{0,0,0},0}/;bv.bv==0

gaExp[0]=1;

Look at an example:

gaExp[{s,{0,0,0},{bx,by,bz},0}]

Green's function in Standard Form

Now we can write Green' s function in Mathematica as

greens[r_,k_]:=(1/(4π))gaPower[gaMag[r],-3]∘ gaExp[I k gaMag[r]]∘(r+I k gaMag[r]∘(gaMag[r]-r))

r

greens[r,k]//MatrixForm

Green's integral

One element of a numerical approximation to Green's integral, associated with an area element

da={0,{dax,day,daz},{0,0,0},0}

is

greens[r,k]∘da∘F//MatrixForm

Coefficients of E and B in Green's integral

With numerical values for position r, wavevector magnitude k, and surface area element a,  such as

r1={0,{3.1,4.2,5.3},{0,0,0},0}

k1=0.001

da1={0,{0.11,0.12,0.13},{0,0,0},0}

and leaving the electromagnetic field F symbolic -

F

it is more apparent that this expression is a linear function of F,

greens[r1,k1]∘da1∘F//MatrixForm

To extract the coefficients, we need a Flattened structure,

Flatten[greens[r1,k1]∘da1∘F]//MatrixForm

Now we can use CoefficientArrays to calculate a matrix of coefficients of the field

F

by writing

Normal[CoefficientArrays[#,{Ex,Ey,Ez,Bx,By,Bz}][[2]]]& /@ Flatten[greens[r1,k1]∘da1∘F]//MatrixForm

A Simple Illustration:  The field of an oscillating current loop

Position vectors and current vectors of the elements of a current loop

We will do all calculations in a List representation, but it is convenient to enter data using symbolic basis vector names; so we define

basisVectorsToLists={e1->{0,{1,0,0},{0,0,0},0},e2->{0,{0,1,0},{0,0,0},0},e3->{0,{0,0,1},{0,0,0},0}};

R is a Rotor that rotates any multivector (including just a vector) by angle ϕ in a plane containing  e1 and e2:

R[ϕ_]:=gaExp[(ϕ/2) e1∘e2]

rs is the position r of a point-like piece of the source current, which we take to be a ring of 8 points:

rs[ϕ_]:=R[-ϕ]∘e1∘R[ϕ]/.basisVectorsToLists

Spot check:

rs[0]

rs[0.1]

rs[π/2]

rs[π]

rs[2π]

We take the source at each of these points to be a current of unit magnitude.  This is represented by a vector equal to the position rs, rotated by π/2 in the plane containing e1 and e2:

current[ϕ_]:=rs[ϕ]∘e1∘e2  /.basisVectorsToLists

current[0.1]

current[π]

Greens integral

Now we can write a simple expression to calculate the Green's integral over the current loop to give the resulting field F :

currentLoopField[r_,k_]:=Sum[  greens[r-rs[ϕ],k]∘current[ϕ],
{ϕ,(1/16)2π,(15/16)2π,(1/8)2π}]

We calculate a table of field values at points {x, 0, 0} on the x axis:

currentLoopFieldTableA=Table[currentLoopField[{0,{x,0,0},{0,0,0},0},1]//N,{x,0,10,1/3}]

Plot of resulting fields

Now we Plot these field values.   We need to load VectorFieldPlots:

Needs["VectorFieldPlots`"]

We will use one expression from this package:

?ListVectorFieldPlot3D

The list of points (just the vector parts) is

pointListA=Table[{x,0,0},{x,0,10,1/3}]

The list of electric field values is

electricTableA=currentLoopFieldTableA[[All,2]]

The list of magnetic field values is

magneticTableA=currentLoopFieldTableA[[All,3]]

We'll use a very short cylinder to indicate the location of the current loop:

?Cylinder

Now we finally Plot these electric and magnetic field values due to the loop of oscillating current:

EPlot=ListVectorFieldPlot3D[Transpose[{pointListA,10^1 Re[electricTableA]}],VectorHeads->True,ColorFunction->(Hue[2/3]&)];
BPlot=ListVectorFieldPlot3D[Transpose[{pointListA,10^1 Re[magneticTableA]}],VectorHeads->True,ColorFunction->(Hue[1/3]&)];
Show[EPlot,BPlot,Graphics3D[{FaceForm[Opacity[0]],Cylinder[{-{0,0,1},{0,0,1}}/100]}]]

An animated version is below.

Animate[
EPlot=ListVectorFieldPlot3D[Transpose[{pointListA,10^1 Re[Exp[-2π I u]electricTableA]}],VectorHeads->True,ColorFunction->(Hue[2/3]&)];
BPlot=ListVectorFieldPlot3D[Transpose[{pointListA,10^1 Re[Exp[-2π I u]magneticTableA]}],VectorHeads->True,ColorFunction->(Hue[1/3]&)];
Show[EPlot,BPlot,Graphics3D[{FaceForm[Opacity[0]],Cylinder[{-{0,0,1},{0,0,1}}/100]}],PlotRange->{{-1,10},{-4,4},{-7,7}}],
{u,0,15/16,1/16},AnimationRate->2]

Examples of Computational Results in Electrodynamics

Hollow Magnetically Permeable Sphere in an Applied Oscillating Magnetic Field

These Plots show electric and magnetic field vectors around Helmholz coils with oscillating current and a small hollow magnetically weakly permeable sphere.

The first Plot shows the Helmholz coils and permeable sphere from a distance of a few coil diameters.

The second Plot shows the field near one short section of the coils.

The next Plot shows the magnetic field along the x axis, inside and outside the permeable spherical shell.

The Plot below shows the magnetic field and scaled electric field on the x axis, near the surface of the sphere.

The Plot below shows the field at large distances:  the coils and magnetic sphere are near the left edge of the Plotted vectors, and are so small as to be imperceptible.

Electromagnetic Field and Poynting Vector Near a Spherical Microlens

In the Plot below, the black arrows are Poynting vectors.  The large sphere is a small glass lens, with radius equal to 1 wavelength.  A plane wave is incident from the left.  The small light gray-violet sphere is positioned at the focal point as calculated by a paraxial approximation, and has a diameter equal to a diffraction-limited estimated spot size.

The plane of the Plot below contains the focal point of the lens that is in the Plot above, and is perpendicular to the incident radiation (the average direction of the black Poynting vectors).  The green arrows represent the electric field and the blue arrows represent the magnetic field.

Poynting Vector Near a Convex Microlens

This Plot shows a thin slice of a 3D simulation of a plane wave incident on a convex glass microlens.  The radius of curvature of the microlens is 1 wavelength.  The black arrows are the Poynting vectors.