Journal of Approximation Theory 143 (2006) 15 – 25

www.elsevier.com/locate/jat

Bivariate Lagrange interpolation at the Padua points:

The generating curve approach
Len Bosa , Marco Caliarib , Stefano De Marchic,∗ , Marco Vianellob ,
Yuan Xud
a Department of Mathematics and Statistics, University of Calgary, Canada
b Department of Pure and Applied Mathematics, University of Padua, Italy
c Department of Computer Science, University of Verona, Italy
d Department of Mathematics, University of Oregon, USA

Received 20 October 2005; accepted 3 March 2006

Communicated by Juan Manuel Peña
Available online 15 May 2006

Abstract
We give a simple, geometric and explicit construction of bivariate interpolation at certain points in a
square (called Padua points), giving compact formulas for their fundamental Lagrange polynomials. We
show that the associated norms of the interpolation operator, i.e., the Lebesgue constants, have minimal
order of growth of O((log n)2 ). To the best of our knowledge this is the first complete, explicit example of
near optimal bivariate interpolation points.
© 2006 Elsevier Inc. All rights reserved.

1. Introduction

Suppose that K ⊂ Rd is a compact set with non-empty interior. Let
n (Rd ) be the space
of polynomials of degree
n in d variables, and
n (K) be the restrictions of this space to K.
Then
n (K) is a vector space of some dimension, say, N := dim(
n (K)). Given N points
A := {Ak }Nk=1 ⊂ K, the polynomial interpolation problem associated to
n (K) and A is to find
for each f ∈ C(K) (the space of continuous functions on K) a polynomial p ∈
n (K) such that

p(Ak ) = f (Ak ), k = 1, . . . , N.

∗ Corresponding author. Fax: +39 0458027928, +39 0458027068.

E-mail address: demarchi@sci.univr.it (S. De Marchi).

0021-9045/$ - see front matter © 2006 Elsevier Inc. All rights reserved.
doi:10.1016/j.jat.2006.03.008
16 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25

If this is always possible the problem is said to be unisolvent, and if this is indeed the case we
may construct the so-called fundamental Lagrange polynomials
j (x) with the property that

j (Ak ) =
j,k ,

the Kronecker delta. Further, the interpolant itself may be written as

N
(Lf )(x) = f (xk )
k (x).
k=1

The mapping f  → Lf may be regarded as an operator from C(K) (equipped with the uniform
norm) to itself, and as such has an operator norm L. Classically, when K = [−1, 1], this norm
is known as the Lebesgue constant and it is known that then L C log n and that this minimal
order of growth is attained, for example, by the Chebyshev points (see e.g., [4]).
In the multivariate case much less is known. Sündermann [6] has shown that for K = B d , the
unit ball in Rd , d 2, the Lebesgue constant has a minimal rate of growth of at least O(n
d
(d−1)/2 ).

In the tensor product case, when K = [−1, 1] and with the polynomial space k=1
1n ,
d

LC(log n)d and this minimal rate of growth is attained for the tensor product of the univariate
Chebyshev points. However, even for the cube and the polynomials of total degree n, i.e., for K =
[−1, 1]d , was not known until now whether the minimal rate of growth of O((log n)d ) could be
attained. In [2] it was shown that Xu’s interpolation scheme [8,1] in the square (which uses a
space of polynomials intermediate to
n−1 (R2 ) and
n (R2 )) the Lebesgue constant does grow
like O((log n)2 ). In this paper we show that there is also such an interpolation scheme for the
case d = 2, K = [−1, 1]2 , and the total degree polynomial space
n (R2 ).
We give an explicit, simple, geometric construction of a point set in K = [−1, 1]2 , nicknamed
“Padua points” (cf. [5]), for which we may construct compact formulas for the fundamental
Lagrange polynomials. These compact formulas allow us to easily analyze the growth rate of the
associated Lebesgue constants.
We point out that this establishes a remarkable
√ dichotomy between the cases of the disk where
the optimal growth would be of at least O( n) and the square where we now know that it is
O((log n)2 ).

2. The Padua points and their generating curve

For the remainder of this paper we take K = [−1, 1]2 for which
n (K) =
n (R2 ) and has
dimension
 
n+2
dim(
n (K)) = .
2

For degree n = 0, 1, 2, . . . , consider the parametric curve

