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."
This is more or less the same sentiment, as applied to cubic
interpolation.
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. 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 on lines and within equilateral triangles. The gambit is to formulate cubic interpolation, squeeze it as hard as possible, and then not do it. The symmetries of the equilateral triangle will be used along with the homogeneous coordinates, (e,n,w) s.t. e+n+w=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=-1/ e=0/ \w=0 \w=-1
/ \ / / \ \
/ \ n=0-----ww------w-------e------ee
--w=1---------e=1-- \ / \ /
/ n=0 \ w=2\ / \ /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,'wn' and 'nw', are meant
to suggest west of north and north of west, for
example. For a function f(e,n,w), its values,
f(e)=f(0,1,0),
f(ww)=f(-1,0,2), etc., are to be used to
construct a cubic
polynomial for interpolations in and on the boundary of
the interior triangle.
We begin with the linear interpolation,
fo(e,n,w) = ef(e) + nf(n) + wf(w),and we seek to correct that with polynomial approximations that are zero at the points 'e', 'n' and 'w'. (This approach will be modified later to make it more efficient.) Some polynomials with
P(e)=P(n)=P(w)=P(/)=P(\)=P(-)=P(|)=0are:
ew eew eww
0 0 0 0 0 0
\ / \ / \ /
-1 n -1 1 n -1 -1 n 1
/ \ / \ /
-2 w e -2 2 w e -4 -4 w e 2
/ \ / \ / \
0 1 0 0 1 0 0 1 0
enw en(1-n)w ee+nn+ww-1
0 0 0 0 4 4
\ / \ /
-1 n -1 0---n---0 2 0 2
/ \ / \
0---w---e---0 0---w---e---0 4 0 0 4
/ \ / \
0 -1 0 0 -2 0 4 2 4
From these and two rotated sets
of diagrams, we are to construct a basis
of nine independent polynomial functions, and the basis must include
three independent quadratics and the four independent cubics.
The easy choice, not among the required ones, is the quartic(s)
Q(e,n,w)=en(1-n)w/2
No algebra is required to see that
Q(e,n,w)+Q(n,w,e)+Q(w,e,n)=enw
Thus the three Q's are not independent, nontrivial quartics,
but they are independent among quartics plus cubics.
Algebraically, or from the diagrams, symmetric and antisymmetric
combinations of the cubic polynomials are
ew(e+n+w)=ew ew(e-w)
0 0 0 0 0
\ / \ | /
\ / \ | /
\ / \|/
-1 n -1 2 n -2
/ \ /|\
/ \ / | \
/ \ / | \
-2 w e -2 4 w | e -4
/ \ / | \
/ \ / | \
/ \ / | \
0 1 0 0 0 0
Finally, it follows from (e+n+w)(e+w+n)=1, or from the diagrams,
that
en+nw+we=(1-ee-nn-ww)/2
Now we drop the references to physical directions altogether but keep the coordinate names 'enw'. Then there are just three basis functions, eew, eww and en(1-n)w, and they are used in the six orientations indicated below:
o o + o o -
\ / / \
- n + + e + + w -
/ \ / \
- w e + o---n---w---o o---e---n---o
/ \ / \
o + o o - - + + o
o + o - - o o + +
\ / / \
+ e w - o---w---n---o o---n---e---o
\ / / \
+ n - + e + - w +
/ \ / \
o o o + - o
Given are the coordinates, e,n & w, and the values of f at the
markeded positions on the diagrams above.
Various types of corrections are added to the flawed, linear
approximations
to compensate for their systematic errors. We begin with
some polynomial approximations.
Setup: f=n*f(n)+e*f(e)+w*f(w); (5N flops)
One pass (3 times):
Cew=e*w/4; (2 flops)
f=f-(f(ww)+f(ee)-f(w)-f(e))*Cew; (5N flops)
Although it is only the average of f(ww) and f(ee) that is
used, this is still better than
no correction at all. The operation count of 6 + 20N is pretty
good, but the coefficient of N can be improved by moving the
terms proportional to f(w)+f(e) from the 3N-times repeated calculation
to the setup. Then we have
Setup: Cew=e*w/4; Cne=n*e/4; Cwn=w*n/4; (6 flops)
Cn=n+Cwn+Cne; Cw=w+Cew+Cwn; Ce=e+Cne+Cew; (6 flops)
f=Cn*f(n)+Cw*f(w)+Ce*f(e); (5N flops)
One pass(3 times, and permute the pointers to Cew, Cne and Cwn too):
f=f-Cew*(f(ww)+f(ee)); (3N flops)
The coefficient of N in the cubic and related algorithms can be minimized by adopting a different point of view. Consider the missing diagrams,
P=n P=w P=e
2 2 0 -1 -1 0
1 1 1 1 0 -1 -1 0 1
0 0 0 0 2 1 0 -1 -1 0 1 2
-1 -1 -1 2 1 0 0 1 2
The elements of these that are not at the vertices, n, w and e,
are the coefficients of f(n), f(w) and
f(e) that appear in the linear
extrapolations of ef(e)+nf(n)+wf(w) that are needed if the cubic
fit is treated as a correction of the linear interpolation.
It is far more efficient to construct a basis that does all of that
without regard to a specific function f. The hitch is that no
new polynomials can be introduced; i.e. the only quartic
polynomials allowed are the rotations of ewn(n-1).
Instead of algebra, consider
eww eew ew(2w+e)
0 0 0 0 0 0
-1 0 1 1 0 -1 -1 0 1
-4 0 0 2 2 0 0 -4 -6 0 0 0
0 1 0 0 1 0 0 3 0
ewn(n-1)/2 ewn(3(1-n)+e-w)/2 ew(n-1)(1+2n+w)/6
0 0 0 0 0 0
0 0 0 1 0 -1 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0
0 1 0 0 -3 0 0 0 0
The first and third of the second row of diagrams above, plus
reflections and rotations, are what is needed to deal with all of
the exterior vertices. The result for the interior vertices is
in the diagrams,
n(1+ew(2-e-w)/2) nn(e+w)/2 n(1+n)(2-n+ew)/2
2 2 -2 -2 0 0
0 1 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0
-1 -1 -1 1 1 1 0 0 0
With these basis functions the cubic+2 algorithm is something like
Setup: C1=e*w*n; C2=C1+C1; C3=C2+C1: (4 flops)
Cn=(1+n)/2; Cn=Cn*(n*(2-n)+C1); (6 flops)
Cw=(1+w)/2; Cw=Cw*(w*(2-w)+C1); (6 flops)
Ce=(1+e)/2; Ce=Ce*(e*(2-e)+C1); (6 flops)
f=Cn*f(n)+Ce*f(e)+Cw*f(w); (5N flops)
One pass (3 times):
Cs=(n-1)/6; Cew=Cs*e*w; Ct=Cew+Cs*C2; (6 flops)
Cww=Ct+w*Cew; (2 flops)
Cee=Ct+e*Cew; (2 flops)
Cs=Cs*C3; (1 flop)
f=f+Cww*f(ww)+Cee*f(ee)+Cs*f(s); (6N flops)
The result, 55 + 23N flops, has the lowest coefficient of N that
can accomodate twelve degrees of freedom. Algorithms with ten and with
nine degrees of freedom will be discussed later.
After quadratic interpolation, which might be good enough for the task at hand, the next level of not doing cubic interpolation follows from the observation that the polynomials like eww are only used within the interior triangle. The element eww can also be regarded as the cubic element whose slope is 1 at 'w' and zero at 'e', and estimating a slope of f at 'w' by doing cubic interpolation is not the only possibility. 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 o | o e=(1-w)/2
/ o o \
o|
/ o o \
o | o
/ o o \
n=0 ------ w -------- | -------- 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 ewn(1-n)
0 0 0 0
\ / \ /
n n
/ \ / \
1/4 / | \ 0 1/8 / | \ 1/8
o o o o
/ o| \ / o| \
/ o o \ / o o \
1 --w-----|-----e-- 0 0 --w-----m-----e-- 0
/ \ / \
0 -3/8 0 0 1/4 0
Next, it follows without much effort that polynomials that are zero
at 'e', 'w', and 'n' and have only one nonzero slope at a corner
or at the midpoint of an edge of the interior triangle are:
4ewn(1-2n) ew(w+n(2e-3n+1/2))
0 0 0 0
\ / \ /
n n
/ \ / \
0 / | \ 0 0 / | \ 0
o o o o
/ o| \ / o| \
/ o o \ / o o \
0 --w-----m-----e-- 0 1 --w-----|-----e-- 0
/ \ / \
0 1 0 0 0 0
On the line 'ww'-'w'-'e' the values of corrections of the linear approximation of f(e,n,w) are f(ww)+f(e)-2f(w) at 'ww' and 0 at 'w' and 'e'. The estimate of the associated slopes at 'w' and 'e' are
S(ww)=a*(2f(w)-f(e)-f(ww)).
S(ee)=a*(2f(e)-f(w)-f(ee)).
The similar results on the line 's'-'m'-'n' are f(s)+f(n)-f(e)-f(w) at 's' and
0 at 'm' and 'n'. The estimate of the associated slope at 'm' is
S(s)=b*(f(e)+f(w)-f(n)-f(s)).
In the symmetric case the six a's are the same and the three b's are the same. Notable values of a and b are: 0 for no correction, 1/2 for central differences and 1 (say) for an estimate that overcompensates for numerical dissipation. Some approximations that appear to be worth the implementations and comparisons of results are all sums of three terms of the forms,
a=b=0: nf(n)
a=1/2: nf(n)+(f(w)-(f(e)+f(ww))/2)*ew*(w+n*(2e-3n+1/2))
+(f(e)-(f(w)+f(ee))/2)*ew*(e+n*(2w-3n+1/2))
b=1/2: +2(f(e)+f(w)-f(n)-f(s))*ewn*(1-2n)
a=1: nf(n)+(2f(w)-f(e)-f(ww))*ew*(w+n*(2e-3n+1/2))
+(2f(e)-f(w)-f(ee))*ew*(e+n*(2w-3n+1/2))
b=1: +4(f(e)+f(w)-f(n)-f(s))*ewn*(1-2n)
The last item in this contribution to the annals of not doing cubic interpolation is the suggestion that cubic (eg. eww) and quartic (eg. ewn(1-2n)) elements can be replaced by piecewise linear or even piecewise constant elements that are not qualitatively different from the cubic and quartic ones. The rougher elements tend to further overcompensate for curvature. Here are two examples of algorithms for such elements:
(1) for eww: do t=(1-n)/3; Ceww=(t>e)?e:(t>w)?0:w-t; (2) for 4ewn(1-2n) do Cew=(e>w)?w:e; Cewn=(n>Cew)?Cew:n; Cs=(1-2n)*Cewn;