# Differential Geometry with GNU/Linux Maxima and Geomview: The Monkey Saddle Surface with Embedded Curve

## Introductory Spiel

Geometry impresses most people as being a highly visual subject, but opening any textbook on differential geometry (DG) one will find only an endless thicket of abstract mathmatical symbols, usually inconsistent from text to text, with only a small smattering of visual graphics if any at all.  Even actual concrete examples, i.e. practical problems, are rare to find. Mathematical rigor is, beyond doubt, of prime importance but the true beauty, I believe, of DG can only be conveyed through practical examples and, above all, vivid imagery. The difficulty presumably lies in the extremely tedious and complex nature of problem solving in DG. This is likely why the only practical examples that are given are limited to simple, and nearly useless, cases like the Euclidian plane (anti-example?) or the cylinder and sphere.

The rise of computer algebra software (CAS), however, has now made this tedium a thing of the past.  Using Maxima CAS, a product of FOSS, as I will attempt to show, can allow differential geometry to be explored in a way that reveals its true power and beauty by nearly effortlessy plowing through the complex calculations and graph production.

This demo is not intended to be a primer of any sort on DG.  The subject is far too vast to be covered in a few web pages.  I only hope to convey an appreciation for the use of Maxima or other CAS in making DG more intuitive by introducing practical examples and graphics. The reader should already be somewhat versed in the rudimets of DG as I will devote little or no effort in explaining any concepts.

The reader is encouraged to download the actual Maxima workbook for this demo and to execute it on his machine.  Only in this way can the graphics be best appreciated.

## The Monkey Saddle and Curves Thereon

We first initialize the Maxima environment:

kill(all)$

load(vect)$ load(draw)$ load(eigen)$ load(ctensor)$

set_draw_defaults(head_angle=10, head_length=0.03, proportional_axes=xyz, nticks=800, xyplane=0, surface_hide=true, line_width=2, view=[62, 140])$

### Define the surface, the monkey saddle

As usual in DG, we consider a mapping from the u,v plane to the x,y,z space.  Geometrically, the monkey saddle, a relatively complex surface, has the following mapping:

\[ f(u, v) = [x, y, z] = [u, v, u^3-3*u*v^2] \]

In these Maxima demos, we will always use \(\mathbf f\) as our surface mapping and \(\mathbf f\) as our curve function, with derivatives indicated by a following letter:

E.g.:
\[ \frac{\partial \mathbf{f}}{\partial u} \to fu\] 
\[ \frac{\mathrm{d}\mathbf{c}}{\mathrm{d}t} \to ct \]

f: [u,v,u^3-3*u*v^2];
fu:diff(f,u);
fv:diff(f,v);

Check if our surface is regular by taking the cross product of the partial derivatives: \[ \frac{\partial \mathbf{f}}{\partial u} \times \frac{\partial \mathbf{f}}{\partial v} \]

ratsimp(express(fu~fv));

This is a regular surface because the cross product is nowhere 0.

### Define a curve on the surface parameterized with t

We now map a linear function, \( u = v^2 \), onto the surface.  We also compute several derivatives with repect to t as they will be needed later.   The first, second, and third derivatives are indicated, respctively, by ct, ctt, cttt:

/* define u,t in terms of t for later use */
tu:t$ tv:t^2$
c:ev(f, u=tu, v=tv);
ct:diff(c,t);
ctt:diff(c,t,2);
cttt:diff(c,t,3);

The curve "velocity," or the length of the first derivative, will also be needed later:

vel:ratsimp(sqrt(ct.ct));

Now comes a plot of the curve in the surface:

f[3];

wxdraw3d(color=blue, line_width=1, grid=true, xrange=[-2,2], yrange=[-2,2], 
    parametric_surface(f[1], f[2], f[3], u,-1,1,v,-1,1),
    color=red, line_width=2, parametric(c[1], c[2], c[3], t,-1,1))$

To better appreciate this surface and curve, a narrated Geomview animation is provided.  Click on the above image or click here (opens in new window).

### Curvature parameters of the curve on the surface

Before we determine the basic curvature parameters of this curve, we need to find the arc length to establish an arc length parameterization:

integrate(sqrt(ct.ct),t);

Maxima cannot solve this integral as the solution is not expressible with standard math functions.  Even with a solution, function inversion to produce the arc-length parameterization would be problematic.  We therefore use the formulas for abitrary speed curves to determine the curvature parameters.

### The Darboux moving frame, the unit tangent, T, and unit normal, N, vectors.

Appropriate to a curve embedded within a surface is the Darboux frame which expresses curve properties in terms of surface parameters, namely the geodesic and normal curvatures and the geodesic torsion.
Because the Darboux frame is essentially a rotated version of the Frenet-Serret frame, we need to determine also the Frenet unit vectors.