n (t) := (cos(nt), cos((n + 1)t)), 0
t
. (1)

Obviously, this curve is contained in the unit square [−1, 1]2 . We also remark that it is an algebraic
curve given by Tn+1 (x) = Tn (y) where Tn denotes the classical Chebyshev polynomial of the
first kind.
 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25 17

The curve γn (t) and associated points for n = 5

1

0.8

0.6

0.4

0.2

0.2

0.4

0.6

0.8

1
1 0.5 0 0.5 1

Fig. 1.

Now, on n (t) consider the point set of equally spaced points (in the parameter t)
   
k
An := n  : k = 0, 1, . . . , n(n + 1) . (2)
n(n + 1)

Fig. 1 shows n (t) for n = 5 as well as the associated point set A5 . Note that despite the fact
that we use n(n+ 1) + 1 = 31 parameter values to generate A5 one may count that there are
5+2
precisely 21 = = dim(
5 (K)) different points in A5 , i.e., exactly the correct number
2
for polynomial interpolation.
Notice that each of the interior points lies on a self-intersection point of the curve. Indeed, it is
easy to see that
   
j n + m(n + 1) j n − m(n + 1)
n  = n 
n(n + 1) n(n + 1)

for any integers j, m. Hence we may naturally describe the set of distinct points An as
   
j n + m(n + 1)
An = Aj m := n  : j, m 0 and j + m
n . (3)
n(n + 1)

More specifically, we may classify the points of An by

18 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25

Interior points:
      
j n + m(n + 1) jn m(n + 1)
A j m = n  = (−1)m cos  , (−1)j cos  ,
n(n + 1) n+1 n
j, m > 0 and j + m
n.

Vertex points:

A00 = (1, 1) and A0,n = ((−1)n , (−1)(n+1) ).

Vertical edge points:

  
m(n + 1)
A0m = (−1)m , cos  , m = 1, 2, . . . , n.
n

Horizontal edge points:

   
jn
Aj 0 = cos  , (−1)j , j = 1, 2, . . . , n.
n+1
 
n+2
It follows then that the number of distinct points in An is indeed .
2

3. The fundamental Lagrange polynomials

Consider the weighted integral


1 1 1
I (f ) := 2 f (x, y) √  dx dy
 [−1,1]2 1 − x2 1 − y2

and the associated inner product


1 1 1
f, g := 2 f (x, y)g(x, y) √  dx dy. (4)
 [−1,1]2 1 − x 2 1 − y2

Note that the polynomials {Tj (x)Tm (y) : j + m
n} form an orthogonal basis for
n (R2 )
with respect to this inner product. Further, the finite dimensional space
n (R2 ) has a reproducing
kernel with respect to this inner product, which we denote by Kn (·; ·). Specifically, for each
A, B ∈ [−1, 1]2 , Kn (A; ·) ∈
n (R2 ), Kn (·; B) ∈
n (R2 ) and

p, Kn (A; ·) = p(A)

for all p ∈
n (R2 ).

We will show that there is a remarkable quadrature formula for I (f ) based on the Padua points
(Theorem 1).
Because of this quadrature formula, it is perhaps not unexpected that there are formulas for
the fundamental Lagrange polynomials of the Padua points in terms of Kn . We begin with the
following two lemmas.
 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25 19

Lemma 1. For all p ∈
2n (R2 ) we have

 
1 1 1 1 
p(x, y) √  dx dy = p(n (t)) dt.
2 [−1,1]2 1 − x2 1 − y2  0

Proof. We need only show that this holds for

p(x, y) = Tj (x)Tm (y), j + m
2n.

But, when j = m = 0, p(x, y) = 1, and clearly both sides are equal to 1. Otherwise, if
(j, m) = (0, 0), the left side equals 0 while the right side equals
 
1  1 
Tj (cos(nt))Tm (cos((n + 1)t)) dt = cos(j nt) cos(m(n + 1)t) dt
 0  0
=0
provided j n = m(n + 1). But since n and n + 1 are relatively prime, j n = m(n + 1) implies that
j = (n + 1) and m = n for some positive integer . Hence j + m = (n + (n + 1)) > 2n. 

Theorem 1. Consider the following weights associated to the Padua points A ∈ An :

