Minimal energy surfaces on Powell–Sabin type triangulations

Applied Numerical Mathematics 58 (2008) 635–645 www.elsevier.com/locate/apnum

Minimal energy surfaces on Powell–Sabin type triangulations ✩ D. Barrera, M.A. Fortes ∗ , P. González, M. Pasadas Departamento de Matemática Aplicada, Universidad de Granada, Edificio Politécnico, C/Severo Ochoa, s/n 18071 Granada, Spain Available online 9 February 2007

Abstract In this paper we present a method to obtain an explicit surface on a polygonal domain D which approximates a Lagrangian data set and minimizes a certain “energy functional”. The minimization space is the C 1 -quadratic spline space constructed from an α-triangulation over D and its Powell–Sabin subtriangulation, i.e., we obtain a C 1 -polynomial with the minimal possible degree. A convergence result is established and some numerical and graphical examples are analyzed. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: α-Triangulation; Powell–Sabin element; Variational spline; Minimal energy; Approximation; Smoothing

1. Introduction In the past few years, variational methods in CAD, CAGD and Earth Sciences have received considerable attention, due to their efficiency and usefulness in the fitting and design of curves and surfaces. The basic idea of these methods is to minimize a functional which typically contains two terms: the first indicates how well the curve or surface approximates a given data set, while the second controls the degree of smoothness or fairness of the curve or surface. A wide range of minimization functionals have been proposed, derived from physical considerations (e.g. stretch energy, bending energy) or geometric entities (e.g. curve length, surface area, curvature). See, for example [4,6–8] and references therein. Discrete smoothing D m -splines [1,2] and discrete smoothing variational splines [7] provide specific examples of variational curves and surfaces. These splines minimize, in a finite element space, some quadratic functionals that contain terms associated with Sobolev seminorms. In [8] a functional of the above type is minimized in a parametric space of bicubic splines. Moreover, in all cases, the splines obtained approximate a Lagrangian data set (that is, a set of values of a known or unknown function) or a Hermite data set (that is, a set of values of a known or unknown function as well as the values of some of their partial derivatives). In this work we present a method to obtain a C 1 -quadratic spline surface on a polygonal domain D ⊂ R2 which approximates a Lagrangian data set and minimizes a certain “energy functional” given by a linear combination of the usual semi-norms | · |j , j = 1, 2, on the Sobolev space H 2 (D). The minimization space is the C 1 quadratic spline ✩

This work has been supported by the Junta de Andalucia (Research Group FQM/191) and by the Ministerio de Educación y Ciencia of Spain (Research Projects MTM2005-01331 and MTM2005-01403). * Corresponding author. E-mail address: [email protected] (M.A. Fortes). 0168-9274/$30.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2007.02.001

636

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

space constructed from an α-triangulation Tα over D and its Powell–Sabin associated subtriangulation Tα (cf. [11]). Thus, the solution of the problem is a polynomial spline of degree two and class C 1 . The remainder of this paper is organized as follows: In Section 2, we recall some preliminary notations and results. Section 3 is devoted to formulate the problem and to present a method to solve it, while in Section 4 we briefly describe the method to obtain the basis functions with local support over the unit reference triangle. In Section 5 a convergence result is proved. In the last section, we give some numerical and graphical examples for the Franke’s and the Nielson’s test functions. 2. Notation and preliminaries Let D ⊂ R2 be a polygonal domain and let us consider the Sobolev space H 2 (D), whose elements are (classes of) functions u defined in D such that their partial derivatives (in the distribution sense) ∂ β u belong to L2 (D), with β = (β1 , β2 ) ∈ N2 and |β| = β1 + β2  2. In this space we consider the usual norm    1 2 u = ∂ β u(x)2 dx , |β|2 D

the semi-norms  1 2 β 2 |u|j = ∂ u(x) dx ,

j = 1, 2;

|β|=j D

and the corresponding inner semi-products  (u, v)j = ∂ β u(x)∂ β v(x) dx,

j = 1, 2.

|β|=j D

In what follows, Pn (X) denotes the space of bivariate polynomials of total degree at most n over X. Given α  1, let Tα be an α-triangulation of D, i.e., a triangulation that satisfies the condition 1

RT α 2rT

