AN INTERPOLATION ERROR ESTIMATE ON ... - Semantic Scholar

SIAM J. NUMER. ANAL. Vol. 45, No. 6, pp. 2368–2391

c 2007 Society for Industrial and Applied Mathematics 

AN INTERPOLATION ERROR ESTIMATE ON ANISOTROPIC MESHES IN Rn AND OPTIMAL METRICS FOR MESH REFINEMENT∗ WEIMING CAO† Abstract. In this paper, we extend the work in [W. Cao, Math. Comp., to appear] to functions of n dimensions. We measure the anisotropic behavior of higher-order derivative tensors by the “largest” (in certain sense) ellipse/ellipsoid contained in the level curve/surface of the polynomial for directional derivatives. Given the anisotropic measure for the interpolated functions, we derive an error estimate for piecewise polynomial interpolations on meshes that are quasi-uniform under a given metric. By using the inertia properties for matrix eigenvalues [R. C. Thompson, J. Math. Anal. Appl., 58 (1977), pp. 572–577] and H¨ older’s inequality, we can identify the optimal mesh metrics leading to the smallest error bound in various norms. Furthermore, we develop a dimensional reduction method to find the anisotropic measure approximately. We present two numerical examples for linear and quadratic interpolation on various anisotropic meshes generated with the optimal mesh metrics developed in this paper. Numerical results show that the smallest interpolation error is attained exactly on meshes optimal for the corresponding error norm as predicted. Key words. interpolation error, anisotropic mesh, anisotropic measure, aspect ratio, mesh alignment, mesh metric AMS subject classifications. 65D05, 65L50, 65N15, 65N50 DOI. 10.1137/060667992

1. Introduction. It is well known in mesh generation and finite element analysis communities that in order to minimize linear interpolation errors, the eigenvalues and eigenvectors of Hessian matrices can be used to determine the element aspect ratio and mesh alignment direction for anisotropic mesh generation or refinement; see, e.g., [3, 4, 11, 13, 15, 19, 21, 23, 24]. In the case of quadratic or higher-order interpolations, the error is determined by the third- or higher-order partial derivatives of the interpolated functions. Measuring their anisotropic behavior is the key for anisotropic mesh design and refinement [2]. In particular, in order to determine an ideal element orientation and aspect ratio, one needs to define the “principal direction” and the “strength” to characterize the anisotropic behavior of the derivative tensors. In a previous paper [6], the author developed a method to measure the orientation and anisotropic ratio of the higher-order derivative tensors for two-dimensional functions. The technique is based on decomposing the homogeneous polynomials for directional derivatives into the product of linear and nonnegative quadratic polynomials. Then the anisotropic measure is defined by the directions of the lines and ellipses corresponding to those factors. An interpolation error estimate is further derived on anisotropic meshes that are quasi-uniform under given metrics. Optimal mesh metrics can be identified to minimize the error bound in various norms. In this paper, we extend the work in [6] to functions of n dimensions. The difficulty in making such an extension is that it is generally impossible to factor a homogeneous polynomial in three or higher dimensions into the product of linear and nonnegative ∗ Received

by the editors August 22, 2006; accepted for publication (in revised form) April 6, 2007; published electronically November 21, 2007. This work was supported in part by NSF grant DMS-0209313. http://www.siam.org/journals/sinum/45-6/66799.html † Department of Mathematics, University of Texas at San Antonio, San Antonio, TX 78249 ([email protected]). 2368

INTERPOLATION ERROR ON ANISOTROPIC MESHES

2369

quadratic functions. However, the idea in [6], which characterizes the anisotropic behavior of derivative tensors by the largest ellipse contained in the level curve of the polynomial for directional derivatives (see Remark 2.1 in [6]), is still valid in higher dimensions. Algebraically, for any positive integer k, let ∇k+1 u(x) = ∇(∇k u)(x) be the (k + 1)th-order tensor for the partial derivatives of function u at a point x. We may characterize its anisotropic behavior by a suitable n × n positive definite matrix Q satisfying the following constraint: |(ξ · ∇)k+1 u(x)| ≤ |ξ · Qξ|

k+1 2

∀ξ ∈ Rn .

Based on the anisotropic measure of the interpolated functions, we derive an error estimate for the higher-order polynomial interpolation on meshes that are quasi-uniform under a given mesh metric. The error bound involves the same convergence rate with respect to the number of elements as in classical results on quasi-uniform meshes under Euclidean metrics. However, the constant in the error bound depends on an interplay between the anisotropic mesh (through the mesh metric) and the anisotropic behavior of ∇k+1 u, which can be much smaller than that in the classical error bound. Furthermore, by using the inertia properties of matrix eigenvalues proved by Thompson [26] and H¨ older’s inequality, we show that the smallest bound for the W m,p -error of kth-order interpolations is attained when the metric is defined as Mk+1,m,p = c(λmax (Q)) (k+1−m)p+n | det(Q)|− (k+1−m)p+n · Q, mp

1

where λmax (Q) is the largest eigenvalue of the matrix Q measuring the anisotropic behavior of ∇k+1 u. In the case of linear interpolation (k = 1), Q can be chosen as Q = [(∇2 u)2 ]1/2 + δ · In , where δ is a user-specified small parameter and In is the identity matrix. Then the above optimal metric for minimizing the Lp -error (i.e., with m = 0) is identical to that presented in Chen, Sun, and Xu [7]. In the case of higher-order interpolations, the optimal choice of Q relies on a constrained minimization. We develop a dimensional reduction method to find an approximate Q to measure the anisotropic behavior of ∇k+1 u. To test the optimality of the proposed mesh metrics, we generated various anisotropic meshes using the black-box anisotropic mesh generator bamg [3, 16] supplied with the metric Mk+1,m,p . We compare various error norms for linear and quadratic interpolations in two dimensions. It is found in all the cases that the smallest error norm is attained exactly with meshes based on the corresponding optimal mesh metric. We also present a three-dimensional example to demonstrate that the smallest error norms are attained at the best aspect ratios, as expected. An outline of this paper is as follows. In section 2 we present the error estimate for interpolations on anisotropic meshes that are quasi-uniform under a given metric. The metrics leading to the smallest error bounds (in various norms) are identified. In section 3 we introduce the notion of anisotropic measure of higher-order derivative tensors, and present a dimension reduction algorithm to determine the measure approximately. In section 4 we present two numerical examples demonstrating the optimality of the proposed mesh metrics. We conclude the paper with some discussion. 2. Error estimates on anisotropic meshes. 2.1. Basic assumptions and lemmas. We introduce in this subsection some assumptions regarding the anisotropic behaviors of the meshes and the higher-order derivative tensors of the interpolated functions, and list several lemmas needed for

2370

WEIMING CAO

deriving the interpolation error estimates. Denote by {TN } a family of triangulations of a polygonal domain Ω ∈ Rn into simplicial elements (see [9]). Here N represents the total number of elements. We study the error estimates for piecewise polynomial interpolations over a class of anisotropic meshes {TN } that are quasi-uniform under a given metric. More specifically, let M be a Riemannian metric on Ω. For each element τ ∈ TN , define Mτ to be the average of M over τ , i.e.,  1 Mτ = M (x)dx. |τ | τ Since Mτ is an n × n symmetric positive definite matrix, we may express it in its eigen-decomposition form, (1)

Mτ = Tτ · Dτ · Tτ ,

where Dτ = diag(d1 , d2 , . . . , dn ) is composed of all the eigenvalues of Mτ , and Tτ is the orthogonal matrix composed of all the eigenvectors. Define (2)

−1

Fτ = Tτ Dτ 2 .

Clearly (Mτ )−1 = Fτ ·Fτ . Let xc be the center of element τ . Define an affine mapping ˜ = Fτ−1 (x − xc ). Then τ is transformed into a simplex element τ˜ with its center at x the origin. We call a family of triangulations {TN } quasi-uniform under metric M if there exist positive constants c1 and c2 independent of N such that (i) for all τ ∈ TN , the smallest internal angle of τ˜ = Fτ−1 (τ − xc ) is bounded from below by c1 , i.e., τ˜ is shape regular; and (ii) maxτ ∈TN |˜ τ | ≤ c2 minτ ∈TN |˜ τ |. Note that if {TN } is quasi-uniform under metric M , then the geometric features of the mesh are determined by M . Specifically, the element size (area/volume) √ is −1/n proportional to | det(M )| , the element aspect ratio is proportional to 1/ d1 : √ √ 1/ d2 : · · · : 1/ dn , and the element orientation (or mesh alignment) is determined by the directions of the eigenvectors of M . Indeed, some anisotropic mesh generators, such as bamg developed by Borouchaki et al. [3] and Hecht [16], take user-specified metrics to control directly the anisotropic behaviors of the meshes. We also would like to point out that not all anisotropic meshes can be considered as quasi-uniform under suitable metrics. For instance, the mesh in Figure 1 is anisotropic. But it cannot be quasi-uniform under any metric, since for a quasiuniform mesh there is at most a fixed number of neighboring elements next to each vertex; otherwise the metric would be singular at certain points. For general meshes in n dimensions, a necessary condition for a family of triangulations being quasiuniform under a metric is that there is an upper limit for the number of neighboring elements at each element interface (i.e., vertex, edge, and face). We call it the limited connectivity condition for anisotropic meshes. In order to derive the continuous form of the interpolation error estimate, we make an assumption on the smoothness of the mesh metric M . Assumption A. Let M be a given Riemannian metric on Ω. There exists δ > 0 such that for any neighborhood Nz of any z ∈ Ω with radius (under metric M ) less than δ, the following is true for all x ∈ Nz : (3)

