From 3fb87bb8f664dc60e86894cf140ed0529ebf0060 Mon Sep 17 00:00:00 2001 From: jthorn Date: Tue, 15 Apr 2003 17:27:45 +0000 Subject: 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 --- doc/documentation.tex | 159 ++++++++++++++++++++++++++++++++++---------------- 1 file 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. -- cgit v1.2.3