⎧
1 ⎨ 1/2 if A ∈ An is a vertex point,
wA := 1 if A ∈ An is an edge point,
n(n + 1) ⎩ 2 if A ∈ An is an interior point.

Then, if p(x, y) ∈
2n (R2 ) is such that p(x, y), T2n (y) = 0, we have


1 1 1
p(x, y) √  dx dy = wA p(A).
 2
[−1,1]2 1 − x2 1 − y2 A∈An

In particular, since T2n (y) = 2(Tn (y))2 − 1, this quadrature formula holds for p = f g with
f, g ∈
n (R2 ) and either f (x, y), Tn (y) = 0 or g(x, y), Tn (y) = 0.

Proof. First, note that the quadrature formula

1
 
1 1
k  1
N−1
q(t) dt = q(0) + q  + q() (5)
 0 N 2 N 2
k=1

holds for all even trigonometric polynomials q(t) (i.e., has an expansion in cosines only) for which
the expansion of q(t) does not contain terms with frequency an even multiple of N , i.e., such that
 
q(t) cos(2Nt) dt = 0,  = 1, 2, . . . .
0
Now, by Lemma 1, we have
 
1 1 1 1 
p(x, y) √  dx dy = p(n (t)) dt
2 [−1,1]2 1 − x2 1 − y2  0

where

p(n (t)) = p(cos(nt), cos((n + 1)t)) =: q(t)

20 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25

is an even trigonometric polynomial (of degree at most n(n + 1)) to which we may apply the
quadrature
 formula (5) with N = n(n + 1) provided, under the assumptions placed on p, we have

q(t) cos(2N t) dt = 0. However, if we write
0

2n 2n−j
p(x, y) =: aj m Tj (x)Tm (y)
j =0 m=0

then



2n 2n−j
p(n (t)) = aj m cos(j nt) cos(m(n + 1)t)
j =0 m=0

aj m
2n 2n−j
= {cos((j n + m(n + 1))t) + cos((j n − m(n + 1))t)}
2
j =0 m=0

so that a0,2n , the coefficient of cos(2n(n + 1)t) in the expansion of p(n (t)), is exactly the
coefficient of T2n (y) in the expansion of p(x, y). By assumption this is zero. Hence we apply
Lemma 1 and (5) (with N = n(n + 1)) to obtain

1 1 1
p(x, y) √  dx dy
2
[−1,1]2 1−x 2 1 − y2

1 
= p(n (t)) dt
 0
1 1
  k  1
N−1
= p(n (0)) + p n  + p(n ())
n(n + 1) 2 N 2
k=1


= wA p(A)
A∈An

since n (0) and n () are the two vertex points and the interior points double up. 

The quadrature just proved has precision 2n − 1 as it is exact for all polynomials of degree
< 2n. From the general theory on quadrature formulas and polynomial ideals, it follows that the
set of Padua points is a variety of a polynomial ideal generated by quasi-orthogonal polynomials
with respect to ·, · , from which a compact formula for the Lagrange interpolation polynomial
can be derived. This approach based on ideal theory will be developed fully in [3]. Here, we give
a proof of the compact interpolation formula.
For A ∈ An we will let yA denote the y-coordinate of A.

Theorem 2. For B ∈ An , set

LB (x, y) := wB [Kn (B; (x, y)) − Tn (y)Tn (yB )], B ∈ An . (6)

Then we have the interpolation formula

LB (A)p(B) = p(A) ∀p ∈
n (R2 ) ∀A ∈ An . (7)
B∈An
 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25 21

Proof. Since, as it is easy to see, Tn (y), Tn (y) = 1

2 and Kn (·; A), Tn (y) = Tn (yA ), we have
Kn (·; A) − 2Tn (yA )Tn (y), Tn (y) = Tn (yA ) − 2Tn (yA )Tn (y), Tn (y)
= Tn (yA ) − 2Tn (yA ) 21
= 0.
Hence, we may apply the quadrature formula of Lemma 1 to (Kn (·; A) − 2Tn (yA )Tn (y)) × p for
any p ∈
n (R2 ) and obtain
p(A) − 2Tn (yA )Tn (y), p
= Kn (·; A) − 2Tn (yA )Tn (y), p


= wB (Kn (B; A) − 2Tn (yA )Tn (yB ))p(B)
B∈An