Unit tangent to curve, \( \mathbf{T} = \frac {c'(t)}{\lvert c'(t) \rvert} \)

T:(fullratsimp(ct/sqrt(ct.ct)));

Check

ratsimp(T.T);

The B, or binormal, vector: \( \frac{c'(t) \times c''(t)}{\lvert c'(t) \times c''(t) \rvert} \)

/* define tmp1 = ct x ctt  for later use*/
tmp1:fullratsimp(express(ct~ctt))$
B:fullratsimp(tmp1/sqrt(tmp1.tmp1));

Check:

ratsimp(B.B);

Check if T and B are perpendicular:

ratsimp(B.T);

The unit normal: \( N = B \times T )\

N:ratsimp(factor(express(B~T)));

Check if the Frenet frame vectors are mutually perpendicular:

ratsimp(N.N);
ratsimp(N.T);
ratsimp(N.B);

We now can determine the Darboux Frame 

/* the unit surface normal and its expression along the curve*/
n:express(fu~fv)$
n:fullratsimp(n/sqrt(n.n));
nt:ev(n, u=tu, v=tv);
/* check */
ratsimp(n.n);
ratsimp(nt.nt);

U:ratsimp(factor(express(nt~T)));
ratsimp(U.U);

/* cos(θ) using two methods*/
radcan(N.U);
radcan(B.nt);

### Curvature, \( \kappa \), and the geodesic, \( \kappa g \), and normal, \( \kappa n \), curvatures 

\( \kappa = \frac {\lvert c'(t) \times c''(t) \rvert}{\lvert c'(t) \rvert} \)

κ:(factor(sqrt(tmp1.tmp1))/(vel^3));

/* geodesic curvature */
κg:(ratsimp(ct.(radcan(express(ctt~nt)))))/vel^3;

/* geodesic curvature second method */

/* normal curvature */
/* project curvature,  \( \kappa \cdot N \), onto the surface normal, \( n \) */
κn:ratsimp((κ*N).nt);

A check of the equality: κ^2 = κn^2 + κg^2

fullratsimp(factor(κ^2));
fullratsimp(κn^2+κg^2);

Now comes a plot of both the geodesic and normal curvature.  Note that of the "action" occurs in the interval [-1, 1] with both curvatures falling to zero beyond.
The normal curvature expresses the rotation of the Darboux frame about the U vector while the geodesic curvature expresses the rotation about the normal, N, vector. 

wxdraw2d(line_width=1, background_color=light_yellow,
    xlabel="Parameter, t", ylabel="Geodesic Curvature, κg", title="Geodesic Curvature",
    proportional_axes=none, grid=true,
    explicit(κg, t, -2, 2))$

wxdraw2d(line_width=1, background_color=light_yellow,
    xlabel="Parameter, t", ylabel="Normal Curvature, κn", title="Normal Curvature",
    proportional_axes=none, grid=true,
    explicit(κn, t, -2, 2))$

### Torsion, τ, and the geodesic torsion, τg  

Although not strictly necessary for our purposes, we determine the torsion, τ, of the free curve.

