Lecture 4:
Object Oriented Mathematica

John Armstrong

Data Types

Circle[ {x,y}, r]

The Circle function doesn’t do anything!

We’re just tagging a list with the name “Circle” to indicate how the data should be interpreted

The hard graft is done by the Graphics function which knows how to interpret the data

As the Graphics function demonstrates, you can write sophisticated code using data types

More Data Types

FullForm[ 2 +  3 I]
FullForm[ {1,2,3}]
FullForm[ Rectangle[{0,1},{2,3}]]

In fact, you developed three data types last week: a[i], b[i] and c[i,j] representing lines on a cubic surface

A Differential Geometry Library

There’s no coding needed to define a DataType, you just start using it.

torusChart = chart[
{(2+Cos[u])Cos[v], (2+Cos[u])Sin[v], Sin[u]},
{{u,0,2Pi},{v,0,2Pi}}]

This represents the inverse of a chart if one recalls the standard definition of a manifold.

The first parameter defines a mapping from (R^2) to (R^3)

The second parameter lists the variables we’re using and their ranges.

Working with a DataType

getMapping[c_chart]:= c[[1]]
getVar1[ c_chart]:= c[[2,1,1]]
getVar2[ c_chart]:= c[[2,2,1]]
getVars[ c_chart ] := {getVar1[c], getVar2[c]}
getRange1[ c_chart]:= c[[2,1]]
getRange2[ c_chart]:= c[[2,2]]
getScale[c_chart]:= Min[getRange1[c][[3]]-getRange1[c][[2]],
getRange2[c][[3]]-getRange2[c][[2]]]
getRangeMatrix[ c_chart]:= {c[[2,1,2;;3]],c[[2,2,2;;3]]}
plot[ c_chart ] := ParametricPlot3D[ Evaluate[getMapping[c]], 
Evaluate[getRange1[ c]],
Evaluate[getRange2[c]]]

* Note that the pattern c_chart matches only items with the Head chart
* The ;; syntax for creating sublists is new
* Note that we need to use Evaluate when calling ParametricPlot3D. It behaves very like Table.

Example

getVar1[ torusChart]
getVar2[ torusChart]
Show[plot[ torusChart], ImageSize->Large]

Creating some more charts

Let’s create 3 charts of the sphere representing spherical polar coordinates, the Mercator chart and the Lambert chart.

All 3 charts are “cylindrical charts” of the sphere so the main ingredient required is a height function.

sphericalPolarsHeight[y_] := -Cos[y]
mercatorHeight[ y_ ] := (-1 + Exp[2 y ])/(1 + Exp[2 y] ) 
lambertHeight[ y_ ] := y

Let’s create a cylindricalSpec data type to contain everything we need

mercatorSpec = cylindricalSpec[mercatorHeight, {-3,3}, thetaM, yM];
sphericalSpec = cylindricalSpec[sphericalPolarsHeight, {0,Pi}, thetaS, yS];
lambertSpec = cylindricalSpec[lambertHeight, {-1,1}, thetaL, yL];

So this contains the range of y values and the name for the coordinates. It’s a good idea to use different variable names in different charts.

Working with cylindricalSpec objects

We write some functions that make cylindricalSpec objects easy to use:

getHeightFunction[ s_cylindricalSpec ]:= s[[1]];
getRange[ s_cylindricalSpec ]:= s[[2]];
getThetaCoord[ s_cylindricalSpec ]:= s[[3]];
getYCoord[ s_cylindricalSpec ]:= s[[4]];
cylindricalProjection[ s_cylindricalSpec ]:= Module[{hf,r,theta,y},
hf = getHeightFunction[s];
r = getRange[s];
theta = getThetaCoord[s];
y = getYCoord[s];
chart[
{Cos[theta] PowerExpand[ Simplify[Sqrt[ 1 - hf[y]^2 ]]],
 Sin[theta] PowerExpand[ Simplify[Sqrt[ 1 - hf[y]^2 ]]],
hf[y]},
{{theta,0,2Pi}, {y,r[[1]], r[[2]]}}]
]

Creating and displaying our charts

sphericalPolars := cylindricalProjection[  sphericalSpec]
lambertProjection := cylindricalProjection[ lambertSpec]
mercatorProjection := cylindricalProjection[ mercatorSpec]
GraphicsGrid[{{plot[ sphericalPolars ],
plot[ lambertProjection ],
plot[ mercatorProjection ]}}, ImageSize->Full]

CountryData

