Interpolation on Edges
Before specific algorithms that use interpolation comes a somewhat more general treatment of the processes.
In the general case of interpolation on an edge, the vertices, at which points data are 'known', are not
uniformly separated. The case where vertices are evenly spaced, to be treated here, is not difficult
to generalize, and the more general case is not needed for interpolation on the global grid
for this project.
To make this more like what goes on in triangular patches, and to anticipate selective grid bisections,
the familiar idea of a basic coordinate 'x' with evenly spaced points, xk, will be replaced with a
local "skeleton" made up of equally spaced vertices, 'ww,w,e,ee', and local pairs of coordinates, (e,w).
The relations between these and the global 'x' are indicated below:
| vertices: | xk-1 | | xk |
| xk+1 | | xk+2 |
| separations: | | dx |
| dx | | dx |
| F(x): | Fk-1 | | Fk |
| Fk+1 | | Fk+2 |
| |
| vertices: | 'ww' | | 'w' |
| 'e' | | 'ee' |
| separations: | | 1 |
| 1 | | 1 |
| F(e,w): | Fww | | Fw |
| Fe | | Fee |
| (e,w) | (-1,2) | | (0,1) | |
(1,0) | | (2,-1) |
| |
|
e = (x - xk)/dx and w = (xk+1 - x)/dx |
(1) |
|---|
From eq(1) it follows that e + w = 1, so the coordinates (e,w) are like the triangular
coordinates of the global grid, as they occur on patches at their bases, where .n = .s = 0.
(Similar ones also apply on the slanted sides where
.e = 0 or .w = 0.) The point of doing things this way is to
make the 'interpolators' universal ; i.e. independent of what is to be interpolated
(one or many functions),
independent of where on the global grid it takes place, and independent of the level of
grid refinement (dx → dx/2 →
dx/4 . . . ).
The cases to
be considered here are: Linear interpolation, Hermite (cubic)
interpolation, a quadratic approximation, and simulations of Hermite interpolation
that will be called Babylonian.
Linear interpolation
This is quickest and easiest because it doesn't use Fww (west of west) or
Fee (east of east). Unfortunately, its error is systematic, and it needs improvement.
Nonetheless, it is a good, uncluttered example of the use of a basis of polynomial
functions for interpolation __ an idea that gets ever better as the number of
functions to be interpolated increases, and better yet when the move
from edges to patches is made.
Now, since e and w are linear functions of x,
F1(e,w) = eFe + wFw (2)
is the linear interpolation of F(x) that is Fw at xk and Fe
at xk+1.
The next polynomial that might be added to the two linear ones could be e2,
w2 or ew, and the preferred choice is the symmetric quadratic interpolation,
F2(x) = eFe + wFw + kew ,(3)
where k, to be discussed later, is a measure of the curvature of F. For present purposes
it suffices to point out that |k| is large, and k changes its sign in the neighborhood
of a steep front (e.g. a cold front). In an attempt to follow the motion of such a front with a
semi-Lagrangian algorithm (which uses interpolation), linear interpolation systematically underestimates
|k|, and the result is a numerical degradation of steepness that
has nothing to do with the physics of the situation. A rough diagram of it, for a cold front moving
eastward, is this:
a << 0 a < 0
sooner oooooooooo later o ooo
o o
o
==> o
o
o o
oooooooooo ooo o
a >> 0 a > 0
Cubic interpolation
Consider first the quartic polynomial, P = (e+1)ew(w+1), which is zero at
the four vertices 'ww', 'w, 'e' and 'ee' (listed here in the same order as the factors in P).
If we remove one of the factors from P and divide
the result by its value at the associated vertex, the resulting cubic polynomials are:
Pw = (e+1)(w+1)w/2 = (ew +2)w/2 and Pww = -(w+1)ew/6
Pe = (e+1)(w+1)e/2 = (ew +2)e/2 and Pee = -(e+1)ew/6 .
(4)
(Note, the second form of Pw follows from e + w = 1, and the second line follows
from the first by symmetry.) These are the Lagrange polynomials; their values are one
at the vertices indicated by the subscripts, and, zero at the other three. Thus, the
cubic function that is correct at the four vertices is
F3(x) = PwwFww + PwFw
+ PeFe + PeeFee
= (ew +2)(eFe + wFw)/2 -
ew((e+1)Fee + (w+1)Fww)/6 .(5)
If there is just one function to be interpolated, the two forms of eq(5) involve
pretty much the same arithmetic, but even for two functions it is more efficient to
evaluate the Lagrange polynomials first, and then to use the first form of eq(5) for the several functions.
The polynomials in eq(4) constitute a basis for cubic interpolation, but it
is not the one that is suitable for entities like shock-waves or cold fronts, where F(x) varies
rather slowly with x except at the fronts, where F looks more like the diagram at the left
between eqs(3) and (4). Because the approximation is cubic over three intervals, and correct at four
vertices, there are alternate maxima and minima between the vertices, as roughly indicated below:
o
o
o-- --w-- --e-- o-- --w-- o--
o
o o o
o
--o --e-- --o --w-- --e-- --o
o
o
These overshoots are not physical, and the use of the Lagrange polynomial basis introduces
numerical waves, with wavelengths that are near 2dx to 4dx. The phenomenon
is somewhat like
the Gibbs phenomenon, which occurs in the
Fourier analysis of discontinuities, but fixing it is a little easier here. (See
Gibbs fix for a relatively easy approach
to the problem with Fourier series.)
Hermite interpolation
The key is to give up the requirement that the cubic approximation shall have correct values
at 'ww' and 'ee'. After all, the linear interpolation needs correcting only in the 'we'-interval,
so there is no reason not to use Fww and Fee to estimate the slopes
of F at 'e' and 'w'.
The basis is constructed in two stages, beginning with the assumption that
we know the slopes,
F'w = dF/de at e = 0 and
F'e = dF/dw at w = 0 .(6)
Consider first the function P(e,w) = eFe + wFw + F'eww. From e + w = 1 (and
dw/de = de/dw = -1) it follows that
P(0,1) = Fw,
P'(0,1) = F'w + Fe - Fw and
P(1,0) = P'(1,0) = 0. Hence the Hermite interpolation
of F(x) on the 'ew'-interval is
F3(e,w) = eFe + wFw +
eew(F'e + Fe - Fw) +
eww(F'w + Fw - Fe)
= (e + eew - eww)Fe + (w + eww - eew)Fw + eewF'e
+ ewwF'w ,(7)
and the second line is an evaluation of the polynomial basis.
The second stage is to introduce approximate Hermite interpolations by providing
estimates of F'e and F'w. The most direct approach
is to begin with the values of a few functions:
| F \ V | 'ww' | 'w' | 'e' | 'ee' |
| F1 = eFe + wFw | 2Fw - Fe
| Fw | Fe | 2Fe - Fw |
| F - F1 | Fww + Fe - 2Fw | 0 |
0 | Fee + Fw - 2Fe |
| eww | -4 | 0 | 0 | 2 |
| eew | 2 | 0 | 0 | -4 |
Different estimates of (F-F1)'w and (F-F1)'e
now give different estimates of the curvature
k = wkw + eke.
The one-parameter family of results like eq(3) is
F3(e,w)
= eFe + wFw + aew(w(2Fw - Fe - Fww)
+ e(2Fe - Fw - Fee))
= (e + a(2e-w)ew)Fe +
(w + a(2w-e)ew)Fw
- aew(eFee + wFww) .(8)
The special cases are:
| no correction | __ | __ | a = 0 |
| collocation | F3ww = Fww |
F3ee = Fee | a = 1/4 |
| central difference | F'w = (Fe - Fww)/2 |
F'e = (Fw - Fee)/2 | a = 1/2 |
| extrapolation | F'w = Fw - Fww |
F'e = Fe - Fee | a = 1 |
The choice of a will be left open for now, but it is unlikely that any value less than 1/2 will be used.
The quadratic approximation
Evaluations at e = w = 1/2 of the coefficients of ew in eq(8) give the symmetric quadratic corrections,
F2(e,w) = (eFe + wFw)
+ aew(Fe - Fee + Fw - Fww)/2 .(9)
This, if it works well enough, is a bargain.
Babylonian approximations
This is about the use of linear approximations of F(e,w) to simulate the cubic
(Hermite) interpolations in eq(8). The
linear functions are the two extrapolations and the interpolation,
Fxw(e) = Fw + eF'w ,
F1(e,w) = wFw + eFe and
Fxe(w) = Fe + wF'e .(10)
Only the values a = 1 and a = 1/2 will be treated, and there are no multiplications
or divisions in the evaluations,
F'w = Fw - Fww or (Fe - Fww)/2
and
F'e = Fe - Fee or (Fw - Fee)/2 .(11)
(Note: /2, a bit shift, isn't really a division.) Likewise, there are no multiplications in the tests for three cases:
| case | test | action |
| positive | Fw + F'w > Fe
and
Fe + F'e > Fw |
Fx is the lesser of Fxe and Fxw
|
| negative | Fw + F'w < Fe
and
Fe + F'e < Fw |
Fx is the greater of Fxe and Fxw
|
| inflection | neither of the above | F1 = eFe + wFw |
The argument here is that the linear interpolation, F1, is not systematically wrong in an interval where F has an inflection point.
Just for fun, the case where a = 1/2 will be called superlinear interpolation,
and the case where a = 1, ultralinear interpolation. They will be used in some examples of simple transport equations,
where the piecewise linear profiles make it easier to understand how the semi-Lagrangian algorithms manage to capture and propagate
shock-waves.
Fearful Asymmetry
As helpful as the symmetry of the formulas in this section may be, when an algorithm is to be used
tens of thousands of times per time-step of a simulation, it is worth the effort to hack it,
even if that only trades one multiplication for an add or two. Hence,
F1(e,w) = eFe + wFw
→
F1(e,w) = Fw + e(Fe - Fw)(2a)
trades 2'x'1'+' for 1'x'2'+'. It's not a very large improvement here in one dimension, but it becomes more
important in two.