\documentclass[mathserif,9pt]{beamer}
\parindent=0pt \parskip=8pt

\beamertemplatenavigationsymbolsempty
\beamertemplateballtoc

\usepackage{beamerthemeCambridgeUS} 
\usepackage{pgf,pgfarrows,pgfnodes,pgfautomata,pgfheaps,pgfshade}
\usepackage{amsmath,amssymb}
\usepackage[mathscr]{eucal}
\usepackage[latin1]{inputenc}
%\usepackage{colortbl}
\usepackage[english]{babel}
\setlength{\parindent}{0cm}
\usepackage{psfrag}

%\usepackage{times}

\input{slidesyms}
\input{syms}

%% ================================================================
\title{Numerical differential geometry in Matlab}
\author[J.\ Kerl]{John Kerl}
\institute[Arizona]{Department of Mathematics, University of Arizona \\
	\vspace*{0.2in}
	Graduate Student Colloquium}
\date{January 16, 2008}

\begin{document}

\begin{frame}
\titlepage
\end{frame}


%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Topics}

\begin{itemize}
\item Why bother?
\vspace*{0.1in}
\item 1D and 2D meshes and numerical differentiation
\vspace*{0.1in}
\item Coordinates and parameterizations
\vspace*{0.1in}
\item Surfaces embedded in $\R^3$
\vspace*{0.1in}
\item Numerical estimation of tangent and unit-normal vectors
\vspace*{0.1in}
\item Metric coefficients
\vspace*{0.1in}
\item Christoffel symbols, curvature tensors, scalar curvature \ldots
\end{itemize}

I care not so much about how far we get today (I can finish in a part-two talk,
or you can read the rest of the Matlab code yourself), and more about
identifying the essential concepts you need to translate between pure math and
numerical methods.  \end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Why?}

\begin{center}
\emph{
The soul cannot think without a picture.
}
\\--- Aristotle\index{Aristotle} (384-322 B.C.).
\end{center}

Level-curve problem (Jan. 2001 geometry qual problem \#3):

\emph{
Let $g:\mathbb{S}^2 \to \R: (x,y,z) \mapsto y^2 - z$.  Determine the critical
points, the critical values, and qualitatively describe the level sets.
}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Solution}

Use the method of Lagrange multipliers with
constraint curve $f(x,y,z)=x^2+y^2+z^2=1$.
There are four critical points:
$$
	\colvecthree{0}{0}{1},
	\colvecthree{0}{0}{-1},
	\colvecthree{0}{\sqrt{3}/2}{-1/2}, \qand
	\colvecthree{0}{-\sqrt{3}/2}{-1/2},
$$
The critical values are found by $y^2-z$.  We have
$$
	-1,
	1,
	5/4, \qand
	5/4,
$$
respectively.  Level sets take a bit more work.
(See 
\texttt{http://math.arizona.edu/{\textasciitilde}kerl/doc/prolrevqual.pdf}.)

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}

Here's a pair of Matlab plots of $g$ on the sphere.  (See \texttt{y2z.m} in the
directory where you found this file.) The first plot is a top view (positive
$z$ is up); the second plot is a bottom view (negative $z$ is up).  The level
curves are simply isotherms; you can read off their topologies at a glance:
point, loop, figure eight, pair of loops, pair of points.  This doesn't replace
the explicit computations \ldots but it certainly adds another perspective
which I want to be able to see.

\begin{center}
\psfragscanon
\includegraphics[scale=0.5]{figures/y2z.eps}
\includegraphics[scale=0.5]{figures/y2zb.eps}
\end{center}
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{1D mesh}

\begin{verbatim}
dx = 0.1; x = [0: dx: 10];
y = sin(x);
y = sin(x)  + cos(x); plot(x,y)
y = sin(x) .* cos(x); plot(x,y)
\end{verbatim}

%Key points:  Multiple input and output points are contained in \emph{vectors}.
%(This is not the way we use vectors in linear algebra.) Matlab supports vector
%operations, e.g. \texttt{y=sin(x)}.

Key points:
\begin{itemize}
\item Multiple input and output points are contained in
\emph{vectors}.  (This is not the way we use vectors in linear algebra.)

\item Matlab supports vector operations, e.g.
\texttt{y=sin(x)}.