earthDegrees = CountryData["World",{"Shape","Equirectangular"}];
grid = {Lighter[Lighter[Blue]],
Rectangle[{0,0},{2Pi,Pi}],
White,
Thick,
Table[ Line[{{i,0},{i,Pi}}],{i,0,2Pi, Pi/12}],
Table[ Line[{{0,i},{2Pi,i}}],{i,0,Pi, Pi/12}],
Thin};
Graphics /@ {earthDegrees, grid }

Combining the images

CountryData produces images with coordinates in degrees.
We use Scale and Translate to convert to radians.
We use PlotRangePadding of zero to eliminate unwanted margins

earth = Graphics[{grid,
Scale[Translate[earthDegrees[[1]],{180,90}],2Pi/360, {0,0}]
},{ImageSize->Large, PlotRangePadding->0}]

A bit of texture

plot[ c_chart, texture_ ] := ParametricPlot3D[ Evaluate[getMapping[c]],
Evaluate[ getRange1[c] ],
Evaluate[ getRange2[c] ],
PlotStyle->Texture[texture], Mesh->None, Lighting->"Neutral", Axes->None, Boxed->False, ImageSize->Large
]

Show[ plot[ sphericalPolars, earth ], ImageSize->Medium]

Trying to find Transition Functions

Solve[ mercatorHeight[y1] == sphericalPolarsHeight[y2], {y1}, Reals]

Solve[ mercatorHeight[y1] == sphericalPolarsHeight[y2] && y2>0 && y2<Pi, {y2},Reals]

Simplify[Solve[ mercatorHeight[y1] == sphericalPolarsHeight[y2] && y2>0 && y2<Pi, {y2},Reals]] /. {ConditionalExpression[x_, _] -> x}

Resulting function

On[Assert]
cylindricalTransition[ s1_cylindricalSpec,
					s2_cylindricalSpec ]:= Module[{ soln, tf, y1, y2},
y1 = getYCoord[s1];
y2 = getYCoord[s2];
soln = Simplify[ Solve[ getHeightFunction[s1][y1]==getHeightFunction[s2][y2]
	   && y2> getRange[s2][[1]] && y2<getRange[s2][[2]], y2, Reals]
];
soln = soln /. {ConditionalExpression[x_,_]->x};
Assert[ Length[soln ] == 1];
tf = {getThetaCoord[s1], soln[[1,1,2]]};
(* store the result *)
transitionFunction[ cylindricalProjection[s1],
cylindricalProjection[s2]] = tf;
tf
]

Compute and store all the transition functions we need

cylindricalTransition[ sphericalSpec, mercatorSpec ]

projectionSpecs= {sphericalSpec, mercatorSpec, lambertSpec };
Table[ If[i==j,{},cylindricalTransition[ projectionSpecs[[i]], projectionSpecs[[j]] ]],
{i,1,Length[ projectionSpecs] }, {j,1,Length[ projectionSpecs] }];

Definition[ transitionFunction]

Default behaviour

transitionFunction[ from_chart, to_chart]:= Throw[{ "No transition function defined between the charts", from, to}]
transitionFunction[ from_chart, from_chart]:= getVars[from]

So even though we can’t define a general transition function, we’ve created a useful lookup table.

transitionFunction[ sphericalPolars, mercatorProjection]

transitionFunction[ sphericalPolars, torusChart]

Using transition functions

transformImage[ 
image_,
fromChart_,
toChart_ ]:= Module[{function,vars},
vars = getVars[toChart];
function[{s_,t_}]:= 
transitionFunction[toChart,fromChart] /. {vars[[1]]->s, vars[[2]]->t};
ImageTransformation[
image, function ,
{DataRange->getRangeMatrix[fromChart],PlotRange->getRangeMatrix[toChart]}
]
];

So function is a function which performs the mapping encoded in the transition function. This function is just what the built in ImageTransformation function needs

Testing our code

transformImage[ earth, sphericalPolars, mercatorProjection]

The Lambert Projection

transformImage[ earth, sphericalPolars, lambertProjection]

We’ve put in quite a bit of effort here, but the end result is quite nice. We’ve got a transformImage function that will work with any image and any charts for which we’ve defined a transition function. It will get nicer and nicer...

earth3dAgain = plot[ lambertProjection, transformImage[ earth, sphericalPolars, lambertProjection]];
Show[earth3dAgain, ImageSize->Small]

Vector Fields

We’re following the usual order of definitions in differential geometry. Charts, then transition functions. So vector fields must be next on our agenda.

