2. The Bernstein frame

10 December 2013

Multivariate Bernstein operators and redundant systems

Tan Duc Do and Shayne Waldron

Department of Mathematics, University of Auckland, Private Bag 92019, Auckland, New Zealand e–mail: [email protected] (http://www.math.auckland.ac.nz/˜waldron)

ABSTRACT

The Bernstein operator B_n for a simplex in IR^d is naturally defined via the Bernstein basis obtained from the barycentric coordinates given by its vertices. Here we consider a generalisation of this basis and the Bernstein operator, which is obtained from generalised barycentric coordinates that are given for more general configurations of points in IR^d. We call the associated polynomials a Bernstein frame, as they span the polynomials of degree ≤ n, but may not be a basis. By using this redundant system we are able to give geometrically motivated proofs of some basic properties of the corresponding generalised Bernstein operator, such as the fact it is degree reducing and converges for all polynomials.

We also consider the conditions for polynomials in this Bernstein form to join smoothly.

Key Words: multivariate Bernstein operator, barycentric coordinates, frames, redundant system Stirling numbers, blossoming, de Casteljau algorithm

AMS (MOS) Subject Classifications: primary 41A10, 41A36 65D17, secondary 15A18, 42C15,

1. Introduction

The Bernstein operator [L53] and its variants [AC94] have important applications, e.g., B´ezier curves and surfaces [F02]. Here we consider a generalisation of this operator, which is based on a redundant “Bernstein basis”.

The Bernstein operator Bn for a simplex in IR^d is defined via the Bernstein basis for Πn(IR^d) (thed–variate polynomials of degree≤n). This basis is obtained by taking powers of the barycentric coordinates given by the vertices of the simplex. In the next section, we outline the basic properties of generalised barycentric coordinates which are given by more general configurations of points in IR^d, e.g., the vertices of a convex polygon. These lead naturally to an analogue of the Bernstein basis, a set of polynomials of degree n which span Πn(IR^d). These are not a basis if they are given by more than d+ 1 points, and so we refer to this possibly redundant system as a Bernstein frame (cf. [C03]).

In Section 3, we define the generalised Bernstein operator given by a Bernstein frame.

We give geometrically motivated proofs of some basic properties of it. These include show- ing that it is degree reducing and converges for all polynomials, that it reproduces the linear polynomials, and more generally has the same spectral structure as the classical Bernstein operator. Similar arguments in terms of a basis would be far more cumbersome. Finally, we explore some applications of our generalised Bernstein operator. These include a de Casteljau algorithm, shape preservation properties (Section 4), and smoothness conditions in terms of the control points of the associated B´ezier surfaces (Section 5).

2. The Bernstein frame

Let V consist of d+ 1 affinely independent points in IR^d, i.e., be the vertices of a d–

simplex. Thebarycentric coordinates(cf. [B87],[LS07]) of a pointx∈IR^d with respect to V are the unique coefficients (ξv(x))v∈V ∈ IR^V for which x can be written as an affine combination of the points in V, i.e.,

x= X

v∈V

ξ_v(x)v, X

v∈V

ξ_v(x) = 1. (2.1)

We follow [B87] and index the barycentric coordinates by the points v ∈ V that they correspond to, and use standard multiindex notation. It follows, from (2.1), that the ξv

are linear polynomials which are a basis for Π1(IR^d). More generally, for any n ≥ 1, the polynomials

B_α :=

|α|

α

ξ^α, |α|=n (α ∈ZZ^V₊) are a basis for Πn(IR^d). Here |α|=P

vαv, ⁿ_α

= ^n!_α!, and ξ^α =Q

vξ_v^αv.

From now on, letV be a sequence (or multiset) ofm=|V|points with affine hull IR^d, so that each point x∈IR^d can be written as an affine combination

x= X

v∈V

avv, X

v∈V

av = 1, (2.2)

where the coefficients a = (av)v∈V are unique if and only if V consists of d+ 1 points.

Following [W11₂], we call the unique minimal ℓ₂–norm coefficientsa ∈IR^V satisfying (2.2) the (generalised barycentric)coordinatesgiven byV, and denote them by (ξv(x))v∈V. By construction, they satisfy (2.1). Each ξv is a linear polynomial, and they span Π₁(IR^d), since (2.1) gives

1 = X

v∈V

ξ_v(x), x_j = X

v∈V

ξ_v(x)vj, j = 1, . . . , d, (2.3)

which is equivalent to the following reproduction formula for affine functions f = X

v∈V

f(v)ξv, ∀f ∈Π1(IR^d). (2.4) From the formula for ξv given in [W11₂] it is easy to see:

• The coordinates of the barycentre c:= _m¹ P

v∈V v of V are ξ_v(c) = _m¹ , ∀v.

• ξ_v is constant (equal to _m¹) if and only if v is the barycentre c.

• ξ_v =ξ_w if and only if v=w.

• The ξ_v are continuous functions of the points v ∈V (with affine hull IR^d).

These imply that the set of points where the coordinates are nonnegative

NV :={x∈IR^d :ξv(x)≥0,∀v∈V} (2.5) is a convex polytope, with the barycentre ofV as an interior point. We callN_V theregion of nonnegativity for the coordinates given by V (See Fig. 1).

Fig 1. The region of nonnegativityN_V for V given by the vertices of a triangle, square and pentagon.

We write V \w for the sequence (or multiset) obtained by removing the point w from V (once), and Aff(V) for the affine hull of the points inV. We recall:

Proposition 2.6 ([W11₂]). The generalised barycentric coordinates satisfy (a) _m¹ < ξ_v(v)≤1.

(b) ξ_v(w) =ξ_w(v).

(c) ξv(v) = 1if and only if v6∈Aff(V \v), in which case ξv = 0 on Aff(V \v).

(d) P

vξ_v(v) =d+ 1.

Expanding the monomial basis for Πn(IR^d) in terms of ξ using (2.3), shows that the polynomials

Bα :=

|α|

α

ξ^α, |α|=n (2.7)

span Πn(IR^d). By a dimension count, these ^n+m−1_m−1

polynomials are basis if and only if V consists of d+ 1 affinely independent points. Thus, we refer to {Bα : |α| = n} as the Bernstein frame given by the points V. This is a partition of unity, since applying the multinomial theorem to (2.1) gives

X

|α|=n

Bα = X

|α|=n

n α

ξ^α = X

v∈V

ξv

n

= 1. (2.8)

A Bernstein frame is nonnegative onN_V, the region of nonnegativity given by (2.5). There have been studies of the approximation properties of linear operators given by partitions of unity which may take negative values on the region of interest, see, e.g., [BD85].

Example 2.9. (See Fig. 2) For V ={0,¹₂,1} ⊂IR the generalised barycentric coordinates are

ξ₀(x) = 5

6 −x, ξ¹

2(x) = 1

3, ξ₁(x) =x− 1 6.

Here, some polynomials in the Bernstein frame have degree < n. This is the case if and only if the barycentre of V is a point of V. The coordinates for V ={0,¹₃,²₃,1} are

ξ₀(x) = 7 10 − 9

10x, ξ¹

3(x) = 2 5 − 3

10x, ξ²

3(x) = 3

10x+ 1

10, ξ₁(x) = 9 10x− 1

5.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−0.2

−0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−0.2

−0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

Fig. 2. The Bernstein frame for Π₂(IR) for V ={0,¹₂,1} and V ={0,¹₃,²₃,1}.

We say a polynomial p∈Πn(IR^d) is in (Bernstein-B´ezier) B-formif

p= X

|α|=n

c_αB_α, (2.10)

The mesh function c : α 7→ c_α is unique if and only if V consists of d+ 1 points. The mesh function with minimalℓ₂–norm gives a canonical B-form, i.e., what[W11₁] calls the canonical coordinates of p with respect to {Bα}_|α|=n.

Many familiar formulas for the Bernstein basis extend to a Bernstein frame. Here are a couple of examples (also see Section 4). Let e_v be the multiindex given by

e_v(w) :=

1, v=w;

0, and define B_α := 0 if α6≥0.

Proposition 2.11. The Bernstein frame {Bα}_|α|=n can be calculated recursively via Bα = X

v∈V

ξvB_α−ev, B₀ = 1, (2.12)

and expressed in terms of the Bernstein frame for polynomials of degree n+ 1 via Bα = X

v∈V

α_v+ 1

|α|+ 1B_α+ev. (2.13)

Proof: We calculate X

v∈V

ξ_vB_α−ev = X

v∈V

|α| −1 α−e_v

ξ^α = X

v∈V

α_v

|α|

|α|

α

ξ^α =B_α, and, using P

vξ_v = 1, that Bα =Bα

X

v∈V

ξv = X

v∈V

|α|!

α! ξ^αξv = X

v∈V

αv+ 1

|α|+ 1

|α+ev|!

(α+e_v)!ξ^α+ev = X

v∈V

αv + 1

|α|+ 1B_α+ev.

Proposition 2.14 (Differentiation). For u, v, w ∈V, we have D_v−wξu =ξu(v)−ξu(w) =ξv(u)−ξw(u).

Thus the Bernstein frame satisfies

D_v−wB_α =|α|X

u∈V

ξ_u(v)−ξ_u(w)

B_α−eu. Proof: Sinceξ_u is affine

(D_v−wξu)(x) = lim

t→0

ξu(x+t(v−w))−ξu(x) t

= lim

t→0

ξ_u(x) +tξ_u(v)−tξ_u(w)−ξ_u(x)

t =ξ_u(v)−ξ_u(w).

By the product and chain rules, we have D_v−wBα = |α|!

α! D_v−w Y

u∈V

ξ_u^αuξ^α−αueu = |α|!

α!

X

u∈V

αuξ_u^αu−1 ξu(v)−ξu(w)

ξ^α−αueu

=|α|X

u∈V

ξu(v)−ξu(w)(|α| −1)!

(α−e_u)!ξ^α−eu =|α|X

u∈V

ξu(v)−ξu(w)

B_α−eu.

3. The generalised Bernstein operator

For a Bernstein frame (2.7) given by points V in IR^d, we define a (generalised) Bernstein operator Bn =Bn,V of degree n≥1 by the usual formula

B_n(f) := X

|α|=n

B_αf(vα), v_α := X

v∈V

α_v

|α|v, (3.1)

which is equivalent to

B_n(f) = X

v1∈V

· · · X

vn∈V

fv₁+· · ·+v_n n

ξ_v1· · ·ξ_vn.

Fig. 3. The points {vα}_|α|=n used in the definition of Bn,Vf, where n= 7 and V is the vertices of a triangle, square, pentagon and hexagon (respectively).

This maps functionsf which are nonnegative at the points (vα)_|α|=k (which are contained in T = conv(V), the convex hull of the points in V) to polynomials of degree ≤n which are nonnegative on the convex polytope (region of nonnegativity) N_V given by (2.5), so that

f ≥0 on T =⇒ B_nf ≥0 on N_V, (3.2)

and reproduces the linear polynomials (cf. Theorem 3.19).

We now show that the generalised Bernstein operator is degree reducing, i.e., B_n(f)∈Πk(IR^d), ∀f ∈Πk (k = 0,1, . . .).

Define the univariate and multivariate (falling) shifted factorials by [x]ⁿ:=x(x−1)· · ·(x−n+ 1), [α]^β := Y

v∈V

[αv]^βv, and the multivariate Stirling numbers of the second kind by

S(τ, β) := Y

v∈V

S(τv, βv),

where S(τv, β_v) are the Stirling numbers of the second kind. We note that

S(τ, β) = 0, β 6≤τ, (3.3)

and define

|α|

α

:= 0, α 6≥0. (3.4)

These are related by

α^τ = X

β≤τ

S(τ, β)[α]^β. (3.5)

Lemma 3.6. For any τ and n, we have X

|α|=n

α^τ |α|

α

ξ^α = X

β≤τ

S(τ, β)[n]^|β|ξ^β. (3.7)

Proof: Since [|α|]^{|β| |α−β|}_α−β

= ^|α|_α

[α]^β (without restriction on α and β), (3.5) gives

X

β≤τ

S(τ, β)[|α|]^|β|

|α−β|

α−β

= |α|

α X

β≤τ

S(τ, β)[α]^β = |α|

α

α^τ. Thus, we calculate

X

|α|=n

α^τ |α|

α

ξ^α = X

|α|=n

X

β≤τ

S(τ, β)[|α|]^|β|

|α−β| α−β

ξ^α

= X

β≤τ

S(τ, β)[n]^|β|ξ^β X

|α|=n α≥β

|α−β|

α−β

ξ^α−β = X

β≤τ

S(τ, β)[n]^|β|ξ^β,

with the last equality given by the multinomial identity.

Theorem 3.8 (Degree reducing). The generalised Bernstein operator B_n is degree reducing. More precisely,

B_n(ξ^β) = [n]^|β|

n^|β| ξ^β + X

0<|γ|<|β|

[n]^|γ|

n^|β| a(γ, β)ξ^γ, (3.9) where w₁, . . . , wm is the sequence of points in V, and

a(γ, β) := X

|τ1|=βw1

· · · X

|τm|=β_wm

βw1

τ₁

ξ^τ1(w₁)· · · βwm

τ_m

ξ^τm(w_m)S(τ₁+· · ·+τ_m, γ).

(3.10) Proof: Since each ξw is an affine function, and ξw(v) =ξv(w), we have

ξw(vα) =ξw

X

v∈V

α_v

|α|v

= X

v∈V

α_v

|α|ξw(v) = X

v∈V

α_v

|α|ξv(w),

and the multinomial identity gives (ξ^β)(v_α) = Y

w∈V

X

v∈V

αv

|α|ξ_v(w)βw

= Y

w∈V

X

|τ|=βw

βw

τ

α^τ

|α|^βwξ^τ(w)

= X

|τ1|=βw1

· · · X

|τm|=βwm

βw1

τ₁

· · · βwm

τ_m

ξ^τ1(w₁)· · ·ξ^τm(wm)α^{τ1+···+τm}

|α|^|β| . Thus, by rearranging (3.1) and Lemma 3.6, we have

B_n(ξ^β) = X

|τ1|=βw1

· · · X

|τm|=β_wm

β_w1 τ₁

ξ^τ1(w1)· · · β_wm

τm

ξ^τm(wm) X

|α|=n

α^{τ1+···+τm} n^|β|

n α

ξ^α

= X

|τ1|=βw1

· · · X

|τm|=β_wm

β_w1 τ₁

ξ^τ1(w₁)· · · β_wm

τ_m

ξ^τm(w_m)

× X

γ≤τ1+···+τm

[n]^|γ|

n^|β| S(τ₁+· · ·+τm, γ)ξ^γ.

Here B_n(ξ^β) is written as a polynomial in ξ of degree≤ |β|, so thatB_n is degree reducing.

The terms of degree |β| can be simplified using the multinomial identity, ξv(wj) =ξwj(v), and (2.4), as follows

X

|τ1|=βw1

· · · X

|τm|=β_wm

βw1

τ₁

ξ^τ1(w₁)· · · βwm

τ_m

ξ^τm(w_m)[n]^|β|

n^|β| ξ^{τ1+···+τm}

= [n]^|β|

n^|β|

Ym

j=1

X

|τj|=β_wj

β_wj τj

ξ^τj(wj)ξ^τj

= [n]^|β|

n^|β|

Ym

j=1

X

v∈V

ξ_v(wj)ξv

β_wj

= [n]^|β|

n^|β|

Ym

j=1

X

v∈V

ξwj(v)ξv

β_wj

= [n]^|β|

n^|β|

Ym

j=1

ξw^βwjj = [n]^|β|

n^|β| ξ^β.

By collecting the terms of degree < |β|, we obtain (3.9). Here, (3.3) allows us to remove the restrictionγ ≤τ₁+· · ·+τ_m, and there are no terms of degree 0 since S(1,0) = 0.

Remark 3.11. If V consists of d+ 1 affinely independent points, then a(γ, β) =

S(β, γ), γ ≤β;

0, γ 6≤β,

and (3.9) simplifies to

Bn(ξ^β) = [n]^|β|

n^|β| ξ^β+ X

γ<β

[n]^|γ|

n^|β| S(β, γ)ξ^γ.

This was proved in [CW02](Lemma 2.1) for the case when β_v0 = 0 for some v₀ ∈V. Example 3.12 (Linear reproduction). For |β|= 1, we have

Bn(ξv) =ξv, ∀v∈V, (3.13)

i.e., B_n reproduces the linear polynomials Π1(IR^d) = span{ξv}v∈V. This is equivalent to

x = X

|α|=n

B_α(x)v_α, X

|α|=n

B_α(x) = 1, x ∈IR^d. (3.14)

Example 3.15 (Quadratics). For |β|= 2, we recall S(1,0) = 0,S(2,1) = 1, so that a(eu,2ew) =ξ_u²(v) =ξ²_v(u), a(eu, e_v +e_w) =ξ_u(v)ξu(w) = (ξvξ_w)(u), v 6=w, (3.16) and we obtain

B_n(ξ^β) = 1− 1

n

ξ^β + 1 n

X

u∈V

ξ^β(u)ξ_u →ξ^β, as n→ ∞, |β|= 2. (3.17)

SinceB_n is not a positive operator in general, the application of the Korovkin theory is more involved (see Theorem 3.31). We easily obtain the following encouraging result.

Corollary 3.18 (Convergence). For all polynomialsf, Bn(f)→f, as n→ ∞.

Proof: It suffices to consider f =ξ^β. For n≥ |β|, (3.9) gives B_n(ξ^β)−ξ^β =[n]^|β|

n^|β| −1

ξ^β − X

|γ|<|β|

[n]^|γ|

n^|β| a(γ, β)ξ^γ, where ^[n]_n|β|^|β| −1,^[n]_n|β|^|γ| =O(¹_n), as n→ ∞.

The remaining eigenstructure of B_n is as follows.

Theorem 3.19 (Diagonalisation). The generalised Bernstein operator B_n given by points V ⊂IR^d is diagonalisable, with eigenvalues

λ⁽ⁿ⁾_k := [n]^k

n^k , k= 1, . . . , n, 1 =λ⁽ⁿ⁾₁ > λ⁽ⁿ⁾₂ >· · ·> λ⁽ⁿ⁾_n >0.

Let P_k,V⁽ⁿ⁾ denote the λ⁽ⁿ⁾_k –eigenspace. Then

P_1,V⁽ⁿ⁾ = Π1(IR^d), ∀n. (3.20)

For k >1, P_k,V⁽ⁿ⁾ consists of polynomials of exact degree k, and is spanned by p⁽ⁿ⁾_ξβ =ξ^β+ X

0<|α|<|β|

c(α, β, n)ξ^α, |β|=k, (3.21)

where the coefficients can be calculated using (3.10) and the recurrence formula c(α, β, n) := a(α, β)

1− |β|, |α|=|β| −1, c(α, β, n) := [n]^|α|

λ⁽ⁿ⁾_|β| −λ⁽ⁿ⁾_|α|

a(α, β)

n^|β| + X

|α|<|γ|<|β|

c(γ, β, n)a(α, γ) n^|γ|

, |α|<|β| −1.

(3.22) Proof: By Example 3.12, theλ⁽ⁿ⁾₁ = 1 eigenspace P_1,V⁽ⁿ⁾ contains Π₁(IR^d). Recall, from (3.9), that B_n(ξ^β) has the form

B_n(ξ^β) =λ⁽ⁿ⁾_|β|ξ^β+ X

0<|γ|<|β|

[n]^|γ|

n^|β| a(γ, β)ξ^γ. (3.23) Motivated by this, we seek λ⁽ⁿ⁾_k –eigenfunctions of the form

f =ξ^β+ X

0<|α|<|β|

c(α, β, n)ξ^α, |β|=k >1.

We observe that for such an eigenfunction the coefficients c(α, β, n) are not unique– even when V consists of d+ 1 points. Expanding B_n(f) =λ⁽ⁿ⁾_k f using (3.23) gives

Bn(f) =Bn(ξ^β) + X

0<|γ|<|β|

c(γ, β, n)Bn(ξ^γ)

=λ⁽ⁿ⁾_|β|ξ^β+ X

0<|α|<|β|

[n]^|α|

n^|β| a(α, β)ξ^α + X

0<|γ|<|β|

c(γ, β, n)

λ⁽ⁿ⁾_|γ|ξ^γ+ X

0<|α|<|γ|

[n]^|α|

n^|γ| a(α, γ)ξ^α

=λ⁽ⁿ⁾_|β|ξ^β+ X

0<|α|<|β|

λ⁽ⁿ⁾_|β|c(α, β, n)ξ^α.

(3.24) Equating the ξ^α, 0<|α|<|β| coefficients gives

λ⁽ⁿ⁾_|β|c(α, β, n) = [n]^|α|

n^|β| a(α, β) +c(α, β, n)λ⁽ⁿ⁾_|α| + X

|α|<|γ|<|β|

c(γ, β, n)[n]^|α|

n^|γ| a(α, γ). (3.25) Since λ⁽ⁿ⁾_|α| > λ⁽ⁿ⁾_|β|, this is equivalent to

c(α, β, n) = 1 λ⁽ⁿ⁾_|β| −λ⁽ⁿ⁾_|α|

[n]^|α|

n^|β| a(α, β) + X

|α|<|γ|<|β|

c(γ, β, n)[n]^|α|

n^|γ| a(α, γ) . From this we can define suitable c(α, β, n) recursively, as in (3.22), starting from

c(α, β, n) := 1 λ⁽ⁿ⁾_|β| −λ⁽ⁿ⁾_|α|

[n]^|α|

n^|β| a(α, β) = a(α, β)

1− |β|, |α|=|β| −1.

A simple dimension count shows that the eigenfunction {p_ξ^β}_|β|=k, so defined, span a space P_k,V⁽ⁿ⁾ of dimension ^k+d−1_d−1

. Again, by dimension counting, we conclude that B_n is diagonalisable, withP_1,V⁽ⁿ⁾ = Π₁(IR^d).

Example 3.26 (Quadratic eigenfunctions). Using (3.16), we have p⁽ⁿ⁾_ξβ =ξ^β −X

u∈V

ξ^β(u)ξu, |β|= 2. (3.27)

In general, p⁽ⁿ⁾_ξα does depend on n (cf. [CW00]).

Despite the fact the coefficients in (3.21) are not unique, we can take their limit as n→ ∞. This indicates that the redundant expansion (3.21) is natural.

Corollary 3.28 (Limits of the eigenfunctions). For 0<|α|<|β|,

n→∞lim c(α, β, n) =c^∗(α, β), where

c^∗(α, β) := a(α, β)

1− |β|, |α|=|β| −1,

c^∗(α, β) := 2

(|β| − |α|)(−|α| − |β|+ 1) X

|γ|=|α|+1

c^∗(γ, β)a(α, γ), |α|<|β| −1. (3.29)

Thus, the eigenfunctions of (3.21) satisfy p⁽ⁿ⁾_ξβ →p^∗_ξβ :=ξ^β + X

0<|α|<|β|

c^∗(α, β)ξ^α, as n→ ∞. (3.30)

Proof: Fix β. We use strong induction on j = |β| − |α| = 1, . . . ,|β| to prove the limit exists. For |α| = |β| −1 (j = 1) the limit is clear. Suppose the limit of c(γ, β, n) exists for all γ with |α|<|γ|<|β|. Then taking the limit of (3.22) gives

n→∞lim c(α, β, n) = 2

(|β| − |α|)(|β| − |α| −2|β|+ 1) X

|γ|=|α|+1

c^∗(γ, β)a(α, γ).

This follows from the calculations λ⁽ⁿ⁾_|β| −λ⁽ⁿ⁾_|α| = [n]^|α|[n− |α|]^|β|−|α|

n^|β| − [n]^|α|

n^|α| = [n]^|α|

n^|β| ([n− |α|]^|β|−|α|−n^|β|−|α|), [n− |α|]^|β|−|α|−n^|β|−|α| = 1

2(|β| − |α|)(−|α| − |β|+ 1)n^{|β|−|α|−1}+ lower order powers of n.

Since p⁽ⁿ⁾_ξβ , p^∗_ξβ ⊂Π_|β|(IR^d), we have the convergence asserted in (3.30).

We now prove the strongest Korovkin theorem that the restricted positivity property (3.2) allows. Since this requires a modification of the usual argument, which is not stated in the literature, we give a self contained proof. This result supercedes Corollary 3.18.

Theorem 3.31 (Korovkin). Let T =T_V be the convex hull of V, and N_V be the region of nonnegativity. For f ∈C(T), Bnf →f uniformly on NV.

Proof: Let ε > 0, and M be the maximum of f over T. Since f is uniformly continuous on the compact set T, there is aδ >0 such that |f(s)−f(t)|< ε, ∀ks−tk< δ.

Thus, we obtain the estimate

|f(s)−f(t)| ≤ε+ 2Mks−tk²

δ² , ∀s, t∈T, For fixed t, we have

−ε− 2M

δ² k · −tk² ≤f−f(t)≤ε+ 2M

δ² k · −tk² onT , and so applying B_n (which reproduces constants) and using (3.2) gives

−ε− 2M

δ² B_n k · −tk²

≤B_nf −f(t)≤ε+ 2M

δ² B_n k · −tk²

on N_V. This last step is the main difference in argument. For t ∈N_V, evaluating at t gives

|Bnf(t)−f(t)| ≤ε+ 2M

δ² B_n k · −tk²

(t), ∀t∈N_V. (3.32) We now estimate B_n(k · −tk²)(t). From (2.1), we obtain

k · −tk² =hX

v

{ξ_vv−ξ_vt},X

w

{ξ_ww−ξ_wt}i=X

v

X

w

ξ_vξ_whv−t, w−ti.

Thus, (3.17) gives Bn k · −tk²

=X

v

X

w

n 1− 1

n

ξvξw+ 1 n

X

u∈V

ξv(u)ξw(u)ξu

o

hv−t, w−ti

= 1− 1

n

k · −tk²+ 1 n

X

v

X

w

X

u∈V

ξv(u)ξw(u)ξuhv−t, w−ti

= 1− 1

n

k · −tk²+ 1 n

X

u∈V

ξuhu−t, u−ti, so that

B_n k · −tk²

(t) = 1 n

X

u∈V

ku−tk²ξ_u(t) = 1 n

X

u∈V

kuk²ξ_u(t)− ktk²

≤ 1 nD², where D:= diam(T) is the diameter of T. Thus, we obtain the uniform estimate

|Bnf(t)−f(t)| ≤ε+ 1 n

2M D²

δ² , t ∈NV, and conclude that Bnf →f uniformly on NV.

4. Applications to CAGD

We now consider how our generalised Bernstein operator B_n = B_n,V might be used in CAGD (computer aided geometric design) to describe polynomials defined on convex polyhedra which are not simplices. We first observe that the Korovkin theory does not allow a Bernstein type operator which is positive on the entire region T = conv(V) if the points of V are not the vertices of a simplex or a cube.

Theorem 4.1 ([MRS94]). LetT be a convex polygon with five or more vertices. There is no positive linear operator L:C(T)→C(T)which reproduces Π₁(IR²) and maps Π₂(IR²) to itself other than the identity.

We recall, that our B_n reproduces Π1(IR²) and maps Π2(IR²), but has the restricted positivity property (3.2), i.e.,

f ≥0 on T =⇒ Bnf ≥0 on NV,

where NV is the region of nonnegativity (2.5). In the multivariate case, the description of other shape preserving properties is involved (cf. [S91]). We mention some which do generalise easily. These suggest that the target region for representing polynomials on should perhaps be the convex polyhedron N_V, rather than T = conv(V).

Proposition 4.2 (Shape preservation). Let T =TV be the convex hull of V, andNV

be the region of nonnegativity. If f is convex onT, then

B_nf ≥B_n+1f ≥ · · · ≥f on N_V. (4.3) Proof: Fixx∈N_V. By (2.1), and a calculation, we have the convex combinations

x= X

|α|=n

Bα(x)vα, vβ = X

v∈V

β_v

|β|v_β−ev. Hence for f convex, Jensen’s inequality gives

Bnf(x) = X

|α|=n

Bα(x)f(vα)≥f(x), X

v∈V

βv

|β|f(v_β−ev)≥f(vβ). (4.4) By the degree raising formula (2.13), we have

B_nf = X

|α|=n

B_αf(v_α) = X

v∈V

X

|α|=n

αv + 1

|α|+ 1B_α+evf(v_α), so that

B_nf−B_n+1f = X

|β|=n+1

c_βB_β, c_β := X

v∈V

β_v

|β|f(vβ−ev)−f(vβ).

By (4.4), cβ ≥0, so that Bnf ≥B_n+1f on NV.

A polynomial in B-form can be calculated via the de Casteljau algorithm. We present this in terms of the blossom (polar form) of a polynomial p∈Πn(IR^d), which we recall (cf. [R89], [DLG91], [DMS92]) is the unique symmetric n–affine function ^ωp with

ωp(t, . . . , t) =p(t), ∀t ∈IR^d. (4.5) Proposition 4.6 (de Casteljau algorithm). Suppose that p∈Πn(IR^d) has the B-form

p= X

|α|=n

c_αB_α. (4.7)

Then the blossom of p at t₁, . . . , t_n can be calculated from (cα)_|α|=n via for j = 1 to n do

cα := X

v∈V

ξv(tj)c_α+ev, |α|=n−j end for

with c₀ =^ωp(t₁, . . . , tn). In particular, taking t₁ =· · ·=tn =t gives c₀ =p(t).

Proof: Suppose p∈Πn(IR^d) is in the B-form (4.7), i.e., equivalently p=X

v1

· · ·X

vn

c_ev

1+···+e_vnξ_v1· · ·ξ_vn. (4.8) Then its blossom is given by

P(t₁, . . . , tn) :=X

v1

· · ·X

vn

cev1+···+e_vnξv1(t₁)· · ·ξvn(tn), (4.9) since P clearly defines an n–affine function, which, by (4.8), satisfies (4.5). The n steps of the algorithm are the averagings given by the n sums in (4.9).

In the terminology of [DMS92], this is a symmetric simplicial algorithm (here the points V need not be the vertices of a simplex).

The coefficients c_α in the B-form (4.7) of p are not unique, unlessV is the vertices of a simplex. Using (2.1) to expand, we have

p(t) =^ωp(t, . . . , t) =^ωp(X

v1

ξ_v1(t)v₁, . . . ,X

vn

ξ_vn(t)v_n)

=X

v1

· · ·X

vn

ξ_v1(t)· · ·ξ_vn(t)^ωp(v1, . . . , v_n), and so the coefficients c_α in (4.7) can be chosen to be the blossoms

b_α =^ωp(α·V), α·V := (. . . , v, . . . , v

| {z }

αv times

, . . .), (4.10)

which we call the blossoming coefficients.

Proposition. For any p of the form (4.7), the blossoming coefficients b = (bα)_|α|=n are given by the matrix multiplication

b=Qc, Q= [qαβ]|α|=n,|β|=n, q_αβ :=B^ω_β(α·V).

For n >1, our calculations show that the blossoming coefficients are not the ℓ₂–norm minimising choice (which is given by multiplication by an orthogonal projection matrix).

Remark 4.11. The previous discussion extends to p = (p₁, . . . , ps) : IR^d → IR^s, where p_j ∈Πn(IR^d) and c_α ∈IR^s, by considering coordinates. Since vector valued p are used in practice, e.g., B´ezier curves in IR³, we henceforth state our results in this setting.

We define the control points of the curve (d = 1), surface (d = 2), etc, given by t7→p(t)∈IR^s, where p has B-form (4.7), to be

nv_α c_α

:|α|=no

⊂IR^d+s.

Since the B-form is not unique when V has more thand+ 1 points, a given curve (surface, etc) may be given by different sets of control points. This redundancy has advantages from the point of view of design, as a given surface could be arrived at by a number of different choices of control points. Equivalently, each control point has less influence on the the surface, and so moving a control point causes a small change of shape making the fine tuning of a surface easier. Once a suitable surface has been obtained it can be presented in terms of a set of canonical control points – if so desired.

By combining (3.14) and (4.7), we have x

p(x)

= X

|α|=n

cαBα(x), cα :=

vα

c_α

. (4.12)

Thus a point x, p(x)

on the curve (surface, etc) can be calculated by the de Casteljau algorithm applied to the vectors (cα)_|α|=n ⊂IR^d+s. Many basic properties of B´ezier curves and surfaces follow from (4.12). We now list these using the terminology of [F02].

Convex hull property: Since P

|α|=nB_α(x) = 1 and B_α(x) ≥ 0, x ∈ N_V, the point x, p(x)

is an affine combination (also called a barycentric combination in CAGD) of the control points (cα)_|α|=n, which is a convex combination if x∈N_V.

Affine invariance: For A : IR^s →IR^ℓ an affine map, we have A p(x)

= X

|α|=n

(Ac_α)B_α(x).

Invariance under affine parameter transformations: ForA: IR^d →IR^d an invertible affine map, B_α(Ax) =A(B_α(x)), so that

p(Ax) = X

|α|=n

c_αA B_α(x) .

Symmetry: Suppose that A is an affine map which maps (the multiset) V to V. Then ξv(x) =ξAv(Ax) (see [W11₂]), so that

B_α(x) =B_α◦A⁻¹(Ax),

where α◦A⁻¹ :V →ZZ₊ denotes the multiindex v7→α_A⁻¹_v. Using this we obtain p(x) = X

|α|=n

c_αB_α(x) = X

|α|=n

c_αB_α◦A⁻¹(Ax) = X

|β|=n

c_β◦AB_β(Ax).

Invariance under barycentric combinations: The affine combination of two curves (surfaces, etc) is given by corresponding affine combination of the control points, i.e.,

λ X

|α|=n

b_αB_α(x) + (1−λ) X

|α|=n

c_αB_α(x) = X

|α|=n

λb_α + (1−λ)c_α B_α(x), λ∈IR.

Endpoint interpolation: If v∈V satisfies v6∈Aff(V \v), then Prop. 2.6 gives (Bnf)(v) =f(v).

The remaining property given in[F02]ispseudolocal control, i.e., the fact thatB_α is peaked at the point vα, and so moving the control point cα has the most influence on the curve (surface, etc) near the point v_α. Numerical calculations indicate some degree of localisation of the Bα near vα, but we do not make any quantitative statements here.

To summarise, the main features of our Bernstein polynomial approximations on a nonsimplicial convex polytope are:

• A surface may be determined by several different choices of control points.

• Moving the control points has less influence on the shape if there are many of them.

• On N_V the surface is a convex combination of the control points.

• Convex functions are monotonely approximated from above by B_n on N_V.

Typical applications, that require a single polynomial defined on a convex polytope, include the description and construction of finite elements on regular polygons and the design of hexagonal lenses (cf. [D87]).

5. Smoothness and multivariate splines

We now consider the possible application of our results to multivariate spline theory, where surfaces are constructed by joining polynomial pieces together smoothly. There is a highly developed theory (and associated software) that involves polynomials defined on simplices, based on the description of the smoothness conditions in terms of the control

points (cf. [B87], [LS07]). We assume that the reader is familiar with this, and we inves- tigate how it extends to partitions involving nonsimplical cells. There was work in this direction (see [L89], [CL90]) on mixed grid partitionswhere in the bivariate case the cells were trianglesandparallelograms, and in the trivariate case they weretetrahedrons,prisms and parallelopipeds. The B-form developed for nonsimplicial cells was for polynomials of coordinate degree n (as were the spline spaces), rather than total degree n, though the positioning of the control points is the same as we propose. Bivariate C¹–quadratics on a polygonal partition were considered in [Wh90].

We will give our results in terms of the blossoming coefficients. The following is adapted from the[B87],[R89]and[L91](cf.[W13]). If (cα) are the blossoming coefficients, then we refer to the (cα) of (4.12) as the blossoming control points of p = P

cαBα. Let

r·x:=x, . . . , x

| {z }

r times

,

Theorem 5.1.. LetW be a set of points inIR^d, with L= Aff(W). Then the polynomials f, g ∈Πn(IR^d) have all derivatives of order ≤r onL equal if and only if

ω

f(w1, . . . , w_n−r, r·x) =^ωg(w1, . . . , w_n−r, r·x), w₁, . . . , w_n−r ∈W, x∈IR^d. (5.2) Let

α·V := (. . . , v, . . . , v

| {z }

αv times

, . . .), α∈ZZ⁺_V.

Lemma 5.3. Let cbe the blossoming coefficients of (4.10) for p∈Πn(IR^d). Then

ωp(α·V, r·x) = X

|β|=r

c_α+βBβ(x), |α|=n−r, r = 0, . . . , n.

Proof: Since the blossom^ωp is n–affine, we have

ωp(α·V, r·x) =^ωp(α·V, X

v1∈V

ξ_v1(x)v1, . . . , X

vr∈V

ξ_vr(x)vr)

= X

v1∈V

· · · X

vr∈V

ξ_v1(x)· · ·ξ_vr(x)^ωp(α·V, v₁, . . . , v_r)

= X

|β|=r

c_α+βB_β(x).

Combining Theorem 5.1 and Lemma 5.3 gives smoothness conditions for the C^r– joining of polynomials in terms of their blossoming B–form coefficients.

We now illustrate this for C¹–quadratics on a region consisting of a triangle and a quadrilateral with a common edge.

Example 5.4 (C¹–quadratics on triangle and quadrilateral). Without loss of gen- erality, suppose the vertices of the triangle and quadrilateral are

V ={v1, . . . , v₃}={0, e1, a}, a₂ 6= 0, V˜ ={˜v₁, . . . ,v˜₄}={0, e1, b,−e2}, so the common edge has endpoints v₁ = ˜v₁ = 0. and v₂ = ˜v₂ =e₁ = (1,0). We index the generalised barycentric coordinatesξ and ˜ξ by both their order and the points themselves, e.g., ξ₃ and ξ_a =ξ_(a1,a2), whichever is the most convenient. We have

ξ₁(x, y) = (1−x) + a₁−1

a₂ y, ξ₂(x, y) =x− a₁

a₂y, ξ₃(x, y) = 1 a₂y, ξ˜_(0,0)(x, y) = −(b²₂+b₁b₂+b₁+ 1)x+ (b²₁+b₁b₂−b₁+ 1)y+b²₁+b²₂+ 1

2(b²₁+b²₂−b₁b₂+b₂−b₁+ 1)

ξ˜_(1,0)(x, y) = (2b²₂−b₁b₂+ 2b₂−b₁+ 2)x+ (b²₁−2b₁b₂−b₁)y+b²₁−b₁b₂−b₁ 2(b²₂+b₂+ 1−b₁b₂−b₁+b²₁) , ξ˜_(b1,b2)(x, y) = (2b1−b₂−1)x+ (2b2−b₁+ 1)y+b₂−b₁+ 1

2(b²₂+b₂+ 1−b₁b₂−b₁+b²₁) ,

ξ˜_(0,−1)(x, y) = (−b²₂+ 2b₁b₂−b₂)x+ (−2b²₁+b₁b₂+ 2b₁−b₂−2)y+b²₂−b₁b₂+b₂ 2(b²₂+b₂+ 1−b₁b₂−b₁+b²₁) , Suppose that f, g ∈Π₂(IR²) are quadratics with blossoming coeffient B–forms

f = X

|ga|=2 α∈ZZV

+

c_αB_α, g= X

|ga|=2 α∈ZZV˜ +

˜ c_αB˜_α.

In Theorem 5.1, we take n= 2, and

W ={v₁, v₂}={˜v₁,v˜₂}={0, e₁}.

For r = 0, the condition forf and g to have a continuous join on the line L= Aff(W) is given by the two element sequences from W ={v₁, v₂}, i.e.,

ω

f(v₁, v₁) =^ωg(v₁, v₁), f^ω(v₁, v₂) =^ωg(v₁, v₂), ^ωf(v₂, v₂) =^ωg(v₂, v₂).

which by Lemma 5.3 is

c_(2,0,0) = ˜c_(2,0,0,0), c_(1,1,0) = ˜c_(1,1,0,0), c_(0,2,0) = ˜c_(0,2,0,0).

For r = 1, the condition for f and g to have a C¹–join on L is given by the one element sequences from W ={v₁, v₂}, i.e.,

ωf(v₁, x) =^ωg(v₁, x), f^ω(v₂, x) =^ωg(v₂, x).

By Lemma 5.3, these equalities of linear polynomials in x can be written

c_(2,0,0)ξ₁+c_(1,1,0)ξ₂+c_(1,0,1)ξ₃ = ˜c_(2,0,0,0)ξ˜₁+ ˜c_(1,1,0,0)ξ˜₂+ ˜c_(1,0,1,0)ξ˜₃+ ˜c_(1,0,0,1)ξ˜₄, c_(1,1,0)ξ₁+c_(0,2,0)ξ₂+c_(0,1,1)ξ₃ = ˜c_(1,1,0,0)ξ˜₁+ ˜c_(0,2,0,0)ξ˜₂+ ˜c_(0,1,1,0)ξ˜₃+ ˜c_(0,1,0,1)ξ˜₄. Since the blossoming control points for the quadrilateral are simply those for corresponding triangles, it follows that these smoothness conditions have the usual geometric interpreta- tion: that all the control points involved (3 and 4 respectively, see Fig. 4) lie in a common plane.

Fig. 4. The C¹–smoothness conditions of Example 5.4.

6. Future work

For applications where it is desirable to have generalised barycentric coordinates which are nonnegative on the convex hull of V, one could modify the definition of (ξv(x))_v∈V for x∈conv(V) to be the unique minimal ℓ₂–norm coefficients a∈IR^V satisfying

x= X

v∈V

avv, X

v∈V

av = 1, av ≥0. (6.1)

With this definition ξv is a continuous piecewise linear polynomial, which coincides with with the original on N_V.

References

[AC94] F. Altomare and M. Campiti, “Korovkin-type approximation theory and its applica- tions”, Walter de Gruyter, Berlin, 1994.

[B87] C. de Boor, B-form basics, in “Geometric Modeling: Algorithms and New Trends”

(G. E. Farin Ed.), pp. 131–148, SIAM Publications, Philadelphia, 1987.

[BD85] C. de Boor and R. DeVore, Partitions of unity and approximation,Proc. Amer. Math.

Soc. 93 (1985), 705–709.

[C03] O. Christensen, “An introduction to frames and Riesz bases”, Birkh¨auser, Boston, 2003.

[CL90] C. K. Chui and M.-J. Lai, Multivariate vertex splines and finite elements,J. Approx.

Theory 60 (1990), 245–343.

[CW00] S. Cooper and S. Waldron, The eigenstructure of the Bernstein operator, J. Approx.

Theory 105 (2000), 133-165.

[CW02] S. Cooper and S. Waldron, The diagonalisation of the multivariate Bernstein operator, J. Approx. Theory 117 (2002), 103–131.

[DMS92] W. Dahmen, C. A. Micchelli, and H.-P. Seidel, Blossoming begetsB-spline bases built better by B-patches, Math. Comp. 59(1992), 97–115.

[DLG91] T. DeRose, M. Lounsbery, and R. Goldman, A tutorial introduction to blossoming, Comput. Graph. Systems Appl. (1991), 267–286.

[D11] Tan Duc Do, Generalised Bernstein operators on polytopes, Honours dissertation, 2011.

[D87] C. F. Dunkl, Orthogonal polynomials on the hexagon, SIAM J. Appl. Math. 47 (1987), 343–351.

[F02] G. Farin, “Curves and surfaces for computer-aided geometric design: A practical guide (5th ed.)”, Academic Press, San Diego, 2002.

[L89] M.-J Lai, On construction of bivariate and trivariate vertex splines on arbitrary mixed grid partitions, Dissertation (Texas A & M Univ.), 1989.

[L91] Ming Jun Lai, A characterization theorem of multivariate splines in blossoming form, Comput. Aided Geom. Design 8 (1991), 513-521.

[LS07] M.-J. Lai and L. L. Schumaker, “Spline functions on triangulations”, Cambridge Uni- versity Press, Cambridge, 2007.

[L53] G. G. Lorentz, “Bernstein Polynomials”, Toronto Press, Toronto, 1953.

[MRS94] F.-J. Mu˜noz-Delgado, V. Ram´irez Gonz´alez, and T. Sauer, Domains for Bernstein polynomials, Appl. Math. Lett. 7 (1994), 7–9.

[S91] T. Sauer, Multivariate Bernstein polynomials and convexity, Comput. Aided Geom.

Design 8 (1991), 465–478.

[R89] L. Ramshaw, Blossoms are polar forms, Comput. Aided Geom. Design 6 (1989), 323- 358.

[W11₁] S Waldron, Frames for vector spaces and affine spaces, Linear Algebra Appl. 435 (2011), 77–94.

[W11₂] S Waldron, Affine generalised barycentric coordinates, Jaen J. Approx. 3 (2011), 209–226.

[W13] S Waldron, Blossoming and smoothly joining polynomials on affine subspaces, Preprint, 2013.

[Wh90] W. Whiteley, The geometry of bivariate C₂¹-splines, Preprint, 1990.