\item Arrays are added componentwise, but the \texttt{*}
symbol means matrix multiplication (which isn't even defined if the dimensions
aren't compatible).  Componentwise multiplication, exponentiation, etc.
are done using \texttt{.*}, \texttt{.{\textasciicircum}}, etc.
\end{itemize}
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{1D numerical differentiation}

\begin{verbatim}
y = sin(x);
plot(x, y, '-o', 'LineWidth', 2);

yp = gradient(y, dx);
hold on
plot(x, yp, '-o', 'LineWidth', 1);
\end{verbatim}

\begin{center}
\psfragscanon
\includegraphics[scale=0.7]{figures/yyp.eps}
\end{center}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{1D numerical differentiation (continued)}

How would you numerically estimate a derivative?
\begin{itemize}
\item Forward difference:   $\frac{f(x+h) - f(x)  }{h}$
\item Backward difference:  $\frac{f(x)   - f(x-h)}{h}$
\item Centered difference:  $\frac{f(x+h) - f(x-h)}{2h}$
\end{itemize}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{1D numerical differentiation (continued)}

Matlab does:
\begin{itemize}
\item forward  differences on the left edge,
\item backward differences on the right edge,
\item centered differences in the middle.
\end{itemize}

The numerical estimates do not exactly match the analytical results (indeed,
they match rather poorly here since $h=1$):
\begin{verbatim}
>> x=[1:5]
x   =  1     2     3     4     5
>> y = x.^2
y   =  1     4     9    16    25
>> gradient(y, 1)
ans =  3     4     6     8     9
\end{verbatim}
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{1D numerical differentiation (continued)}

What we learn (or teach) in calculus class: accuracy of numerical estimates
increases with smaller mesh size ($h$ or \texttt{dx}) \emph{if} arbitrary
precision is available.

What we don't learn until numerical analysis \emph{if} we ever take/teach it at
all (!):  subtraction of nearly equal numbers --- e.g.  $f(x+h)-f(x)$ for
continuous $f$ --- results in loss of significant digits.  You \emph{cannot}
just try to ``take the limit'' as $h$ goes to zero.  You have to stop sooner.

There is a trade-off:  smaller $h$ increases ideal accuracy but worsens
cancellation error.

Example: base-10 arithmetic with 9 significant digits (cf. base-2 with 53
significant digits for double precision).  Inputs have 9 significant digits;
output has only one:
$$
	1.23456789 -
	1.23456788 =
	0.00000001
$$
Numerical analysis is a non-trivial subject, but a rule of thumb (not important
for this talk since I'm using very coarse meshes) is to keep $h$ greater than
square root of machine epsilon ($10^{-9}$ for double precision).
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{2D meshgrid}

\begin{verbatim}
>> [x,y]=meshgrid(1:4,5:9)
x =  1     2     3     4  |  y = 5     5     5     5
     1     2     3     4  |      6     6     6     6
     1     2     3     4  |      7     7     7     7
     1     2     3     4  |      8     8     8     8
     1     2     3     4  |      9     9     9     9
                          |
>> z = x + y              | >> z = x.^2 .* y
z =  6     7     8     9  | z =  5    20    45    80
     7     8     9    10  |      6    24    54    96
     8     9    10    11  |      7    28    63   112
     9    10    11    12  |      8    32    72   128
    10    11    12    13  |      9    36    81   144
\end{verbatim}
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{Surface plots}

\begin{verbatim}
dx = 0.1; dy = 0.1;
[x,y] = meshgrid(-2:dx:2,-2:dy:2);
z = x.^2 - y.^2;
surf(x,y,z);
\end{verbatim}

\begin{center}
\psfragscanon
\includegraphics[scale=0.8]{figures/zx2y2.eps}
\end{center}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{2D numerical differentiation}

\begin{eqnarray*}
	\frac{\partial f}{\partial x} &\approx& \frac{f(x+h, y) - f(x-h, y)}{2h} \\
	\frac{\partial f}{\partial y} &\approx& \frac{f(x, y+h) - f(x, y-h)}{2h}.
\end{eqnarray*}

\begin{verbatim}
>> [dzdx,dzdy]=gradient(z,dx,dy)
dzdx = 15  20  30  35  |  dzdy =  1   4   9  16
       18  24  36  42  |          1   4   9  16
       21  28  42  49  |          1   4   9  16
       24  32  48  56  |          1   4   9  16
       27  36  54  63  |          1   4   9  16
\end{verbatim}

Again, numerical estimates don't exactly match analytical results.

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Dual roles for vectors and matrices}

Purely, we tend to think of matrices as implementing linear transformations.

Here, a matrix is an $M \times N$ container of sample points --- just as
in the 1D case, a vector served as a container of sample points.

Given $z = f(x,y)$, e.g. $z = x^2 - y^2$, the $x$, $y$, and $z$ are 2D arrays.
Think of this $z$ as a \emph{discretely sampled scalar field} on the plane.

Likewise
$$
	u = f(x,y), \quad v = g(x,y), \quad w = h(x,y)
$$
would be a discretely sampled vector field.

In Matlab, you can loop over array indices, or you can have the looping done
for you.  Furthermore, that implicit looping can be done on the entire array
(e.g. \texttt{z=x+y}), or on slices (using the colon).

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]