exampleVectorField = vectorField[{Sin[thetaS]Sin[yS], Sin[yS]},
sphericalPolars];

getComponents[ vectorField[ components_, chart_ ]]:= components;getChart[ vectorField[ components_, chart_ ]]:= chart;

Notice that a vector field is a pair consisting of a set of components and a chart.

The definition of a vector field on a manifold is as an equivalence class of such pairs (you may have seen other more algebraic or axiomatic definitions, but this is the most “intuitive” definition)

Plotting a Vector Field

plot[ vf_vectorField]:= Module[{field,chart, v1, v2,min1,max1,min2,max2, arrows, plotRange, factor, epsilon},
(* Some aesthetic parameters *)
epsilon = getScale[ getChart[vf]]/10.0;
factor = epsilon;

field = getComponents[vf];
chart = getChart[vf];
{v1,min1,max1}= getRange1[chart];
{v2,min2,max2}= getRange2[chart];
arrows = Table[ Evaluate[Arrow[ {{v1,v2}, {v1,v2}+factor*field }]], Evaluate[{v1,min1+epsilon/2,max1,epsilon}],
Evaluate[{v2, min2+epsilon/2, max2, epsilon}]];
plotRange = {{min1,max1},{min2,max2}};
Graphics[ {Darker[Red],Arrowheads[Medium],arrows}, PlotRange->plotRange,PlotRangeClipping->True, PlotRangePadding->0]
]

Testing the plot function

plot[ exampleVectorField]

Plotting a vector field in 3D

plot3D[ twoD_]:= Module[{texture,chart},
texture = plot[twoD];
chart=getChart[twoD];
plot[chart,texture]
]
plot3D[ exampleVectorField ]

Vector operations

You can do more with vectors than just plot them of course. For example you can add them

vectorField[ c1_, chart_ ] + vectorField[ c2_, chart_]:=
vectorField[ c1 + c2, chart ]
lambda * vectorField[ c1_, chart_ ]:=
vectorField[ c1 lambda, chart ] /; Element[ lambda, Reals]

Rats!

You now have two options.

* You can Unprotect[Plus] and Unprotect[Times]
* You can associated the definition with vectorField rather than with Plus and Times

vectorField /: vectorField[ c1_, chart_ ] + vectorField[ c2_, chart_]:=
vectorField[ c1 + c2, chart ]
vectorField /: lambda * vectorField[ c1_, chart_ ]:=
vectorField[ c1 lambda, chart ] /; Element[ lambda, Reals]

2 * exampleVectorField

When looking for rules to apply to an expression, Mathematica usually looks for rules associated with the Head of that expression. However, it will also look down a level for rules associated with parameters. Rules defined in this way are called “up values”. The types of rules we’ve used before are called “down values”

Up or Down?

It doesn’t make much difference logically whether you define a rule as an up value or a down value. It makes a performance difference, however. If our rule affects Plus then our rule will be briefly considered every time we add two objects, even just two numbers. However, if our rule affects only vector fields, then it will be ignored by Mathematica except when processing vector fields.

This means that if we were to use Unprotect and then associate our rule with Plus and Times we will slow down Mathematica’s general processing speed. This would be bad and is one reason why Plus and Times are protected in the first place.

As a general rule of thumb, it is better to associate a definition with the least frequently used part of the definition. Since vectorFields come up less often than addition in general, it’s better to associate the definition with vectorField rather than with Plus

Code reuse in action

torusVectorField1 = vectorField[{1/2,0}, torusChart];
torusVectorField2 = vectorField[{0,1/2}, torusChart];

GraphicsGrid[ {{plot3D[ torusVectorField1], 
plot3D[ torusVectorField2],
plot3D[ torusVectorField1 + torusVectorField2]}}, ImageSize->Full]

transform[ vf_vectorField, toChart_ ] := Module[{yAsFunctionOfX,xAsFunctionOfY,
xVars, fromChart, dim, fromComponents, toComponents, dyByDx, toComponentsAsFunctionsOfX, toComponentsAsFunctionsOfY,varChange},
fromComponents = getComponents[ vf ];
fromChart = getChart[ vf];
xVars = getVars[ fromChart];
yAsFunctionOfX = transitionFunction[ fromChart, toChart ];
xAsFunctionOfY = transitionFunction[ toChart, fromChart];
dim = Length[ xVars];
dyByDx = Table[ D[ yAsFunctionOfX[[i]], xVars[[j]] ], {i,1,dim},{j,1,dim}];
toComponentsAsFunctionsOfX = dyByDx . fromComponents;
varChange = Table[ xVars[[i]]->xAsFunctionOfY[[i]],{i,1,dim}];
toComponentsAsFunctionsOfY = toComponentsAsFunctionsOfX /. varChange;
Simplify[ vectorField[toComponentsAsFunctionsOfY, toChart ]]
];

