document new tensor-product and maximum-degree Lagrange interpolation

operators


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@146 df1f8a13-aa1d-4dd4-9681-27ded5b42416

1 files changed, 110 insertions, 49 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index 53e3cf0..a2b97d9 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -96,10 +96,12 @@
 \setlength{\unitlength}{1mm}	% length unit for latex picture environment
 \def\csmash#1{\hbox to 0em{\hss{#1}\hss}}
 
+\def\defn#1{{\bf #1}}
 \def\thorn#1{{\bf #1}}
 
 \def\Cplusplus{\hbox{C\hspace{-.05em}\raisebox{.4ex}{\tiny\bf ++}}}
 \def\cf{cf.\hbox{}}
+\def\Ie{I.e.\hbox{}}
 \def\ie{i.e.\hbox{}}
 \def\eg{e.g.\hbox{}}
 
@@ -133,7 +135,7 @@ interpolators in the future):
 	This interpolator was written in C (plus Maple to generate
 	the coefficients) by Jonathan Thornburg in winter 2001-2002.
 	It supports the \verb|CCTK_InterpLocalUniform()| API.
-	It provides Lagrange polynomial interpolation in
+	It provides two variants of Lagrange polynomial interpolation in
 	1~dimension (orders~1-6) and~2 and 3~dimensions (orders~1-4),
 	and Hermite polynomial interpolation in
 	1~dimension (orders~2-4) and~2 and 3~dimensions (orders~2 and~3).
@@ -146,7 +148,7 @@ interpolators in the future):
 \label{localinterp/sect-basic-terminology}
 
 Within Cactus, each interpolator has a character-string name;
-this is mapped to a Cactus {\bf interpolator handle} by
+this is mapped to a Cactus \defn{interpolator handle} by
 \verb|CCTK_InterpHandle()|.  For any given interpolator handle,
 there may be a separate interpolator defined for each of the
 interpolation APIs (both the processor-local ones provided by this
@@ -154,14 +156,14 @@ thorn, and the global ones provided by driver-specific interpolator
 thorns such as \thorn{PUGHInterp}).
 
 Terminology for interpolation seems to differ a bit from one author
-to another.  Here we refer to the {\bf point-centered} interpolation
-of a set of {\bf input arrays} (defining data on a uniformly or
-nonuniformly spaced {\bf grid} of {\bf data points}) to a set of
-{\bf interpolation points} (specified by a corresponding set of
-{\bf interpolation coordinates}), with the results being stored
-in a corresponding set of {\bf output arrays}.  Alternatively,
-we may refer to {\bf cell-centered} interpolation, using a grid
-of {\bf data cells} and a set of {\bf interpolation cells}.
+to another.  Here we refer to the \defn{point-centered} interpolation
+of a set of \defn{input arrays} (defining data on a uniformly or
+nonuniformly spaced \defn{grid} of \defn{data points}) to a set of
+\defn{interpolation points} (specified by a corresponding set of
+\defn{interpolation coordinates}), with the results being stored
+in a corresponding set of \defn{output arrays}.  Alternatively,
+we may refer to \defn{cell-centered} interpolation, using a grid
+of \defn{data cells} and a set of \defn{interpolation cells}.
 (This latter terminology is common in hydrodynamics interpolators.)
 
 At present all the interpolators do polynomial interpolation, and
@@ -185,11 +187,13 @@ In the future, we may add other interpolators where the choice of
 molecule is data-dependent (and thus may vary from one input array to
 the next), and/or where the entire input grid is used in each interpolation.
 
-We define the {\bf order} of the interpolation to be the order of
+We define the \defn{order} of the interpolation to be the order of
 the fitted polynomial.  That is, in our terminology linear interpolation
 is order~1, quadratic is order~2, cubic is order~3, etc.  An order~$n$
 interpolator thus has $O\big((\Delta x)^{n+1})\big)$ interpolation errors
 for generic smooth input data.