Examples:

\begin{verbatim}
>> x
x =   1     2     3     4
      1     2     3     4
      1     2     3     4
      1     2     3     4
      1     2     3     4

>> x(1,1)=400; x(2,:)=77; x(:,3)=999
x = 400     2   999     4
     77    77   999    77
      1     2   999     4
      1     2   999     4
      1     2   999     4
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Coordinates and parameterizations}

\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{figures/coordparam.eps}
\end{center}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Surfaces embedded in $\R^3$}

Parameterization:
\begin{eqnarray*}
	\colvecthree{x}{y}{z} &=& \colvecthree{x(s,t)}{y(s,t)}{z(s,t)}
\end{eqnarray*}
Coordinates:
\begin{itemize}
\item Graph coordinates: $x$ and $y$.
\item Surface of revolution: $\theta$ and $y$.
\item Cylinder:  $r$ and $z$.
\item Spherical/toroidal: $\phi$ and $\theta$.
\end{itemize}
Generically (in \texttt{gdgproj.m}) think of these as $s$ and $t$.

For the hyperboloid of one sheet,
use $s=z$ and $t=\theta$ coordinates.  Find the hyperbola in $r$ and $z$;
form a surface of revolution by angle $\theta$ about the $z$ axis.
\begin{align*}
	x^2 + y^2 - z^2 &= 1 \\
	r^2 - z^2 &= 1 \\
	r &= \sqrt{1 + z^2}
\end{align*}

% One sheet  (this file):     r = sqrt(1 + z.^2)

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{Doing that in Matlab}

\begin{verbatim}
zmin     =  -5.0; zmax     =   5.0;
nz       =  50;   dz       = (zmax-zmin)/nz;

thetamin =   0.0; thetamax =   2.0*pi;
ntheta   =  50;   dtheta   = (thetamax-thetamin)/ntheta;

[paramz,theta] = meshgrid(zmin:dz:zmax, ...
    thetamin:dtheta:thetamax+dtheta);

% Apply the parameterization to obtain R^3 coordinates
r = sqrt(1 + paramz.^2);
x = r .* cos(theta);
y = r .* sin(theta);
z = paramz;
\end{verbatim}
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{Doing that in Matlab (continued)}

\begin{verbatim}
surf(x,y,z)
\end{verbatim}

\begin{center}
\psfragscanon
\includegraphics[scale=0.7]{figures/hyp1surfaceplot.eps}
\end{center}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Numerical estimation of tangent vectors}

\begin{eqnarray*}
	\frac{\Dx}{\Ds}\pickpipe{\colvecthree{x}{y}{z}}
		&\approx& \frac{x(s+h, t) - x(s-h,t)}{2h} \\
	\frac{\Dy}{\Ds}\pickpipe{\colvecthree{x}{y}{z}}
		&\approx& \frac{y(s+h, t) - y(s-h,t)}{2h} \\
	\frac{\Dz}{\Ds}\pickpipe{\colvecthree{x}{y}{z}}
		&\approx& \frac{z(s+h, t) - z(s-h,t)}{2h}
\end{eqnarray*}

So,
\begin{eqnarray*}
	\frac{\D}{\Ds}\pickpipe{\colvecthree{x}{y}{z}} &=&
		\colvecthree {\Dx/\Ds} {\Dy/\Ds} {\Dz/\Ds}
			\pickpipe{\colvecthree{x}{y}{z}}
		\\
	&\approx&
	\frac{1}{2h}
	\left[
		\colvecthree {x(s+h,t)} {y(s+h,t)} {z(s+h,t)} -
		\colvecthree {x(s,t)} {y(s,t)} {z(s,t)}
	\right]