for all triangle T ∈ Tα , RT and rT being the radius of the circumscribed and inscribed circles of T , respectively (see [11]). The associated Powell–Sabin triangulation Tα is obtained by joining the centre ΩT of the inscribed circle of each interior T ∈ Tα to the vertices of T and to the centres ΩT  of the inscribed circles of the neighbouring triangles T  ∈ Tα . When T has a side that is on the boundary of D, the point ΩT is joined to the mid-point of this side, to the vertices of T and to the centres ΩT  of the inscribed circles of the neighbouring triangles T  ∈ Tα . It is well known (cf. [10]) that given the values of a function f (defined on D) and its first partial derivatives at all the vertices of Tα , there exists a unique S in   P12 (D, Tα ) = S ∈ C 1 (D): S|T ∈ P2 (T ) for all T ∈ Tα such that the values of S and its first partial derivatives at the vertices of Tα coincide with those of f . For the sake of completeness, we give two more definitions that will be used through the work: Let (V ,  · ) be a normed vector space. A bilinear form f : V × V → R is called V -elliptic if there exists β > 0 such that βv2  f (v, v) for all v ∈ V (see [3]). On the other hand, a subset {a0 , a1 , a2 } ⊂ D is called P1 -unisolvent if given real scalars {β0 , β1 , β2 }, there exists a unique function p ∈ P1 (D) such that p(ai ) = βi for i = 0, 1, 2. 3. Formulation of the problem Given a finite set of points D r in D and a set of values Z r = {g(a)}a∈D r (for r ∈ N∗ = N − {0}), we are looking for a C 1 -quadratic surface that approximates the points {(a, g(a))}a∈D r ⊂ R3 and minimizes a “functional energy” that is described below:

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

637

Let k = k(r) = card(D r ) and let us denote · k (respectively ·,· k ) the usual Euclidean norm (respectively Euclidean inner product) in Rk . Let ρ r be the evaluation operator ρ r : H 2 (D) → Rk defined by ρ r (v) = (v(a))a∈D r and let us suppose that  Ker ρ r ∩ P1 (D) = {0} for all r ∈ N∗ . (1) Given τ = (τ1 , τ2 ), where τ1 ∈ [0, +∞) and τ2 ∈ (0, +∞), let us consider the functional defined on H 2 (D) by 2

2  τj |v|2j . Jr (v) = ρ r (v − g) k + j =1

Observe that the first term of Jr measures how well v approximates the values {g(a)}a∈D r over the set of points D r (in the least squares sense), while the second one represents the “minimal energy condition”. This second term is included in the functional Jr in order to dispose of certain fairness control. Recall that the semi-norms | · |1 and | · |2 can be seen as a simplified measure of the surface area and of the bending energy, respectively. The parameters τ1 and τ2 weight the importance given to each of the two semi-norms. The minimization problem we want to solve is: Given an α-triangulation Tα of D and its associated Powell–Sabin triangulation Tα , we look for an element σr ∈ P12 (D, Tα ) such that Jr (σr )  Jr (v)

for all v ∈ P12 (D, Tα ).

(2)

Proposition 1. Problem (2) has a unique solution that is also the unique solution of the following variational problem: Find σr ∈ P12 (D, Tα ) such that (3) ρ r (σr ), ρ r (v) k + τ1 (σr , v)1 + τ2 (σr , v)2 = Z r , ρ r (v) k , ∀v ∈ P12 (D, Tα ). Proof. Condition (1) allows us to assure that the application

1 2 2

r 2  2 v → [[v]] = ρ (v) k + τj |v|j j =1

is a norm on P12 (D, Tα ) equivalent to  · . As a consequence, it is easy to check that the symmetric continuous bilinear form a : P12 (D, Tα ) × P12 (D, Tα ) → R defined by 2  

τj (u, v)j a(u, v) = ρ r (u), ρ r (v) k + j =1

is P12 (D, Tα )-elliptic. On the other hand, the application ϕ : P12 (D, Tα ) → R defined by ϕ(v) = Z r , ρ r (v) k is a linear and continuous form and thus, Lax–Milgram Lemma allows us to assure that the functional 12 a(v, v) − ϕ(v) has a unique minimum that is the solution of (3). Finally, since   Z r 2k 1 , Jr (v) = 2 a(v, v) − ϕ(v) + 2 2 then the solution of (3) is also the unique minimum of Jr .

2