= wB (Kn (B; A) − Tn (yA )Tn (yB ))p(B) − Tn (yA ) wB Tn (yB )p(B).
B∈An B∈An

It follows that


wB (Kn (B; A) − Tn (yA )Tn (yB ))p(B)
B∈An
⎧ ⎫
⎨
⎬
= p(A) − Tn (yA ) 2Tn (y), p − wB Tn (yB )p(B) . (8)
⎩ ⎭
B∈An

But, writing p = p̃ + Tn (y) where  := 2p, Tn (y) and p̃, Tn (y) = 0, we have



wB Tn (yB )p(B) = wB Tn (yB )[p̃(B) + Tn (yB )]
B∈An B∈An


= Tn (y), p̃ +  wB Tn2 (yB )
B∈An


=0+ wB × 1
B∈An
=×1
= 2p, Tn (y) ,
where we have used the fact that Tn2 (yB ) = 1 for all B ∈ An and that the sum of the weights wB
is one.
Therefore, by (8), we have


wB (Kn (B; A) − Tn (yA )Tn (yB ))p(B) = p(A)
B∈An

for all p ∈
n (R2 ) and all points A ∈ An , which completes the proof. 

The polynomials LB are indeed the fundamental Lagrange polynomials, i.e., they satisfy

LB (A) =
A,B , A, B ∈ An . (9)

This can be easily verified using the following compact formula for the reproducing kernel Kn
proved in [7].
22 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25

Lemma 2. For A, B ∈ [−1, 1]2 write A = (cos(1 ), cos(2 )), B = (cos(1 ), cos(2 )). Then
Kn (A; B) = Dn (1 + 1 , 2 + 2 ) + Dn (1 + 1 , 2 − 2 )
+Dn (1 − 1 , 2 + 2 ) + Dn (1 − 1 , 2 − 2 )
where
1 cos((n + 1/2)) cos(/2) − cos((n + 1/2)) cos(/2)
Dn (, ) = .
2 cos() − cos()

We remark that (9) does not follow from the interpolation formula (7) (as a simple example,
consider the formula (f (A) + f (B))/2 = f (C) valid for all f ∈
0 and arbitrary A, B, C). A
direct derivation of the compact formulas in (6) and (7) will be given in [3].

4. The growth of the Lebesgue function

We first give a more convenient expression for Dn , which is obtained by simple trigonometric
manipulations.

Lemma 3 (cf. Lemma 1 of [2]). The function Dn (, ) can be written as

⎧    
⎪ + −
⎪ sin n
1⎨
sin n
2 2
Dn (, ) =    
4⎪⎪  +   − 
⎩ sin sin
2 2
   ⎫
+ − ⎪
sin(n + 1) sin(n + 1) ⎪
⎬
2 2
+     . (10)
+ − ⎪
⎪
sin sin ⎭
2 2

Now, note that the norm of the polynomial projection

Lf := f (A)LA
A∈An

considered as a map from C(K) to C(K), equipped with the uniform norm, is given by the
maximum of the so-called Lebesgue function


n (x, y) := |LA (x, y)|. (11)
A∈An

Specifically,
L = max n (x, y).
(x,y)∈[−1,1]2

We now proceed to give a bound on n (x, y).

Theorem 3. There is a constant C > 0 such that the Lebesgue function is bounded,
n (x, y)
C(log n)2 , n2, (x, y) ∈ [−1, 1]2 .
 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25 23

Proof. First, recall that the Padua points are given by

 
j n + m(n + 1)
Aj m = n  , j, m0, j + m
n.
n(n + 1)
Some simple manipulations show that, in fact,
    m 
j
Aj m = (−1)j +m cos  , cos  . (12)
n+1 n
Now,


n (x, y) = |LA (x, y)|
A∈An


= wA |Kn (A; (x, y)) − Tn (y)Tn (yA )|
A∈An




wA |Kn (A; (x, y))| + wA |Tn (y)| |Tn (yA )|
A∈An A∈An



wA |Kn (A; (x, y))| + 1 (13)
A∈An

since |Tn (y)|
1 and the sum of the weights wA is 1.

To bound this latter sum (13) we use the compact formula for Kn given in Lemma 2 and the
alternate expression for Dn given in Lemma 3. Indeed it is clear that our result will follow if we
find an upper bound for the sum
    
  + A A − A 
 sin n A sin n 

 2 2 