\end{eqnarray*}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{Unit-normal vectors}

\begin{align*}
	\vecn &= \frac{\D}{\Ds} \times \frac{\D}{\Dt} \\
	&\\
	\vecnhat &= \frac{\vecn}{\|\vecn\|}
\end{align*}

Matlab:
\begin{verbatim}
N1    = cross(DDs1, DDt1);
N1hat = N1 / norm(N1);
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{One problem}

I had the $x$, $y$, and $z$ coordinates of a surface as three separate 2D
arrays.  Likewise, the partials $\Dx/\Ds$, $\Dy/\Ds$, $\Dz/\Ds$, $\Dx/\Dt$,
$\Dy/\Dt$, and $\Dz/\Dt$ are six separate 2D arrays.

\begin{center}
\psfragscanon
\includegraphics[scale=0.4]{figures/stack.eps}
\end{center}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{Solution}

Snippet from \texttt{gdgproj.m}:

\begin{verbatim}
[DxDs,DxDt] = gradient(x,ds,dt); % Three scalar fields
[DyDs,DyDt] = gradient(y,ds,dt);
[DzDs,DzDt] = gradient(z,ds,dt);
DDs1 = [0 0 0]; DDt1 = [0 0 0]; % Single vectors

% Allocate space for vector fields.
% Put the s and t indices first so that
% size(DDs(:,:,3) is ns x nt.
DDs  = zeros(ns, nt, 3);
DDt  = zeros(ns, nt, 3);
Nhat = zeros(ns, nt, 3);
% continued on next slide ...
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]

\begin{verbatim}
% ... continued from next slide.
% Here is the one place in this file where I loop over
% s and t meshgrid indices.
for ss = 1:ns
    for tt = 1:nt
        DDs1(1) = DxDs(ss,tt); DDt1(1) = DxDt(ss,tt);
        DDs1(2) = DyDs(ss,tt); DDt1(2) = DyDt(ss,tt);
        DDs1(3) = DzDs(ss,tt); DDt1(3) = DzDt(ss,tt);

        N1 = cross(DDs1, DDt1); N1hat  = N1 / norm(N1);

        % Store the tangent basis and unit normals for later.
        DDs (ss,tt,:) = DDs1; % Now a vector field
        DDt (ss,tt,:) = DDt1;
        Nhat(ss,tt,:) = N1hat;
    end
end
\end{verbatim}

%At this point, run \texttt{gdgproj hyp1}.  Point the audience to
%
%\texttt{http://math.arizona.edu/{\textasciitilde}kerl/gdg/gdgproj.m}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{Quiver plots}

The first three arguments to Matlab's \texttt{quiver3} are the
foot coordinates; the last three are the vector coordinates.

\begin{verbatim}
% Draw the manifold for background
surf(x,y,z,zeros(ns,nt)); shading faceted; hold on;

% Draw the vector field
quiver3(x, y, z, DxDs, DyDs DzDs));
title 'D/Ds';
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\begin{center}
\psfragscanon
\includegraphics[scale=1.0]{figures/hyp1DDs.eps}
\end{center}
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\begin{center}
\psfragscanon
\includegraphics[scale=1.0]{figures/hyp1DDt.eps}
\end{center}
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\begin{center}
\psfragscanon
\includegraphics[scale=1.0]{figures/hyp1Nhat.eps}
\end{center}
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Metric coefficients}

The Riemannian metric is a twice-covariant tensor field
which is symmetric and positive-definite at each point.
(Think ``field of dot products'').  Due to bilinearity,
it is specified by the \emph{metric coefficients}:
\begin{eqnarray*}
	g_{ij} &=& g\left(\frac{\D}{\Ds_i}, \frac{\D}{\Ds_j}\right) \\
	%g_{ij} &=& g(\D/\Ds_i, \D/\Ds_j) \\
\end{eqnarray*}
For surfaces:
\begin{align*}
	g_{11} &= g\left(\frac{\D}{\Ds}, \frac{\D}{\Ds}\right) &
	g_{12} &= g\left(\frac{\D}{\Ds}, \frac{\D}{\Dt}\right) \\
	g_{21} &= g\left(\frac{\D}{\Dt}, \frac{\D}{\Ds}\right) &
	g_{22} &= g\left(\frac{\D}{\Dt}, \frac{\D}{\Dt}\right) \\
