This is a somewhat easier, one-dimensional version of the
sort of simulated cubic interpolation that is described
in On not doing cubic interpolation
The context is the same in both cases - 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. Each of the eight faces of the octahedron shares its
three edges with three other faces, and data in the rows of nodes closest
to the edges can be used to find external boundary data for
the neighboring faces. The
required interpolations are described in
Stitching the patches.
The generic case for interpolation of f(x) is:
Given data here: o o o o o Find approximations here: x x x xFor present purposes, the o's can be taken to be uniformly spaced, with values of f(o) specified at integers in [0,N-1] (divided by N-1). Then the generic problem is to find interpolations of f(x) at integers in [0,N+1] (divided by N+1). Again, the example with N=5 on the edge that connects the atlantic and east poles is:
common domains: [0,1] = [0,4]/4 = [0,5]/5 = [0,6]/6 interpolated: n x x x x x n one row north: o n o o o n o on the equator: o A o o o o E o one row south: o s o o o s o interpolated: s x x x x x s(Note how the diagram is slightly bent.)
If f( [-1 0 1 2] ) are used for the interpolation on the interval [0,1] a cubic polynomial can be fitted, and that will certainly take care of the systematic errors in linear interpolation. Indeed that will be done first, and then simulations of cubic interpolation that behave qualitatively like the real thing will be introduced.
It is convenient (and efficient when several functions are to be interpolated) to introduce homogeneous coordinates, e and w with e+w=1, for each interval [k,k+1]. The naming of nodes is then
e=-1 e=0 e=1 e=2
| | | |
ww - - - - - - w - - - - - - e - - - - - - ee
| | | |
w=2 w=1 w=0 w=-1
and f(ww), f(w), f(e) and f(ee) will be used for the values of f at the marked nodes.
The linear interpolation of f is then
L[f] = ef(e)+wf(w)and the deviations from linear are
| | | | f(ww)+f(e)-2f(w) - - - 0 - - - - - - 0 - - - f(ee)+f(w)-2f(e) | | | |The values of two cubic functions that are zero at e=1 and w=1 are
| | | | | | | | eew: 2 0 0 -4 and eww: -4 0 0 2 | | | | | | | |and the values of e and w are
| | | | | | | |
e: -1 0 1 2 and w: 2 1 0 -1
| | | | | | | |
Personally, I prefer using the diagrams over doing the algebra to find the
Lagrange polynomials,
P(e,w) = -(eww+2eew)/6 = -ew(e+1)/6 ( 0 0 0 1 ) and Q(e,w) = e - 2P(e,w) + P(w,e) = e - eew/2 ( 0 0 1 0 )That gives the cubic interpolation formula,
C[f] = f(ww)P(w,e) + f(w)Q(w,e) + f(e)Q(e,w) + f(ee)P(e,w) .
In the application of this to atmospheric circulation there are quite a few functions to be interpolated (certainly no less than five), and the use of the Lagrange polynomials is more efficient than other schemes that might be somewhat better for a single interpolation. For K functions the number of floating point operations is S + 7K, where S is the number needed to evaluate the P's and Q's, and 7K flops is the least number that can accomodate the four degrees of freedom. Given e and w, the best result I have found so far for S is 11 flops.
The situation in one dimension is different from that in higher dimensions: here any symmetric correction of linear interpolation values of e or w in [0,1] involves the four values of f(e) or f(w) at the very least. However, that does not imply that the use of the Lagrange polynomials, P and Q, is somehow the best procedure among those that use data for e or w in [-1,2] only.
The first of two ways to abandon the Lagrange polynomials is to consider an approximate Hermite interpolation in which the slopes of the corrections to linear interpolation are estimated from data for e in [-1,1] on the left (west) side and data for w in [-1,1] on the right (east) side. Again, the values of deviations from the linear interpolation, written D[f] = f-L[f], are
D[ww] = f(ww)+f(e)-2f(w), D[ee] = f(ee)+f(w)-2f(e) and D[e]=D[w]=0,and a one-parameter family of corrected interpolations is
C[f] = ef(e) + wf(w) - αew(eD[ee]+wD[ww]).
Notable values of α are: 0 for no correction, 1/2 for the central difference
approximations of the slopes,
and 1 (say) for estimates that overcompensate for numerical dissipation.
The result of taking α to be greater than 1/2 in numerical algorithms for
problems like
f(x,t) : df/dt ≡ ∂f/∂t + ∂f/∂x = 0 <=> f(x,t+dt) = f(x-dt,t)is that weak fronts that span many intervals of dx (where de=1 and dw=1) steepen with time until they span just a few intervals of dx (2 to 3). It is somewhat like what might be expected if negative dissipation were introduced, but these algorithms are stable if α is not greater than one. (I have not investigated cases where α > 1.)
When the coefficients of the values of f in C[f] are collected, the result is
C[f] = (e+2Ce-Cw)f(e) + (w+2Cw-Ce)f(w) - f(ee)Ce - f(ww)Cw
where Ce = αeew and Ce = αeww.
Again the number of floating point operations is S + 7K, and the
best result I have found so far for S is 10 flops. The saving of a flop or
two in the evaluation of S is of no consequence: it is that α>1/2
can be used to reverse the effects of numerical dissipation that may
make this kind of interpolation more attractive than the use of Lagrange
polynomials in some situations.
The only reduction of the coefficient of K in these one-dimensional formulae comes from the replacement of the separate slopes at 'e' and 'w' by the average of the two. Then the more roughly corrected result is
C[f] = (e+Cew)f(e) + (w+Cew)f(w) - Cew(f(ee)+f(ww))
where Cew = αew/2.
This, with 5 + 6K flops, appears to be about the least that can possibly be done
to correct the (3K flop) linear interpolation.
The last level of not doing cubic interpolation comes with the replacement of eew and eww by piecewise linear or even piecewise constant functions that are not qualitatively different from the polynomials. Piecewise linear examples are:
eew ==> (e<1/3)?e:(w<1/3)?0:w-1/3
eww ==> (w<1/3)?w:(e<1/3)?0:e-1/3
(The notation, from C, means (if this is true)?(do this):(otherwise do this).)
Again, the saving of a flop or
two in the evaluation of S is of no consequence, but the corners in the
results make it easier to follow the progress of a discontinuity of f(x,t)
as it passes through
an interval where it cannot be resolved because dx is fixed and greater than zero.
(The information is in the cubic
polynomials too, but it is harder to see it there.)
The other feature of these approximations is that they also tend to overcompensate
for numerical dissipation.