Let us denote N = dim(P12 (D, Tα )). If {v1 , . . . , vN } is a basis of the finite element space P12 (D, Tα ) whose elements  have local support, and σr = N i=1 αi vi , then, Problem (3) gives rise to the system of linear equations CX = B, where B = (Z r , ρ r (v1 ) k , . . . , Z r , ρ r (vN ) k )t , X = (α1 , . . . , αN )t and C = (cij )N i,j =1 , with

r  cij = ρ (vi ), ρ r (vj ) k + τ1 (vi , vj )1 + τ2 (vi , vj )2 .

(4)

638

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

Lemma 2. C is a symmetric, positive definite and banded matrix. Proof. From the expression of the elements cij , we easily deduce that C is symmetric. Let p = (p1 , . . . , pN ) ∈ RN , then

pCp = t

N 

pi cij pj =

i,j =1

 = ρ

N 

r

= ρr

pi

N 



pi v i , ρ

N  i=1

r

 2   τs (vi , vj )s pj ρ (vi ), ρ (vj ) k +



i,j =1

r

r

s=1 N 



2 pi v i k

+

pj v j

j =1

i=1



k

2 

τs

s=1

N  i=1

pi v i ,

N 

 pj v j

j =1

s

 N 2 2      + τs  pi vi   0.   s=1

i=1

s

 On the other hand, let p = (p1 , . . . , pN ) ∈ RN be such that pCp t = 0 and let us consider w = N i=1 pi vi . Then, 2 2 r 2 ρ (w) k + s=1 τs |w|s = 0 and hence, [[w]] = 0. So, w = 0 and as a consequence we get p = 0. Finally, C is banded because each vi has local support. 2 4. Computation We recall that if A is a triangle in R2 with vertices a1 , a2 and a3 and x ∈ R2 , then the barycentric coordinates of x with respect to A are defined as unique numbers (λ1 , λ2 , λ3 ) such that λ1 + λ2 + λ3 = 1 and x = λ1 a1 + λ2 a2 + λ3 a3 . Let us consider the reference triangle T0 with vertices A1 = (0, 1), A2 = (0, 0) and A3 = (1, 0) and the linear ∂f functionals given by Li (f ) = f (Ai ) for i = 1, 2, 3; Li (f ) = ∂f ∂x (Ai−3 ) for i = 4, 5, 6, and Li (f ) = ∂y (Ai−6 ) for i = 7, 8, 9. In order to determine a basis of functions {v1 , . . . , vN } with local support of the finite element vector space P12 (D, Tα ), and as a consequence the matrices B and C and the solution of the system of linear equations (4), we have considered a basis of functions {w1 , . . . , w9 } over T0 that verify Li (wj ) = δij . Let {T1 , . . . , T6 } be the microtriangles of the Powell–Sabin triangulation of T0 (see Fig. 1).  d λi λj λk Over each triangle Td every polynomial p of total degree two can be expressed as p(x) = i,j,k=0,1,2 cij k 1 2 3 i+j +k=2