\end{align*}
%By symmetry of the metric, $g_{12}=g_{21}$.  However, I find it pleasing to go
%ahead and compute all four, and see the numbers doing what the theory says they
%would do (approximately, at least).

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Dot products}

Remember that we multiply vectors componentwise using \texttt{.*}.
For dot products, we want plain old \texttt{*}, along with transpose.
The dot product of $\vecu$ and $\vecv$, if $\vecu$ and $\vecv$ are
written as row vectors, is the matrix product
$$
	\vecu \cdot \vecv = \vecu \; \vecv^t
$$
e.g.
\begin{align*}
	\rowvecthree{1}{2}{3}
	\cdot
	\rowvecthree{4}{5}{6}
	&=
	\rowvecthree{1}{2}{3}
	\;
	\colvecthree{4}{5}{6}
	\\
	&= 1\cdot 4+2\cdot 5+3\cdot 6
	&= 32.
\end{align*}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{Metric coefficients in Matlab}

\begin{verbatim}
g = zeros(ns,nt,2,2); invg = zeros(ns,nt,2,2);
for ss = 1:ns
  for tt = 1:nt
    % The above DDs1 & DDt1 computations were here.
    a = DDs1 * DDs1'; b = DDs1 * DDt1';
    c = DDt1 * DDs1'; d = DDt1 * DDt1';
    g(ss,tt,1,1) = a; g(ss,tt,1,2) = b;
    g(ss,tt,2,1) = c; g(ss,tt,2,2) = d;
    gdet = a*d - b*c; % Inverse metric coefficients
    invg(ss,tt,1,1) =  d/gdet; invg(ss,tt,1,2) = -b/gdet;
    invg(ss,tt,2,1) = -c/gdet; invg(ss,tt,2,2) =  a/gdet;
  end
end
\end{verbatim}
For the hyperboloid of one sheet, I used $s=z$, $t=\theta$.

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}

\begin{center}
\includegraphics[scale=0.55]{figures/hyp1g11.eps}
\quad
\includegraphics[scale=0.55]{figures/hyp1g12.eps}

\includegraphics[scale=0.55]{figures/hyp1g21.eps}
\quad
\includegraphics[scale=0.55]{figures/hyp1g22.eps}
\end{center}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\begin{center}
\frametitle{Christoffel symbols}
\end{center}

For surfaces, there are eight Christoffel symbols, with $i,j,k=1,2$:
$$
	\Gamma_{ij}^k = \frac{1}{2}\, g^{km} \,
		\left( \frac{\Dg_{jm}}{\Di} + \frac{\Dg_{im}}{\Dj} - \frac{\Dg_{ij}}{\Dm} \right).
$$
Here the Einstein summation convention (which Einstein, tongue-in-cheek, called
his ``one great contribution to mathematics'') is used.  That is, we mean
$$
	\Gamma_{ij}^k = \sum_{m=1}^2 \frac{1}{2}\, g^{km} \,
		\left( \frac{\Dg_{jm}}{\Di} + \frac{\Dg_{im}}{\Dj} - \frac{\Dg_{ij}}{\Dm} \right).
$$

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\begin{center}
\frametitle{In Matlab}
\end{center}

First compute the $\frac{\D g_{jk}}{\Di}$'s; Christoffel symbols on the next slide.

\begin{verbatim}
DD_g = zeros(ns,nt,2,2,2);

for j = 1:2
  for k = 1:2
    [DD_g(:,:,1,j,k), DD_g(:,:,2,j,k)] =
      gradient(g(:,:,j,k),ds,dt);
  end
end
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]

\begin{verbatim}
Gamma = zeros(ns,nt,2,2,2);
for i = 1:2
  for j = 1:2
    for k = 1:2
      for m = 1:2
        Gamma(:,:,i,j,k) = Gamma(:,:,i,j,k) ...
          + 0.5 * invg(:,:,k,m) .* ( ...
          DD_g(:,:,i,j,m) + ...
          DD_g(:,:,j,i,m) - ...
          DD_g(:,:,m,i,j));
      end
    end
  end
end
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}

\begin{center}
\includegraphics[scale=0.5]{figures/hyp1G111.eps}
\quad
\includegraphics[scale=0.5]{figures/hyp1G121.eps}