wA       , (14)
A∈An  sin A + A A − A 
 sin 
2 2
where we write (x, y) =: (cos(1 ), cos(2 )) and
j m
A = Aj m , A :=  + 1 , A :=  + 2 .
n+1 n
The analysis of the other terms that result from expanding Kn by means of Lemmas 2 and 3 is
entirely analogous. We should note, however, that the (−1)j +m factor in (12) has no effect on the
sum (14) due to the fact that | sin(k( + ))| = | sin(k)| for any integer k.
We now proceed to bound (14). Using the definitions of wA , A and A , we have
    
 A + A A −  A 
 sin n sin n 

 2 2 
wA      
A∈An  sin A + A A − A 
 sin 
2 2
  
 1 + 2 j  m 
 sin n + + 
2


n n−j
2 n+12 n 2 

  
n(n + 1)  1 + 2 j  m  
j =0 m=0  sin + + 
2 n+12 n2
24 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25

  
 1 − 2 j  m 
 sin n + − 
 2 n+12 n 2 
×   
 sin 1 − 2 + j  − m  
 
2 n+12 n2
  
 1 + 2 j  m 

n
n  sin n + + 
2  2 n+12 n 2 

  
n(n + 1)  1 + 2 j  m  
j =0 m=0  sin + + 
2 n+12 n2
  
 1 − 2 j  m 
 sin n + − 
 2 n+12 n 2 
×   (15)
.
 sin 1 − 2 + j  − m  
 
2 n+12 n2
Now this last expression (15) is a Riemann sum for a multiple of the (improper) integral
     
 1 + 2    1 − 2 
 /2  /2  sin n ++   sin n +− 
 2   2 
        d d. (16)
 1 + 2 
0 0 
+  +    sin 1 − 2 +  −  
 sin  
2 2

The integrand of (16) is bounded by

1 1
      .
   
sin 1 + 2 +  +   sin 1 − 2 +  −  
 2   2 

Note that the denominator of this upper bound is zero at most along two lines in the domain
(, ) ∈ [0, /2]2 , which corresponds to at most Cn points in the corresponding sum (15), for some
constant C > 0. At each of these “singularities”, we may use the fact that | sin(n)/ sin()|
n to
reduce that portion of the sum (15) that involves the singular points to a univariate sum similar to
  
 k 
n  sin n

+  
2  n 
   
n+1  k

k=0  sin  +  

n

which a classical analysis (cf. Lemma 4 of [2]) shows to be bounded by C log n for some constant
C > 0.
Thus, by the monotonicity and periodicity of csc(), it follows that the sum (15) is bounded by
a constant multiple of
 /2  /2
| csc() csc()| d d. (17)
/(n+1) /(n+1)

This in turn is easily evaluated and seen to be bounded by C(log n)2 . The result follows. 
 L. Bos et al. / Journal of Approximation Theory 143 (2006) 15 – 25 25

Acknowledgment

Work supported by the ex-60% funds of the Universities of Padua and Verona, and by the
GNCS-INdAM.

References

[1] L. Bos, M. Caliari, S. De Marchi, M. Vianello, A numerical study of the Xu polynomial interpolation formula in two
variables, Computing 76 (3–4) (2005) 311–324.
[2] L. Bos, S. De Marchi, M. Vianello, On the Lebesgue constant for the Xu interpolation formula, J. of Approx. Theory,
to appear.
[3] L. Bos, S. De Marchi, M. Vianello, Y. Xu, Bivariate Lagrange interpolation at the Padua points: the ideal theory
approach, submitted.
[4] L. Brutman, Lebesgue functions for polynomial interpolation—a survey, Ann. Numer. Math. 4 (1997) 111–127.
[5] M. Caliari, S. De Marchi, M. Vianello, Bivariate polynomial interpolation on the square at new nodal sets, Appl. Math.
Comput. 165 (2005) 261–274.
[6] B. Sündermann, On projection constants of polynomials spaces on the unit ball in several variables, Math. Z. 188
(1984) 111–117.
[7] Y. Xu, Christoffel functions and Fourier series for multivariate orthogonal polynomials, J. Approx. Theory 82 (1995)
205–239.
[8] Y. Xu, Lagrange interpolation on Chebyshev points of two variables, J. Approx. Theory 87 (1996) 220–238.