where (λ1 , λ2 , λ3 ) is the barycentric coordinate vector of x with respect to Td , for all x ∈ Td . By applying the relations (see [5]) that must verify the B-coefficients of a given function f in order to be of class C 1 , we determine the Bcoefficients of the basis functions {wi }9i=1 . If we denote E0 the common vertex of the six microtriangles, then the d } of the functions w over every triangle T are as it is shown in Fig. 2. B-coefficients {cij i d k In Fig. 3 we give the B-coefficients of a multiple of the basis functions {wi }9i=1 . It is also possible to determine the basis functions {wi }9i=1 by using truncated powers (see [9]).

Fig. 1. Powell–Sabin triangulation of T0 .

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

639

Fig. 2. Representation of the B-coefficients over Td .

Fig. 3. B-coefficients of a multiple of the basis functions {wi }9i=1 .

5. Convergence Let A0 = {a10 , a20 , a30 } be a P1 -unisolvent subset of D and let us suppose that   1 , r → +∞. sup minr x − a 2 = O a∈D r x∈D

(5)

Then, there exists C > 0 and r1 ∈ N∗ such that for all r  r1 there exist {a1r , a2r , a3r } ⊂ D r verifying

ai0 − air

 2



C r

for all i = 1, 2, 3 and for all r  r1 .

(6)

Lemma 3. Let us suppose that (5) is verified and let Ar = {a1r , a2r , a3r } ⊂ D r be any subset verifying (6). Then, there exists r ∗ ∈ N such that for each r  r ∗ the application [[·]]r defined by

3 1 2   2 r r 2 v ai + |v|2 [[v]] = i=1

is a norm on

H 2 (D)

uniformly equivalent with respect to r to the norm  · .

640

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

Proof. Because of A0 is a P1 -unisolvent subset of D, there exists r0 ∈ N such that Ar is also a P1 -unisolvent subset of D for all r  r0 and as a consequence, [[·]]r is a norm on H 2 (D) for all r  r0 . From the continuous injection of H 2 (D) into C 0 ( D ), we obtain that there exists C1 > 0 such that maxx∈D |v(x)|  C1 v for all v ∈ H 2 (D). In particular, |v(air )|  C1 v for all i = 1, 2, 3; r ∈ N and v ∈ H 2 (D). Hence, there exists C2 > 0 such that 3   2 v air + |v|22  C2 v2 , i=1

and as a consequence, 1/2

[[v]]r  C2 v for all r  r0 and v ∈ H 2 (D).

(7)

On the other hand, for every r ∈ N it holds  2   r 2 1   0 2   0 v ai − v air v ai  + v ai , 2 3

3

3

i=1

i=1

i=1

and by applying again the continuous injection of H 2 (D) into C 0 ( D ), we obtain 

  2 2 1   0 2 ai0 − air 2 v2 + v ai + |v|22  v air + |v|22 2 3

3

3

i=1

i=1

i=1

for all r ∈ N.

Since A0 is a P1 -unisolvent subset, we have that the application

3 1 2 1   0 2 2 v → v ai + |v|2 2 i=1

is a norm on H 2 (D) that, besides, is equivalent to  ·  (see Proposition I-2.2 of [2]), hence, there exists C3 > 0 such that C3 v2 

3 

3   2 2 v air + |v|22 ai0 − air 2 v2 +

i=1

for all r ∈ N.

i=1

Moreover, by (6) we have that   2 3C 2 v2 + v air + |v|22 2 r 3

C3 v2 

for all r  r1 .

i=1

Let r2  r1 be such that

3C 2 r22

< C3 . Then, C4 = C3 −

  2  3C 2 C4 v2  C3 − 2 v2  [[v]]r r and as a consequence, 1/2

C4 v  [[v]]r

3C 2 r22

satisfies

for all r  r2 ,

for all r  r2 .

Finally, from (7) and (8) we obtain that the norms

(8) [[·]]r

and  ·  are uniformly equivalent with respect to r.

2

Let g ∈ C 3 (D), α  1, H ⊂ (0, +∞) a subset that admits 0 as accumulation point and, for each h ∈ H, let Tα be an α-triangulation such that maxT ∈Tα 2RT = h and sh ∈ P12 (D, Tα ) be the unique function that interpolates the values of g and its first partial derivatives at all the vertices of Tα . Lemma 4. There exist positive constants C, C1 and C2 , depending only on α and g, such that for all h ∈ H it holds   max(sh − g)(x)  Ch3 and (9) x∈D

|sh − g|j  Cj h3−j

for j = 1, 2.

(10)

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

641

Proof. From Theorem of Section 5 in [11] we deduce that there exist constants C, C1 and C2 , depending on α and g, but not on h and T , such that   max(sh − g)(x)  Ch3 and x∈T   max ∂ α (sh − g)(x)  Cj h3−j for j = 1, 2 |α|=j x∈T

for all T ∈ Tα , from which we conclude (9). Moreover    2   2 2 |sh − g|j  max∂ α (sh − g)(x) dx  nj · meas(D) · (Cj )2 h3−j , |α|=j D

x∈D

nj being the number of partial derivatives of order j for all j = 1, 2, and therefore we obtain (10).

2

Theorem 5. Let us suppose that, in addition to hypotheses (5), it holds: there exist C > 0 and r0 ∈ N − {0} such that k(r)  Cr 2 for all r  r0 ;  τ2 = o r 2 , r → +∞;

(12)

τ1 = o(τ2 ),

(13)

r → +∞;

r 2 h6 = o(1), τ2

r → +∞.

(11)

(14)

Let σr be the unique solution of problem (2) for the triangulation Tα fixed before. Then, lim g − σr  = 0.

(15)

r→+∞

Remark 6. Before going to the proof of this theorem, we are going to sketch the main ideas of it: From (12) and (14) we can deduce that h → 0, which together with (11) and (13) allow us to assure that |σr |2 is bounded as r → +∞. By using Lemma 3 we will show that σr  is also bounded, and hence the sequence {σr } has a convergent subsequence to g, from which we can prove that limr→+∞ g − σr  = 0. Proof of Theorem 5. Since σr is the solution of problem (2), it holds that Jr (σr )  Jr (sh ), that is: 2 2 2 

2  ρ r (σr − g) k + τj |σr |2j  ρ r (sh − g) k + τj |sh |2j .



j =1

(16)

j =1

By (9) we know that there exists C > 0 such that (sh (a) − g(a))2  Ch6 for all a ∈ D r , and by using (10) we can assure that 2  2 2 τ 1 1 r 1 τ1  ρ (sh − g) k + |sh |21 + |sh |22  k(r)Ch6 + C1 h2 + |g|1 + C2 h + |g|2 , τ2 τ2 τ2 τ2 and by using (11) we obtain that there exist C > 0 and r0 ∈ N∗ such that 2  2 2 τ 1 Cr 2 h6 τ1  1 r ρ (sh − g) k + |sh |21 + |sh |22  C1 h2 + |g|1 + C2 h + |g|2 + τ2 τ2 τ2 τ2

(17)

for all r  r0 . Moreover, from (12) and (14) we obtain h = o(1),

r → +∞,

(18)

and taking into account (13) and (14) we can assure that there exist C > 0 and r1 ∈ N∗ such that 2  2 Cr 2 h6 τ1  C1 h2 + |g|1 + C2 h + |g|2  C + τ2 τ2

for all r  r1 .

(19)

642

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

By using (16), (17) and (19) we obtain that |σr |22 

2  2 Cr 2 h6 τ1  C1 h2 + |g|1 + C2 h + |g|2  C + τ2 τ2

and that 2

r ρ (σr − g) k  Cτ2

for all r  r1

for all r  r1 .

(20)

(21)

Let A0 = {a10 , a20 , a30 } ⊂ D be a P1 -unisolvent subset and let r1 and C be the constants given in Eq. (6). Let r1  r1 be such that   C B aj0 ,  ⊂ D. r1 From (5) and (6) we know that    C 0 C ⊂ B aj ,  − r1 r



a∈D r ∩B(aj0 , C ) r1

C B a, r



for all r  r1 and j = 1, 2, 3.

If we denote Nj = card(D r ∩ B(aj0 , rC )), then we have that 

C C − r1 r

2

 2 C  Nj r

and as a consequence,   1 1 2 2 Nj   − r r1 r

1

for all r  r1 and j = 1, 2, 3,

for all r  r1 and j = 1, 2, 3.

From (12) and (21) we deduce that     σr (a) − g(a)2 = o r 2 ,

r → +∞, j = 1, 2, 3.

(22)

(23)

a∈D r ∩B(aj0 , C ) r1

Moreover, for all r  r1 and j = 1, 2, 3, there exists at least one point ajr ∈ D r ∩ B(aj0 , rC ) such that 1   r     σ r a − g a r  =   min σr (a) − g(a) . j j a∈D r ∩B(aj0 , C )

(24)

r1

From (22), (23) and (24) we have that   r   σr a − g a r  = o(1), r → +∞. j

j

(25)

= {a1r , a2r , a3r } ⊂ D r . By applying Lemma 3 to the subset Ar = B r , we obtain from (20) and (25) > 0 and r ∗ ∈ N such that σr   C for all r  r ∗ , what implies that the sequence {σr }rr ∗ is bounded

Br

Let us consider that there exist C in H 2 (D). As a consequence, there exists a subsequence {σrk }k∈N of {σr }rr ∗ , with τk = τ (rk ), limk→+∞ rk = +∞ and limk→+∞ hk = 0, and an element g ∗ ∈ H 2 (D) such that σrk

converges weakly to g ∗ in H 2 (D).

(26)

Let us suppose that g ∗ = g. Then, by using the compact injection of H 2 (D) into C 0 ( D ), we obtain that there exist θ > 0 and an open set D0 ⊂ D such that |g(x) − g ∗ (x)| > θ for all x ∈ D0 , and from (26) we obtain k0 ∈ N∗ such that |g ∗ (x) − σrk (x)|  θ2 for all k  k0 and x ∈ D0 . As a consequence, we have       g(x) − σr (x)  g(x) − g ∗ (x) − g ∗ (x) − σr (x) > θ , ∀k  k0 and ∀x ∈ D0 . (27) k k 2 Nevertheless, from (5) we deduce that for k large enough there exists a point a rk ∈ D rk ∩ D0 such that |σrk (a rk ) − g(a rk )| = o(1) as k → +∞, what contradicts (27). Therefore g ∗ = g.

(28)

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

643

For all k ∈ N we have that |σrk − g|22 = |σrk |22 + |g|22 − 2(σrk , g)2 , and then, by using (26), (28) and (20) we obtain that lim |σrk − g|2 = 0.

(29)

k→+∞

Besides, since H 2 (D) is compactly injected in H 1 (D), by using (26) and (28) we obtain that g = lim σrk k→+∞

in H 1 (D),

that jointly with (29) allow us to assure that limk→+∞ σrk − g = 0. Finally, let us suppose that (15) is not satisfied. Then, there exist α > 0 and sequences {τk }k∈N , {hk }k∈N and {rk }k∈N with τk = τ (rk ), limk→+∞ rk = +∞ and limk→+∞ hk = 0 such that σrk − g > α

for all k ∈ N.

(30)

As the sequence {σrk }k∈N is bounded in H 2 (D), a similar argument to that used for {σrk }k∈N would prove that there exists a subsequence of {σrk }k∈N that converges to g in H 2 (D), being in contradiction with (30); therefore (15) holds. 6. Numerical and graphical examples We have considered, for different values of r, sets D r consisting of k points arbitrarily distributed over the domain i n−1 }i=0 of the interval [0, 1] into D = [0, 1] × [0, 1]. We have taken, for different values of n, uniform partitions {ti = n−1 n − 1 subintervals, from which we obtain uniform partitions of D whose elements are {[ti , ti+1 ] × [tj , tj +1 ]}n−1 i,j =0 . By dividing each rectangle {[ti , ti+1 ] × [tj , tj +1 ]} by the diagonal joining the points (ti , tj +1 ) and (ti+1 , tj ), we obtain an ( 12 + √1 )-triangulation Tn of D formed by n2 vertices, from which we consider its associate Powell–Sabin’s 2 triangulation Tn . Example 1. Fig. 4 shows the graphic of Franke’s function in D, defined by f (x, y) = 0.75e−

(9x−2)2 +(9y−2)2 4

+ 0.75e−

(9x+1) (9y+1)2 10 − 49

+ 0.5e−

(9y−7)2 +(9x−3)2 4

− 0.2e−(9y−4)

2 −(9x−7)2

,

and Fig. 5 shows the approximating surfaces obtained for different values of k, τ1 and τ2 , and different triangulations Tn . The values of the parameters τ1 and τ2 have been chosen in such a way to obtain approximating surfaces with different degrees of “energy” and “closerness” to the graphic of the original function. In the middle and right approximating surfaces in Fig. 5, we consider small values of the parameters τ , in such a way that the summand of the functional Jr related to the approximating in the least square sense to the data point set has more weight, as seen in the figures. Meanwhile, the left approximating surface in Fig. 5 has been obtained with bigger values of the parameters τ giving as a result a surface which is smoother than the other ones but not reflecting as well the approximation to the data set point. The same criteria has been followed to obtain the different approximating surfaces of Nielson’s functions given in Fig. 7.

Fig. 4. Franke’s function in [0, 1] × [0, 1].

644

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

Fig. 5. Approximating spline surfaces for Franke’s function.

Fig. 6. Nielson’s function in [0, 1] × [0, 1].

Fig. 7. Approximating spline surfaces for Nielson’s function.

The error estimations have been computed by using the following relative error formula:  2500 2 1/2 ν=1 (f − σ )(aν ) E= , 2500 2 ν=1 f (aν ) where {a1 , . . . , a2500 } are arbitrary points in [0, 1]2 , f is Franke’s (or Nielson’s) function and σ = σr . Table 1 shows the error estimations obtained for the case k = 5000, different values of τ and different triangulations Tn . Example 2. In this case, we have considered Nielson’s function, defined by f (x, y) = y2 cos(4(x 2 + y − 1))4 (Fig. 6), and the approximating surfaces obtained for different values of k, τ1 and τ2 , and different triangulations Tn (Fig. 7). The errors committed in the preceding approximations are E = 2.21 × 10−2 , E = 6.68 × 10−3 and E = 2.28 × −3 10 , respectively. In order to pay a special attention to the hollow of Franke’s function, we have also generated approximating spline surfaces constructed on a non-uniform Powell–Sabin triangulation T∗ of D. Such a hollow appears in Franke’s 6 9 3 6 , 10 ] × [ 10 , 10 ] of D. Triangulation T∗ is the function graphic when it is represented on the subdomain D  = [ 10 refinement of T11 obtained by substituting T11 |D  for the uniform triangulation of D  into seventy two triangles, and

D. Barrera et al. / Applied Numerical Mathematics 58 (2008) 635–645

645

Fig. 8. Triangulation T∗ of D and its Powell–Sabin associated triangulation. Table 1 τ1

τ2

T5

T10

T20

T30

T40

10−2 10−2 10−5 0

10−2 10−5 10−5 10−8

3.37 × 10−2 2.53 × 10−2 2.53 × 10−2 2.53 × 10−2

1.61 × 10−2 2.48 × 10−3 2.48 × 10−3 2.49 × 10−3

1.3 × 10−2 3 × 10−4 2.54 × 10−4 2.49 × 10−4

1.35 × 10−2 3.85 × 10−4 1.26 × 10−4 1.17 × 10−4

1.36 × 10−2 3.71 × 10−4 1.49 × 10−4 9.83 × 10−5

Table 2

Table 3

τ1

τ2

T∗

T11

τ1

τ2

T∗

T11

10−2

10−2

1.25 × 10−1

1.43 × 10−1

10−2

10−2

1.66 × 10−4

2.37 × 10−2 1.43 × 10−3 1.35 × 10−3 1.34 × 10−3

10−2 10−5 0

10−5 10−5 10−8

5.58 × 10−2 1.26 × 10−2 1.18 × 10−2

5.91 × 10−2 1.38 × 10−2 1.25 × 10−2

10−2 10−5 0

10−5 10−5 10−8

1.12 × 10−3 1.04 × 10−3 1.21 × 10−3

replacing appropriately the triangles surrounding D  in order to assure that the final scheme becomes a triangulation, as shown on the left of Fig. 8. On the right of Fig. 8 we show the Powell–Sabin subtriangulation associated to T∗ . We want to remark that T∗ is also a ( 12 + √1 )-triangulation. Tables 2 and 3 show the error estimations in the cases 2 k = 100 and k = 2000, respectively, for different values of τ by using the triangulations T∗ and T11 . References [1] R. Arcangeli, B. Ycart, Almost sure convergence of smoothing D m -splines for noisy data, Numer. Math. 66 (1993) 281–294. [2] R. Arcangèli, M.C. López de Silanes, J.J. Torrens, Multidimensional Minimizing Splines, Kluwer Academic, Dordrecht, 2004. [3] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, SIAM – Classics in Applied Mathematics, vol. 40, North-Holland, Amsterdam, 2002. [4] J. Duchon, Interpolation des fonctions de deux variables suivant le principle de la flexion des plaques minces, RAIRO 10 (12) (1976) 5–12. [5] G. Farin, Bézier polynomials over triangles and the construction of piecewise C r -polynomials, Technical Report TR/91, Brunel University, Uxbridge, England, 1980. [6] G. Greiner, Surface construction based on variational principles, in: P.J. Laurent, A. Le Méhauté, L.L. Schumaker (Eds.), Wavelets, Images and Surface Fitting, A.K. Peters, Wellesley, 1994, pp. 277–286. [7] A. Kouibia, M. Pasadas, Approximation by discrete variational splines, J. Comput. Appl. Math. 116 (2000) 145–156. [8] A. Kouibia, M. Pasadas, Approximation of surfaces by fairness bicubic splines, Adv. Comput. Math. 20 (2004) 87–103. [9] P. Lancaster, K. Šalkauskas, Curve and Surface Fitting, Academic Press, London, 1986. [10] M.J.D. Powell, M.A. Sabin, Piecewise quadratic approximations on triangles, ACM Trans. Math. Software 3 (4) (1977) 316–325. [11] P. Sablonnière, Error bounds for Hermite interpolation by quadratic splines on an α-triangulation, IMA Journal of Numerical Analysis 7 (4) (1987) 495–508.