\includegraphics[scale=0.5]{figures/hyp1G211.eps}
\quad
\includegraphics[scale=0.5]{figures/hyp1G221.eps}
\end{center}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}

\begin{center}
\includegraphics[scale=0.5]{figures/hyp1G112.eps}
\includegraphics[scale=0.5]{figures/hyp1G122.eps}

\includegraphics[scale=0.5]{figures/hyp1G212.eps}
\includegraphics[scale=0.5]{figures/hyp1G222.eps}
\end{center}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}
\frametitle{Curvature tensors}

The 1,3 curvature tensor:
$$
	R_{ijk}^l =
		  \frac{\D}{\Di} \Gamma_{jk}^l + \Gamma_{ik}^m \Gamma_{jm}^l
		- \frac{\D}{\Dj} \Gamma_{ik}^l - \Gamma_{jk}^m \Gamma_{im}^l.
$$
The 0,4 curvature tensor:
$$
	R_{ijkl} = R_{ijk}{}^m g_{lm}.
$$
The 0,2 curvature tensor:
$$
	R_{ij} = g^{kl} R_{kijl}.
$$
Scalar curvature:
$$
	S = g^{ij} R_{ij}
$$
which in turn is
$$
	g^{kl} g^{ij} R_{kijl}.
$$
\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]
\frametitle{In Matlab}

Compute $\frac{\D}{\Di} \Gamma_{jk}^l$'s.
\begin{verbatim}
DD_Gamma = zeros(ns,nt,2,2,2,2);
for j = 1:2
  for k = 1:2
    for l = 1:2
      [DD_Gamma(:,:,1,j,k,l), DD_Gamma(:,:,2,j,k,l)] = ...
        gradient(Gamma(:,:,j,k,l),ds,dt);
    end
  end
end
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]

\begin{verbatim}
R13 = zeros(ns,nt,2,2,2,2);
for i = 1:2
  for j = 1:2
    for k = 1:2
      for l = 1:2
        for m = 1:2
          R13(:,:,i,j,k,l) = R13(:,:,i,j,k,l) ...
            + DD_Gamma(:,:,i,j,k,l) - DD_Gamma(:,:,j,i,k,l) ...
            + Gamma(:,:,i,k,m) .* Gamma(:,:,j,m,l) ...
            - Gamma(:,:,j,k,m) .* Gamma(:,:,i,m,l);
        end
      end
    end
  end
end
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]

Matlab doesn't support index contraction by a built-in function (as far as I
know), but you certainly can do it by summing.
\begin{verbatim}
R04 = zeros(ns,nt,2,2,2,2);
for i = 1:2
  for j = 1:2
    for k = 1:2
      for l = 1:2
        for m = 1:2
          R04(:,:,i,j,k,l) = R04(:,:,i,j,k,l) ...
            + R13(:,:,i,j,k,m) .* g(:,:,l,m);
        end
      end
    end
  end
end
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]

\begin{verbatim}
R02 = zeros(ns,nt,2,2);
for i = 1:2
  for j = 1:2
    for k = 1:2
      for l = 1:2
        R02(:,:,i,j) = R02(:,:,i,j) ...
          + invg(:,:,k,l) .* R04(:,:,k,i,j,l);
      end
    end
  end
end
\end{verbatim}

\end{frame}

%% ----------------------------------------------------------------
\begin{frame}[fragile]

\begin{verbatim}
S = zeros(ns,nt);
for i = 1:2
  for j = 1:2
    S(:,:) = S(:,:) + invg(:,:,i,j) .* R02(:,:,i,j);
  end
end
\end{verbatim}

\begin{center}
\includegraphics[scale=0.7]{figures/hyp1scalarcurvature.eps}
\end{center}

\end{frame}


%% ----------------------------------------------------------------
\begin{frame}
\frametitle{More information}

The level-curve problem (January 2001 \#3):

\texttt{http://math.arizona.edu/{\textasciitilde}kerl/doc/prolrevqual.pdf}

Matlab code:

\texttt{http://math.arizona.edu/{\textasciitilde}kerl/gdg/gdgproj.m}

\texttt{http://math.arizona.edu/{\textasciitilde}kerl/gdg/y2z.m}

This file:

\texttt{http://math.arizona.edu/{\textasciitilde}kerl/gdg/gdgprojnotes.pdf}

\end{frame}
%% ================================================================
\end{document}