Testing the transform

version1 = plot[ transform[ exampleVectorField, mercatorProjection ]];
version2 = transformImage[ plot[ exampleVectorField ], sphericalPolars, mercatorProjection ];
GraphicsGrid[ {{version1, version2 }}, ImageSize->Full]

The angle looks right...

One Forms

getChart[ oneForm[ components_, chart_ ]]:= chart;getComponents[ oneForm[ components_, chart_ ]]:= components;

r2Chart = chart[{x,y,0},{{x,-3,3},{y,-3,3}}];
exampleOneForm = oneForm[{-y,x},r2Chart];

oneForm /: oneForm[ c1_, chart_] + oneForm[c2_, chart_ ]:=
	oneForm[ c1 + c2, chart ];
oneForm /: lambda_ * oneForm[ c1_, chart_]:=
	oneForm[ lambda * c1, chart ] /; Element[lambda, Reals]

Today’s most complex function...

formDisc[ form_, centre_, alpha_, epsilon_  ]:=Module[{kernelVector, orthogonalVector, offsetVector, numLines,intersectionPoints, lines, colors,colorList, colouredLines},
colouredLines = If[ Norm[form]>0,
kernelVector=NullSpace[{form}][[1]];
orthogonalVector = form ;
offsetVector = alpha  * orthogonalVector/ (orthogonalVector . form);
numLines = Min[Ceiling[epsilon / Norm[offsetVector]],100];
intersectionPoints[ radius_, kernelVector_, point_ ]:= Module[
{a, solns, lhs},
lhs = point  + a kernelVector;
solns = Solve[ lhs[[1]]^2 + lhs[[2]]^2== radius^2 , {a}, Reals];
(lhs /. # &) /@ solns
];
lines = Table[ Line[intersectionPoints[epsilon,kernelVector,
offsetVector*i]], {i,-numLines,numLines}];
colors = Darker /@ {Blue, Cyan, Green };
colorList = Table[ colors[[Mod[i,3]+1]], {i,-numLines,numLines}];
Riffle[colorList,lines],
{}];
 Translate[{ EdgeForm[Darker[Blue]],White,Disk[{0,0}, epsilon], colouredLines}, centre]
]

Drawing covectors

Graphics[{formDisc[ {20.0,10.0}, {0,0}, 1, 1.0 ],
Thick,Red,Arrow[{{0,0},{1/2,2/3}}]}]

The red arrow is a vector. To compute the effect of the form on the vector, count the number of lines the vector crosses in the direction Blue/Cyan/Green. If it crosses the other way, you get a negative number.

plot[  f_oneForm]:=Module[{field,chart, v1, v2,min1,max1,min2,max2, discs, plotRange,graphics, epsilon, alpha, proportion},
(* Aesthetic parameters *)
epsilon = getScale[ getChart[f] ]/10.0;
alpha = epsilon/6.0;
proportion = 9.0/10.0;
(* Main work *)
field = getComponents[f];
chart = getChart[f];
{v1,min1,max1}= getRange1[chart];
{v2,min2,max2}= getRange2[chart];
discs = Table[ formDisc[ field,{v1,v2}, alpha, proportion epsilon/2 ], Evaluate[{v1,min1+epsilon/2,max1,epsilon}],
Evaluate[{v2, min2+epsilon/2, max2, epsilon}]];
(* Fit everything neatly to plot range *)
plotRange = {{min1,max1},{min2,max2}};
Graphics[ {discs},  PlotRange->plotRange,ImageSize->Large,PlotRangeClipping->True, PlotRangePadding->0]
]

Psychedelia

plot[ exampleOneForm ]

Polymorphism

We’ve written plot so it works with vector fields and forms
We’ve written getChart so it works with vector fields and forms
We’ve written plot3D so it only uses plot and getChart
We deduce that plot3D will work just as well for forms

plot3D[ exampleOneForm]

Enjoying ourselves

torusForm1 = oneForm[{1,0}, torusChart];
torusForm2 = oneForm[{0,1}, torusChart];
GraphicsGrid[ {{plot3D[ torusForm1], plot3D[ torusForm2]}},
ImageSize->Full]

Improving plot and plot3D to display multiple images

plot[{items__}]:= Show[{plot /@ {items}}]
plot3D[ {items__ }]:= Module[{texture,chart},
chart = getChart[ {items}[[1]] ];
texture = plot[ {items}];
plot[chart,texture]
];
plot3D[ {torusForm1,torusVectorField1}]

Curves

We’re on a roll now. Let’s define a curve data type and plot function

getPoint[ curve[ point_, range_, chart_ ] ]:= point;
getRange[ curve[ point_, range_, chart_ ] ]:= range;
getChart[ curve[ point_, range_, chart_ ] ]:= chart;

circle = curve[ {3/2Sin[theta], 3/2Cos[theta]},{theta, 0, 2Pi}, r2Chart];

plot[ c_curve ]:= Module[{epsilon},
epsilon = getScale[ getChart[ c] ]/500;
Graphics[ {Darker[Red], Thick,Arrow[ Table[ Evaluate[getPoint[c]], Evaluate[ Join[getRange[c], {epsilon}]] ]]},
PlotRangePadding->0, PlotRangeClipping->True]
]

Integration

Show[  plot[ exampleOneForm],plot[circle] ]

De Rham cohomology

torusCurve1 = curve[{t, 17 Pi/10},{t,0,2 Pi},torusChart];
torusCurve2 = curve[{ Pi/10,t},{t,0,2Pi},torusChart];
GraphicsGrid[ {{plot3D[ { torusForm1, torusCurve1}],
plot3D[ { torusForm2,torusCurve2}]}}, ImageSize->Full]

Penrose Stairs

Penrose Stairs provide a visual proof that (H^1)((S^1))= 0

Summary

The notion of a “data type” allows us to build sophisticated programs that manipulate them

Defining “polymorphic” data types that reuse functions with the same name to produce different effects are especially valuable.

The programming paradigm where you pass round polymorphic data types is called “object oriented programming”

Exercises. (You probably only need to do one to get the point)

Define a 0-form data type including vector operations and plot function.
Define a function “pair” that takes a vector field and a one form and computes the 0-form resulting from pairing them.
Define a function d that acts on 0-forms to produce a one form.
Plot the following 0-form and its exterior derivative: sin(u) cos(v)
Extend transform to 1-forms
Etc. etc.

Exercise

Create your own Prog Rock album cover

Rules for your talks

Ten minutes per team member

You must work in teams of two or three.

Talks must use some computing. Visualization is quite sufficient.

It takes longer than you might expect to get software to work

You can use any computer language you like.

The talks should all be reasonably different from each other.

Extended Exercises

Write your own 3D graphics library (using graphics 3D is cheating)

Write an elliptic curve addition library

Write your own solver for Diophantine equations of the following form (using Solve or Reduce is cheating)

a x^2+b y^2+c z^2+d x y+e y z+f z x == 0

Compute the Riemann curvature tensor of a manifold.

3D Graphics Library

A viewpoint consists of a point representing your eye and a plane representing the canvas. Given another point P construct the line from your eye to P and find where it intersects the canvas. You should be able to write a function “project” that takes a 3-d point and a viewpoint and returns 2-d coordinates for the intersection.

You can now plot a point. It shouldn’t be long till you can plot a line segment and then a cube.

Writing an elliptic Curve mutiplication Library

GIven an elliptic curve and an arbitrary point 0 on that curve, write a function that performs ellitpic curve addition

Write a function that performs multiplication by an integer - i.e. given a point A computes nA

Can you get Mathematica to prove that addition is associative?

Write your own Diophantine equation solver

The more you know about Diophantine equations, the more you can solve. Solving the Diophantine equation a x + b y = c should be within everyone’s capabilities. Solving x^2 + b y^2 = c z^2 is much more challenging, but you can look up algorithms for doing this.

A basic ingredient that is interesting is quadratic reciprocity. You can use the Jacobi Symbol to write an efficient “quadratic reciprocity tester”. See https://brilliant.org/wiki/jacobi-symbol/ for a description of the algorithm.

This is a building block for writing an algorithm to test whether certain quadratic diophantine equations have solutions. See https://math.stackexchange.com/questions/25640/existence-of-solutions-to-diophantine-quadratic-form for a discussion and some references.

Computing the Riemann curvature tensor

By making heavy use of D, Table and . you should be able to
compute the Christoffel Symbols and then the curvature
tensor.

Make sure that you don’t discard the Christoffel symbols too
soon, because you want your code to be fast.

Can you use this to find some interesting surfaces in R3? For
example what surfaces of revolution have constant curvature?