¯ ξ ≤ ξ · M (x)ξ ≤ c2 ξ · M ¯ξ c1 ξ · M

¯ is the average of M over Nz . where M

∀ξ ∈ Rn ,

INTERPOLATION ERROR ON ANISOTROPIC MESHES

2371

Fig. 1. An example of anisotropic meshes that are not quasi-uniform under any metric.

This assumption means that the metric M is equivalent to its local average over a sufficiently small neighborhood of each point. It basically requires the continuity of M . Indeed, if M is uniformly continuous over Ω, the above assumption holds obviously. Lemma 2.1. Let M be a metric on Ω satisfying Assumption A. Let F (x)F  (x) = [M (x)]−1 be the decomposition of its inverse as given in (2). For any neighborhood Nz of z ∈ Ω with radius less than δ, we have for all x ∈ Nz that √ √ (4) F¯ −1 F (x)≤ 1/ c1 , F −1 (x)F¯  ≤ c2 , and that (5)

¯ )| ≤ det(M (x)) ≤ (c2 )n | det(M ¯ )|, (c1 )n | det(M

¯ −1 is the decomposition of M ¯ −1 , ¯ is the average of M over Nz , F¯ F¯  = M where M and  ·  stands for the 2-norm of matrices. Proof. It is easy to see that  1/2 F −1 (x)F¯ ξ (F¯ ξ) · M (x)(F¯ ξ) = maxn F −1 (x)F¯  = maxn ξ∈R ξ∈R ξ ξ·ξ 1/2 1/2   √ ξ · M (x)ξ ξ · M (x)ξ = maxn ¯ −1 = max ≤ c2 , −1 n ¯ ¯ ξ∈R (F ξ∈R ξ) · (F ξ) ξ · Mξ ¯ξ = 1 and similarly for the second inequality in (4). To show (5), we note that ci ξ · M n (i = 1, 2) and ξ · M (x)ξ = 1 are three ellipsoids in R , with their volumes being n n n π[ 2 ] ¯ )|−1/2 , i = 1, 2, and π[ 2 n] | det(M (x))|−1/2 , respectively. · (ci )− 2 · | det(M Γ(1+ n ) Γ(1+ 2 2) Here [ n2 ] represents the integer part of n2 , and Γ is the Gamma function. Assumption A ¯ ξ = 1 (i = 1, 2). means that ellipsoid ξ · M (x)ξ = 1 lies in between ellipsoids ci ξ · M Thus (5) follows easily from the ordering of their volumes. Because the error for polynomial interpolations of degree k is determined by the (k + 1)th derivative tensor ∇k+1 u of the interpolated functions u, we make an assumption on its anisotropic behavior.

2372

WEIMING CAO

Assumption B. For each x ∈ Ω, there exists a positive definite matrix Q(x) such that (6)

|(ξ · ∇)k+1 u(x)| ≤ |ξ · Q(x)ξ|

k+1 2

∀ξ ∈ Rn .

Note that pk+1 (ξ) ≡ (ξ · ∇)k+1 u(x) is a homogeneous polynomial of ξ of degree k + 1. For ξ = 1, pk+1 (ξ) is simply the (k + 1)th-order directional derivative at x along ξ. The right-hand side is also a homogeneous function of ξ. Thus, geometrically Assumption B is equivalent to saying that at each x, the ellipse/ellipsoid ξ · Qξ = 1 is contained in the level curve/surface |pk+1 (ξ)| = 1 for directional derivatives. Some choices of Q can be made readily. For instance, if one is only interested in isotropic mesh refinement, then we may choose Q(x) = |||∇k+1 u(x)||| · In , where |||∇k+1 u(x)||| = max |(ξ · ∇)k+1 u| ξ=1

is the largest (k + 1)th-order directional derivative at x, and In is the n × n identity matrix. This choice corresponds to defining Q by the largest circle/sphere contained in the level curve/surface. In the case of k = 1 (i.e., linear interpolation), Q can be chosen as Q = abs(∇2 u) + δ · In , where abs(A) = [AT A]1/2 is the symmetric matrix of the same eigenvectors as matrix A, but of eigenvalues equal to the absolute eigenvalues of A. δ is a small positive constant to avoid Q from being degenerate in case ∇2 u(x) becomes singular. Such a Q is called the majorization matrix for the Hessian ∇2 u in Chen, Sun, and Xu [7]. In order to reflect accurately the anisotropic behavior of ∇k+1 u at x, the matrix Q(x) in (6) should be chosen as “small” (in certain sense) as possible. Also, its choice should be invariant under translation and rotation transforms of the coordinates, since the anisotropic behavior is so. The best choice will depend on how Assumption B is used. For minimizing the interpolation error in various norms, we present in section 3.1 an ideal choice of Q. We also present in section 3.2 an algorithm to find an approximate Q that characterizes roughly the anisotropic behavior of ∇k+1 u. Next, we list two lemmas regarding the anisotropic behavior of ∇k+1 u under affine transforms.  Lemma 2.2. Let |α| = αi  for any multi-index α = (α1 , . . . , αn ) of nonnegative integers. For any integer k ≥ 0, |α|=k+1 |∂ α u| is equivalent to |||∇k+1 u|||. Proof. Define ⎧ ⎫

⎨ ⎬

Cα (ξ1 )α1 · · · (ξn )αn

∀Cα ∈ R . P¯k+1 = v(ξ) = ⎩ ⎭ |α|=k+1

The P¯k+1 is a finite-dimensional linear space. It is easy to verify that vP¯k+1 = max |v(ξ)| ξ=1

is a norm on P¯k+1 . Indeed, if vP¯k+1 = 0, then for any ξ = 0, |v(ξ)| = ξk+1 · |v(ξ/ξ)| = 0, thus v ≡ 0 and  · P¯k+1 is a norm.  On the other hand, |α|=k+1 |Cα | is clearly also a norm on P¯k+1 . Hence it is equivalent to vP¯k+1 . Now, for pk+1 (ξ) = (ξ · ∇)k+1 u ∈ P¯k+1 , since its coefficients

2373

INTERPOLATION ERROR ON ANISOTROPIC MESHES

 α Cα are multiples of ∂ α u with fixed positive constants, therefore |α|=k+1 |∂ u| is  equivalent to |α|=k+1 |Cα |, and to maxξ=1 |pk+1 (ξ)| = |||∇k+1 u|||, too. Lemma 2.3. Let x = Fτ (˜ x) + xc be the affine mapping from τ˜ to τ, and define ˜ the gradient operator with respect to x ˜ . Then x) + xc ). Denote by ∇ u ˜(˜ x) = u(Fτ (˜ under Assumption B we have ˜ k+1 u |||∇ ˜(˜ x)||| ≤ Fτ QFτ 

k+1 2

.

˜ k+1 u Proof. Let p˜k+1 (ξ) = (ξ · ∇) ˜(˜ x). By the fact that Fτ is constant, it follows ˜ = F  ∇ and p˜k+1 (ξ) = [ξ · (F  ∇)]k+1 u(x) = from the chain rule for derivatives that ∇ τ τ pk+1 (Fτ ξ). Hence Assumption B implies that ˜ k+1 u |||∇ ˜(˜ x)||| = max |pk+1 (Fτ ξ)| ≤ max |(Fτ ξ) Q(Fτ ξ)| ξ=1

k+1 2

ξ=1

= Fτ QFτ 

k+1 2

.

Finally, in order to minimize the interpolation errors and to select the optimal mesh metrics, we need the following inertial properties for matrix eigenvalues established in [26] by Thompson and an elementary inequality. Lemma 2.4. Let A be an n × n symmetric matrix with eigenvalues α1 ≥ · · · ≥ αn ≥ 0, and let S be a nonsingular matrix with singular values s1 ≥ · · · ≥ sn . Denote the eigenvalues of B = S ∗ AS by β1 ≥ · · · ≥ βn . Then βi+j−n ≥ s2i αj for all 1 ≤ i, j ≤ n with i + j > n. In particular, λmax (B) = β1 ≥ max1≤i≤n (si )2 αn+1−i . Lemma 2.5. Let 0 < λ1 ≤ λ2 ≤ · · · ≤ λn and 0 ≤ α < 1. For t = (t1 , t2 , . . . , tn ), let λi . 1≤i≤n ti

f (t) = |tn |α · max Then f attains its infimum on the set

K = { t ∈ Rn | 0 < t1 ≤ t2 ≤ · · · ≤ tn , and Πni=1 ti = 1 } if and only if λtii = const. for all i. Proof. First we show that f attains its infimum on K. Let {t(k) } be a sequence in K such that lim f (t(k) ) = inf t∈K f (t) < ∞. Clearly {t(k) } must be bounded; (k )



otherwise there exists a subsequence {t(k ) } whose nth components tn → ∞, while (k ) its ith components ti → 0 for some i, which would imply f (t(k) ) → ∞. Therefore, {t(k) } is bounded and has a cluster point t∗ ∈ K with f (t∗ ) = inf t∈K f (t). Next we show that λi λn ≤ ∗ t∗i tn

(7)

∀i.

Suppose otherwise; then there exists m ≤ n − 1 such that λi max ∗ 1≤i≤n ti

=

λm λj > ∗ t∗m tj

∀j ≥ m + 1.

For each j ≥ m + 1, since λm ≤ λj , we must have t∗m < t∗j . Let β > 1 and define t˜ = (t˜1 , t˜2 , . . . , t˜n ) with  1/m ∗ ti if i ≤ m, β ˜ ti = β −1/(n−m) t∗i if i ≥ m + 1.

2374

WEIMING CAO

For β close enough to 1, we have t˜ ∈ K and f (t˜) = β −α/(n−m)   λi λi · max β −1/m · (t∗n )α max ∗ , β −1/(n−m) · (t∗n )α max ∗ 1≤i≤m ti m+1≤i≤n ti < f (t∗ ). Finally we prove that λi = const. t∗i

(8)

∀i.

Note that Πni=1 t∗i = 1. It follows from (7) that  n λn n Πi=1 λi ≤ , t∗n which implies t∗n ≤

 Πni=1

λn λi

1/n .

In particular, “=” in the above inequality holds if and only if (8) is satisfied. Thus ∗

inf f (t) = f (t ) = t∈K

λn (t∗n )α−1

 ≥ λn

λn Πni=1 λi

 α−1 n ,

and “=” holds if and only if (8) is satisfied. 2.2. Error estimate and optimal mesh metrics. We first recall some classical results for the interpolation error estimates under Euclidean metrics. Let k be a positive integer. Denote by Pk the set of all the polynomials of x ∈ Rn of total degree less than or equal to k. Let Πk be an interpolation operator whose restriction on each element preserves Pk . It is well known that on any shape regular element τ , and for any 0 ≤ m ≤ k and p, q ∈ [1, ∞] |u − Πk u|m,p,τ ≤ c|τ |(k+1−m)/n+1/p−1/q |u|k+1,q,τ ,

(9) provided that (10)

W k+1,q (τ ) → C s (τ )

and

W k+1,q (τ ) → W m,p (τ ),

where s is the highest degree of derivatives used in defining the interpolation Πk . See, e.g., Theorem 3.1.5 of [9]. If we further assume that {TN } is a family of quasi-uniform triangulations, i.e., all τ ∈ TN , for all N, are shape regular and max |τ | ≤ c min |τ |,

τ ∈TN

τ ∈TN

then we have globally 1/p  p (11) |u − Πk u|m,p,τ ≤ cN −(k+1−m)/n−1/p+1/q |u|k+1,q,Ω . ∀τ ∈TN

2375

INTERPOLATION ERROR ON ANISOTROPIC MESHES

Now we present the main theorem of this paper on interpolation error estimates for anisotropic meshes. Theorem 2.1. Let M be a Riemannian metric on Ω satisfying Assumption A, and let F (x) F  (x) = (M (x))−1 be the decomposition of its inverse at each x ∈ Ω. Let {TN } be a family of triangulations of Ω that is quasi-uniform under metric M . Let k be a positive integer, and let Πk be an interpolation operator whose restriction on each element preserves Pk . For any function u satisfying Assumption B, any 0 ≤ m ≤ k, and p ∈ [1, ∞] satisfying (10), we have (12) 



1/p |u −

Πk u|pm,p,τ

τ ∈TN

≤ cN

−(k+1−m)/n



(k+1−m)/n  | det(M )|

1 2

F

Ω

−1 mp





F QF 

(k+1)p 2

1/p .

Ω

Furthermore, among all the Riemannian metrics, the optimal bound of the above estimate is attained when M is defined to be Mk+1,m,p ≡ c(λmax (Q)) (k+1−m)p+n | det(Q)|− (k+1−m)p+n · Q. mp

(13)

1

If {TN } is quasi-uniform under Mk+1,m,p , we have  (14)



1/p |u −

Πk u|pm,p,τ

τ ∈TN

≤ cN −(k+1−m)/n  |λmax (Q)| 2 | det(Q)| m

k+1−m 2n

Lnp/[(k+1−m)p+n] (Ω) .

Proof. Consider an element τ ∈ TN . Denote by Mτ the average of M over τ , and let Fτ be defined as in (2). By the fact that (see [9])

|∂ α v(x)| ≤ cFτ−1 m ·

|α|=m



|∂ α v˜(˜ x)|

∀v,

|α|=m

we have  |u − Πk u|pm,p,τ =



|∂ α (u − Πk u)(x)|p d˜ x

τ˜ |α|=m

≤c

|τ | ˜ ku Fτ−1 mp |˜ u−Π ˜|pm,p,˜τ . |˜ τ|

Because τ˜ is shape regular, we have from the classical error estimate (9) that ˜ ku |˜ u−Π ˜|m,p,˜τ ≤ c|˜ τ |(k+1−m)/n |˜ u|k+1,p,˜τ , which implies that |u − Πk u|pm,p,τ ≤ c|τ | |˜ τ |(k+1−m)p/n−1 Fτ−1 mp |˜ u|pk+1,p,˜τ .

2376

WEIMING CAO

It follows from Lemmas 2.2 and 2.3 that k+1 ˜ k+1 u |∂ α u ˜(˜ x)| ≤ c |∇ ˜(˜ x)| ≤ c Fτ Q(x)Fτ  2 . |α|=k+1

Hence

   (k+1)p ˜ ∂x dx |u − Πk u|pm,p,τ ≤ c|τ | |˜ τ |(k+1−m)p/n−1 Fτ−1 mp Fτ Q(x)Fτ  2 · det ∂x τ  (k+1)p ≤ c|˜ τ |(k+1−m)p/n Fτ−1 mp Fτ Q(x)Fτ  2 dx. τ

Now we write Fτ on the right-hand side of the above inequality in terms of M (x) directly. For each x ∈ τ , decompose M (x) into its eigenvalues and eigenvectors as M (x) = T (x) · D(x) · T  (x),

(15)

where D(x) is the diagonal matrix composed of all the eigenvalues d1 (x) ≤ d2 (x) ≤ · · · ≤ dn (x) and T is the orthogonal matrix composed of all the eigenvectors. Define 1 also F (x) = T (x)D(x)− 2 as in (2). Then by Assumption A about the smoothness of M , we have from Lemma 2.1 that Fτ Q(x)Fτ  ≤ (F −1 (x)Fτ ) · (F  (x)Q(x)F (x)) · (F −1 (x)Fτ ) ≤ F −1 (x)Fτ 2 · F  (x)Q(x)F (x) ≤ cF  (x)Q(x)F (x) and Fτ−1  ≤ Fτ−1 F (x) · F −1 (x) ≤ cF −1 (x). Therefore,



F −1 (x)mp F  (x)Q(x)F (x)

τ |(k+1−m)p/n (16) |u − Πk u|pm,p,τ ≤ c|˜

(k+1)p 2

dx.

τ

Summing up the above inequality for all τ ∈ TN , we find that (17) |u − Πk u|pm,p,τ τ ∈TN

 ≤c

(k+1−m)p/n  (k+1)p max |˜ τ| F −1 (x)mp F  (x)Q(x)F (x) 2 dx.

τ ∈TN

Ω

Now we estimate |˜ τ |. Note that τ˜ = Fτ−1 (τ − xc ) and det(Fτ−1 ) = | det(Mτ )|1/2 . By the assumption that {TN } is quasi-uniform under metric M , the sizes of all τ˜’s are of the same order. Hence,    −1 −1 τ | ≤ cN |˜ τ | = cN | det(Fτ−1 )| dx max |˜ τ ∈TN τ ∈TN

τ ∈TN

= cN −1 ≤ cN −1 = cN −1

 τ ∈TN

 τ ∈TN



τ

 τ

| det(Mτ )|1/2 dx | det(M (x))|1/2 dx

τ

| det(M (x))|1/2 dx, Ω

INTERPOLATION ERROR ON ANISOTROPIC MESHES

2377

where in the last inequality we used (5) in Lemma 2.1. Putting the above inequality into (17), we have the error estimate (12). Next, we consider for what metric M the error bound on the right-hand side of (12) is the smallest. We determine M through its eigenvalues and eigenvectors. Since they are independent of each other, we proceed in three steps as follows. First, for any x ∈ Ω and a given set of eigenvalues of M (x), we determine the orthogonal matrix T (x) so that the integrands on the right-hand side of (12) are the smallest possible. Then we determine the ratios among the eigenvalues to further minimize the integrands. Finally, the optimal distribution of det(M ) on Ω is determined such that the error bound in (12) is minimized. First, for fixed x we write Q(x) in its eigen-decomposition form, Q(x) = S(x) · diag(λ1 , λ2 , . . . , λn ) · S  (x), where 0 < λ1 ≤ λ2 ≤ · · · ≤ λn are eigenvalues of Q(x) and S is the orthogonal matrix composed of all its eigenvectors. Applying Lemma 2.4 on the inertia properties of matrix eigenvalues (with A = Q(x), S = F (x), αi = λn+1−i , and si = (di )−1/2 ), we have F  (x)Q(x)F (x) = λmax (F  (x)Q(x)F (x)) ≥

max {λi /di }.

1≤i≤n

Moreover, the equality in the above relation is attained when T (x) = S(x). Thus, we conclude that min F  (x)Q(x)F (x) = max {λi /di }, ∀T

1≤i≤n

and the minimum value is attained at T (x) = S(x). Since F −1 (x) = (dn )1/2 is (k+1)p independent of T , then for the integrand F −1 (x)mp F  (x)Q(x)F (x) 2 , this choice of T (x) is also optimal, with the minimum value (18)

(dn )

mp 2



 (k+1)p 2

max {λi /di }

.

1≤i≤n

Now we consider minimizing (18) with respect to the ratios among eigenvalues of M (x). For fixed det(M ) = Πni=1 di , it follows from Lemma 2.5 that the expression in (18) achieves its minimum if and only if λi /di =

1 μ

∀i = 1, 2, . . . , n,

where μ depends on x and will be determined later. This implies that di = μ · λi for all i, and (19)

M (x) = μ(x) · Q(x).

The minimum value of (18) is (20)

μ−

(k+1−m)p 2

(λn )

mp 2

.

This is the smallest possible value of the integrand in (17) at x for all possible T (x) and different ratios among di (x), i = 1, 2, . . . , n. With this optimal choice, the error

2378

WEIMING CAO

estimate (12) becomes (21) |u − Πk u|pm,p,τ τ ∈TN

≤ cN −(k+1−m)p/n

 n

1

μ 2 | det(Q)| 2

(k+1−m)p/n  (k+1−m)p mp 2 · μ− (λn ) 2 .

Ω

Ω

Finally, we determine the distribution of μ on Ω to minimize the right-hand side of the above error estimate. Let α = (k + 1 − m)p/n, and let n

1

f = μ 2 · | det(Q)| 2 , g = μ−

(k+1−m)p 2

· (λn )

mp 2

.

  Then the integrals on the right-hand side of (21) are [ f ]α [ g]. It follows from H¨ older’s inequality that 

α  1/(α+1)  α/(α+1)  1/(α+1) f · g = f · g α

1

= f α+1 L(α+1)/α · g α+1 Lα+1  α 1 ≥ f α+1 · g α+1 , and the equality holds if and only if f is a constant multiple of g. In this case μ = c(λn ) (k+1−m)p+n | det(Q)|− (k+1−m)p+n , mp

(22)

1

and the metric M becomes M = μ · Q = c(λn ) (k+1−m)p+n | det(Q)|− (k+1−m)p+n · Q. mp

1

With this optimal choice of metric M , we have the smallest error bound  α+1 n 1 p − (k+1−m)p n 2 2 μ | det(Q)| |u − Πk u|m,p,τ ≤ cN (23) τ ∈TN

≤ cN

− (k+1−m)p n

 (λn )

mp 2(α+1)

| det(Q)|

α 2(α+1)

α+1 .

Ω

Let β = form: (24) 

p α+1

=

np (k+1−m)p+n .

Then we may write the above inequality in the following

1/p |u −

Πk u|pm,p,τ

≤ cN

− k+1−m n

 m 2

[ (λn ) | det(Q)|

1/β

α 2p

β

]

Ω

τ ∈TN

= cN −

k+1−m n

This completes the proof of the theorem.

  k+1−m  m  (λn ) 2 | det(Q)| 2n 

Lnp/[(k+1−m)p+n] (Ω)

.

INTERPOLATION ERROR ON ANISOTROPIC MESHES

2379

Remark 2.1. Theorem 2.1 covers the error estimate on isotropic meshes as a 2 special case. Indeed, if we choose Q(x) = c |∇k+1 u(x)| k+1 · In , then the optimal 2p metric Mk+1,m based on Q becomes c |∇k+1 u| (k+1−m)p+n · In , and that {TN } is quasi-uniform under Mk+1,m,p implies it is isotropic. In this case, error estimate (14) is reduced to   p1 k+1−m |u − Πk u|pm,p,τ ≤ c N − n |u|k+1,np/[(k+1−m)p+n],Ω . τ ∈TN

The above error bound is sharper than estimate (11) with q = p, since Lp (Ω) ⊂ np L (k+1−m)p+n due to m ≤ k. Remark 2.2. In the case of k = 1, namely for linear interpolation, we may choose Q = abs(∇2 u)+δ·In , where δ > 0 is a small parameter to avoid Q from being singular. In this case, the optimal metrics and error estimate with m = 0 (i.e., for Lp -error) stated in Theorem 2.1 are identical to those in [7] by Chen, Sun, and Xu. They are also identical to those in [18] by Huang. For the case of k ≥ 2, error estimates and mesh metrics are also derived in [17, 18] based on the sum of the Hessians of (k − 1)th partial derivatives. It is shown in our previous study [5] that metrics based on the sum of the Hessians can be problematic for general anisotropic meshes; see Remark 3 in [5] for details. Remark 2.3. It is easy to see that if {TN } is quasi-uniform under the optimal metric Mk+1,m,p , then estimate (16) for the W m,p (τ )-seminorm of the error is of the same magnitude on each element. In other words, Mk+1,m,p is a matrix which makes the W m,p -error evenly distributed over all elements. Therefore, the optimal metrics and meshes follow also the so-called equidistribution principle. This principle has been used extensively to justify the selection of optimal or nearly optimal meshes; see, e.g., [12, 18, 22]. 3. Measuring the anisotropic behavior of ∇k+1 u(x). 3.1. Definition of anisotropic measure of ∇k+1 u(x). It is seen from Theorem 2.1 that a good interpolation error estimate relies on a proper matrix Q in Assumption B that characterizes the anisotropic behavior of ∇k+1 u at each point. In order to produce as tight as possible an error estimate, we need to measure as accurately as possible the anisotropic behavior of ∇k+1 u. Note that the optimal error k+1−m m bound (14) is determined by the function |λmax (Q(x))| 2 | det(Q(x))| 2n , where Q(x) is a positive definite matrix satisfying Assumption B. Therefore, to produce the tightest error bound, we choose matrix Q(x) in (13) and (14) to be Qk+1,m , the solution of the following minimization problem: (25)

min

Q∈Vk+1 (x)

m

|λmax (Q)| 2 | det(Q)|

k+1−m 2n

,

where Vk+1 (x) is the set of all n × n symmetric positive definite matrices Q that satisfy (26)

|pk+1 (ξ)| = |(ξ · ∇)k+1 u(x)| ≤ (ξ · Qξ)

k+1 2

∀ξ ∈ Rn .

We call Qk+1,m an anisotropic measure of ∇k+1 u(x). It depends not only on k, the interpolation degree, but also on the index m, which is associated with the norm used to measure the error. In the special case of m = 0 (i.e., for Lp -error estimates), the

2380

WEIMING CAO

above problem is reduced to minimizing det(Q) for all Q ∈ Vk+1 (x). Geometrically, for any symmetric positive definite matrix Q, the level curve/surface ξ · Qξ = 1 is an ellipse/ellipsoid of area/volume proportional to | det(Q)|−1/2 , and Q ∈ Vk+1 (x) implies that this level curve/surface is enclosed in that of |pk+1 (ξ)| = 1. Thus the solution Qk+1,m to (25) corresponds to the largest ellipse/ellipsoid (in area/volume) contained in |pk+1 (ξ)| = 1. The above definition of Qk+1,m based on constrained minimization formulation is not quite convenient for theoretical study and practical use. We may change it into an unconstrained problem. For this purpose, we write the eigen-decomposition of each Q ∈ Vk+1 (x) in the form Q = ν · S Λ S.

(27)

1

Without loss of generality, we assume ν = | det(Q)| n and −1 Λ = diag(a1 , a2 , . . . , an−1 , [Πn−1 i=1 ai ] )

n−1 with 0 < a1 ≤ a2 ≤ · · · ≤ an−1 ≤ [Πi=1 ai ]−1 . Clearly, minimization with respect to all Q is equivalent to minimization with respect to all ai , ν > 0, and all n × n orthogonal matrices S. Note that det(Q) = ν n and the largest eigenvalue of Q is −1 λmax (Q) = ν · [Πn−1 . The objective function in (25) is indeed i=1 ai ] m

|λmax (Q)| 2 · | det(Q)|

k+1−m 2n

=ν

k+1 2

−2 · |Πn−1 , i=1 ai | m

and the constraint Q ∈ Vk+1 (x) becomes |pk+1 (ξ)| ≤ [ ν (S  ξ) · Λ(S  ξ) ]

k+1 2

∀ξ ∈ Rn .

This condition is equivalent to |pk+1 (SΛ− 2 η)| ≤ ν 1

k+1 2

· ηk+1

∀η ∈ Rn .

Since pk+1 is a homogeneous polynomial of degree k + 1, it is further equivalent to |pk+1 (SΛ− 2 η)| ≤ ν 1

k+1 2

∀η ∈ Rn with η = 1.

Hence we conclude that the following ν value is optimal:  (28)

ν=

max |pk+1 (SΛ

η=1

− 12

2  k+1

η)|

,

with which the constraint Q ∈ Vk+1 (x) is automatically satisfied. Now the constrained minimization problem (25) is reduced to finding the minimum of (29)

−2 · max |pk+1 (SΛ− 2 η)| (Πn−1 i=1 ai ) m

1

η=1

−1 and all orthogonal with respect to all 0 < a1 ≤ a2 ≤ · · · ≤ an−1 ≤ [Πn−1 i=1 ai ] matrices S. It is possible that the minimization problem (29) does not have a solution, or problem (25) has a solution with some λi = 0. In this case, Q becomes singular, and the optimal metrics Mk+1,m,p based on it will lead to elements of infinitely large aspect

INTERPOLATION ERROR ON ANISOTROPIC MESHES

2381

ratio. In order to avoid such a degenerate mesh metric in practice, we put a cap on the ratios of the eigenvalues of Q. More precisely, we may restrict λmin (Q)/λmax (Q) ≥ δ, where δ ∈ (0, 1] is a user-specified parameter. This requirement is guaranteed when n−1 1 ai ∈ [δ n , δ − n ] for all 1 ≤ i ≤ n − 1. Therefore, we may seek the minimizer to n−1 1 −1 ≤ δ − n and all orthogonal (29) over all δ n ≤ a1 ≤ a2 ≤ · · · ≤ an−1 ≤ [Πn−1 n=1 ] S. Since the set of all these ai ’s and T ’s is compact, and the objective function in (29) is continuous with respect to its variables, a positive definite minimizer Qk+1,m is guaranteed to exist for any specified 0 < δ ≤ 1. Examples. We present two examples to explain the definition of Qk+1,m . First we consider the simplest case, n = 2, k = 1, which corresponds to linear interpolation in R2 . Without loss of generality, suppose   μ1 2 RφT ∇ u = Rφ · μ2 with |μ1 | ≤ |μ2 |, where Rφ is the matrix of rotation by angle φ counterclockwise. Let  Q = ν · Rψ ·

a

 a−1

T Rψ

with 0 < a ≤ a−1 or a ∈ (0, 1]. Then Qk+1,m is defined by the solution (a∗ , ψ∗ ) to the following problem:   m 1 a− 2 · max |p2 (Rψ Λ− 2 η)| . min (30) η=1

a∈(0,1],ψ∈[0,2π]

Note that max |p2 (Rψ Λ− 2 η)| = max 1

η=1

ξ=1

|p2 (Rψ ξ)| 1

Λ 2 ξ2

.

Using the polar coordinates for ξ ∈ R2 , the objective function of (30) can be written m as a− 2 · maxt∈[0,2π] J(a, ψ, t), where (31)

μ cos2 (t − φ + ψ) + μ sin2 (t − φ + ψ)

1

2 J(a, ψ, t) ≡

. 2 2 −1 a · cos t + a · sin t

When μ1 · μ2 ≥ 0, we can show that (32)

max J(a, ψ, t) ≥ max J(a, φ, t)

t∈[0,2π]

t∈[0,2π]

∀ψ ∈ [0, 2π].

Indeed, it is easy to see that

μ cos2 t + μ sin2 t



1 2 −1 max J(a, φ, t) = max

|, |μ2 · a|) 2 = max(|μ1 · a t∈[0,2π] t∈[0,2π] a · cos2 t + a−1 · sin t since J(a, φ, t) is a rational function of cos2 t ∈ [0, 1]. On the other hand,

μ cos2 (φ − ψ) + μ sin2 (φ − ψ)



1 2 max J(a, ψ, t) ≥ J(a, ψ, 0) =

≥ |μ1 · a−1 | a t∈[0,2π]

2382

WEIMING CAO

since |μ1 | ≤ |μ2 |, and similarly  π max J(a, ψ, t) ≥ J a, ψ, φ − ψ + 2 t∈[0,2π]





μ2

≥ |μ2 · a|

=

2 −1 2 a · sin (φ − ψ) + a · cos (φ − ψ)

since a ≤ 1. Therefore, we conclude that   −m − 12 2 min a (33) · max |p2 (Rψ Λ η)| η=1

a∈(0,1],ψ∈[0,2π]

 = min

a∈(0,1]

 m a− 2 · max(|μ1 · a−1 |, |μ2 · a|) ,

where the minimum is attained at ψ = φ. Furthermore, since m ≤ k = 1, the minimizer for the right-hand side of the above equation is 

μ1

a = , μ2 which implies by (28) that ν=

(34) and

 Qk+1,m = Rφ ·

(35)

 |μ1 · μ2 |

|μ1 |

 |μ2 |

RφT = abs(∇2 u).

In the case of μ1 ·μ2 < 0, (32) is not true in general. However, our numerical  calculation shows that in this case the minimum of (30) is still attained at a = |μ1 /{μ2 }| and ψ = φ, which results in the same ν and Q as in (34) and (35). This example confirms from another perspective that the conventional choice is optimal by using the eigenvalues and eigenvectors of the Hessian matrix to characterize its anisotropic behavior in mesh generation and refinement for linear interpolation or linear elements; see [3, 4, 7, 8, 11, 13, 17, 18]. Our second example is for u = 16 xy 2 in R2 and k = 2 (quadratic interpolation). In this example, p3 (ξ) = ξη 2 . Matrix Q can be determined as in (27) with a and ψ being the minimizer to the following problem:   −m − 12 2 min a (36) · max |p3 (Rψ Λ η)| . η=1

a∈(0,1],ψ∈[0,2π]

For each a > 0, it can be shown numerically that max |p3 (Λ− 2 η)| ≤ max |p3 (Rψ Λ− 2 η)| 1

η=1

1

η=1

∀ψ ∈ [0, 2π].

Thus ψ = 0 is an optimal solution, with minimum value a−(m−1)/2 max | cos t·sin2 t| = √ t √ 2 3 −(m−1)/2 2 3 −(m−1)/2 . In this case, problem (36) is reduced to mina∈(0,1] 9 a . For 9 a m = 0, which corresponds to Lp -error estimates, the optimal a is 0; while for m = 1

2383

INTERPOLATION ERROR ON ANISOTROPIC MESHES

(corresponding to H 1 -error estimates), any positive a ≤ 1 is optimal. The optimal ν √ 3 value is 34 . Applying Theorem 2.1 to the quadratic interpolation of u = 16 xy 2 , we see that the 2 L -error bound can be driven to 0 if we increase the aspect ratio of the elements while keeping their area fixed and alignment direction fixed along the x-axis. The H 1 -error k+1 m bound does not change (since |λmax (Q)| 2 | det(Q)| 2 is a constant independent of a), which implies that using highly anisotropic triangles does not help reduce the H 1 -error. This conclusion coincides with the study in [5] based on the exact formula for the quadratic interpolation errors; see Remark 2 in [5]. 3.2. Estimate of the anisotropic measure: A dimension reduction method. The definition of the anisotropic measure for ∇k+1 u in the previous subsection involves a nonlinear minimization (25) with respect to the matrix Q. Solving this problem could be expensive in practice. Here we describe a method to find a suboptimal solution to the minimization problem, and give an approximate Q for the anisotropic measure. Let λ1 , λ2 , . . . , λn be the eigenvalues (in ascending order) of Q, and v 1 , v 2 , . . . , v n be the corresponding eigenvectors. Notice that the objective function in (25), k+1−m m |λmax (Q)| 2 | det(Q)| 2n , is the product of the eigenvalues of Q. We determine an approximate Q by choosing successively (λi , v i ) for i = n, n − 1, . . . , 1, such that each λi is the smallest possible to have constraint (26) hold. More precisely, we first choose λn as 2   k+1 2 max |pk+1 (ξ)| ; λn = |||∇k+1 u||| k+1 = ξ=1

k+1

i.e., |λn | 2 is the largest (k+1)th-order directional derivative of u. The corresponding eigenvector v n is chosen as the unit vector along which the (k + 1)th directional derivative is the largest, i.e., v n = arg max |pk+1 (ξ)|. ξ=1

To determine the rest n − 1 eigenpairs, let v ˜1 , . . . , v ˜n−1 be a set of orthonormal bases for the orthogonal complement of span{v n } in Rn . Then any ξ ∈ Rn can be expressed as ξ=

n−1

ζi v ˜i + zv n

i=1

for some ζ = [ζ1 , . . . , ζn−1 ] ∈ Rn−1 and z ∈ R. Furthermore, let ⎤ ⎡  ⎡ v ˜1 λ1 ⎥ ⎢ .. ⎢ .. Tn−1 = ⎣ . Dn−1 = ⎣ ⎦ · [v 1 , . . . , v n−1 ] , .  v ˜n−1 and  Qn−1 = Tn−1 Dn−1 Tn−1 .

Then ξ · Qξ = ζ · Qn−1 ζ + λn z 2 ,

⎤ ⎥ ⎦, λn−1

2384

WEIMING CAO

and constraint (26) is reduced to 2

n−1  k+1







ζ · Qn−1 ζ ≥ h(ζ, z) ≡ pk+1 ζi v ˜i + zv n

− λn z 2



∀ζ ∈ Rn−1 , ∀z ∈ R.

i=1

This condition can be expressed equivalently as (37)

ζ · Qn−1 ζ ≥ g(ζ) ≡ max h(ζ, z)

∀ζ ∈ Rn−1 .

z∈R

It is easy to verify that g(ζ) is a homogeneous function of ζ. Indeed, for any t ∈ R, g(tζ) = max h(tζ, z) = max h(tζ, tz) = max t2 h(ζ, z) = t2 g(ζ). z∈R

z∈R

z∈R

Therefore, to determine Qn−1 under constraint (37) is similar to the original problem to determine Q under constraint (26), except the former is of one dimension less than the later. We may repeat this process n − 1 times to arrive at a one-dimensional problem, whose solution is ready to work out. In practice, the evaluation of g(ζ) in the constraint (37) can be carried out by checking the critical points of h(ζ, z) in the z direction. Indeed, for given ζ, ∂ v n · ∇)u. pk+1 (ξ) = (k + 1)(ξ · ∇)k (˜ ∂z Thus the critical points (excluding those making pk+1 = 0, which are clearly not among the maxima of h(ζ, z)) satisfy the following equation: 1−k

|pk+1 (ξ)| k+1 (ξ · ∇)k (˜ v n · ∇)u = λn z. It is equivalent to v n · ∇)u (ξ · ∇)k (˜

!k+1

= (λn z)k+1 · (ξ · ∇)k+1 u

!k−1

,

which is a polynomial equation for z of degree k(k + 1). Remark 3.1. For two-dimensional problems (n = 2), the author developed in [5, 6] a method to define a matrix Qk+1 characterizing the anisotropic behavior ∇k+1 u by using the factors of polynomial pk+1 (ξ) for directional derivatives. For the cases of k = 1, 2, it can be shown that Qk+1 given in [5, 6] is equivalent to the matrix Q produced by the dimension reduction algorithm described here. For k ≥ 3, we believe they are still equivalent. However, it is yet to be verified. 4. Numerical results. In this section, we present some numerical results to compare the error in various norms for interpolations based on anisotropic meshes generated with the optimal metric Mk+1,m,p developed in Theorem 2.1. Two-dimensional example. We consider linear and quadratic interpolations of the following function on Ω = [0, 1]2 : (38) u(x, y) = x2 + y 2 + x3 + y 3 + exp(−K (y − d1 (x))2 ) + exp(−K (y − d2 (x))2 ), where K = 10000 and d1 (x) = −x(x − a)/[2(1 − a)] + 1,

with a = 1.25;

d2 (x) = (x − b)(x − c)/[2(1 − b)(1 − c)],

with b = 0.2, c = 1.25.

2385

INTERPOLATION ERROR ON ANISOTROPIC MESHES

Table 1 The error in various norms for linear (in big font) and quadratic (in small font) interpolations of function (38) based on meshes generated under metrics M2,m,p . Ne and Nv represent the total number of elements and nodes, respectively.

Ne ≈ 1, 000

Ne ≈ 4, 000

Ne ≈ 16, 000

Ne

Nv

Metric

eL1

eL2

eL∞

|e|H 1

1030

550

M2,0,1

8.55865e-03

2.25843e-02

5.65950e-01

1.00987e+01

7.00047e-04

3.38558e-03

2.82878e-01

3.12289e+00

1021

542

M2,0,2

1.00819e-02

1.91599e-02

4.45866e-01

8.39336e+00

5.39823e-04

2.22625e-03

1.02203e-01

2.17432e+00

1039

546

M2,0,∞

2.01321e-02

2.80857e-02

2.00845e-01

8.41199e+00

6.22584e-04

1.94362e-03

8.06520e-02

2.00037e+00

1017

535

M2,1,2

2.85871e-02

4.33668e-02

6.68077e-01

7.76100e+00

8.23694e-04

2.50331e-03

1.87979e-01

1.99802e+00

3992

2069 M2,0,1

1.70237e-03

4.38095e-03

1.96297e-01

4.08116e+00

7.02603e-05

3.69513e-04

4.32572e-02

6.95222e-01

4000

2063 M2,0,2

2.10463e-03

3.44617e-03

7.64404e-02

3.41493e+00

5.62782e-05

2.20319e-04

1.82189e-02

5.33361e-01

4040

2068 M2,0,∞

6.66663e-03

8.14892e-03

2.83157e-02

3.08884e+00

8.89956e-05

2.19451e-04

6.35742e-03

4.53298e-01

4002

2047 M2,1,2

1.27443e-02

1.74685e-02

1.58691e-01

2.83627e+00

1.43980e-04

2.67353e-04

1.60196e-02

4.16979e-01

15971 8136 M2,0,1

4.09208e-04

2.80185e-03

4.46053e-01

2.44949e+00

1.91768e-05

1.04155e-03

2.30520e-01

8.87463e-01

16088 8171 M2,0,2

4.91963e-04

8.02090e-04

5.75916e-02

1.54053e+00

8.81760e-06

1.64703e-04

6.68060e-02

2.68500e-01

15994 8088 M2,0,∞

1.79767e-03

2.06070e-03

5.39628e-03

1.37961e+00

1.54309e-05

4.24226e-05

9.44125e-04

1.43448e-01

16004 8088 M2,1,2

4.44656e-03

5.47114e-03

9.68989e-02

1.25693e+00

3.17102e-05

5.80640e-05

9.69865e-03

1.36711e-01

This function u has steep layers around two parabolas y = di (x), i = 1, 2. We calculate exactly all the second and third partial derivatives of u and determine matrix Q for measuring their anisotropic behaviors by using the dimension reduction algorithm described in section 3.2. Then we form the mesh metric Mk+1,m,p according to (13), which is optimal for minimizing an upper bound of the W m,p -norm of the interpolation error e = u − Πk u. The constant multiple c in (13) is used to control the total number of elements. We use the two-dimensional mesh generator bamg (bidimensional anisotropic mesh generator) developed by Borouchaki et al. [3] and Hecht [16] to create the anisotropic meshes. This package accepts a user-defined metric to create and refine an anisotropic mesh that is quasi-uniform under the given metric. We choose the following parameter setting in all our experiments: "-NoRescaling -NbSmooth 5 -hmax 0.02 -hmin 0.0000005 -ratio 0 -nbv 100000 -v 9" In order to make the anisotropic mesh as uniform as possible under metric Mk+1,m,p , we call bamg iteratively with the metrics recalculated over the updated mesh. The final mesh is the one after 20 iterations for all the cases, where there is little change of the mesh and the interpolation error. We consider specifically the linear (k = 1) and quadratic (k = 2) interpolations, and measure the errors in (i) L1 -norm (m = 0, p = 1), (ii) L2 -norm (m = 0, p = 2), (iii) L∞ -norm (m = 0, p = ∞), and (iv) H 1 -seminorm (m = 1, p = 2). These error norms are calculated using numerical quadratures based on 7 and 28 Fekete points for

2386

WEIMING CAO

Table 2 The error in various norms for linear (in small font) and quadratic (in big font) interpolations of function (38) based on meshes generated under metrics M3,m,p . Ne and Nv represent the total number of elements and nodes, respectively.

Ne ≈ 1, 000

Ne

Nv

Metric

eL1

eL2

eL∞

|e|H 1

1026

542

M3,0,1

1.19200e-02

2.02953e-02

4.01216e-01

8.96489e+00

5.47899e-04

2.63100e-03

1.41338e-01

2.59358e+00

1016

534

M3,0,2

1.84082e-02

2.70593e-02

4.59322e-01

8.95022e+00

5.84794e-04

2.15282e-03

1.06775e-01

2.42518e+00

1019 1000 4002 Ne ≈ 4, 000

3963 3991 4013

536 526

M3,0,∞ M3,1,2

2060 M3,0,1 2033 M3,0,2 2041 M3,0,∞ 2051 M3,1,2

16038 8140 M3,0,1 Ne ≈ 16, 000

16070 8136 M3,0,2 16058 8113 M3,0,∞ 15995 8080 M3,1,2

2.92068e-02

4.64891e-02

2.47095e-01

8.85002e+00

8.38699e-04

2.23870e-03

3.84038e-02

2.16885e+00

2.86835e-02

4.47000e-02

3.66967e-01

8.01008e+00

8.41909e-04

2.57482e-03

1.39962e-01

2.16436e+00

2.78079e-03

4.75471e-03

1.63081e-01

4.31546e+00

4.10310e-05

1.86169e-04

1.42471e-02

4.26067e-01

4.95528e-03

6.45569e-03

1.04982e-01

4.07441e+00

4.69086e-05

1.34854e-04

6.66449e-03

3.53092e-01

1.50754e-02

2.12624e-02

6.14308e-02

3.97760e+00

1.46617e-04

2.32106e-04

2.95223e-03

3.50702e-01

1.90713e-02

2.86518e-02

1.09707e-01

3.94508e+00

2.04870e-04

3.38047e-04

1.18358e-02

3.02695e-01

7.17854e-04

1.33199e-03

4.08086e-02

2.46107e+00

4.04626e-06

1.64799e-05

1.14591e-03

8.76529e-02

1.28305e-03

1.70431e-03

2.82953e-02

2.36372e+00

5.01463e-06

1.28541e-05

6.64558e-04

7.34430e-02

4.86946e-03

6.23055e-03

2.06393e-02

2.24020e+00

2.18202e-05

2.97789e-05

2.82230e-04

7.30433e-02

8.42677e-03

1.14039e-02

3.99892e-02

2.21819e+00

4.88933e-05

7.11500e-05

1.58626e-03

6.24718e-02

linear and quadratic interpolations, respectively; see [25]. We list in Tables 1 and 2 the linear and quadratic interpolation errors in four norms, and display in Figures 2 and 3 the anisotropic meshes (of about 4, 000 elements). Several observations are clearly at hand. (a) In all the cases, the smallest W m,p -norm of the interpolation error is obtained when the mesh is generated according to the optimal metric Mk+1,m,p . This indicates not only the optimality of the metric Mk+1,m,p stated in Theorem 2.1, but also that the matrix Q produced by the dimension reduction algorithm characterizes fairly accurately the anisotropic behavior of ∇k+1 u. (b) When the number of elements is quadrupled (i.e., the element length scale is halved), the interpolation error is reduced by a factor 2k+1−m as predicted in the error estimate. (c) The mesh metrics and the anisotropic meshes ideal for linear interpolation are not necessarily good for quadratic interpolation, and vice versa. Three-dimensional example. We present here an example for the quadratic interpolation in three dimensions. Consider the following function on Ω = [0, 1]3 : (39)

u(x, y, z) = x3 + (μ1 y)3 + (μ2 z)3 ,

where μ1 = 10, μ2 = 30. For this function, ∇3 u is constant over all Ω. By using the dimension reduction algorithm, it is easy to find the matrix Q that characterizes the anisotropic behavior of ∇3 u as follows: ⎡ ⎤ 1 ⎦. μ1 Q=⎣ μ2

2387

INTERPOLATION ERROR ON ANISOTROPIC MESHES

1 0.5

0.9

Metric: M

0.8

2,0,1

0.45

0.7 0.4

0.6

0.5 0.35

0.4 0.3

0.3

0.2 0.25

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.25

0.3

0.35

0.4

0.45

0.5

0.25

0.3

0.35

0.4

0.45

0.5

0.25

0.3

0.35

0.4

0.45

0.5

1 0.5

0.9

Metric: M2,0,2

0.8

0.45

0.7 0.4

0.6

0.5 0.35

0.4 0.3

0.3

0.2 0.25

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

1 0.5

0.9

Metric: M2,0,∞

0.8

0.45

0.7 0.4

0.6

0.5 0.35

0.4 0.3

0.3

0.2 0.25

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

1 0.5

0.9

Metric: M2,1,2

0.8

0.45

0.7 0.4

0.6

0.5 0.35

0.4 0.3

0.3

0.2 0.25

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

Fig. 2. Anisotropic meshes that are quasi-uniform under respective metrics M2,m,p and their closeup views.

2388

WEIMING CAO

1 0.5

0.9

Metric: M

0.8

3,0,1

0.45

0.7 0.4

0.6

0.5 0.35

0.4 0.3

0.3

0.2 0.25

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.25

0.3

0.35

0.4

0.45

0.5

0.25

0.3

0.35

0.4

0.45

0.5

0.25

0.3

0.35

0.4

0.45

0.5

1 0.5

0.9

Metric: M3,0,2

0.8

0.45

0.7 0.4

0.6

0.5 0.35

0.4 0.3

0.3

0.2 0.25

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

1 0.5

0.9

Metric: M3,0,∞

0.8

0.45

0.7 0.4

0.6

0.5 0.35

0.4 0.3

0.3

0.2 0.25

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

1 0.5

0.9

Metric: M3,1,2

0.8

0.45

0.7 0.4

0.6

0.5 0.35

0.4 0.3

0.3

0.2 0.25

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

Fig. 3. Anisotropic meshes that are quasi-uniform under respective metrics M3,m,p and their closeup views.

2389

INTERPOLATION ERROR ON ANISOTROPIC MESHES

The optimal mesh metric Mk+1,m,p is a constant multiple of the above matrix. Thus in an anisotropic mesh that is quasi-uniform under Mk+1,m,p all elements must be of about the same volume, and all elements must have approximately the same length scale 1 : μ11 : μ12 in x, y, and z directions. We create meshes of a specified length scale 1 : L1y : L1z in the following way. First we generate a quasi-uniform tetrahedral mesh using a three-dimensional Delaunay generator gmsh developed by Geuzaine and Remacle [14] over a rectangular box [0, 1] × [0, Ly ] × [0, Lz ]. Then it is scaled in y and z directions by factors Ly and Lz to obtain the anisotropic mesh on Ω. When Lx = μ1 and Ly = μ2 , the obtained mesh is quasi-uniform under Mk+1,m,p . This mesh should also be the minimizer for the W m,p -error of interpolation Πk (in this example, Mk+1,m,p is identical for all m, k, and p). We calculate the error norms by using a quadrature formula with 24 points supplied by [10], which is exact for numerical integration of polynomials of degree less than or equal to 6. We vary Ly between 1 ∼ 20 and Lz between 10 ∼ 40 to obtain meshes of different aspect ratios, while keeping the total number ofelements around 40, 000 (by setting the characteristic length lc in gmsh to be 0.05 3 Ly Lz ). We display in Figure 4 the error contour plots against Ly and Lz in the cases of 60

0.0

0.019

0.004

0.335

(a): L1−norm 0.019

0.177 0.604 0.335 0.089

50

60

0.005

0.00 0.019 0.008

0.008

0.029

0.029

0.7 1.13 0.1260.013

50

0.019

0.0 0 0.013 0.0 0.029

0.029

(b): L2−norm 0.029

0.416 0.236

0.029 0.013

0.005

0.063 0.335 40

0

0

40

0.013 0.004

0.008

0.002

0.002 Lz

Lz

0.042 30

30

0.008 0.003 0.002

0.003

0.236

20

20

0.236

0.019 10

0.029

10 0.042 0.0890.177 1.756 4.502

1.756 0

0

60

5

10

15

5.415

0.063 0.1260.236

0.335 0.604 2.851

20 Ly

25

1.048 30

4.502 6.9 10.46 35

40

0

5

10

60

(c): L∞−norm

15

2.488 292.15 181.64 107.52

5.415 1.432 5.415 50

2.665 3.924 5.645 7.954

2.665 0

50

0.416 1.13 1.762

30

35

(d): H1−seminorm

40

2.488 2.488

2.488

2

451.04 59.947

2.488

2.488

1.074

2

0.692 0.636

0.297

5.645 7.9 11 14.96

40

0.015

2.86

25

2.488

0.023

5.415 40

20 Ly

0.7

0.05

30

Lz

Lz

0.023 30

0.0150.012

31.056

0.05

14.709

2.86 0.023

0.674 20

20 1.432

59.947

0.123 0.297

10

2.86 5.41516.988 0

10

0.674 1.432

0

5

73.004 112.51 10

15

20 Ly

46.217 25

30

6.291

35

40

14.709 31.056 107.52 181.64 292.15 451.04 1898.6974.13

59.947

9.787 46.217 11 73.004 73.00428.46 169.5 0

0

5

10

15

20 Ly

25

30

672.72570.74482.5 34 1375.2 1898.6 35

40

Fig. 4. Contour plot of the (a) L1 , (b) L2 , (c) L∞ , and (d) H 1 errors with respect to various element aspect ratios 1 : L1 : L1 . The star-circle mark at (10, 30) indicates the optimal aspect ratio y z predicted by the error estimates.

2390

WEIMING CAO

m = 0, p = 1, 2, ∞, and m = 1, p = 2. It is noted that the smallest interpolation error is achieved when the element aspect ratio is approximately equal to 1 : μ11 : μ12 in all four cases. This indicates the optimality of the metric Mk+1,m,p . It is also noted that the error is relatively insensitive to variation of the length scales in y and z directions when they are close to the optimal values. We believe this is because the anisotropic behavior of ∇3 u in this example is relatively “mild.” For interpolated functions of stronger anisotropic behaviors, the improvement by using the optimal mesh metrics can be more drastic. Due to the lack of a reliable anisotropic mesh generator in three dimensions, we are unable to test the optimality of Mk+1,m,p for u with variable ∇3 u. 5. Conclusion and discussions. In the previous sections, we presented an error estimate for higher-order interpolations over anisotropic meshes in Rn . It involves an interplay of the mesh features controlled by a given mesh metric and the anisotropic measures of the higher-order derivative tensors of interpolated functions. Based on the error estimate, we were able to identify the optimal mesh metrics leading to the smallest error bound for a subset of interpolated functions exhibiting similar anisotropic behaviors. Numerical results indicate that the meshes generated according to the optimal metrics produce the smallest interpolation error in the corresponding norms. A critical component in applying the error estimate for anisotropic mesh generation or refinement is to measure the anisotropic behavior of higher-order derivative tensors. We define such a measure by the “largest” ellipse/ellipsoid contained in the level curve for directional derivatives. To avoid solving the minimization problem for defining the anisotropic measure, we developed a dimension reduction algorithm to produce the measure approximately. The practical application of the results in this paper is always associated with the development of a reliable and efficient anisotropic mesh generator. While there have been many two-dimensional packages available (bamg is clearly among the best of them), general three-dimensional anisotropic meshing packages are yet to be developed and tested. Also, it is natural to apply the results in this paper to quadratic and higher-order finite element methods for solving PDEs. To this end, one needs to recover the higher-order derivatives of the PDEs solution from its numerical approximation. While there have been extensive studies along this direction for isotropic triangulations and linear elements [1, 20, 27, 28], the analysis and application for higher-order elements on anisotropic meshes are yet to be developed. REFERENCES [1] M. Ainsworth and J. T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Wiley Intersciences, New York, 2000. [2] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., Teubner, Stuttgart, 1999. [3] H. Borouchaki, P. L. George, F. Hecht, P. Laug, and E. Saltel, Delaunay mesh generation governed by metric specifications. I. Algorithms, Finite Elem. Anal. Des., 25 (1997), pp. 61–83. [4] W. Cao, On the error of linear interpolation and the orientation, aspect ratio, and internal angles of a triangle, SIAM J. Numer. Anal., 43 (2005), pp. 19–40. [5] W. Cao, Anisotropic measure of third order derivatives and the quadratic interpolation error on triangular elements, SIAM J. Sci. Comput., 29 (2007), pp. 756–781. [6] W. Cao, An interpolation error estimate in R2 based on the anisotropic measures of higher order derivatives, Math. Comp., to appear. [7] L. Chen, P. Sun, and J. Xu, Optimal anisotropic meshes for minimizing interpolation errors in Lp -norm, Math. Comp., 79 (2007), pp. 179–204.

INTERPOLATION ERROR ON ANISOTROPIC MESHES

2391

[8] L. Chen and J. Xu, Optimal Delaunay triangulations, J. Comput. Math., 22 (2004), pp. 299–308. [9] P. G. Ciarlet, The Finite Element Methods for Elliptic Problems, Classics Appl. Math. 40, SIAM, Philadelphia, 2002. [10] R. Cools, An encyclopaedia of cubature formulas, J. Complexity, 19 (2003), pp. 445–453. [11] E. F. D’Azevedo and R. B. Simpson, On optimal triangular meshes for minimizing the gradient error, Numer. Math., 59 (1991), pp. 321–348. [12] C. de Boor, Good approximation by splines with variable knots II, in Conference on the Numerical Solution of Differential Equations, Lecture Notes Math. 363, Springer-Verlag, Berlin, 1974, pp. 12–20. [13] L. Formaggia and S. Perotto, New anisotropic a priori error estimates, Numer. Math., 89 (2001), pp. 641–667. [14] C. Geuzaine and J.-F. Remacle, Gmsh: A Three-Dimensional Finite Element Mesh Generator with Built-in Pre- and Post-Processing Facilities, Technical report, Dept. of Electrical Engineering and Computer Sciences, University of Li`ege, Li` ege, Belgium, 2005. [15] W. G. Habashi, J. Dompierre, Y. Bourgault, D. Ait-Ali-Yahia, M. Fortin, and M. Vallet, Anisotropic mesh adaptation: Towards user-independent, mesh-independent and solver-independent CFD. I. General principles. Internat. J. Numer. Methods Fluids, 32 (2000), pp. 725–744. [16] F. Hecht, Bidimensional Anisotropic Mesh Generator (User’s Manual), INRIA, Rocquencourt, 1997, http://www-rocqq.inria.fr/gamma/cdrom/www/bamg/eng.htm [17] W. Huang, Measuring mesh qualities and application to variational mesh, SIAM J. Sci. Comput., 26 (2005), pp. 1643–1666. [18] W. Huang, Mathematical principles of anisotropic mesh adaptation, Commun. Comput. Phys., 1 (2006), pp. 276–310. [19] E. J. Nadler, Piecewise Linear Approximation on Triangulations of a Planar Region, Ph.D. thesis, Division of Applied Mathematics, Brown University, Providence, RI, 1985. [20] A. Naga and Z. Zhang, A posteriori error estimates based on the polynomial preserving recovery, SIAM J. Numer. Anal., 42 (2004), pp. 1780–1800. [21] S. Rippa, Long and thin triangles can be good for linear interpolation, SIAM J. Numer. Anal., 29 (1992), pp. 257–270. [22] R. D. Russell and J. Christiansen, Adaptive mesh selection strategies for solving boundary value problems, SIAM J. Numer. Anal., 15 (1978), pp. 59–80. [23] J. R. Shewchuk, What Is a Good Linear Finite Element? Interpolation, Conditioning, Anisotropy, and Quality Measure, Preprint, Dept. of Electronic Engineering and Computer Sciences, University of California, Berkeley, CA, 2002. [24] R. B. Simpson, Anisotropic mesh transformations and optimal error control, Appl. Numer. Math., 14 (1994), pp. 183–198. [25] M. A. Taylor, B. A. Wingate, and R. E. Vincent, An algorithm for computing Fekete points in the triangle, SIAM J. Numer. Anal., 38 (2000), pp. 1707–1720. [26] R. C. Thompson, Inertial properties of eigenvalues. II, J. Math. Anal. Appl., 58 (1977), pp. 572–577. [27] Z. Zhang and A. Naga, A new finite element gradient recovery method: Superconvergence property, SIAM J. Sci. Comput., 26 (2005), pp. 1192–1213. [28] O. C. Zienkiewicz and J. Z. Zhu, The superconvergence patch recovery and a posteriori error estimates, Part I: The recovery technique, Internat. J. Numer. Methods Engrg., 33 (1992), pp. 1331–1364.