\( \tau = \frac{(c'(t) \times c''(t)) \cdot c'''(t)}{{\lvert c'(t) \times c''(t) \rvert}^2} \)

τ:ratsimp(tmp1.cttt)/ratsimp(tmp1.tmp1);

The geodesic torsion, τg, expresses the rotation of Darboux frame about the tangent, T, vector.

ntt:ratsimp(factor(diff(nt,t)))$
τg:factor((ratsimp(ct.(radcan(express(ntt~nt)))))/vel^2);

Now the geodesoc torsion is plotted.  Again note that most of the "twisting occurs within the interval [-1, 1].

wxdraw2d(line_width=1, background_color=light_yellow,
    xlabel="Parameter, t", ylabel="Geodeic Torsion, τg", title="Geodesic Torsion",
    proportional_axes=none, grid=true,
    explicit(τg, t, -2, 2))$

Plot the Darboux frame at several points along the curve

Put Geomview graphic and animation here

The following narrated animation, from Geomview, shows the movement of the frame and the rotation of the osculating plane.

For a Geomview graphic click the above image or click here (opens in new window).

### Curvature Parameters of the Surface (Monkey Saddle)

### We now plot the surface unit normal, n, the curve unit normal, N, and the tangent plane, Tp

The general surface tangent plane, \( T_p \), with dummy parameters \(u_p\) and \(v_p\) specifying the point \( p \):

Tp:[f[1]+up*fu[1]+vp*fv[1], f[2]+up*fu[2]+vp*fv[2], f[3]+up*fu[3]+vp*fv[3]];

/* for plotting along the curve*/
ev(Tp, u=tu, v=tv);

A plot the tangent plane, and both surface and curve normals, at a specific point, \( u=0.8, v=0.64 \):

s1:0.2$ s2:s1^2$
Tpp:ev(Tp, u=s1, v=s2)$
wxdraw3d(color=gray30, line_width=1, parametric_surface(f[1], f[2], f[3], u,-1,1,v,-1,1), line_width=2, color=red, parametric(c[1], c[2], c[3], t,-1,1),
color=blue, view=[62, 140], line_width=1, parametric_surface(Tpp[1], Tpp[2], Tpp[3], up, -1/4, 1/4, vp, -1/4, 1/4),
line_width=2, color=purple, vector(ev(c, t:s1), ev(n, u:s1, v:s2)),
color=forest_green, vector(ev(c, t:s1), float(ev(N, t:s1)))
)$

Note that the surface normal, \(\mathbf n \), and the curve normal, \(\mathbf N \), are not the same.

For a narrated Geomview graphic click the above image or click here (opens in new window).

### First differential form (metric), g

We denote the metric matrix as \(g\):

E:fu.fu$
F:fu.fv$
G:fv.fv$
g:matrix([E, F],[F, G]);

ratsimp(determinant(g));

We will also require the metric inverse, \(ginv\) later:

ginv:ratsimp(invert(g));

Check. Should give the identity matrix:

ratsimp(g.ginv);

### Second differential form, \( b \)

Note: The second differential form variables are usually designated L, M, N, but here we use NN to avoid stepping on the curve unit normal variable, N.

L:diff(fu,u).n$
M:diff(fu,v).n$
NN:diff(fv,v).n$
b:matrix([L, M], [M, NN]);

### Shape operator, \( Sp \)

The shape operator is defined as the rate of change of the unit normal, i.e. curvature, in a particular direction on the surface.  It is an extrinsic curvature.

The coefficients can be found by taking the partial derivatives of the surface normal with respect to the surface parameters:

\( Sp_u = \frac{\partial \mathbf{n}}{\partial u} \)

\( Sp_v = \frac{\partial \mathbf{n}}{\partial v} \)

However, the shape operator matrix can also be found as the matrix product of the metric, \(g\), with the second fundamental form, \(b\).
We use this method as it will be useful immediately below.

Sp:factor(ratsimp(expand(ginv.b)));

### Principal curvatures of the surface

The principal curvatures can be found from the eigenvalues of the shape operator matrix and the principal directions from the eigenvectors.
Using Maxima we detrmine both in one fell swoop:

[spvals,spvects]:eigenvectors(Sp)$

In all their horrendous glory, here are the eigenvalues, simplified by Maxima CAS: 

k1:factor(radcan(ratsimp(spvals)))[1][1];
k2:factor(radcan(ratsimp(spvals)))[1][2];

Now comes plot of the principal curvature (PC) values along the embedded curve.  Note that the PC values are equal at the origin (t=0).  This point is thus an umbilical point that is identical to a flat plane in terms of curvature.

The PC values are defined anywhere on the surface but here they are plotted only along the embedded curve.

k1t:ev(k1, u=tu, v=tv)$
k2t:ev(k2, u=tu, v=tv)$
wxdraw2d(line_width=1, background_color=light_yellow,
    xlabel="Parameter, t", ylabel="Principal Curvature: pc1 (green), pc2 (blue)", title="Principal Curvature Values",
    proportional_axes=none, grid=true,
    color=dark_green, explicit(k1t, t, -2, 2),
    color=blue, explicit(k2t, t, -2, 2)
)$

In all their formidable glory, here are the eigenvectors, which contain the components of the principal directions with respect to the tangent basis, i.e.:

\(spvecs[1][1] * fu + spvecs[1][2] * fv\)

\(spvecs[2][1] * fu + spvecs[2][2] * fv\)

spvec1:radcan(spvects[1][1]);
spvec2:radcan(spvects[2][1]);

One can see that without CAS (Maxima) this would not be so easily accomplished.

Now we check if the principle directions are orthogonal, as they should be.  This is done by taking their dot product:

tmp2:(spvec1[1]*fu+spvec1[2]*fv).(spvec2[1]*fu+spvec2[2]*fv);

This mess needs to be simplified:

ratsimp(tmp2);

They are orthogonal, and, once again, the power of CAS (Mavima) in making a horrendous calculation easy is demonstrated.

Since the length of the eigenvectors have no particular geometric significance, we need to first normailze them (make them unit vectors).

tmp3:fullratsimp(spvec1[1]*fu+spvec1[2]*fv)$
pcv1:fullratsimp(tmp3/fullratsimp(sqrt(tmp3.tmp3)));

tmp3:fullratsimp(spvec2[1]*fu+spvec2[2]*fv)$
pcv2:fullratsimp(tmp3/fullratsimp(sqrt(tmp3.tmp3)));

Check if they are both unit vectors and also check if they lie within the tangent plane.

fullratsimp(pcv1.pcv1);
fullratsimp(pcv2.pcv2);
fullratsimp(pcv1.n);
fullratsimp(pcv2.n);

In order to plot, we first express both vectors along the curve, c, in terms of the parameter t.

pcv1t:fullratsimp(ev(pcv1, u=tu, v=tv));
pcv2t:fullratsimp(ev(pcv2, u=tu, v=tv));

All is good, so now the plot.  Here, the PC unit vectors are scaled with the respective PC values.

The plot is not animated but different points along the curve can be plotted by varying the "s1" variable in the plot command.

s1:-0.8$
wxdraw3d(color=gray30, line_width=1,
parametric_surface(f[1], f[2], f[3], u,-1,1,v,-1,1), line_width=2, color=red, parametric(c[1], c[2], c[3], t,-1,1),
line_width=2, color=forest_green, vector(ev(c, t:s1), ev(k1t*pcv1t, t=s1)),
color=blue, vector(ev(c, t:s1), ev(k2t*pcv2t, t=s1))
)$

### The Christoffel Symbols of the Monkey Saddle Surface

In Maxima, the Christoffel symbols are found using the "ctensor" package:

init_ctensor()$

dim:2$ ct_coords:[u,v]$
lg:g$
ug:ginv$

christof(mcs)$

A peculiarity of Maxima is that the indexing of the Christoffel symbols differs from the conventions used in most textbooks.  We now rectify this by defining new varialbles.

The pattern is \( \Gamma_{jk}^i = mcs(j, k, i) \)

Γ111:mcs[1,1,1];
Γ112:mcs[1,2,1];
Γ121:mcs[2,1,1];
Γ122:mcs[2,2,1];
Γ211:mcs[1,1,2];
Γ212:mcs[1,2,2];
Γ221:mcs[2,1,2];
Γ222:mcs[2,2,2];

### Geodesics of the Monkey Saddle Surface

The geodesic curvature of the parabolic curve that is embedded in the Monkey Saddle surface is not everywhere equal to 0 and hence the curve is not a geodesic. But what are the geodesics?

The trivial geodesics are the initial coordinate curves: [0, t, 0], [t,0,t^3] 

gd1:ev(f, u=0, v=t);
gd2:ev(f, u=t, v=0);

These trivial geodesics are now plotted at a specific point:

wxdraw3d(color=gray30, line_width=1, parametric_surface(f[1], f[2], f[3], u,-1,1,v,-1,1), line_width=2, color=red, parametric(gd1[1], gd1[2], gd1[3], t,-1,1),
line_width=2, color=blue, parametric(gd2[1], gd2[2], gd2[3], t,-1,1)
)$

For the general geodesic curve, we need to solve the differential equation system:

\[ u''(t) + \Gamma_{11}^1 u'(t)^2 + 2 \Gamma_{12}^1 u'(t) v'(t) + \Gamma_{22}^1 v'(t)^2 = 0 \]

\[ v''(t) + \Gamma_{11}^2 u'(t)^2 + 2 \Gamma_{12}^2 u'(t) v'(t) + \Gamma_{22}^2 v'(t)^2 = 0 \]

In what follows, the variables \(u\), \(v\) are to be understood as being functions of the parameter \(t\), with subscripts indicating derivatives w/r to \(t\).

depends(u,t)$
depends(v,t)$
diff(u,t,2) + Γ111*diff(u,t)^2 + 2*Γ112*diff(u,t)*diff(v,t) + Γ122*diff(v,t)^2;
diff(v,t,2) + Γ211*diff(u,t)^2 + 2*Γ212*diff(u,t)*diff(v,t) + Γ222*diff(v,t)^2; 

This is a coupled set of second-order, non-linear differential equations that would need to be solved numerically.
The numerical solution, although straightforward and within the capability of Maxima/WxMaxima, would only yield an array of points that would provide no mathematical insight.
However, we will demonstrate finding the numerical solution to the Monkey Saddle geodesics in another, dedicated, web page. 

## Epilogue

Differential geometry is a most beautiful and exciting area of mathematics, yet these inspiring aspects are usually obscured by the dull tedium of the often complex derivations and computations
that are involved. Using a CAS, however, can liberate one from this dull tedium and allow the underlying beauty and excitement to come forth. The free and open source Maxima/WxMaxima package,
as (hopefully) adequately demonstrated in this article, is a perfect tool for exploring differential geometry without getting stalled in the dense and tiresome thicket of manual calculation. 