+Section~\ref{localinterp/sect-multi-dim-interp} explains how the
+interpolating polynomial is defined for multi-dimensional interpolation.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -203,7 +207,13 @@ with $O(\Delta x^n)$ jump discontinuities each time the interpolating
 polynomial changes.  Correspondingly, the interpolation error is
 generically a ``bump function'' which is zero at each grid point
 and rises to a local maximum in each grid cell.  This is the case,
-for example, for Lagrange polynomial interpolation.
+for example, for tensor-product Lagrange polynomial interpolation.
+
+[For maximum-degree Lagrange polynomial interpolation, the
+interpolation error may not be zero at the grid points, and the
+interpolant itself may not be continuous there.  
+Section~\ref{localinterp/sect-multi-dim-interp} explains this
+in detail.]
 
 This thorn also provides Hermite polynomial interpolation, which
 guarantees a smooth ($C^1$) interpolant and (for smooth input data)
@@ -239,7 +249,7 @@ We use \verb|pt|, \verb|in|, and \verb|out| as generic 0-origin integer
 subscripts into the arrays of interpolation points, input arrays, and
 output arrays respectively.  We use \verb|(i,j,k)| as a generic
 \verb|N_dims|-vector of integer subscripts into the input array
-\verb|input_arrays[in]|.  (Thus {\bf {\tt (i,j,k)} space} refers to
+\verb|input_arrays[in]|.  (Thus \defn{{\tt (i,j,k)} space} refers to
 the grid of data points.)  We usually only write array subscripting
 expressions for the 3-D case; the restrictions/generalizations to
 other dimensions should be obvious.
@@ -265,11 +275,11 @@ $x_{\min}$, $x_{\max}$, $y_{\min}$, $y_{\max}$, $z_{\min}$,
 $z_{\max}$, \dots.
 
 When describing Jacobians and domains of dependence, it's useful to
-introduce the notion of {\bf molecule position}, a nominal reference
+introduce the notion of \defn{molecule position}, a nominal reference
 point for the interpolation molecule in \verb|(i,j,k)| coordinate
 space.  (For example, the molecule position might just be the \verb|(i,j,k)|
 coordinates of the molecule's center.)  We also introduce
-{\bf molecule coordinates} \verb|(mi,mj,jk)|, which are just
+\defn{molecule coordinates} \verb|(mi,mj,jk)|, which are just
 \verb|(i,j,k)| coordinates relative to the molecule position.
 We use \verb|m| or as a generic molecule coordinate.  Thus
 (in notation which should be obvious) a generic molecule operation
@@ -295,7 +305,7 @@ to return outputs back to the caller.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{The Uniform Cartesian Interpolator}
-\label{sect-localinterp/uniform-Cartesian-interp}
+\label{localinterp/sect-uniform-Cartesian-interp}
 
 The uniform Cartesian interpolator is the one which used to live
 in the \thorn{PUGHInterp} thorn.  It was written in C by Thomas Radke
@@ -326,7 +336,16 @@ API, and some examples of how to use it.
 Internally, this interpolator always does 3-D interpolation, inserting
 zero coordinates as appropriate for lower dimensionalities.  The
 interpolation is done by successive 1-D interpolations along each
-axis.  See the \verb|README| file in the source code directory
+axis.%%%
+\footnote{%%%
+	 Note that this means that different axes are treated
+	 slightly differently by the interpolator.  In other
+	 words, at the level of finite differencing errors,
+	 interpolation does {\em not\/} commute with permuting
+	 the axes.  However, in practice the differences are
+	 likely to be small, at least for smooth input data.
+	 }%%%
+{}  See the \verb|README| file in the source code directory
 \verb|LocalInterp/src/UniformCartesian/| for further details.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -339,12 +358,22 @@ The generalized polynomial interpolator was written in C
 winter 2001-2002.  It provides:
 \begin{itemize}
 \item	Lagrange polynomial interpolation of orders 1--6 for 1-D arrays,
-	and 1--4 for 2- and 3-D arrays, registered under the operator
-	name \verb|"Lagrange polynomial interpolation"|.%%%
+	and 1--4 for 2- and 3-D arrays.  This interpolator actually provides
+	two slightly different flavors of Lagrange interpolation,
+	\defn{tensor product} and \defn{maximum degree}, registered
+	under the operator names\\
+	\verb|"Lagrange polynomial interpolation (tensor product)"| and\\
+	\verb|"Lagrange polynomial interpolation (maximum degree)"|
+	respectively.  Section~\ref{localinterp/sect-multi-dim-interp}
+	of this thorn guide explains the two variants in detail, but
+	for most purposes the \underline{tensor-product} variant is
+	preferable.  For convenience this is also registered under
+	the operator name\\
+	\verb|"Lagrange polynomial interpolation"|.%%%
 \footnote{%%%
-	 For backwards compatability, the operator name
-	 \hbox{\tt "generalized polynomial interpolation"}
-	 is also allowed.
+	 For backwards compatability it's {\em also\/}
+	 registered under the operator name
+	 \hbox{\tt "generalized polynomial interpolation"}.
 	 }%%%
 \item	Hermite polynomial interpolation of orders 2-4 for 1-D arrays,
 	and 2--3 for 2- and 3-D arrays, registered under the operator
@@ -355,7 +384,7 @@ but at the moment only hypercube-shaped molecules are implemented.
 It would be easy to add additional orders and/or dimensions if desired.
 
 This interpolator supports a number of options specified by a
-{\bf parameter table} (a Cactus key-value table, manipulated by the
+\defn{parameter table} (a Cactus key-value table, manipulated by the
 \verb|Util_Table*()| APIs).  Note that all the interpolator options
 described here apply only to the current interpolator call: there is
 no visible state kept inside the interpolator from one call to another.
@@ -376,6 +405,7 @@ order~2, cubic is order~3, etc.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \subsection{Molecule Size and Centering}
+\label{localinterp/sect-molecule-size+centering}
 
 If no grid boundaries or excised points are nearby, the interpolator
 centers the molecules around the interpolation point as much as possible.
@@ -581,11 +611,11 @@ near boundaries and excised points.
 To describe how these parameters work, it's useful to define two
 regions in the data coordinate space:
 \begin{description}
-\item[The ``valid-data region'']
+\item[The \defn{valid-data region}]
 	is the entire grid minus any ommited-on-boundary
 	or excised points; to make this a region we take the
 	Lego-block bounding box of the non-excised grid points.
-\item[The ``default-centering region'']
+\item[The \defn{default-centering region}]
 	is that region in interpolation-point space where the
 	interpolator can use the default molecule centering
 	(described in table~\ref{localinterp/tab-molecule-size+centering}),
@@ -844,6 +874,57 @@ out-of-range point will be reported.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+\subsection{Multi-Dimensional Interpolation}
+\label{localinterp/sect-multi-dim-interp}
+
+In multiple dimensions, there are two plausible definitions for the
+generic interpolating polynomial of a given order (degree) $n$.  For
+convenience I'll describe these for the 2-D case, but the generalization
+to any number of dimensions should be obvious:
+
+\begin{itemize}
+\item	The generic ``\defn{tensor product}'' polynomial of degree $n$ is
+	\begin{equation}
+	f(x,y) = \sum_{{0 \le i \le n} \atop {0 \le j \le n}} a_{ij} x^i y^j
+					\label{localinterp/eqn-polynomial-2d-TP}
+	\end{equation}
+\item	The generic ``\defn{maximum degree}'' polynomial of degree $n$ is
+	\begin{equation}
+	g(x,y) = \sum_{0 \le i+j \le n} a_{ij} x^i y^j
+					\label{localinterp/eqn-polynomial-2d-MD}
+	\end{equation}
+\end{itemize}
+
+Because it has $(n+1)^2$ coefficients, the tensor-product polynomial $f$
+can (and does) pass through all the $(n+1)^2$ input data points in a
+square molecule of size $n+1$.  This implies that the interpolation
+error vanishes at the input grid points, and that the overall
+interpolating function is continuous (up to floating-point roundoff
+errors).  However, $f$ does have the slightly peculiar property of
+having terms up to $x^n y^n$ despite being of ``degree'' $n$.
+(For example, the ``linear'' interpolant for $n=1$ would have $xy$ terms,
+even though those are formally quadratic in the independent variables
+$x$ and $y$.)
+
+In contrast, the maximum-degree polynomial $g$ has only
+$\frac{1}{2} (n+1)(n+2)$ coefficients, so for generic input data
+(\ie{} input data which isn't actually sampled from a polynomial of
+the form~$(\ref{localinterp/eqn-polynomial-2d-MD})$) $g$ can't pass
+through all the $(n+1)^2$ input data points in a square molecule of
+size $n+1$.  The interpolator actually does a least-squares fit of
+the polymomial $g$ to the input data in the molecule, so in general
+the molecule won't pass through {\em any\/} of the data points!
+Moreover, each time the interpolation point crosses a grid square
+(for odd~$n$) or the center lines of a grid square (for even~$n$),
+the set of points used in the molecule changes 
+(\cf{}~section~\ref{localinterp/sect-molecule-size+centering}),
+so the interpolant generally has a jump discontinuity!
+For these reasons, the tensor-product
+choice~$(\ref{localinterp/eqn-polynomial-2d-TP})$ is generally
+preferable.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 \subsection{Molecule Family}
 \label{localinterp/sect-molecule-family}
 
@@ -1253,7 +1334,7 @@ interp_molecule_positions[0][13] =  9    /*   1.00     [0.8, 0.9, 1.0] */
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\subsection{Smoothing and Multi-dimensional Interpolation}
+\subsection{Smoothing}
 
 The way the generalized polynomial interpolator is implemented it's
 easy to also do Savitzky-Golay smoothing.%%%
@@ -1283,34 +1364,14 @@ of the usual 4-point one for cubic interpolation).  2 would mean to
 enlarge by 2 points, (\eg{} to use a 6-point molecule for cubic
 interpolation).  Etc etc.
 
-Note that in $>1$~dimension, the default hypercube-shaped molecules
-already use more data points than the number of free parameters in
-the interpolating polynomials, \ie{} they already do some Savitzky-Golay
+Note that in $>1$~dimension, the maximum-degree Lagrange interpolation
+already uses more data points than the number of free parameters in
+the interpolating polynomials, \ie{} it already does some Savitzky-Golay
 smoothing.  For example, in 2-D a generic cubic polynomial
 $f(x,y) = \sum_{i+j \le 3} c_{ij} x^i y^j$ has 10 free parameters
 $c_{ij}$, which we least-squares fit to the 16 data points in the
 $4 \times 4$ molecule.
 
-A somewhat suprising consequence of this is that in $>1$~dimension,
-it's {\em not\/} necessarily the case that the interpolant passes through
-the input data points, \ie{} if an interpolation point coincides with
-one of the input grid points, in general the output value (the interpolation
-result) will {\em not\/} necessarily match the input value.  Of course,
-if the input data are actually a polynomial of degree $\le$ the interpolation
-order, then the interpolant will match them exactly (up to floating-point
-roundoff errors).
-
-[An alternative design choice would be to define a generic cubic
-polynomial in 2-D as $g(x,y) = \sum_{i \le 3, j \le 3} c_{ij} x^i y^j$,
-so there would be the same number (16) of free parameters as data points,
-and thus the interpolant would be guaranteed to pass through the input
-data points (up to floating-point roundoff errors) for any choice of
-the input data.  (This is similar to how the older uniform Cartesian
-interpolator described in
-section~\ref{sect-localinterp/uniform-Cartesian-interp} works.)
-No such molecule family is implemented yet for the generalized polynomial
-interpolator, but this might get done in the future if someone wants it\dots]
-
 Savitzky-Golay smoothing is basically free apart from the increase in
 the molecule size, e.g. a \verb|smoothing|=2 cubic interpolation has
 exactly the same cost as any other 6-point--molecule interpolation.