According to Pete Seeger, Frank Proffitt, when asked about the
Earl Scruggs style of banjo picking, responded with something like:
"I'd really like to be able to play like that, and then not do it."
The sentiment here is more or less the same, but the subject
is cubic interpolation. (Incidentally,
if I could play like Scruggs, or Proffitt or Seeger,
I surely would - now and then.)
By way of explanation of the context: This is part of a rough
model of global atmospheric circulations (under construction). The
computational grid, composed of equilateral triangles on the octahedron,
corresponds to spherical triangles with more or less equal separations
of vertices (example).
E.g. with 32,500 vertices, the lengths of the sides are
near 111 km (1 degree). Regardless of the instantaneous resolution of
the physical grid, algorithms for modeling motions of cold or warm
fronts often predict an unphysical smearing of the fronts. In Lagrangian
algorithms, where interpolation of data from an earlier time is used, the
smearing can be controlled and even reversed by compensating for the
numerical dissipation that is induced by linear interpolation.
So the task at hand is to find a reliable way to compensate for the systematic error of linear interpolation in the equilateral triangles of a hexagonal grid. The emphasis is not so much on numerical accuracy, but rather on not getting it qualitatively wrong in a systematic way.
A one-dimensional version of this is also needed for another purpose, and it is in Correcting linear interpolation. In both cases the gambit is to formulate cubic interpolation, squeeze it as hard as possible, and then not do it. Here the symmetries of the equilateral triangle will be used along with the homogeneous coordinates, (e,w,n) s.t. e+w+n=1.
The orientations of the coordinate lines are indicated by:
n=2 ---nw------ne---
/ \ / \
\ / e=-1 \ / w=-1
----n=1---- / \ / \
/ \ n=1 ---wn------n-------en---
/ \ / \ / \ / \
e=0 w=0 / \ e=0 w=0 / \
\ / \ / / \ / \ / \
\ / \ / n=0 ---ww------w-------e-------ee---
--w=1---n=0---e=1-- \ / \ / \ /
/ \ / \ w=2 / w=1 e=1 \ e=2
\ / \ / \ / \ /
n=-1 ---ws------s-------es---
The names that
have been assigned to the vertices are meant to suggest directions,
but they are not quite points of the compass. The
pairs,'ww', 'wn' and 'nw', are meant
to suggest west of west, west of north and north of west, for
example. For a function f(e,w,n), the values,
f(e)=f(1,0,0),
f(ww)=f(-1,2,0), etc., are to be used to
construct a cubic
polynomial for interpolations in and on the boundary of
the interior triangle. Note: there are twelve independent values of f( ) to be
reckoned with, and there are twelve different elements in a basis of
interpolation functions.
We begin with the linear interpolation,
L[f] = ef(e) + wf(w) + nf(n),and we seek to correct that with nine more polynomial approximations that are zero at the points 'e', 'w' and 'n'. Some polynomials with
P(e)=P(w)=P(n)=0are:
ew eww ewn
0 0 0 0 0 0
\ / \ / \ /
-1 0 -1 -1 0 1 -1 0 -1
/ \ / \ / \
-2 0 0 -2 -4 0 0 2 0---0---0---0
/ \ / \ / \
0 1 0 0 1 0 0 -1 0
ewn(1-n) ew+wn+ne ewn+ew+wn+ne
0 0 -2 -2 -2 -2
\ /
0---0---0 -1 0 -1 -2 0 -2
/ \
0---0---0---0 -2 0 0 -2 -2 0 0 -2
/ \
0 2 0 -2 -1 -2 -2 -2 -2
From these and their reflections and rotations
we are to construct a basis
of nine independent polynomial functions, and the basis must include
three independent quadratics and four independent cubics.
The easy choice, not among the required ones, is for the quartic(s)
Q(e,w,n) = ewn(1-n)/2
It follows easily from the diagrams that
Q(e,w,n)+Q(w,n,e)+Q(n,e,w)+ewn = 0
Thus the three Q's are not independent, nontrivial quartics,
but they are independent among quartics plus cubics. (The quartic
terms that are not directly included are the permutations of
eee(e, w or n) and ee(ww or nn).)
By inspection, the three Q's are Lagrange polynomials - zero at eleven of the twelve points on the diagrams. They are necessarily quartic, but no new quartics will be needed for the construction of the other nine Lagrange polynomials. On the line from 'ww' to 'ee', e+w=1, and the function
e(1-e)(2-e) = ew(1+w) = (1-w)w(1+w),
is zero at 'w', 'e', and 'ee'. Its inward slope at 'w' (the
e-derivative) is 2 and its value at 'ww' is -6. The steps
in the construction of its extension to the entire diagram are in the
diagrams,
ewn(1-e)/2 ew(1+w) ew(1+w+n(n-e))/6
0 0 0 0 0 0
\ / \ / \ /
1 0 0 -2 0 0 0 0 0
/ \ / / \ / \
0---0---0---0 -6 0 0 0 -1 0 0 0
/ / \ / \ / \
0 0 0 0 2 0 0 0 0
Reflections and rotations of the first and third diagrams give the nine
exterior Lagrange polynomials, and the interior ones are
the rotations of nC(e,w,n), where C is the invariant cubic polynomial,
0 0
0 1 0
C(e,w,n)=1+(ew+wn+en+ewn)/2
0 1 1 0
0 0 0
Now let us drop the references to physical directions altogether but keep the coordinate names 'e', 'w' and 'n'. Then there are just the four basis functions,
Cs(e,w,n) = ewn(1-n)/2,
Cww(e,w,n) = ew(1+w+n(n-e))/6,
Cee(e,w,n) = ew(1+e+n(n-w))/6,
Cn(e,w,n) = n(1+(ew+wn+en+ewn)/2),
and they are used in the six orientations indicated below:
0 0 + 0 0 -
\ / / \
- n + + e + + w -
/ \ / \
- w e + 0---n---w---0 0---e---n---0
/ \ / \
0 + 0 0 - - + + 0
0 + 0 - - 0 0 + +
\ / / \
+ e w - 0---w---n---0 0---n---e---0
\ / / \
+ n - + e + - w +
/ \ / \
0 0 0 + - 0
Given are the coordinates, e, w & n, and the values of f at the
marked positions in the diagrams above. In numerical algorithms the
locations of data for the several orientations are gotten by permutations
of pointers.
The number of functions to be approximated in the interior triangle will be denoted as K. Even for relatively crude atmospheric circulation models K is of the order, ten or more. In applications where there are several functions to be interpolated it is the coefficient of K in the counting of floating point operations that matters more in the construction of more efficient algorithms. In fact, the coefficient of K in cubic/quartic algorithms is minimized by the use of the twelve Lagrange polynomials. The number of flops involved is then S+23K, where S is the number of flops required to evaluate the twelve C's. The K interpolations are all instances of
C[f] = Cn*f(n)+Cs*f(s)-Cww*f(ww)-Cee*f(ee); (permute pointers)
C[f] = C[f]+Cn*f(n)+Cs*f(s)-Cww*f(ww)-Cee*f(ee); (permute pointers)
C[f] = C[f]+Cn*f(n)+Cs*f(s)-Cww*f(ww)-Cee*f(ee); (permute pointers)
The twelve-dimensional algorithm, with a basis composed of twelve linearly
independent polynomials, also has twelve independent parameters, and
the coefficient of K in the operation count cannot be further reduced.
Symmetric algorithms with less than twelve dimensions can be constructed by discarding some of the parameters, or by replacement of sets of symmetrically disposed ones by their averages. To be discussed are: a ten-dimensional, twelve-parameter algorithm in which the rotations of f(s) are replaced by (f(s)+f(en)+f(wn))/3, a nine-dimensional, nine-parameter one in which the rotations of f(s) are not used, and a six-dimensional, six-parameter one in which only the rotations of f(s) are used. Linear interpolation, it may be noted, is a three-dimensional, three parameter algorithm, and there are a number of four-dimensional algorithms that may also be too crude for any serious consideration.
Before going on to lower order approximations that might be good enough for the task at hand, the first level of not doing cubic interpolation follows from the observation that the polynomials like eww are only used within the interior triangle. On the line, n=0, the element eww can also be regarded as that cubic element whose slope is 1 at 'w' and zero at 'e', and eww can be used in the construction of a basis for Hermite interpolation. The nine slopes associated with the elements are directed along lines, as indicated by:
n=1-e n=1-w
w=0 e=0
\ /
\ /
n
/ \
|
/ \
|
n=(1-e)/2 / \ n=(1-w)/2
w=(1-e)/2 - | - e=(1-w)/2
/ - - \
-
/ - - \
- | -
/ - - \
n=0 ------ w - - - - - m - - - - - e ------ n=0
w=1-e / e=(1-n)/2 \ e=1-w
w=(1-n)/2
/ \
e=0 w=0
w=1-n e=1-n
From these come the slopes (directed inward) at the edges,
eww eww(1-2n)
o o o o
o o
o o o o
o + o o - o
1/4 o + + o 0 o o o o o o o o
o + + + o o + + + o
o + + + + o o + + + + o
1 o + + + + + o 0 1 o + + + + + o 0
o -3/8 o o -5/8 o
Shown above by o o o are straight lines where P(e,w,n) is constant (=0).
At least one directional derivative is zero on such lines, and all
directional derivatives are zero at intersections of two noncoincident ones.
Next, it follows without much effort that polynomials that are zero
at 'e', 'w', and 'n' and have only one nonzero slope at a midpoint or corner
are:
4ewn(1-2n) ew(w+n(1-6n+4e)/2)=ew(1-2n)(w+5n/2)
o o o o
o o
o o o o
o - o o - o
o o o o o o o o o o o o o o o o
o + + + o o + + + o
o + + + + o o + + + + o
o o o o o o o o o 1 o + + + + + o 0
o 1 o o 0 o
If the six slopes at corners and the three at midpoints were prescribed,
the rotations of these two and ew(1-2n)(e+5n/2) would define nine
elements of a basis for corrections of linear interpolation
by Hermite interpolation on a hexagonal grid.
On the line segment, 'ww' to 'e' (not to 'ee'), the values of corrections of the linear approximation of f(e,w,n) are f(ww)+f(e)-2f(w) at 'ww' and 0 at 'w' and 'e'. For approximate Hermite interpolations, estimates of the associated slopes at 'w' and 'e' are
S(ww)=α*(2f(w)-f(e)-f(ww)).
S(ee)=α*(2f(e)-f(w)-f(ee)).
The values of corrections on the line, 's' to 'n', are f(s)+f(n)-f(e)-f(w) at 's' and
0 at 'm' and 'n'. The estimates of the associated slope at 'm' are
S(s)=β*(f(e)+f(w)-f(n)-f(s)).
Notable values of α and β are: 0 for no correction, 1/2 for central
difference approximations of the slopes, and
1 (say) for estimates that overcompensate for numerical dissipation.
When α=β=1/2, the twelve-dimensional algorithm of this kind is
qualitatively simalar to the one based upon Lagrange polynomials, and its operation
count is S+23K. There is little to be gained by by further consideration of it.
Algorithms with nine degrees of freedom are considerably simpler than those with ten or twelve. There is a twelve-parameter one in which the pairs of function values at the corners (e.g. f(nw) and f(ne)) are replaced by their averages, and there is this nine-parameter one in which values of f(s) f(en) and f(wn) are not used. Then the parameter β does not appear in the approximate Hermite interpolations, and the function, eww, already has the slope one in just one of the six ways slopes are to be specified. The reflections and rotations of eww are not satisfactory elements of a basis, however, because the rotations of ew cannot be constructed from them (ew(e+w)=ew(1-n)!). To define a cubic scheme that interpolates quadratic functions exactly, it suffices to introduce the reflection and rotations of ew(w+n/2). Even though the symmetric term, ewn, appears in it, the approximate, almost-cubic algorithm has only three independent cubic polynomials to accompany the three independent quadratic polynomials in its basis. The algorithm consists of a sum of three terms (orientations) of the form
nf(n)+αew(1-2n)(w+5n/2)(2f(w)-f(e)-f(ww))
+αew(1-2n)(e+5n/2)(2f(e)-f(w)-f(ee))
The reflections, (n,e) -> (e,n) and (w,n) -> (n,w) in the coefficients of
f(e) and f(w), serve to collect the coefficient of f(n). Then the sum is of three
terms of the form
nf(n)+αn(e(2n-e+w/2)+w(2n-w+e/2))f(n)
-αew((w+n/2)f(ww)+(e+n/2)f(ee)).
By inspection, this contains the terms proportional to ewn that
are required for exact interpolation of quadratic functions (when α=1/2).
The only cubic function for which it is qualitatively wrong is f=ewn.
The total number of floating point operations is S+17K.
The term, ewn, reappears explicitly in the ten-dimensional approximation where the averages of terms that contain the rotations of ewnn all contribute ewn/3 to the average of corrections that come from the slopes at midpoints. The ten-dimensional algorithm can be derived from the elements
eww ewn ew(w+n/6)
o o o o o o
o o o
o o o o o o
o + o o + o o + o
1/4 o + + o 0 1/4 o + + o 1/4 7/24 o + + o 1/24
o + + + o o + + + o o + + + o
o + + + + o o + + + + o o + + + + o
1 o + + + + + o 0 o o o o o o o o o 1 o + + + + + o 0
o -3/8 o o 1/4 o o -1/3 o
Then the contributions from the corners are rotations of
nf(n)+αn(e(2n-e+w/6)+w(2n-w+e/6))f(n)
-αew((w+n/6)f(ww)+(e+n/6)f(ee)).
Term by term, the average of slopes at midpoints is
zero, so the ten-dimensional algorithm is a sum of three terms of the form
nf(n)+αn(e(2n-e+w/6)+w(2n-w+e/6))f(n)
-αew((w+n/6)f(ww)+(e+n/6)f(ee))
+(4βewn/3)(f(n)-f(s)).
The contributions from rotations of f(s) are
-(4βewn/3)(f(s)+f(en)+f(wn)),
and the number of flops for this one is S+21K. It will be argued presently that
this one is probably not worth doing.
There are a number of different six-dimensional approximations that fit a quadratic correction that depends on some or all of the nine inward slopes. The fastest of them is the one where only the midpoint slopes are used, and the ones at corners are discarded. The element n(1-n)=n(e+w) is zero at 'e', 'w', and 'n', the inward slope at e=w=1/2 is one, and the slopes at e=n=1/2 and w=n=1/2 are both zero (on a line where n(1-n) is constant and a maximum). The interpolation is the sum of three terms of the form
nf(n)+βn(e+w)(f(e)+f(w)-f(n)-f(s)),
and after the coefficient of f(n) has been found by reflections, the result is
three terms of the form
(n+2βew)f(n)-βn(1-n)f(s).
This one should always be
tried first - the operation count is S+11K, it does quadratic functions exactly
when β=1/2, and if that's good enough, it's truly a bargain.
Beyond this there are a number of different approximations in which the only "correction" of linear interpolation is a multiple of the invariant quadratic,
e(1-e)+w(1-w)+n(1-n) = 1-ee-ww-nn = 2(ew+wn+ne).
It only takes one of the corner or midpoint slopes
(nine parameters) to estimate the multiple, so
there is a host of algorithms that employ symmetric averages of various
kinds or asymmetric averages over slopes that are chosen at random.
The fastest algorithm of the lot is: choose one inward slope randomly
- hence in an unbiased way - then use it to estimate the amplitude of the
symmetric, quadratic correction. The operation count, S+7K, is a mere 2K flops more
than the count for uncorrected linear interpolation. For problems of
fluid mechanics,
where the fronts that traverse the inner triangle are generally asymmetric,
the use of the symmetric correction will probably turn out to
be a case of carrying not doing cubic interpolation just
a bit too far.
Next comes the matter of the tenth dimension. The nine-dimensional algorithm can be viewed as a cubic correction of quadratic interpolation that allows for the representation of functions that have inflection points on the slice from 'w' to 'e' and the rotations of that. (In smoothed representations of moving discontinuities, the position of the inflection point is sometimes thought of as the position of the discontinuity - if only we could resolve it.) If we consider the rotations of terms of the form
Aew(w+n/2)+Bew(e+n/2) = (A+B)ew/2+(A-B)ew(w-e)/2,
then the contributions to the quadratic approximation are
in the symmetric terms and the cubic
corrections are in the antisymmetric terms. Now it is possible that the above
quadratic contributions to the symmetric term, proportional to
1-ee-ww-nn, can differ significantly from those of
the other quadratic interpolation, which has the form
Ce(1-e)+Dw(1-w)+En(1-n).
The test is easy: if the sum of the rotations and reflections of f(ww) is four
times the sum of the rotations of f(s), the expressions are identical.
If not so within reason, we probably need the full twelve-dimensional algorithm
or even more; otherwise
we are comparing common multiples of
(1-ee-ww-nn)/3 and 4ewn/3.
The quadratic expression totally dominates the cubic one
near the edges of the interior triangle, and it is 9/2
times larger at the center where e=w=n=1/3. It's an intuitive result,
of course, but it looks to me that the tenth dimension cannot add
any substantial improvement of the nine-dimensional description of moving fronts.
As for what to do if quadratic interpolation is not producing an acceptable resolution of advancing fronts, the choices are to refine the computational grid
o o
/ \ / \
from / \ to o---o
/ \ / \ / \
o-------o o---o---o ,
or to employ the nine-dimensional, almost cubic interpolation.
Either of these remedies can be employed selectively in regions where
there are steep fronts to be resolved. The schematic diagram,
ew -ew ew(e-w)
* * * *
* * * *
---o-----------o-----------o-----------o-----*-----o--- ,
* * * *
* * * *
is a very rough indication that the resolution of the nine-dimensional
algorithm is near twice that of the six-dimensional ones. Even if the two is
an overestimate (make it one and a half then), the computational efforts,
4(S+11K) for grid refinement vs S+17K for
the almost cubic algorithm, provide a strong indication of what to try first.
The next level of not doing cubic interpolation comes from the suggestion that cubic and quadratic elements can be replaced by piecewise linear elements that are not qualitatively different from the cubic and quadratic ones. The rougher elements tend to overcompensate for curvature in ways that differ from setting α>1/2 or β>1/2. Here are some examples:
for n(1-n): t = (n<1/4)?n:(n>3/4)?1-n:1/4;
for ew : t = (1-n)/4; t = (1-n)*((e<t)?e:(w<t)?w:t);
for eww : t = (1-n)/4; t = (1-n)2*((e<t)?e:(w<t)?0:(w-t)/2);
for eew : t = (1-n)/4; t = (1-n)2*((w<t)?w:(e<t)?0:(e-t)/2);
for ewn/2 : t = ((e<w)?((e<n)?e:n):(w<n)?w:n;)/8
The notation, from C, means (if this is true)?(do this):(otherwise do this).
The last item in this contribution to the annals of not doing cubic
interpolation (coming soon)
That's it for now. Algorithms to follow.
A few examples of algorithms for various interpolations: