\documentclass{article} % Your input file must contain these two lines \usepackage{amsfonts} \usepackage{graphics} \newcommand{\Mone}{\ensuremath{M_1}} \newcommand{\Mtwo}{\ensuremath{M_2}} \newcommand{\float}[1]{\ensuremath{\mathit{float}(#1)}} \newcommand{\iprod}[1]{\ensuremath{\left\langle #1 \right\rangle}} \begin{document} % plus the \end{document} command at the end. \title{Numerical Manifold Madness} \author{Scot Free Kennedy} \maketitle \section{Drive-By Differential Geometry} \float{\mathbb{R}} \subsection{Manifolds} \begin{center}{\em ``A manifold is a set which looks locally like a vector space.'' }\end{center} That is, there exist a family of ``charts'' which map sets in an open cover on a manifold $M$ into a vector space. In our case, we'll use the vector space $\mathbb{R}^n$ here\footnote{Actually, we'll use \float{\mathbb{R}} here, but I'm not sure where to go with that...}, \[\varphi_i : U_i \subset M \rightarrow \mathbb{R}^n\] though complex and stranger manifolds exist. In order that we be able to pass from the domain of one chart to that of another with no discontinuity, we require that the charts be well behaved on the intersection of these domains. That is, where defined, the composition of one with the inverse of the other should be a diffeomorphism. \subsection{Tangent Spaces} The tangent space at a point $p \in M$ is written as $T_p M$, and the union of all these is the Tangent Bundle $TM$. We may think of $T_p M$ for now as the set of all tangent vectors at a point, like the ``tangent plane'' in Calculus. The problem with this definition is that it's not ``intrinsic''; it requires stuctures ``outside'' the manifold itself. We resolve this by defining a ``tangent vector'' at a point $p \in M$ to be an equivalence class, namely the set of (differentiable) curves sharing first derivatives at $p$. This elegant definition is entirely intrinsic, though perhaps less intuitive. Note that, in general, there is no trivial relationship between the tangent spaces at different points on $M$. \subsection{Vector Fields and Flows} A vector field is a selection of a tangent vector from $T_p M$ at every $p \in M$. \[ X : M \rightarrow T M \] An integral curve or trajectory of a vector field is a curve $c : [a,b] \rightarrow M$ such that \[ c'(\tau) := \frac{d}{dt} c(t) |_{t=\tau} = X(c(\tau)) \; \forall \tau \in [a,b] \] and the ``union'' of all these is a diffeomorphism of $M$ called the Flow: \[ \Phi : [a,b] \times M \rightarrow M\] which we may sometimes call the ``Time $\tau$ map'' and write: \[ \Phi_{\tau} : M \rightarrow M \] \section{Numerical Simulation of Manifold Flows} \subsection{Basic Ideas} \subsubsection{Nontrivial Vector Field Evaluation and Approximation} Since the tangent spaces at different points on $M$ are different (though isomorphic), we cannot blithely evaluate the vector field where we like. In addition, simply ``moving in a particular direction'' is a fairly subtle concept - we can't add points on a manifold the way we can in a vector space. In general we solve the former problem via parallel transport or a connection, and the latter by something like an exponential map, \[ \mathit{exp}_p : [a,b] \times T_p M \rightarrow M \] which allows us to move along geodesics in a direction of our choice. But these are complicated structures, and vary with each manifold. This makes computation a grueling effort. \subsubsection{Portrait of the Chartist as a Young Man} Our solution is to do everything in charts. We'll map everything to a vector space, simulate the induced flow there, and map everything (well, just the next position actually) back to the manifold. We need to show that the induced flow maps back to the original flow on $M$. That is, given a curve $c(t)$ in $M$, it is an integral curve of the flow $\Phi$ iff \[ c'(t) = X(c(t)) \; \forall t\] similarly, a curve $\tilde{c}(t)$ in $\mathbb{R}^n$ is an intgral curve of the induced vector field $\tilde{X}$ iff \[ \tilde{c}'(t) = \tilde{X}( \tilde{c}(t)) = d\phi \cdot X(\varphi^{-1}(\tilde{c}(t))) \] but this is a direct result of the chain rule and the properties of these various mappings. \subsection{Issues} \subsubsection{It's {\em nice} to work in charts} We can always be sure we're starting on the surface, and moving on it, without hassling with the sphere as a set in $\mathbb{R}^n$. We can use all existing algorithms with no change: they still performa standard vector field evaluation. It just involves a composition with a chart, which they don't need to see. \subsubsection{Interplay of Math with Machine} Manifolds provide a nice epistemological framework for numerics. For instance, think of the chart as a mapping to the computer screen. Encourages well generalized program structure. This can make things tricky though: a decent framework for numerics can happily compute nonsense. In addition, the Tangent Bundle formalism lends itself well to programming: a computer program {\em needs} to think of a tangent vector as having ``base point'' and ``principal part''. Start on charts rather than the manifold. Everything works otherwise, but the meaning is unclear if we start with an initial point not on our manifold. This possibility is precluded by giving initial conditions in terms of chart coordinates. Problems with switching charts. Seemingly good behavior by luck when code is wrong. \subsection{Specifics} Here's the relevant maps for the Southern Hemisphere: The chart: \begin{eqnarray*} \varphi_S :\mathcal{S}^2 \subset \mathbb{R}^3 &\rightarrow& \float{\mathbb{R}^2} \subset \mathbb{R}^2\\ x_i &\mapsto& \frac{x_i}{1-x_{n+1}} := y_i \end{eqnarray*} The Tangent map induced by the chart: \[ d\varphi : T\mathcal{S}^2 \subset T\mathbb{R}^3 \rightarrow T\float{\mathbb{R}^2} \subset \mathbb{R}^2 \] is given by \[ d\varphi = \left( \frac{\partial y_i}{\partial x_j} \right) = \left( \begin{array}{ccc} \frac{1}{1-x_3} & 0 & \frac{x_1}{(1-x_3)^2} \\ 0 & \frac{1}{1-x_3} & \frac{x_2}{(1-x_3)^2} \end{array}\right) \] and a Euler Step in charts: \[y_{n+1} = y_n + \underbrace{h d\varphi \cdot X(\varphi^{-1}(y_n))}_{\mathit{exp}_{y_{n}} ?} \] \subsection{Code Snippets} The chart for the northern hemisphere: \begin{verbatim} int PolarProjS (int m,int n,double **x,double **y) { int i,j; if (m != n+1) {return (1);} /* DIMENSION ERROR!! */ for (i=1;i<=n;i++) { y[i][1] = x[i][1]/(1+x[m][1]) ; } return(0); } \end{verbatim} Here we see the main loop, a veritable snakepit of references and referents: (Note: f is the vector field on the manifold) \begin{verbatim} for (i=1;i<=iterations;i++) { /* MAP BACK TO MANIFOLD FROM CHART */ (*PhiInv) (n,m,yold,x); /* EVAULATE VECTOR FIELD */ (*f) (3,3,x,v); /* PUSH BACK TO CHART */ (*dPhi) (3,2,x,v,w); /* TAKE STEP IN CHARTS */ EulerStep (2,stepsize,w,yold,ynew); if (EucNorm(n,ynew) > 1) {/* SWITCH CHARTS */ ...} /* DATA OUTPUT: MANIFOLD OR CHART COORDS? */ (charts == 1) ? ShowDataVector(yold,n) : ShowDataVector (x,m); /* RESET yold FOR NEXT ITER */ for (j=1;j<=n;j++) {yold[j][1] = ynew[j][1];} } \end{verbatim} In this implementation, an Euler step is particularly elegant: \begin{verbatim} int EulerStep ( int n, double stepsize,double ** v,double ** xold,double ** xnew) { /* v is VF eval'd on manifold then pushed down into charts */ LinearCombo (n,xnew,1,xold,stepsize,v); } \end{verbatim} A sample vector field or two: \begin{verbatim} int HarmonicVertical (int m,int n,double ** x,double ** y) { int i; if (m != n) {return (1);} /* DIMENSION ERROR!! */ y[1][1] = -1* x[3][1]; y[3][1] = x[1][1]; y[2][1] = 0; return (0); } int NeatoFlow (int m,int n,double ** x,double ** y) { float theta,phi; int i; if (m != 3) {return (1);} /* DIMENSION ERROR!! */ theta = atan(x[2][1]/x[1][1]); phi = acos(x[3][1]); y[1][1] = -pow(sin(phi),2)*sin(theta)-pow(cos(phi),3)*cos(theta); y[3][1] = pow(sin(phi),2)*cos(theta)-pow(cos(phi),3)*sin(theta); y[2][1] = sin(phi)*pow(cos(phi),2); return (0); } \end{verbatim} \pagebreak \section{Examples} \subsection{Naivete Foresworn} This is an illustration of what we'd like to fix by simulating in charts. As should be clear from the given system of ODE's, solutions {\em should} remain on the sphere, which they clearly do not. Note that restricting the system to the sphere via projection will have other problems, such as a loss of ``z-constancy''. \begin{eqnarray*} \dot{x} &=& - y \\ \dot{y} &=& x \\ \dot{z} &=& 0 \end{eqnarray*} \includegraphics{offsphere.eps} \pagebreak \subsection{Longitudinal (or Latitude or some such...)} Here we see that while we still have error, it is error {\em on} our manifold. Note that the error ``pushes'' orbits toward the equator. \begin{eqnarray*} \dot{x} &=& - y \\ \dot{y} &=& x \\ \dot{z} &=& 0 \end{eqnarray*} \includegraphics{hsphere.eps} \pagebreak \subsection{Latitudinal (the other one)} This is a trophy: note that we're crossing charts. This was a milestone... Note also that the ``spreading'' due to error (these should be closed orbits) is nonconstant. This seems an artifact of our charts: the error may be magnified differently by different slopes where the projection intersects the sphere. \begin{eqnarray*} \dot{x} &=& - z \\ \dot{y} &=& 0 \\ \dot{z} &=& x \end{eqnarray*} \includegraphics{vsphere.eps} \pagebreak \subsection{Mystery Flow} This was intended to be a nicely behaved vector field, but apparently some fairies got ahold of it while I wasn't looking. It's a changling; I accept no behavior for its bizarre behavior. \begin{eqnarray*} \dot{x} &=& \sin^2\phi \sin \theta - \cos ^3 \phi \cos \theta \\ \dot{y} &=& \sin^2 \phi \cos \theta - \cos^3 \phi sin \theta \ \\ \dot{z} &=& \sin \phi cos^2 \theta \end{eqnarray*} where \begin{eqnarray*} \theta &=& \arctan (y/x)\\ \phi &=& \arccos(z) \end{eqnarray*} \includegraphics{nsphere.eps} \end{document} % The input file ends with this command.