\documentclass{article} % Your input file must contain these two lines \usepackage{amsfonts} \usepackage{graphics} \usepackage{html} \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} This is in no way intended as a self-contained introduction to Differentiable Manifolds. It is simply a very rough overview, and I accept no responsibility for confusion, injury, or death resulting from the misuse of this product. \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: (Note: f is the vector field on the manifold) \begin{verbatim} for (i=1;i<=iterations;i++) { /* TAKE STEP IN CHARTS: */ switch (anum) { case 0 : EulerStep (n,stepsize,yold,ynew,narr);break; case 1 : RK4Step (n,stepsize,yold,ynew,narr);break; } /* SWITCH CHARTS IF NEED BE: */ TendToCharts (&achart,ynew,narr); /* DATA OUTPUT: MANIFOLD OR CHART COORDS? */ (charts == 1) ? ShowDataVector(yold,n) : ShowDataVector (x,m); /* OUTPUT NORM? */ if (narr >= 3) {printf("Norm: %f\n",EucNorm(m,x));} /* RESET yold FOR NEXT ITER: */ for (j=1;j<=n;j++) {yold[j][1] = ynew[j][1];} } \end{verbatim} The code for an Euler step: \begin{verbatim} int EulerStep ( int n, double stepsize,double ** xold,double ** xnew,int narr) { double ** v = newmatrix(n,1); /* TEMP FOR LOCAL VF */ /* VF EVAL IN CHARTS...NOTE HEAVY USE OF GLOBALS...*/ PullBack(xold,v); /* UPDATE ESTIMATE */ LinearCombo (n,xnew,1,xold,stepsize,v); if (narr >= 2) { printf("\nTaking an Euler Step...\n"); printf("Vector Field Pullback: ");ShowDataVector (v,n); printf("new chart coords: ");ShowDataVector (xnew,n); } } \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} \subsection{Notes on the Code} \subsubsection{Implementation} I've made some serious changes to the code since this all started. While I think it makes the code easier to use and generalize, I fear it has been at the cost of additional abstruseness. Much more is being done with pointers, which I think increases efficiency and generality (makes the code faster, smaller, and easier to extend to other manifolds). Let's take a look at out Inverse Chart, \[ \varphi ^{-1}: \mathbb{R}^2 \rightarrow \mathcal{S}^2 \subset \mathbb{R}^2 \] The old version of calling this function looked like: \begin{verbatim} (*PhiInv) (n,m,y,x); \end{verbatim} Simple! \texttt{PhiInv} is a pointer to a function sending $\mathbb{R}^2 \rightarrow \mathbb{R}^3$. So the dereferencing, \texttt{(*PhiInv)} is that function, and we see it being applied. Now, for the new version, we have: \begin{verbatim} (*((*achart).PhiInv)) (n,m,y,x); \end{verbatim} What happened? That's progress? Well, of a sort. The idea is that \texttt{achart} is a pointer to a structure containing all the chart information: \begin{verbatim} typedef int (*funptr) (int,int,double **,double **); typedef struct chart { funptr Phi; funptr PhiInv; int (*dPhi) (int,int,double **,double **,double **); /* Tangent Map */ int (*Psi) (int,double **); } Chart; Chart *achart; \end{verbatim} This allows us to switch and apply the various functions associated with charts without much overhead: before, all the functions had to be changed; now, we simply point \texttt{achart} at a different chart. So let us dissect the above function call (that is what it is, believe it or not...). \texttt{achart} points to a \texttt{struct chart}, so \texttt{*achart} is that structure. Thus, \texttt{(*achart).PhiInv} is simply the "PhiInv" component of the Chart currently being pointed at by \texttt{achart}. This happens to be a function pointer pointing at the desired map. So we dereference it: \texttt{(*((*achart).PhiInv))} is thus a function. \subsubsection{Error and Accuracy} Some more Vague Thoughts on the Numerical Issues involved in all this. How worthwhile is all this? Runge-Kutta in the ambient space seems to perform just fine without all this extra stuff. A real error analysis is needed, but somewhat complicated by all the composed mappings and Pullbacks. So a few brief thoughts (for now): Performing the algorithms in charts should yield the same order of accuracy as doing everything in the ambient space. As long as our charts are more-or-less ``well-behaved,'' we're really not doing anything new. Aside from the few {\em flops} involved in the chartistry, the algorithms are the same. But we are able to do it a dimension lower, which would seem to reduce roundoff error more fundamentally than a few extra {\em flops} will. But that's treading on thin ice: I shouldn't even say things like that until I've done more analysis. How well behaved are these charts? They're fine...more-or-less. Note that we're always dividing by a reasonably sized number. The only problem area is when the height of a point on our manifold (that is, its z-coordinate in $\mathbb{R}^3$ is nearly equal to one of the other coordinates. This will happen midway between the poles and the equator. Note that this seems to agree with my observation in example 2 below. So in the end is it all worth it? Well, it's been a very good exercise in C programming and Differential Geometry, so of course it has. I also think being able to run an intrinsic model would be very useful for certain types of problem. It would also seem that for less simplistic vector fields, we might find that these intrinsic integrations would outperform extrinsic ones. \section{Resources} Here's a list of all the various files involved in this project, and where to find them. First, this file itself may be obtained as a \htmladdnormallinkfoot{web site}{http://www.asis.com/~scotfree/latexlab/latex2html/sphere}, a \LaTeX\ \htmladdnormallinkfoot{file}{http://www.asis.com/~scotfree/latexlab/latex2html/sphere.tex}, or even a PostScript File (not yet). The numerical code is all written in C. The main code is \htmladdnormallinkfoot{here}{http://www.asis.com/~scotfree/public_html/cgi} \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...)} \subsubsection{Euler Step} 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 \subsubsection{4th order Runge-Kutta} Here we see that the Runge-Kutta is just a dandy algorithm: the following is 1000 iterations at stepsize .1.... \includegraphics{rk4harmony.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 changeling; 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.