Polynomial interpolation on interlacing rectangular grids
Michael S. Floater
∗June 1, 2017
Abstract
In this paper we establish the unisolvence of any interlacing pair of rectangular grids of points with respect to a large class of associated polynomial spaces. This includes interpolation on Padua points and the schemes of Morrow and Patterson, Xu, and Erb et al. We use Newton polynomials and divided differences and an apparently new formula for determinants of certain divided difference matrices.
Math Subject Classification: Primary: 65D05, 65F40 Secondary: 65D32
Keywords: Multivariate interpolation, Padua points, divided differences, determi- nants.
1 Introduction
In [12, 14, 3, 1, 9, 8] point sets consisting of various interlacing pairs of rectangular grids in the domain [−1,1]2 have been shown to have some remarkable properties when used for approximation and cubature. Using properties of Chebyshev poly- nomials, these points sets have been shown to be unisolvent for interpolation in appropriate spaces of polynomials, with slowly growing Lebesgue constants, and as- sociated cubature formulas with high degrees of precision with respect to a Chebyshev weighting.
A notable example is the Padua points [3], which are unisolvent for polynomial interpolation of full degree N. The Lebesgue constant grows with minimal order O(log2(N)) and the associated cubature rule has degree of precision 2N −1 with respect to the Chebyshev weighting [1]. These points can be defined as the intersec- tions of the Lissajous curveγ(t) = (cosN t,cos(N+ 1)t),t∈[0, π], with itself and the boundary of [−1,1]2. Figure 1 shows the exampleN = 4. The nine points shown as circles form one grid, the six points shown as squares form another, interlacing grid.
Padua points of even degree are a modification of a similar point set proposed by Morrow and Patterson [12]. Point sets for interpolation in other polynomial spaces
∗Department of Mathematics, University of Oslo, Moltke Moes vei 35, 0851 Oslo, Norway, email:
-1 -0.5 0 0.5 1 -1
-0.5 0 0.5
1
Figure 1: Padua points, degree 4.
were studied in [12], and by Xu [14]. Recently, Erb et al. [9] and Erb [8] have studied large classes of point sets generated by Lissajous curves with respect to both interpolation and cubature.
A natural question arising from these works is whether these points sets continue to be unisolvent when their coordinates are arbitrarily spaced. For the Padua points, these more general ‘Padua-like’ points were studied in [2], [6], and [13]. It was shown that the associated Vandermonde determinant has a factorization into smaller determinants, each of which depends only on either thex coordinates of the points or they coordinates. From this it was argued that the Vandermonde determinant is non-zero.
In this paper we show unisolvence for any interlacing pair of grids with respect to a large class of associated polynomial spaces. This implies, in particular, the unisolvence of all the interpolation schemes in [12, 14, 3, 9, 8]. The basic idea is to use tensor-product interpolation, Newton polynomials, and divided differences to reduce the problem to solving a sequence of smaller linear systems. The elements of the matrices in these linear systems are divided differences. We show that these matrices are non-singular by deriving an apparently new formula for their determinants.
2 The point set
Defining
Bk,l={(i, j) : 0≤i≤k,0≤j≤l}, k, l≥0, the point set is the unionU ∪X of two rectangular grids of points inR2,
U ={(ui, vj) : (i, j)∈Bµ,ν}, X={(xi, yj) : (i, j)∈Bm,n}, that are interlaced:
u0 < x0< u1 < x1 <· · · or x0< u0 < x1 < u1 <· · · , (1) and
v0 < y0 < v1 < y1 <· · · or y0 < v0 < y1 < v1 <· · ·. (2)
It follows that|m−µ| ≤1 and|n−ν| ≤1. Therefore,m≤µ+ 1, and, by reversing the roles of U and X if necessary, we may assume without loss of generality that n≤ν. Three examples of such grids are shown in Figure 2, with black points for U and white points forX. For simplicity the figure shows uniform spacing but we allow arbitrary spacing in the two coordinate directions. For example, the dimensions and interlacing of U and X in Figure 2 (left) are the same as for the Padua points of Figure 1.
Figure 2: Examples of grids U, black points, X, white points.
3 Interpolation
We next associate with the point set an appropriate space of polynomials. LetN0 = {0,1,2, . . .}, and for any two index-pairs (i, j) and (k, l) in N2
0 we write (i, j)≤(k, l) if both i≤kand j ≤l. We call any setL⊂N2
0 a lower set [7] if it has the property that if (k, l)∈L and (i, j) ∈N20 and (i, j) ≤(k, l) then (i, j) ∈L. We can associate with any lower setL the linear space of polynomials
Π(L) := span{xiyj : (i, j)∈L}.
In the special case that
L=Ts:={(i, j) :i+j ≤s}, Π(L) = Πs, the space of bivariate polynomials of degree ≤s.
Next, let K1⊂Bm,n be a lower set, with the condition that
(0, n)∈K1 ifm=µ+ 1. (3)
Then let
K2 ={(i, j) : (m−i, n−j)∈Bm,n\K1},
which is clearly also a lower set: we can think of it as the result of rotating the ‘upper set’Bm,n\K1 through 180 degrees. In the case thatm=µ+ 1 the condition (3) is equivalent to (m,0)6∈K2.
We now define L to be the disjoint union
L=Bµ,ν∪ {(i+µ+ 1, j) : (i, j)∈K1} ∪ {(i, j+ν+ 1) : (i, j)∈K2}, which has the same cardinality asU ∪X, specifically,
(µ+ 1)(ν+ 1) + (m+ 1)(n+ 1).
L is a lower set, since, by construction, the largest j coordinate of K1 is at most ν, and the largest icoordinate ofK2 is at most µ.
We will show
Theorem 1 Given the values of a real function f on U ∪ X, there is a unique polynomialp∈Π(L) such thatp=f onU ∪X.
Let us consider some examples of the theorem. Figure 3 illustrates possible lower sets L corresponding to the three grid pairs of Figure 2. Black points are used for Bµ,ν, white points forK1 shifted to the right ofBµ,ν, and white points forK2 shifted up aboveBµ,ν.
0 0
1 2 3 4
1 2 3 4
0 0
1 2 3
1 2 3 4
0 0
1 2 3 4
1 2
5 3
4
Figure 3: Possible index sets L for the examples in Figure 2.
1. Padua and Morrow-Patterson points. For Padua points [3] of even degree and the Morrow-Patterson points of Sec. 2.2 of [12], (µ, ν) = (s, s) and (m, n) = (s, s−1), for s ≥ 1. Setting K1 = Ts−1 gives K2 = Ts−1 and Π(L) = Π2s. Figures 2 and 3 (left) show the cases= 2.
For Padua points of odd degree, (µ, ν) = (s+ 1, s) and (m, n) = (s, s), s≥1.
SettingK1=Ts−1 givesK2=Ts and Π(L) = Π2s+1.
2. Xu and Morrow-Patterson points. In [14] and in Sec. 2.3 of [12], points were studied in which (µ, ν) = (s−1, s) and (m, n) = (s, s−1), s ≥ 1. With the choice K1 = Ts−1 (which satisfies Condition (3) because (0, s−1)∈ K1), we obtain K2 = Ts−1 and Π2s−1 ⊂ Π(L) ⊂ Π2s. Figures 2 and 3 (middle) show the cases= 2.
In Sec. 2.3 of [12], Morrow and Patterson also considered grids with (µ, ν) = (s, s), (m, n) = (s−1, s−1), s≥1. Choosing K1=Ts−1 gives K2 =Ts−2 and Π2s−1 ⊂Π(L)⊂Π2s.
In Sec. 2.3 of [14], Xu also considered grids with (µ, ν) = (m, n) = (s, s),s≥1.
With the choiceK1=Ts, we obtainK2=Ts−1 and Π2s⊂Π(L)⊂Π2s+1. 3. Lissajou points of Erb et al. Large classes of points sets satisfying the conditions
of Theorem 1 arise from Lissajou curves, as studied by Erb et al. [9] and Erb [8].
Figure 2 (right) shows the type of grid pair generated by the Lissajou curve γ(t) = (sin 2t,sin 3t), t ∈ [0,2π], the first example of Sec. 2 of [9]. Figure 3 (right) shows the index setLused in Sec. 4 of [9]. Here, Π4 ⊂Π(L)⊂Π5.
4 Tensor-product interpolation
Let us start by using tensor-product interpolation on U to reduce the complexity of the problem. Letting
a(x) = (x−u0)(x−u1). . .(x−uµ), b(y) = (y−v0)(y−v1). . .(y−vν), we can express anyp∈Π(L) uniquely as
p(x, y) =p0(x, y) +a(x)p1(x, y) +b(y)p2(x, y),
wherep0∈Π(Bµ,ν), p1 ∈Π(K1), and p2 ∈Π(K2). We note that for Padua interpo- lation, the polynomialsa(x) andb(y) are the same as those used to make a basis for Π(L) in [2, Sec. 2]. So letp0 be the tensor-product interpolant tof on U,
p0(x, y) =
µ
X
k=0 ν
X
l=0
f(uk, vl)αk(x)βl(y),
where
αk(x) =
µ
Y
r6=kr=0
x−ur
uk−ur, βl(y) =
ν
Y
s=0s6=l
y−vs vl−vs,
and we will have p = f on U ∪X if we can find p1 ∈ Π(K1) and p2 ∈ Π(K2) such that
a(xi)p1(xi, yj) +b(yj)p2(xi, yj) =gij, (i, j)∈Bm,n, (4) where
gij :=f(xi, yj)−p0(xi, yj).
5 Newton form
To solve the equations in (4) we representp1 andp2in Newton form [10] with respect to the grid X,
p1(x, y) = X
(k,l)∈K1
cklφk(x)ψl(y), p2(x, y) = X
(k,l)∈K2
dklφk(x)ψl(y),
for coefficients ckl, dkl∈R, where
φ0(x) = 1, φk(x) = (x−x0)(x−x1)· · ·(x−xk−1), k≥1, (5) ψ0(y) = 1, ψl(y) = (y−y0)(y−y1)· · ·(y−yl−1), l≥1. (6) Substituting these expressions into (4) gives
X
(k,l)∈K1
ckl(aφk)(xi)ψl(yj) + X
(k,l)∈K2
dklφk(xi) (bψl)(yj) =gij, (i, j)∈Bm,n.
Next, by taking divided differences of these (m+ 1)(n+ 1) equations with respect to the grid X, we transform them into an equivalent set of equations which are easier
to solve. Applying the divided difference operator [x0, . . . , xi;y0, . . . , yj] to each side of the equation gives
X
(k,l)∈K1
ckl[x0, . . . , xi](aφk)[y0, . . . , yj]ψl+
X
(k,l)∈K2
dkl[x0, . . . , xi]φk[y0, . . . , yj](bψl) =hij, (i, j)∈Bm,n,
where
hij = [x0, . . . , xi;y0, . . . , yj](f −p0).
These equations now simplify because
[x0, . . . , xi]φk =δik, and [y0, . . . , yj]ψl =δjl, and, by the Leibniz rule for divided differences [5],
aki := [x0, . . . , xi](aφk) =
i
X
ρ=0
[xρ, . . . , xi]a[x0, . . . , xρ]φk=
([xk, . . . , xi]a, k≤i;
0, k > i,
and similarly,
blj := [y0, . . . , yj](bψl) =
([yl, . . . , yj]b, l≤j;
0, l > j.
We have thus obtained the equivalent equations X
k:(k,j)∈K1
akickj+ X
l:(i,l)∈K2
bljdil=hij, (i, j)∈Bm,n. (7)
6 Solution
The equations (7) can be solved by block back substitution, as follows. An index pair (i, j) of the lower setK1 is a maximal point ofK1 if there is no pair (k, l)∈K1, distinct from (i, j), such that (i, j)≤(k, l). In general the maximal points ofK1 are (ir, jr),r= 1,2, . . . , s, for somes≥1, with the ordering
0≤i1 <· · ·< is ≤m, n≥j1>· · ·> js≥0.
Let us assume initially that (0, n)∈ K1 and (m,0)6∈K1. Then j1 = nand is < m andK2 also hassmaximal points, (m−ir−1, n−jr+1−1),r=s, s−1, . . . ,1, where js+1:=−1. This is in fact the case in all three examples of Figure 3, wheres= 2.
Figure 4 illustrates, on the left, an example of K1 ⊂ Bm,n, where (m, n) = (20,12). Here, s = 3 and (i1, j1) = (5,12), (i2, j2) = (8,9), (i3, j3) = (16,4). The corresponding set K2 has maximal points (3,12), (11,7), (14,2), and is illustrated on the right.
We now group the equations (7) into 2sblocks. For each r= 1,2, . . . , swe apply the following two steps (we also define i0 =−1 andis+1=m):
8 12
9
4
0 5 16 20
K1
12
0 3 11 14 20
2 7
K2
Figure 4: Example of K1 and K2, both with three maximal points.
(i) For j = jr+1 + 1, . . . , jr, we write equations (7) with i = m−ir, . . . , m as the linear system
ir
X
k=0
akickj =hij−
n−jα−1
X
l=0
bljdil, α= 1, . . . , r, i=m−iα, . . . , m−iα−1−1. (8)
Since the coefficientsdil,i≥m−ir, are known, the right hand side is known and we can solve forc0,j, . . . , cir,j.
(ii) Fori=m−ir+1, . . . , m−ir−1, we write equations (7) withj=jr+1+ 1, . . . , n, as the linear system
n−jr+1−1
X
l=0
bljdil =hij −
iα
X
k=0
akickj, α= 1, . . . , r, j=jα+1+ 1, . . . , jα. (9)
Since the coefficientsckj,j≥jr+1+ 1, are known, the right-hand side is known and we can solve fordi,0, . . . , di,n−jr+1−1.
12 9
4
0 3 11 14 20
1 2
3
5 4
6
Figure 5: The equations are grouped into six blocks.
Figure 5 illustrates how the equations in the example of Figure 4 are grouped together into blocks. There are six blocks, numbered according to the order in which they are solved. Odd-numbered blocks find rows of c-coefficients, one row of coefficients for each horizontal line of equations in the block. Even-numbered blocks find columns of d-coefficients, one column of coefficients for each vertical line of equations in the block.
The algorithm is similar for the remaining cases ofj1andis. Ifj1 =nandis=m, K2 hass−1 maximal points, and the algorithm stops after 2s−1 steps. If j1 < n, the first block of equations findsd-coefficients, the second findsc-coefficients, and so on.
7 Unisolvence
We have now shown that the original problem has a solution if the linear systems in (8) and (9) are solvable. To discuss their solvability consider the two divided difference matrices
A=
[x0, . . . , xm]a · · · [xm−1, xm]a [xm]a [x0, . . . , xm−1]a · · · [xm−1]a 0
... . ..
. .. ...
[x0]a 0 · · · 0
(10)
and
B =
[y0, . . . , yn]b · · · [yn−1, yn]b [yn]b [y0, . . . , yn−1]b · · · [yn−1]b 0
... . ..
. .. ... [y0]b 0 · · ·
. (11)
For each r = 0,1, . . . , m, let Ar be the submatrix of A consisting of its first r+ 1 rows and columns (the upper left corner ofA). Similarly, for eachr = 0,1, . . . , n, let Br be the submatrix of B consisting of its first r+ 1 rows and columns.
By ordering the equations in (8) by decreasing i and the unknowns ck,j by in- creasing k, we see that the linear system (8) has a unique solution if the matrix Air is non-singular. Similarly, by ordering the equations in (9) by decreasing j and the unknowns di,l by increasing l, the linear system (9) has a unique solution if the matrixBn−jr+1−1 is non-singular.
Thus the proof of Theorem 1 will be complete if we can show that all the subma- tricesAr,r = 0,1, . . . , m, and Br,r= 0,1, . . . , n, are non-singular. It is sufficient to studyA. Consider first A0, which consists of the single element
[x0, . . . , xm]a=
m
X
k=0
qk,
where
qk=a(xk)
m
Y
j=0 j6=k
(xk−xj)−1. (12)
Since x0, x1, . . . , xm are increasing, the sign of the product in (12) is (−1)m−k. By the interlacing ofX withU, the values a(x0), . . . , a(xm) alternate in sign,
sgna(xk) = (−1)m−ka(xm),
and it follows that sgnqk= sgna(xm), and so
sgn ([x0, . . . , xm]a) = sgna(xm)6= 0, and A0 is non-singular.
In order to extend this argument to general r, we derive a formula for the deter- minant ofAr.
Lemma 1 Forr = 0,1, . . . , m, and with qk as in (12),
detAr= X
0≤k0<k1<···<kr≤m
Y
0≤i<j≤r
(xkj −xki)2qk0qk1· · ·qkr. (13)
Here, the product over an empty set is defined as 1 (corresponding to the caser= 0).
From this lemma, the proof of Theorem 1 is complete because the sign of every term in the sum in (13) is sgn (a(xm))r, and soAr is non-singular. The proof that theBr are non-singular is similar.
As an example, for m= 2, Lemma 1 shows that the determinant ofA1 is
[x0, x1, x2]a [x1, x2]a [x0, x1]a [x1]a
= (x1−x0)2q0q1+ (x2−x0)2q0q2+ (x2−x1)2q1q2,
where
q0 = a(x0)
(x0−x1)(x0−x2), q1= a(x1)
(x1−x0)(x1−x2), q2 = a(x2)
(x2−x0)(x2−x1). This can be verified by direct calculation using the observation that
[x0, x1, x2]a=q0+q1+q2,
[x0, x1]a= (x0−x2)q0+ (x1−x2)q1, [x1, x2]a= (x1−x0)q1+ (x2−x0)q2,
[x1]a= (x1−x2)(x1−x0)q1.
Sincex0< x1< x2, and assuming thata(x0)>0,a(x1)<0,a(x2)>0, we find that q0, q1, q2>0 and so detA1>0.
To prove Lemma 1 observe that the elements of A in (10) are Aij = aj,m−i for i, j = 0, . . . , m. We start by expressing these elements as linear combinations of q0, . . . , qm.
Lemma 2 Fori, j= 0, . . . , m,
Aij =
m
X
k=0
φˆi(xk)φj(xk)qk, (14)
where qk is given by (12),φj is given by (5), and
φˆ0(x) = 1, φˆi(x) = (x−xm)(x−xm−1)· · ·(x−xm−i+1), i≥1.
Proof. Consider first the case that i+j≤m. Then
Aij = [xj, . . . , xm−i]a=
m−i
X
k=j
a(xk)
m−i
Y
s=j s6=k
(xk−xs)−1 =
m−i
X
k=j
φˆi(xk)φj(xk)qk.
Then (14) follows from the fact that φj(xk) = 0 whenever k < j and ˆφi(xk) = 0 whenever k > m−i.
Otherwise, i+j > m. Then (14) follows from the fact that ˆφi(xk)φj(xk) = 0 for
all k= 0, . . . , m. ✷
Proof of Lemma 1. By Lemma 2, we can writeAr as the matrix productAr=CrDr where
Cr= [ ˆφi(xk)]i=0,...,r,k=0,...,m, Dr= [φj(xk)qk]k=0,...,m,j=0,...,r. By the Cauchy-Binet theorem [11, p. 17] applied to this product,
|Ar|= X
0≤k0<k1<···<kr≤m
|φˆi(xkj)|i,j=0,...,r|φj(xki)qki|i,j=0,...,r.
Since
|φj(xki)qki|i,j=0,...,r =qk0qk1. . . qkr|φj(xki)|i,j=0,...,r, the proof is completed by observing that
|φˆi(xkj)|i,j=0,...,r =|φj(xki)|i,j=0,...,r =|xikj|i,j=0,...,r = Y
0≤i<j≤r
(xkj−xki).
8 Final remark
While the purpose of this work was to prove unisolvence, the method of proof provides an algorithm for finding the polynomial interpolant p in Theorem 1. However, this algorithm will be subject to the numerical sensitivity of divided differences, and the conditioning of the matrices Ar and Br, and it might be better, numerically, to represent p in some other way, such as, for example, using a Chebyshev basis, as proposed for Padua points in [4]. This is a topic for future research.
Acknowledgement. I wish to thank the referees for their careful reading of the manuscript and their valuable comments which have helped to improve the paper.
References
[1] L. Bos, M. Caliari, S. De Marchi, M. Vianello, and Y. Xu, Bivariate Lagrange interpolation at the Padua points: the generating curve approach, J. Approx.
Theory 143 (2006), 15–25.
[2] L. Bos, S. De Marchi, and S. Waldron, On the Vandermonde determinant of Padua-like points, Dolomites Research Notes on Approximation 2(2009), 1–15.
[3] M. Caliari, S. De Marchi, and M. Vianello,Bivariate Lagrange interpolation on the square at new nodal sets, Appl. Math. Comp.165 (2005), 261–274.
[4] , Bivariate Lagrange interpolation at the Padua points: computational aspects, J. Comp. Appl. Math. 221 (2008), 284–292.
[5] C. de Boor, Divided differences, Surveys in Approximation Theory 1 (2005), 46–69.
[6] S. De Marchi and K. Usevich, On certain multivariate Vandermonde determi- nants whose variables separate, Lin. Alg. Appl. 449 (2014), 17–27.
[7] N. Dyn and M. S. Floater,Multivariate polynomial interpolation on lower sets, J. Approx. Theory 177 (2014), 34–42.
[8] W. Erb,Bivariate Lagrange interpolation at the node points of Lissajous curves - the degenerate case, Appl. Math. Comp. 289(2016), 409–425.
[9] W. Erb, C. Kaethner, M. Ahlborg, and T. M. Buzug,Bivariate Lagrange inter- polation at the node points of non-degenerate Lissajous curves, Numer. Math.
133 (2016), 685–705.
[10] E. Isaacson and H. B. Keller,Analysis of numerical methods, Dover, 1994.
[11] S. Karlin,Total positivity, Vol. I, Stanford University Press, Stanford, Califor- nia, 1968.
[12] C. R. Morrow and T. N. L. Patterson,Construction of algebraic cubature rules using polynomial ideal theory, SIAM J. Numer. Anal. 15(1978), 953–976.
[13] A. Pierro de Carmago and S. De Marchi, A few remarks on “On certain Van- dermonde determinants whose variables separate”, Dolomites Research Notes on Approximation8 (2015), 1–11.
[14] Yuan Xu, Lagrange interpolation on Chebyshev points of two variables, J. Ap- prox. Theory 87(1996), 220–238.
