SIAM J. SCI. COMPUT. Vol. 29, No. 2, pp. 756–781
c 2007 Society for Industrial and Applied Mathematics
ANISOTROPIC MEASURES OF THIRD ORDER DERIVATIVES AND THE QUADRATIC INTERPOLATION ERROR ON TRIANGULAR ELEMENTS∗ WEIMING CAO† Abstract. The main purpose of this paper is to present a closer look at how the H 1 - and L2 -errors for quadratic interpolation on a triangle are determined by the triangle geometry and the anisotropic behavior of the third order derivatives of interpolated functions. We characterize quantitatively the anisotropic behavior of a third order derivative tensor by its orientation and anisotropic ratio. Both exact error formulas and numerical experiments are presented for model problems of interpolating a cubic function u at the vertices and the midpoints of three sides of a triangle. Based on the study on model problems, we conclude that when an element is aligned with the orientation of ∇3 u, the aspect ratio leading to nearly the smallest H 1 - and L2 -norms of the interpolation error is approximately equal to the anisotropic ratio of ∇3 u. With this alignment and aspect ratio taken, the H 1 -seminorm of the error is proportional to the reciprocal of the anisotropic ratio of ∇3 u, the L2 -norm of the error is proportional to the − 32 th power of the anisotropic ratio, and both of them are insensitive to the internal angles of the element. Key words. anisotropic mesh, quadratic interpolation, anisotropic ratio, aspect ratio, mesh alignment AMS subject classifications. 65D05, 65L50, 65N15, 65N50 DOI. 10.1137/050634700
1. Introduction. In anisotropic mesh generation and refinement, elements or cells can be of different length scales in different directions in order to better fit the behavior of the interpolated functions or the properties of the PDEs to be solved. There are two anisotropic characters of a triangular element. One is the orientation, which is roughly the direction of its longest side. The other is the aspect ratio, which measures how long and thin the triangle is. For linear interpolation over triangular meshes, the error is determined by the Hessian ∇2 u of the interpolated function u. There have been extensive studies on the optimal shape and orientation for minimizing the linear interpolation error in various norms. For instance, Nadler [16] studied the model problem of linear interpolation of quadratic functions and concluded that the optimal triangle for minimizing the L2 -error is an equilateral under the Hessian metric. Rippa [17] concluded that to minimize the Lp -error for linear interpolation, the triangular elements should be long and thin in the direction where the second order directional derivative of u is small. D’Azevedo and Simpson [11, 12] made similar conclusions for the maximum error and maximum derivative error in linear interpolation. They also developed several strategies for the anisotropic mesh generation and refinement based on ∇2 u [11, 19]. Shewchuk [18] summarized various formulas for the linear interpolation error on a triangle in terms of its area and edge lengths. In [6] the author expresses the H 1 - and L2 -errors for the linear interpolation of a quadratic function u on a triangle K in terms of the size, orientation, and aspect ratio of K, and ∗ Received
by the editors June 28, 2005; accepted for publication (in revised form) October 25, 2006; published electronically April 10, 2007. This work was supported in part by the NSF through grant DMS-0209313. http://www.siam.org/journals/sisc/29-2/63470.html † Department of Mathematics, University of Texas at San Antonio, San Antonio, TX 78249 (
[email protected]). 756
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
757
the eigenvalues and eigenvectors of ∇2 u. It is identified that among all possible triangles of a fixed area, the smallest H 1 -seminorm of the interpolation error is attained with the acute isosceles triangle that is aligned with the eigenvector (associated with the smaller eigenvalue) of ∇2 u, and has an aspect ratio 0.8|λ2 /λ1 |, where |λ1 | ≤ |λ2 | are the eigenvalues of ∇2 u. For general triangular elements aligned with ∇2 u without the restriction on the maximum internal angles, the best aspect ratio for the H 1 -error is about |λ2 /λ1 |. This is also the best aspect ratio for minimizing the L2 -error. Meanwhile, there have been quite a few papers on generalizing the classical interpolation error estimates on quasi-uniform meshes to those on anisotropic meshes. For instance, Formaggia and Perotto [13], Huang [14], Huang and Sun [15], Chen, Sun, and Xu [8], and Chen and Xu [9] derived various estimates for the W m,p -norm of the error for linear interpolation of general functions. They conclude similarly that the mesh leading to the smallest error bound is the one where all elements are aligned with the eigenvector of ∇2 u, and their aspect ratios are approximately equal to the square root of the ratio of the eigenvalues of ∇2 u. Chen, Sun, and Xu [8] further proved that their estimate of the Lp -error is optimal for the class of functions with positive Hessians. In the case of quadratic interpolation, the understanding of the interplay between the geometric feature of elements and the anisotropic feature of interpolated functions has been mostly qualitative and crude [18]. As pointed by Apel in [1], a more precise understanding of their effects on the interpolation errors relies on “sufficiently fine description of the properties” of the interpolated functions. We introduced in [7] two quantities to measure the anisotropic characters of third order and higher derivative tensors. One is the orientation, along which the absolute value of the higher order directional derivative is about the smallest, while along its perpendicular direction it is about the largest. Another is the anisotropic ratio, which measures the strength of the anisotropic behavior, or how differently the directional derivative behaves in different directions. These measures are invariant under rotation and translation of the coordinates. We also derived an interpolation error estimate in terms of these anisotropic measures and identified the optimal alignment direction and aspect ratio which minimize an upper bounds of the error estimate; see [7] for details. In this paper, we take a closer look at how the H 1 - and L2 -errors for quadratic interpolation on a triangle are determined by the triangle geometry and the anisotropic features of ∇3 u. Since any smooth function can be approximated locally by a cubic function to a higher order accuracy than by a quadratic function, we simplify our local error analysis over an element by assuming the interpolated functions u are cubics. Unfortunately, unlike the case of linear interpolation studied in [6], a closed formula for the quadratic interpolation error in terms of the magnitude, orientation, and anisotropic ratio of ∇3 u turns out to be too complicated. Nevertheless, for isosceles triangles that are aligned with ∇3 u, we are able to get some handy formulas for the H 1 - and L2 -norms of the error. It is found from these formulas that when a triangle is aligned with the orientation of ∇3 u, the optimal aspect ratio for minimizing both H 1 - and L2 -norms of the quadratic interpolation error is approximately equal to the anisotropic ratio of ∇3 u. This is true regardless of whether the element is an acute or an obtuse isosceles triangle. Therefore, in contrast to the case of linear interpolations, there seems to be no “superaccuracy” in favor of acute isosceles triangles for quadratic interpolations of general cubic functions; see [18]. For general triangle shape and orientations, we study numerically the effects of element alignment, aspect ratio, and internal angles on the H 1 - and L2 -norms of the interpolation errors. It is found that
758
WEIMING CAO
(i) the interpolation error is nearly the smallest when an element is aligned with the orientation of ∇3 u; (ii) when an element is aligned with ∇3 u, the optimal aspect ratio is about equal to the anisotropic ratio of ∇3 u. With this aspect ratio, the H 1 seminorm of the interpolation error is proportional to the magnitude of ∇3 u and to the reciprocal of the anisotropic ratio of ∇3 u. The L2 -norm of the error is proportional to the magnitude of ∇3 u and to the − 32 th power of the anisotropic ratio of ∇3 u. (iii) When proper alignment and proper aspect ratio are maintained, the H 1 -error is insensitive to the internal angles of the element. Otherwise, it increases as the largest internal angle of the element increases. In this situation the maximum angle condition is essential. The L2 -error is always insensitive to the internal angles. An outline of this paper is as follows: We first explain in sections 2 and 3 how to measure the orientation and aspect ratio of a triangle, and how to define the orientation and anisotropic ratio of ∇3 u. Then we derive in section 4 the exact formulas for the H 1 - and L2 -norms of the error for quadratic interpolation of cubic functions on a triangle. In section 5 we study the special case with isosceles triangles that are aligned with the orientation of ∇3 u. In this case the H 1 - and L2 -errors can be expressed in simple formulas, and optimal aspect ratios can be identified. For general triangular elements we present a numerical study in section 6 on the effects of the element shape and orientation. We conclude this paper with a few discussions. 2. Anisotropic character of a triangular element. We first recall the definiˆ be tions of orientation and aspect ratio for a triangular element described in [6]. Let K √ √ 3 3 1 1 the standard element with vertices ξ 1 = (0, 1), ξ 2 = (− 2 , − 2 ), ξ 3 = ( 2 , − 2 ). Let K be a general triangular element with vertices x1 , x2 , and x3 . Let xc = 13 (x1 +x2 +x3 ) ˆ to K can be expressed as x = F ξ +xc be the center of K. The affine mapping from K with 1 F = √ (x1 − x3 ), x2 − xc . 3 Denote by F = U ΣV ∗ the singular value decomposition (SVD) of matrix F . Without loss of generality, we assume Σ = diag(σ1 , σ2 ) with σ1 ≥ σ2 > 0, and U = Rφu and V = Rφv are the matrices for rotation transform by angles φu and φv (counterclockwise), respectively. ˆ into a triangle centered at the origin. It The mapping x = F ξ transforms K also transforms the unit circle on the ξη-plane into an ellipse on the xy-plane, with σ1 and σ2 being its semimajor and semiminor axes; see Figure 1. It is easy to see that φu determines the direction of the semimajor axis of the ellipse. We define this direction as the orientation of triangle K. We define the aspect ratio of K as r = σσ12 . These definitions are equivalent to the commonly used definition of orientation based on the longest side of K and the definition of aspect ratio based on the diameter of K and the radius of the inscribed circle of K. The angle φv in the SVD of mapping F determines the internal angles of K. In particular, K is an obtuse isosceles triangle when φv = 0 and an acute isosceles triangle when φv = π/2. If K is an obtuse/acute isosceles triangle, then its orientation is parallel/perpendicular to its base. 3. Anisotropic measures of third order derivatives. By Taylor’s expansion, any smooth function u can be approximated locally by lower order polynomials, u(x, y) =
m 1 [(x − xc ) · ∇]j u(xc ) + O(x − xc m+1 ). j! j=0
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
(a)
759
(b) φ
v
1
1
(c)
(d)
σ
2
σ
φ
u
1
ˆ (b) K ˆ under R∗ ; (c) K ˆ under ΣR∗ ; (d) K ˆ under F = Fig. 1. (a) standard element K; φv φv ∗ . Rφu ΣRφ v
Let Πm−1 be an interpolation operator that preserves the set of polynomials of total degree ≤ m − 1. Then the leading term of the interpolation error u − Πm−1 u is equal 1 to the interpolation error for polynomial m! ((x − xc ) · ∇)m u(xc ), which is controlled m by the partial derivative tensor ∇ u(xc ). To understand the local behavior of the interpolation error around xc , we need to measure the anisotropic behavior of ∇m u. Recall that there are three quantities, the area, orientation, and aspect ratio, to characterize the geometric features of a triangular element. We introduced in [7] three quantities to characterize the anisotropic behavior of ∇m u. The first one is its magnitude, which is equivalent to the Frobenius norm of ∇m u. The second is the orientation, along which direction the absolute value of the mth order directional derivative at xc is about the smallest, while along its perpendicular direction it is about the largest. The third quantity is the anisotropic ratio, which measures the strength of the anisotropic behavior of ∇m u at xc . These quantities are determined by the distribution of the mth order directional derivatives. They are invariant under the rotation and translation transforms of the independent variables of u. We describe below in detail these definitions in the cases of m = 2, 3. 3.1. Anisotropic measure of ∇2 u. We start from the case of ∇2 u, which will be needed for introducing the anisotropic measures of ∇3 u. Note that ∇2 u(xc ) is a 2 × 2 symmetric matrix; it admits the eigendecomposition λ1 2 ∇ u = Rφ RφT , (3.1) λ2 where |λ1 | ≤ |λ2 |, and Rφ is the matrix for rotation with angle φ ∈ [− π2 , π2 ] counterclockwise. The two eigenvectors of ∇2 u are v 1 = [cos φ, sin φ]T and v 2 = v ⊥ 1 , which are in the directions of angle φ and φ + π2 from the x-axis. When |λ1 /λ2 | is small, the absolute value of the second order directional derivative is about the smallest along v 1 , and it is the largest along its perpendicular direction. We define v 1 as the orientation of ∇2 u at xc . Furthermore, we define the anisotropic ratio of ∇2 u as λ2 2 S(∇ u) ≡ S2 ≡ λ1
760
WEIMING CAO
and define the magnitude of ∇2 u as |∇2 u| ≡ D2 ≡
1 (|λ1 | + |λ2 |). 2
Clearly, S2 ≥ 1. The larger S2 is, the stronger the anisotropic behavior of ∇2 u. Since 1 2 2 |λ2 | ≤ D2 ≤ |λ2 |, the magnitude |∇ u| is equivalent to the 2-norm, and hence the 2 Frobenius norm, of ∇ u. The above definition of the orientation, anisotropic ratio, and magnitude for ∇2 u can be interpreted geometrically in terms of the level curves of the following polynomial: p2 (ξ, η) =
1 (ξ · ∇)2 u(xc ), 2
where ξ = [ξ, η]T . For ξ = 1, p2 (ξ, η) is one half of the second order directional derivative along ξ for u at xc . We first consider the case λ1 λ2 < 0 (assume, more specifically, λ1 ≥ 0 and λ2 ≤ 0.) Let (r, t) be the polar coordinates of (ξ, η). Then by the eigendecomposition (3.1) for ∇2 u, λ1 (RφT ξ) p2 (ξ, η) = 12 (RφT ξ) · λ 2 (3.2) = 12 r2 [λ1 cos2 (t − φ) + λ2 sin2 (t − φ)].
Let α = tan−1 ( | λλ12 |). Then α ∈ [0, π4 ] for |λ1 | ≤ |λ2 |, and sin2 α =
λ1 , |λ1 | + |λ2 |
cos2 α = −
λ2 . |λ1 | + |λ2 |
Therefore, p2 (ξ, η) = 12 r2 (|λ1 | + |λ2 |) [sin2 α cos2 (t − φ) − cos2 α sin2 (t − φ)] = 12 r2 (|λ1 | + |λ2 |) sin(α + t − φ) · sin(α − t + φ). Clearly, p2 (ξ, η) = 0 along the lines of angles φ + α + kπ and φ − α + kπ from the ξ-axis. Let β1 ≤ β2 be the two angles among them that lie in (−π/2, π/2]. Then we can express p2 as p2 (ξ, η) = ± 12 r2 (|λ1 | + |λ2 |) sin(t − β1 ) sin(t − β2 ) (3.3)
= ± 12 r2 (|λ1 | + |λ2 |) (cos t sin β1 − sin t cos β1 ) · (cos t sin β2 − sin t cos β2 ) = ±D2 1 (ξ, η) 2 (ξ, η),
where D2 = 12 (|λ1 | + |λ2 |), and (3.4)
i (ξ, η) = ξ sin βi − η cos βi ,
i = 1, 2.
We are ready to see that the orientation of ∇2 u is just the bisection direction of the smaller sector formed by the level-0 lines of p2 (ξ, η); see Figure 2(a). The angle φ in (3.1) for this direction can be expressed as φ=
1 ¯ (β + β), 2
761
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
case (a)
case (b) l1(ξ,η)=0
l*2(ξ,η)=0 l*1(ξ,η)=0
l2(ξ,η)=0
α α
α
φ β2
β1
β1
α β2
φ
Fig. 2. Contour lines of p2 (ξ, η) and the orientation (bold line) of ∇2 u. (a) p2 (ξ, η) has two linear factors. Here β = β2 and β¯ = β1 + π. (b) p2 (ξ, η) does not change sign for all (ξ, η). ∗1 and ∗2 are the linear factors of the associated polynomial p∗2 . Here β = β1 and β¯ = β2 .
¯ is the smaller interval between [β1 , β2 ] and [β2 , β1 + π]. Furthermore, the where [β, β]
anisotropic ratio S(∇2 u) = | λλ21 | can be expressed in terms of the opening β¯−β = 2α of the smaller sector as S(∇2 u) = cot α = cot((β¯ − β)/2). The magnitude |∇2 u| of ∇2 u is just the coefficient D2 in the factoring (3.3) of p2 . Algebraically, λ1 λ2 ≤ 0 corresponds to the case where p2 (1, z) has two real roots zi = tan(βi ). Now we consider the case λ1 λ2 > 0. In this case, p2 (1, z) has two conjugate complex roots, and p2 (ξ, η) does not change sign for all (ξ, η). In this case, we introduce matrix ˆ = Rφ λ 1 RφT (3.5) H −λ2 ˆ see Figure 2(b). Clearly, p∗2 and define an associated polynomial p∗2 (ξ, η) = 12 (ξ · Hξ); can be factored into the product of two linear functions ∗1 and ∗2 as in (3.4). Since ˆ the magnitude, orientation, and anisotropic ratio of ∇2 u are identical to those of H, ∗ we can determine them equivalently by the level-0 lines for p2 as in the case with λ1 λ2 < 0. In summary, if p2 (ξ, η) can be factored into the product of two linear functions, then the orientation and anisotropic ratio of ∇2 u can be determined by the level-0 lines of p2 (ξ, η). Otherwise, they can be determined equivalently by the level-0 lines ˆ in (3.5). In particular, the level-0 lines of p2 or p∗ are of H of p∗2 (ξ, η) defined by 2 −1 angles ±α = ± tan ( |λ1 /λ2 |) from the eigenvector v 1 of ∇2 u. 3.2. Anisotropic measure of ∇3 u . To characterize the anisotropic behavior of ∇3 u at a point xc , we study p3 (ξ) =
1 (ξ · ∇)3 u(xc ). 3!
For ξ = 1, p3 (ξ) is one sixth of the third order directional derivative along ξ for u at xc . Since p3 (ξ) is an homogeneous polynomial, we may write p3 (ξ, η) = ξ 3 p3 (1, η/ξ),
762
WEIMING CAO
with p3 (1, z) being a real polynomial of degree ≤ 3. p3 (1, z) can be factored into the product of linear polynomials and/or a positive quadratic polynomial of z. Therefore, we may factor p3 (ξ, η) into either the product of three linear functions or the product of a linear function and a positive, homogeneous, quadratic function of ξ, η. We discuss the former case first. Assume (3.6)
p3 (ξ, η) = D3 1 (ξ, η) 2 (ξ, η) 3 (ξ, η),
where D3 ≥ 0, and each i is a linear function of the form (3.4) with −
π π < β1 ≤ β2 ≤ β3 ≤ . 2 2
These factors can be easily determined by the roots of p3 (1, z) = 0. Indeed, tan βi ’s are just the roots of p3 (1, z) = 0. Clearly, i (ξ, η) = 0, i = 1, 2, 3, are the level-0 lines of p3 . We define the magnitude of ∇3 u as |∇3 u| = D3 . This is an upper bound of all the third order directional derivatives at xc . It is not difficult to show that |∇3 u| is equivalent to the Frobenius norm of ∇3 u [7]. Moreover, when all three level0 lines of p3 are close to a common direction, the absolute value of the third order directional derivative of u is about the smallest along this direction, while it is about the largest along its perpendicular direction. In this case, ∇3 u is strongly anisotropic. To determine the common direction to which the level-0 lines are all close, we consider the three sectors formed by every other line of i = 0 for i = 1, 2, 3; see Figure 3(a). The ranges of the central angle of these sectors are, respectively, (3.7)
[β1 , β3 ],
[β2 , β1 + π],
[β3 , β2 + π].
Let [β, β] be the shortest interval among the three of (3.7). Clearly, β − β ≤ 2π 3 . If 1 β−β 1, all the level-0 lines of p3 (ξ, η) are close to the direction of angle φ = 2 (β+β) from the ξ-axis. We define this direction as the orientation of ∇3 u. Furthermore, the smaller β − β is, the stronger the anisotropic behavior of ∇3 u. Hence, we define the anisotropic ratio of ∇3 u as S3 ≡ S(∇3 u) ≡ max(1, cot((β¯ − β)/2)). The artificial cutoff of S(∇3 u) at the value 1 reflects the fact that we consider only stretching an element in two perpendicular directions in anisotropic mesh refinements. When β¯ − β ≥ π2 , rescaling an element in two perpendicular directions (while keeping the area fixed) has little advantage in reducing the H 1 -norm of the interpolation error; see Remark 2 in section 5 below. Therefore, we simply consider in this case that ∇3 u is isotropic and its anisotropic ratio S(∇3 u) = 1. Next we study the case where p3 (ξ, η) is the product of a linear function and a positive, homogeneous, quadratic function. Assume p3 (ξ, η) = 3 (ξ, η) q(ξ, η), where 3 (ξ, η) is in the form of (3.4), and q(ξ, η) > 0 for all ξ = 0. We write q = 12 ξ T Qξ, where Q is the Hessian of q. Decompose Q into (3.8)
Q=R
λ1
λ2
RT ,
763
QUADRATIC INTERPOLATION ERROR ON TRIANGLE l3(ξ,η)=0
case (a)
l3(ξ,η)=0
case (b)
l*2(ξ,η)=0
l2(ξ,η)=0
α
α β2
β3
α
φ
α β2
β1
β3
φ
β1
l*1(ξ,η)=0
l1(ξ,η)=0
l*2(ξ,η)=0
case (c)
l*1(ξ,η)=0
α
α β1
β2
φ
β3
l3(ξ,η)=0
1 Fig. 3. Contour lines of p3 = 3! (ξ · ∇)3 u and the orientation (bold line) of ∇3 u. (a) p3 has three linear (real) factors. (b) p3 has a linear factor 3 (ξ, η) and a positive quadratic factor q(ξ, η). 3 = 0 lies “inside” the level-0 lines of q ∗ (ξ, η). (c) p3 has a linear factor 3 (ξ, η) and a positive quadratic factor q(ξ, η). 3 = 0 lies “outside” the level-0 lines of q ∗ (ξ, η).
where R is an orthogonal matrix, and 0 < λ1 ≤ λ2 . As explained in the previous ˆ subsection, the anisotropic feature of q(ξ, η) is identical to that of q ∗ (ξ, η) = 12 ξ · Qξ, where λ1 ˆ RT . (3.9) Q=R −λ2 Therefore, we define the magnitude, orientation, and anisotropic ratio of ∇3 u based on (3.10)
p∗3 (ξ, η) = 3 (ξ, η) q ∗ (ξ, η);
see Figure 3(b) and 3(c). Since p∗3 (ξ, η) can be factored into the product of three linear functions ∗1 , ∗2 , and 3 , these anisotropic measures can be determined by the position of the three level-0 lines of p∗3 . An important feature about the above definitions of magnitude, anisotropic ratio, and orientation for ∇3 u is that they are invariant under translation and rotation transforms of the independent variables of u. Remark 1. The above definitions of the anisotropic measures can be generalized 1 to the case of ∇m u for m > 3; see [7]. Indeed, since pm (ξ, η) = m! (ξ · ∇)m u is
764
WEIMING CAO
an homogeneous polynomial, we may write pm (ξ, η) = ξ m pm (1, η/ξ), with pm (1, z) being a real polynomial of degree ≤ m. By the fundamental theorem of algebra, pm (1, z) can be factored into the product of linear functions and/or positive quadratic functions of z. Therefore, pm (ξ, η) can be factored into the product of linear functions and/or positive, homogeneous, quadratic functions of ξ, η. We may convert each positive quadratic function q(ξ, η) into q ∗ (ξ, η). Then the magnitude, orientation, and anisotropic ratio of ∇m u can be determined by the positions of the level-0 lines of p∗m (ξ, η). 3.3. Examples for anisotropic measures. We give several examples to demonstrate how to calculate the orientation and anisotropic ratio of ∇3 u. Assume is a small positive number. Example 1. u(x, y) = y 3 + ( x)3 . Here p3 (ξ, η) =
1 (ξ · ∇)3 u = η 3 + ( ξ)3 = (η + ξ)(η 2 − ξη + 2 ξ 2 ), 3!
which is the product of a linear function and a positive quadratic function. Let q(ξ, η) = η 2 − ξη + 2 ξ 2 . Then 2 2 ∇ q= − 2
− . 2
The eigenvalues of ∇2 q are √ λ1 = ( 2 + 1) − 4 − 2 + 1 ≈ ( 2 + 1) − (1 − 2 /2) = 3 2 /2, √ λ2 = ( 2 + 1) + 4 − 2 + 1 ≈ ( 2 + 1) + (1 − 2 /2) = 2 + 2 /2. The √ eigenvector of ∇2 q, associated with the smaller eigenvalue λ1 , is v 1 = [1, ( 2 − 1 + 4 − 2 + 1)/ ]T ≈ [1, /2]T . Note that the level-0 lines of q ∗ (ξ, η) are along the angles φ ± α from the x-axis, where φ ≈ tan−1 ( 2 ) is the angle of v 1 and α =
tan−1 | λλ21 |. Thus in this case the two level-0 lines of q ∗ (ξ, η) are along the angles
√ 3 32 tan−1 ( 2 ) ± tan−1 ( 4+ 2 ) ≈ 2 ± 2 . Combining the direction of linear factor y + x = 0, we have the three level-0 lines of p∗3 (ξ, η) in approximately the directions of √ √ √ 1− 3 1+ 3 1+ 3 ¯ angle − , 2 and 2 , which implies β ≈ − and β ≈ − 2 . Therefore, the orientation of ∇3 u is √ 1 −1 + 3 ¯ φ = (β + β) ≈ , 2 4 and the anisotropic ratio is S(∇3 u) = cot( 12 (β¯ − β)) ≈
4√ 1 . 3+ 3
Example 2. u(x, y) = (x + y)3 + (x3 + y 3 ). Here p3 (ξ, η) = (ξ + η)3 + (ξ 3 + η 3 ) = (ξ + η) · [(1 + )ξ 2 + (2 − )ξη + (1 + )η 2 ],
765
QUADRATIC INTERPOLATION ERROR ON TRIANGLE 0
10
−1
10
[S3]−1 −2
1
10
−3
10
φ+π/4
3
−4
10
−5
10 −10 10
−8
−6
10
10
−4
ε
10
−2
0
10
Fig. 4. Example 3: The anisotropic ratio S3 and orientation φ +
π 4
10
versus parameter .
which is the product of a linear function and a positive quadratic function. Let q(ξ, η) = (1 + )ξ 2 + (2 − )ξη + (1 + )η 2 . Then 2 + 2 2 − ∇ q= . 2 − 2 + 2 2
Its two eigenvalues are 3 and 4 + , and the eigenvector associated with the smaller ∗ eigenvalue is [1, −1].
Therefore, the two level-0 lines of q (ξ, η) are along angles
tan−1 (−1) ± tan−1
3 4+ .
Together with the level-0 line x + y = 0 due to the linear
factor, we are ready to see that the orientation of ∇3 u is along angle φ = tan−1 (−1) = − π4 , and the anisotropic ratio is S(∇3 u) = cot [ tan−1 (
3 4+ ) ]
≈
√2 . 3
Example 3. u(x, y) = (x+y)3 + (x3 −y 3 ). In this case, p3 (ξ, η) = (ξ +η)3 + (ξ 3 − η ), which is again the product of a linear function and a positive quadratic function. However, the asymptotic analysis for the level-0 lines of p∗ (ξ, η) is rather complicated. Instead we use a numerical method to determine the angles of the orientation and the anisotropic ratio; see Figure 4. It turns out that the orientation
is still close to the 3
line y + x = 0, but the anisotropic ratio is of the magnitude O( 3
1 ).
4. Quadratic interpolation error. In this section, we study the model problem of interpolating a cubic function by quadratic polynomials at the three vertices and the midpoints of three sides of a triangular element K. We derive the exact formulas for the H 1 -seminorm and L2 -norm of the interpolation error. Let x1 , x2 , and x3 be the vertices of K (numbered counterclockwise), and let 1 = x3 − x2 , 2 = x1 − x3 , 3 = x2 − x1 . Let ci (x, y), i = 1, 2, 3, be the barycentric coordinates of a point (x, y) in K. We call a triple index (i, j, k) in cyclic order if it is one of the
766
WEIMING CAO
following three: (1, 2, 3), (2, 3, 1), (3, 1, 2). For (i, j, k) in cyclic order, define (4.1)
b0 = c1 c2 c3
bi = cj ck (cj − ck ),
and
i = 1, 2, 3.
We are ready to see that bi ’s are linearly independent. Let u be a cubic function, and let Π2 u be its quadratic interpolation at the vertices and side midpoints of K. Then the error e = u − Π2 u can be expressed as e=
(4.2)
3
vi bi (x, y).
i=0
To determine the constants vi , 0 ≤ i ≤ 3, we first list a few basic formulas. Their proofs are elementary. Lemma 4.1. Let τ i = i /|i |, and let ni be the unit outward normal on the side i . Then we have ∇ci = −
(4.3)
(4.4)
⎧ 0 ⎪ ⎪ ⎪ ⎪ ⎨ −1 τ i · ∇cj = | i | ⎪ ⎪ ⎪ ⎪ +1 ⎩ | i |
|i | ni , 2|K|
if i = j, if (i, j) = (1, 2), (2, 3), (3, 1), if (i, j) = (2, 1), (3, 2), (1, 3).
Lemma 4.2.
(4.5)
⎧ 3 ∂ bi 12 ⎪ ⎪ ⎨ ∂τ 3 = |i |3 i
i = 1, 2, 3,
3 ⎪ ⎪ ⎩ ∂ b3i = 0 ∂τ j
∀i = j.
Note that Π2 u is of degree 2, and we get from the above two lemmas that vi =
(4.6)
|i |3 ∂ 3 u , 12 ∂τ 3i
i = 1, 2, 3.
Next we determine the coefficient v0 . It is natural to consider the symmetric 3 differentiation ∂ τ 1∂τ 2 τ 3 on both sides of (4.2). However, since (4.7)
∂ 3 b0 =0 ∂τ 1 τ 2 τ 3
and
∂ 3 bi 4 , = ∂τ 1 τ 2 τ 3 |1 | · |2 | · |3 |
i = 1, 2, 3,
this operation leads only to an identity ∂3u ∂3u 1 = | i |3 3 . ∂τ 1 τ 2 τ 3 3|1 | |2 | |3 | i=1 ∂τi 3
(4.8)
Therefore, we consider taking the following operation on both sides of (4.2): (4.9)
L=
∂3 ∂3 ∂3 + + . ∂τ 21 τ 2 ∂τ 22 τ 3 ∂τ 23 τ 1
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
Lemma 4.3. Let ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩
767
(i, j, k) be in cyclic order. Then ∂ 3 b0 = −2 , ∂τ 2i τ j |i |2 · |j |
∂ 3 b0 = +2 , ∂τ 2i τ k |i |2 · |k |
∂ 3 bi = −6 , ∂τ 2i τ j |i |2 · |j |
∂ 3 bi = +2 , ∂τ i τ 2j |i | · |j |2
∂ 3 bi = −6 , ∂τ 2i τ k |i |2 · |k |
∂ 3 bi = +2 , ∂τ i τ 2k |i | · |k |2
∂ 3 bi = −2 , ∂τ 2j τ k |j |2 · |k |
∂ 3 bi = −2 . ∂τ j τ 2k |j | · |k |2
We have from the above formulas that Lb0 =
(4.10)
−2 −2 −2 + + 2 2 |1 | |2 | |2 | |3 | |3 |2 |1 |
and (4.11)
Lbi =
−6 −2 +2 + + , |i |2 |j | |j |2 |k | |k |2 |i |
i = 1, 2, 3.
Now the coefficient v0 can be solved from Lu = v0 Lb0 +
3
vi Lbi ,
i=1
which yields v0 =
(4.12)
Lu −
3
vi Lbi
/Lb0 .
i=1
Now we derive the formulas for the H 1 -seminorm and L2 -norm of the interpolation error. Clearly, (4.13) |∇e|2 dxdy = v · B1 v, K
where v = [v0 , v1 , v2 , v3 ]T , and B1 = K ∇bi · ∇bj dxdy . By the identity [10] α1 !α2 !α3 ! (c1 )α1 (c2 )α2 (c3 )α3 dxdy = 2|K| , (α + α2 + α3 + 2)! 1 K we have (4.14) ⎡ 2 |1 | + |2 |2 + |3 |2 ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
B1 = 2(|3 |2 − |2 |2 )
1 720|K| 2(|1 |2 − |3 |2 )
⎤
⎥ ⎥ ⎥ ⎥. 2 ⎥ 4(2|2 | − 31 · 3 ) −42 · 3 ⎦ 4(2|3 |2 − 31 · 2 )
4(2|1 |2 − 32 · 3 ) −41 · 2
sym.
2(|2 |2 − |1 |2 )
−41 · 3
768
WEIMING CAO
For the L2 -norm of the interpolation error, we similarly have (4.15) |e|2 dxdy = v · B0 v, K
where (4.16)
B0 =
⎡
bi bj dxdy
=
K
|K| ⎢ ⎢ 7! ⎣
2
0 6
0 −2 6
sym.
⎤ 0 −2⎥ ⎥. −2⎦ 6
5. Optimal aspect ratios: A case study. Since error formulas (4.13) and (4.15) involve complicated interaction between the geometric features of K and the orientation and anisotropic ratio of ∇3 u, it is rather difficult, if not impossible, to decipher from them some quantitative guidelines for anisotropic mesh generation. In particular, among triangles that are of unit area and aligned with the orientation of ∇3 u, we would like to identify those which produce the smallest quadratic interpolation errors. We study this issue in the case when the triangle is an acute or obtuse isosceles triangle. Because the error is invariant under the rotation transform of the coordinates, we may assume the orientation of ∇3 u is simply the x-axis. Furthermore, since we are mainly concerned with the local behavior of the error, we assume u is an homogeneous cubic function. In this case, p3 (ξ) = u(ξ) is either the product of three linear functions or the product of a linear function and a positive quadratic function. Suppose the anisotropic ratio of ∇3 u is 1 . By the definition of the orientation and anisotropic ratio, u must be a multiple of one of the following three types of functions. Case (i). u is the product of three linear functions. In this case u must be a multiple of (5.1)
u1 (x, y) =
1 (y − δx)(y 2 − 2 x2 ), 3!
where |δ| ≤ ; see Figure 3(a). Case (ii). u is the product of a linear function 3 and a positive quadratic function q, with the linear factor “lying inside” the quadratic factor (see Figure 3(b)). In this case, the orientation of ∇3 u is determined by the directions of ∗1 and ∗2 in q ∗ . Hence, u is a multiple of (5.2)
u2 (x, y) =
1 (y − δx)(y 2 + 2 x2 ), 3!
where |δ| ≤ . Case (iii). u is the product of a linear function 3 and a positive quadratic function q, with the linear factor “lying outside” the quadratic factor (see Figure 3(c)). In this case, the orientation of ∇3 u is determined by the directions of 3 and one of the factors
∗1 or ∗2 in q ∗ . Therefore, u must be a multiple of (5.3) where
u3 (x, y) =
1 (y − x)q(x, y), 3!
2 sin γ q(x, y) = (Rγ x) · 0 T
0 cos2 γ
(RγT x),
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
769
( ). where Rγ is the matrix of rotation by angle γ, and |γ| ≤ tan−1√ 2 The magnitude of the first two types of functions is (1 + ) 1 + δ 2 , while for the √ 2 third type it is 1 + . Under the assumption 1, we have for all three cases the magnitude being approximately 1. We shall simplify the error formulas (4.13) and (4.15) for the first two types of functions when the element is an acute or obtuse isosceles triangle aligned with the ∇3 u. In this case, optimal aspect ratios leading to the smallest interpolation errors can be identified. The error formulas for the third type of function are still too complicated. For this case we study the optimal aspect ratio numerically in the next section. 5.1. H 1 -seminorm of the interpolation error. First, we present the formula for the H 1 -seminorm of the interpolation error for function (5.1). We first consider element K being an acute isosceles triangle aligned with the orientation of ∇3 u, i.e., φu = 0, φv = π2 . Denote by |K| the area of K, and denote by r its aspect ratio. With the help of Maple, we have in this case √ √ 4|K| |K| |K| − 3r 0 3r , 3 = √ 1 = √ , , 2 = √ 1 1 3r −1 3r 3r and
(5.4)
⎡ √ ⎢ ⎢ 3⎢ B1 = 540 ⎢ ⎢ ⎣
3 2 (r
+ 1r )
3 2 (r
0 5r +
9 r
− 1r )
2 r
6r +
8 r
sym.
⎤ − 1r ) ⎥ 2 ⎥ r ⎥ ⎥. 1 ⎥ 3r − r ⎦ 3 2 (r
6r +
8 r
By using formulas (4.6)–(4.12) in section 4, we have the coefficients vi , 0 ≤ i ≤ 3, as v0 = 0, √ 4
v1 = − 2 9 3 |K|3/2 r−3/2 , 33/4 3/2 (3r2 2 108 |K|
√
3) r−3/2 , √ 3/4 v3 = − 3108 |K|3/2 (3r2 2 − 1)(3δr + 3) r−3/2 , v2 =
− 1)(3δr −
and the H 1 -seminorm of the interpolation error e = u − Π2 u is ∇e2L2 (K) =
|K|3 2 4 12960 [r (9
+ 54δ 2 ) + (7 4 + 12δ 2 2 − 6 2 + 16δ 2 )
+ r−2 (33 + 6δ 2 + 6 2 ) + 15r−4 ]. Now consider the optimal choice of the aspect ratio r for a given anisotropy ratio −1 . It is easy to verify that the minimum of ∇eL2 (K) is attained when r solves the equation (3 + 18δ 2 )r6 4 − (11 + 2δ 2 + 2 2 )r2 − 10 = 0. Use the formula for roots of cubic polynomials, we have the only positive solution of the above equation, 11 + 2δ 2 + 2 2 2 θ 2 r = 2 , cos 2 3 3 1 + 6δ
770
WEIMING CAO
where
θ = cos
−1
When |δ| ≤ 1, we have θ ≈ ratio r∗ is approximately
2
45 π 2
r∗ ≈
1 + 6δ 2 (11 + 2δ 2 + 2 2 )3
and r2 ≈ 4
11 −2 . 3
.
Therefore, the optimal aspect
11 −1 ≈ 1.383 −1 . 3
With this aspect ratio, the interpolation error is ∇e∗ 2L2 (K) ≈ (0.002196 2 + 0.001235δ 2 )|K|3 . Next we consider the case where K is an obtuse isosceles triangle aligned with the orientation of ∇3 u, i.e., φu = 0, φv = 0. Simplifying error formula (4.13), we end up with ∇e2L2 (K) =
|K|3 4 2 4 12960 [15r δ
+ r2 (33δ 2 4 + 6 2 δ 2 + 6 4 )
+ (−6 2 δ 2 + 16 4 + 12 2 + 7δ 2 ) + r−2 (54 + 9δ 2 )]. The critical points of the right-hand side function are the solutions of 10r6 δ 2 4 + r4 (11δ 2 4 + 2 2 δ 2 + 2 4 ) − (18 + 3δ 2 ) = 0. Let z = r−2 ; then z satisfies the cubic equation z3 −
11δ 2 4 + 2 2 δ 2 + 2 4 10δ 2 4 z = . 18 + 3δ 2 18 + 3δ 2
By studying the positive solution of the above equation for all |δ| ≤ 1, we conclude that the optimal value of the aspect ratio is in the range 1.047 −1 ≤ r∗ ≤ 1.732 −1 , and the optimal interpolation error is in the range 0.003703|K|3 2 ≤ ∇e∗ 2L2 (K) ≤ 0.007673|K|3 2 . We may study the H 1 -error for the second type of function (5.2) in a similar way. The error formulas and the optimal aspect ratios are summarized in Table 5.1. A noticeable observation is that the optimal aspect ratio for both types of functions is proportional to the anisotropic ratio −1 of ∇3 u, and the H 1 -seminorm of the error is proportional to . This is also true for the third type of function (5.3); see the numerical study of this case in the next section. Thus we conclude that when an element is aligned with the orientation of ∇3 u, the optimal aspect ratio that leads to the smallest H 1 -seminorm of the quadratic interpolation errors is about cS3 , with c ∈ (1.047, 1.732). With this direction and aspect ratio taken, the interpolation error ∇e2L2 (K) is of the magnitude c(D3 )2 (S3 )−2 |K|3 for c ∈ (0.001852, 0.007673). Here D3 and S3 are the magnitude and anisotropic ratio of ∇3 u defined in section 3.
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
771
Table 5.1 H 1 -seminorm of the quadratic interpolation error for functions (5.1) and (5.2) when element K is an isosceles triangle aligned with ∇3 u. |K| and r are the area and aspect ratio of the triangle, respectively.
u=
1 (y 3!
− δx)(y 2 − 2 x2 ):
K = acute isosceles ∇e2L2 (K) =
|K|3 [r2 4 (9 12960
+ 54δ 2 ) + (74 + 12δ 2 2 − 62 + 16δ 2 ) + r−2 (33 + 6δ 2 + 62 ) + 15r−4 ],
r∗ ≈ 1.383 −1 , ∇e∗ 2L2 (K) ≈ (0.0021962 + 0.001235δ 2 )|K|3 . K = obtuse isosceles ∇e2L2 (K) =
|K|3 [15r 4 δ 2 4 12960
+ r2 (33δ 2 4 + 62 δ 2 + 64 ) + (−62 δ 2 + 164 + 122 + 7δ 2 ) + r −2 (54 + 9δ 2 )],
r∗ ≈ c−1 , with c ∈ (1.047, 1.732), ∇e∗ 2L2 (K) ≈ c|K|3 2 , with c ∈ (0.003703, 0.007673).
u=
1 (y 3!
− δx)(y 2 + 2 x2 ):
K = acute isosceles ∇e2L2 (K) =
|K|3 [r2 4 (9 12960
+ 54δ 2 ) + (74 − 12δ 2 2 + 62 + 16δ 2 ) + r−2 (33 + 6δ 2 − 62 ) + 15r−4 ],
r∗ ≈ 1.383 −1 , ∇e∗ 2L2 (K) ≈ (0.0031222 + 0.001235δ 2 )|K|3 . K = obtuse isosceles ∇e2L2 (K) =
|K|3 [15r 4 δ 2 4 12960
+ r2 (33δ 2 4 − 62 δ 2 + 64 ) + (62 δ 2 + 164 − 122 + 7δ 2 ) + r −2 (54 + 9δ 2 )],
r∗ ≈ c−1 , with c ∈ (1.103, 1.743), ∇e∗ 2L2 (K) ≈ c|K|3 2 , with c ∈ (0.001852, 0.004752).
Another important observation is that the optimal aspect ratio is about the same for both acute and obtuse isosceles triangles. This is in contrast to the case of linear interpolations, where a much higher accuracy can be achieved for acute isosceles triangles with a much larger aspect ratio than that for obtuse isosceles triangles. This phenomenon is termed as “superaccuracy” by Shewchuk in [18]. The study in this subsection indicates that there is no “superaccuracy” in favor of acute isosceles triangles for quadratic interpolation of general cubic functions. 5.2. L2 -norm of the interpolation error. For the L2 -norm of the interpolation error, we may study similarly the optimal choice of the aspect ratio r for the first two types of functions (5.1) and (5.2). The results are summarized in Table 5.2 for
772
WEIMING CAO
Table 5.2 L2 -norm of the quadratic interpolation error for functions (5.1) and (5.2) when element K is an isosceles triangle aligned with ∇3 u. |K| and r are the area and aspect ratio of the triangle, respectively.
u=
1 (y 3!
− δx)(y 2 − 2 x2 ):
K = acute isosceles e2L2 (K) =
√ 3|K|4 [274 δ 2 r 3 272160
+ 3(2 + 2δ 2 )2 r + (11δ 2 − 102 )r −1 + 19r −3 ],
r∗ ≈ c−1 , with c ∈ (0.9265, 1.732), e∗ 2L2 (K) ≈ c|K|4 3 , with c ∈ (0.00001959, 0.0003486). K = obtuse isosceles e2L2 (K) =
√ 3|K|4 [194 δ 2 r 3 272160
+ (112 − 10δ 2 )2 r + 3(δ 2 + 22 )r −1 + 27r −3 ],
r∗ ≈ c−1 , with c ∈ (1.079, 1.732), e∗ 2L2 (K) ≈ c|K|4 3 , with c ∈ (0.0001764, 0.0003486).
u=
1 (y 3!
− δx)(y 2 + 2 x2 ):
K = acute isosceles e2L2 (K) =
√ 3|K|4 [274 δ 2 r 3 272160
+ 3(2 − 2δ 2 )2 r + (11δ 2 + 102 )r −1 + 19r −3 ],
r∗ ≈ c−1 , with c ∈ (1.000, 2.517), e∗ 2L2 (K) ≈ c|K|4 3 , with c ∈ (0.00008092, 0.0004073). K = obtuse isosceles e2L2 (K) =
√ 3|K|4 [194 δ 2 r 3 272160
+ (112 + 10δ 2 )2 r + 3(δ 2 − 22 )r −1 + 27r −3 ],
r∗ ≈ c−1 , with c ∈ (1.000, 1.567), e∗ 2L2 (K) ≈ c|K|4 3 , with c ∈ (0.0001300, 0.0004073).
the cases φu = 0, φv = π2 (aligned acute isosceles triangle) and φu = 0, φv = 0 (aligned obtuse isosceles triangle) separately. When the triangle is aligned with the orientation of ∇3 u, the aspect ratio leading to the smallest L2 -error is approximately equal to cS3 , with c ∈ (0.9264, 2.517). With this aspect ratio, e2L2 (K) is of the magnitude c(D3 )2 (S3 )−3 |K|4 for c ∈ (1.959 × 10−5 , 4.072 × 10−4 ). It is interesting to note that when a triangle is aligned with ∇3 u, the optimal aspect ratios for both the H 1 -seminorm and L2 -norm are close to the anisotropic ratio of ∇3 u. In practice, an aspect ratio r ≈ S(∇3 u) can be used generally. Remark 2. We emphasize the fact that the optimal aspect ratios listed in Tables 5.1 and 5.2 are derived under the assumption that K is aligned with the orientation of ∇3 u defined in section 3. This orientation direction is only an estimate
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
773
of the direction along which the absolute value of third order directional derivative is about the smallest, while along its perpendicular direction it is about the largest. It is close to, but may not be exactly, the triangle orientation leading to the minimum interpolation error among all possible choices. Therefore, for elements aligned with the orientation of ∇3 u we may expect only a nearly optimal interpolation error. Unlike the case of linear interpolation studied in [6], we do not know precisely in what alignment direction the quadratic interpolation error would be minimum (among all possible triangle shapes and orientations). Furthermore, it seems that the optimal choice of this direction is dependent on the norm used to measure the error. For instance, if we consider the L2√-norm of the interpolation error for √ 11 3 3 3 function u(x, y) = 61 xy 2 , then e2L2 (K) = 272160 |K|4 r−1 and 272160 |K|4 r−1 on an acute/obtuse isosceles triangle K aligned with the x-axis, respectively. In this case, the x-axis and a large r are good choices for the alignment and aspect ratio. On the other hand, if we consider the H 1 -seminorm of the error for the above function, we 1 1 have ∇e2L2 (K) = 12960 |K|3 (16 + 6r−2 ) and 12960 |K|3 (7 + 9r−2 ) on an acute/obtuse isosceles triangle aligned with the x-axis, respectively. In this case, using triangles of high aspect ratios helps little to reduce the H 1 -error. We may consider in this case ∇3 u isotropic with no preferred alignment direction. Remark 3. In a recent work studying mesh quality measures and variational mesh generation methods [14], Huang suggested use the following matrix to determine the element orientation and aspect ratio in the case of quadratic interpolation: H(∇u) = |∇2 (∂x u)| + |∇2 (∂y u)|, where |A| = (AT A)1/2 for a matrix A. More precisely, let (λ1 , v 1 ) and (λ2 , v 2 ) √ be the λi in eigenpairs of H(∇u). He suggested choosing the element length scale to be 1/ the v i direction, i = 1, 2, i.e., to use v 1 as element orientation and λ2 /λ1 as element aspect ratio (assume λ1 ≤ λ2 ). We believe H(∇u) cannot accurately characterize the anisotropic behavior of ∇3 u, and the guidelines based on H(∇u) may not be helpful. While there is no numerical result in [14] for quadratic interpolation, we tested this choice for the function u = y 3 + ( x)3 , i.e., our Example 1 in section 3.3. In this case 3 6 0 H(∇u) = . 0 6 Its eigenvalues are λ1 = 6 3 , λ2 = 6, and the corresponding eigenvectors are along the x-axis and y-axis. Suppose K is an isosceles triangle aligned with the x-axis and of aspect ratio r; then we may derive the exact formula of the interpolation error. Indeed, it is found that |K|3 ∇e2L2 (K) = 12960 (33r−2 + 15r−4 + 54 6 r2 ), √ 3|K|4 2 eL2 (K) = 272160 (19r−3 + 27 6 r3 ) for acute isosceles triangle K and that |K|3 ∇e2L2 (K) = 12960 [54r−2 + 6 (15r4 + 33r2 )], √ 3|K|4 e2L2 (K) = 272160 (27r−3 + 19 6 r3 )
774
WEIMING CAO
for obtuse isosceles triangle K. If r is chosen as λ2 /λ1 = −3/2 as determined in [14], then ∇e2L2 (K) = O(1) · |K|3 for obtuse isosceles triangle K. This H 1 -error is in the same magnitude as that over a shape regular triangle of the same area. In other words, using obtuse anisotropic elements with such an aspect ratio helps little to reduce the H 1 -error. For the L2 -error, r = −3/2 leads to e2L2 (K) = O( 3/2 ) · |K|4 for both acute and obtuse isosceles triangles K. They are a factor of O( −3/2 ) larger than with the choice suggested in this paper. 6. Effects of element orientation, aspect ratio, and internal angles: A numerical study. We present in this section a numerical study on the effects of the element alignment, aspect ratio, and internal angles in quadratic interpolations. Let K be a triangle of unit area but in general shape and orientation. Since, under rotation of the coordinates, any cubic function can be written as a multiple of the three types of functions (5.1)–(5.3), we study only the interpolation errors over K for these functions. Note that they are all of anisotropic ratio −1 and approximately of magnitude 1. We test with various values of in (0, 1). We also experiment with different δ values for functions (5.1) and (5.2) and with various γ for |γ| ≤ tan−1 ( ) for function (5.3). We calculate exactly the H 1 -seminorm and L2 -norm of the error u−Π2 u over K by a numerical quadrature for triangles based on 28 Fekete points described in [20]. This quadrature formula is exact for integrands of total degree up to 6, which is sufficient to evaluate precisely the H 1 - and L2 -error for interpolation of cubic functions. We study mainly (1) the effects of the triangle aspect ratios when the element is aligned with ∇3 u; (2) the effects of the element alignments; and (3) the effects of the element internal angles. Because the numerical results for functions (5.1) and (5.2) are very similar, we report only those for (5.1). In particular, we report the results for (5.1) with δ = 0 and δ = . The errors with general |δ| ≤ are in between these two cases. Our main observations are the following: (1) Optimal r: If the triangle is aligned with the orientation of ∇3 u, then the optimal aspect ratios to produce the smallest interpolation error in both L2 - and H 2 norms are close to the anisotropic ratio of ∇3 u. We plot in Figure 5 the H 1 -seminorm and L2 -, L1 -, and L∞ -norms of the error versus the aspect ratio of the element for test functions of different anisotropic ratios. The linear dependency of the optimal r on the anisotropic ratio −1 is clearly observed. In addition, with the optimal aspect ratio, the H 1 -seminorm of the interpolation error is proportional to , and the L2 -, L1 -, and L∞ -norms are proportional to 1.5 . (2) Optimal φu : We study the effects of the element alignment by varying the angle φu in the SVD of the affine mapping for it; see section 2. We plot in Figure 6 the L2 -norm and H 1 -seminorm of the error for test functions of various anisotropic ratios. In all the cases, the smallest interpolation error occurs when φu is close to 0, i.e., when the element is oriented in about the same direction of ∇3 u. Furthermore, this conclusion is true regardless of whether the element is an acute or obtuse isosceles triangle. It indicates that mesh alignment is crucial for minimizing the interpolation error. (3) Optimal φv : The internal angles of a triangular element are controlled by φv in the SVD of its affine mapping. When φv = mπ/3, the triangle is an obtuse isosceles triangle, and when φv = (2m + 1)π/6, it is an acute isosceles triangle. We plot in Figures 7 and 8 the L2 - and H 1 -errors versus φv for elements with φu = 0 and r = −1 , 10 −1 . When the aspect ratio r = −1 , the L2 - and H 1 -errors do not depend much on the value of φv , and thus the element internal angles have little
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
775
effect on the interpolation error. When the aspect ratio r = 10 −1 , which is far away from the optimal values, the H 1 -seminorm of the error is about the smallest when φv = π/6, π/2, 5π/6, which correspond to acute isosceles triangles, and the largest error occurs at φv = 0, π/3, 2π/3, π, which correspond to obtuse isosceles triangles. The L2 -norm is still insensitive to the internal angles. Hence, we conclude that when proper aspect ratio r ≈ S(∇3 u) is chosen and proper alignment is maintained, the H 1 -seminorm of the interpolation error is insensitive to the element internal angles. If a much larger aspect ratio is used, then the H 1 -error depends on the internal angles. The larger the maximum internal angle is, the larger the H 1 -error. In this case, the maximum angle condition is essential [2]. The L2 -norm of the interpolation error is always insensitive to the internal angles. 7. Conclusion and discussions. In the previous sections we described how to measure the anisotropic behavior of the third order derivative ∇3 u at a fixed point based on the distribution of the third order directional derivatives at the point. In particular, the orientation of ∇3 u is defined as the direction along which the third order directional derivative is about the smallest, while along the perpendicular direction it is about the largest. The anisotropic ratio S(∇3 u) measures the strength of the anisotropic behavior of ∇3 u. These quantities are invariant under translation and rotation of the independent variables. Then we considered the model problem of quadratic interpolation of cubic functions over a triangular element. In particular, we studied the optimal choice of the aspect ratio for elements aligned with the orientation of ∇3 u. It turns out that (1) when an element is aligned with the orientation of ∇3 u, both the H 1 -seminorm and L2 -norm of the error are about the smallest if its aspect ratio is approximately equal to S(∇3 u). This is true regardless of whether the element is an acute or an obtuse isosceles triangle. In practice, we may use the orientation of ∇3 u for the mesh alignment direction and the anisotropic ratio S(∇3 u) for the element aspect ratio. Our study also indicates that unlike the case of linear interpolation [18, 6], there is no “superaccuracy” in favor of the acute isosceles triangles for quadratic interpolation in general. While there is no rigorous proof of it yet, we believe the reason is because ∇3 u has richer geometric features than a triangle can capture. (Indeed, ∇3 u has four degree of freedoms, while a triangle has only three: area, aspect ratio, and orientation.) (2) When an element is aligned with ∇3 u and is of aspect ratio S(∇3 u), the H 1 -seminorm of the error is proportional to [S(∇3 u)]−1 , and the L2 -norm of the error is proportional to [S(∇3 u)]−3/2 . Both of them are insensitive to the internal angles of the element. However, if the element aspect ratio is much larger than the anisotropic ratio of ∇3 u, the H 1 -seminorm of the error increases greatly as the largest internal angle increases. The maximum angle condition is essential in this situation. The L2 -error is always insensitive to the internal angles. The definitions of the orientation and the anisotropic ratio for ∇3 u can be generalized to the case of ∇m u for m > 3. It is also possible to derive more precise and coordinate independent error estimates for higher order interpolation based on these anisotropic measures; see [7]. These results can be used to better understand and improve the anisotropic mesh generation and local refinement. They can also be used for global mesh refinement strategies like the r-refinement or the moving mesh method. In these techniques, one typically needs to define a mesh metric (or monitor function) M (x) and looks for a mesh that is quasi-uniform under the metric [3, 4, 5, 8, 9, 15]. For quadratic interpolation, one may choose the metric in the following form based on the anisotropic behavior of the third order derivatives of interpolated function u:
776
WEIMING CAO
2
−1
2 2
2
function u=y(y −ε x ) and φu=0, φv=π/2
10
2 2
function u=y(y −ε x ) and φu=0, φv=0
ε=2−1
−1
10 ε=2−3
−7
ε=2
−9
1
ε=2 −4
H1−seminorms
10
1
maximum norms L1−norms
−5
2
10
L −norms
ε=2−3
−2
10
−5
ε=2 −3
10
−7
ε=2
ε=2−9
1
ε=2 −3
10
1.5
various norms of interpolation error
ε=2−1 −5
various norms of interpolation error
−2
10
−4
H1−seminorms
10
1
maximum norms L1−norms
1
L2−norms
1.5
−5
10
1
−6
10
−6
10 0
1
10
2
10
3
10 aspect ratio r
4
10
0
10
1
10
0
3
10 aspect ratio r
4
10
10
0
10
10
function u=(y−ε x)(y2−ε2x2) and φu=0, φv=π/2
function u=(y−ε x)(y2−ε2x2) and φu=0, φv=0 −1 ε=2
−1
−1
10
−3
ε=2
−2
10
ε=2−5 −3
ε=2−7
10
−9
1
ε=2 −4
H1−seminorms
10
1
maximum norms
various norms of interpolation error
ε=2
ε=2−3 −5
ε=2 −2
10
ε=2−7 −9
ε=2
−3
10
1
−1
10
various norms of interpolation error
2
10
−4
H1−seminorms
10
1
maximum norms
1
1
L −norms
L2−norms
1.5
1.5
L −norms −5
10
−5
10
1
1
−6
10
−6
0
1
10
2
10
3
10 aspect ratio r
10
4
10
0
10
1
10
0
2
10
3
10 aspect ratio r
10
10
function (5.3) and φu=0, φv=π/2
function (5.3) and φu=0, φv=0
−1
−1
10
ε=2−1
ε=2−1 various norms of interpolation error
−3
ε=2
−2
10
ε=2−5 ε=2−7
−3
10
1
ε=2−9 −4
H1−seminorms
10
1
maximum norms
−3
ε=2
−2
10
ε=2−5 −3
10
ε=2−7 ε=2−9
1
10
−4
H1−seminorms
10
1
maximum norms
L1−norms
L1−norms −5
L −norms
2
1.5
2
1.5
−5
10
10
1
L −norms 1
−6
10
4
10
0
10
various norms of interpolation error
L2−norms
−6
0
10
1
10
2
10 aspect ratio r
3
10
4
10
10
0
10
1
10
2
10 aspect ratio r
3
10
4
10
Fig. 5. The H 1 -seminorm and L1 -, L2 -, and maximum norms of the quadratic interpolation errors on an aligned element versus the element aspect ratio r. (i) Upper panel: for function u = y(y 2 − x2 ); (ii) middle panel: for function u = (y − x)(y 2 − 2 x2 ); (iii) bottom panel: for function (5.3) with various γ ∈ [0, arctan()]. The left panel is for the case φu = 0 and φv = π/2 (i.e., aligned acute isosceles triangle), and the right panel is for the case φu = 0 and φv = 0 (i.e., aligned obtuse isosceles triangle).
777
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
4
6
10
10 function u=y(y2−ε2x2): L2−error v.s. alignment
3
10
ε=2
10
10 H1−seminorm of interpolation error
L2−norm of interpolation error
4
ε=2−7
1
10
ε=2−5
0
10
ε=2−3
−1
10
−1
ε=2 −2
10
−3
10
−4
ε=2−9
3
10
ε=2−7
2
10
−5
ε=2 1
10
0
ε=2−3
−1
ε=2
10
−1
10
−2
10
10
−5
−3
10
10
−6
−4
10 −0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 alignment angle φu/π
0.2
0.3
0.4
10 −0.5
0.5
4
−0.4
−0.3
−0.1 0 0.1 alignment angle φu/π
0.2
0.3
0.4
0.5
0.3
0.4
0.5
0.3
0.4
0.5
10 function u=(y−ε x)(y2−ε2x2): L2−error v.s. alignment
3
10
ε=2
10
4
10 H1−seminorm of interpolation error
ε=2−7
1
10
ε=2−5
0
10
ε=2−3
−1
functions u=(y−ε x)(y2−ε2x2): H1−error v.s. alignment
5
10
−9
2
10
−1
ε=2 −2
10
−3
10
−4
ε=2−9
3
10
ε=2−7
2
10
ε=2−5 1
10
0
ε=2−3
−1
ε=2−1
10
10
−2
10
10
−5
−3
10
10
−6
−4
10 −0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 alignment angle φu/π
0.2
0.3
0.4
10 −0.5
0.5
4
−0.4
−0.3
−0.2
−0.1 0 0.1 alignment angle φu/π
0.2
6
10
10 function (5.3): L2−error v.s. alignment
function (5.3): H1−error v.s. alignment
−9
2
−0.2
6
10
L2−norm of interpolation error
function u=y(y2−ε2x2): H1−error v.s. alignment
5
10
−9
2
ε=2
4
10
10
ε=2−9
H1−seminorm of interpolation error
L2−norm of interpolation error
ε=2−7 −5
ε=2 0
10
ε=2−3 −1
ε=2 −2
10
−4
10
ε=2−5
ε=2−3 0
10
ε=2−1
−2
10
10
−6
−4
10
−0.5
ε=2−7 2
10 −0.4
−0.3
−0.2
−0.1 0 0.1 alignment angle φu/π
0.2
0.3
0.4
0.5
−0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 alignment angle φu/π
0.2
Fig. 6. The L2 -norm (left) and H 1 -seminorm (right) of the quadratic interpolation errors on an obtuse isosceles element (i.e., φv = 0) versus its alignment angle φu /π for the case r = −1 . (i) Upper panel: for function u = y(y 2 −x2 ); (ii) middle panel: for function u = (y −x)(y 2 −2 x2 ); (iii) bottom panel: for function (5.3) with various γ ∈ [0, arctan()].
778
WEIMING CAO
2
2 2
2
function u=y(y −ε x ): L −error
function u=y(y2−ε2x2): H1−error
−1
10
−2
10
H1−seminorm of interpolation error
ε=2−1
−3
10
ε=2−3
−4
10
−5
ε=2
2
L −norm of interpolation error
ε=2−1
−5
10
−2
10
ε=2−3
ε=2−5 −3
10
ε=2−7
ε=2−7
−9
ε=2
−4
10
ε=2−9
−6
10
0
0.1
0.2
0.3
0.4 0.5 0.6 internal angle φv/π
2
2 2
0.7
0.8
0.9
1
0
0.3
0.4 0.5 0.6 internal angle φv/π
0.7
0.8
0.9
1
0.7
0.8
0.9
1
0.7
0.8
0.9
1
function u=(y−ε x)(y2−ε2x2): H1−error
−1
10
−1
10
ε=2
−1
H1−seminorm of interpolation error
ε=2
ε=2−3
−3
10
−5
ε=2
−4
10
2
L −norm of interpolation error
0.2
2
function u=(y−ε x)(y −ε x ): L −error −2
0.1
ε=2−7
ε=2−3
−2
10
−5
ε=2
−3
10
ε=2−7
−5
10
ε=2−9
−9
ε=2
−4
10
−6
10
0
0.1
0.2
0.3
0.4 0.5 0.6 internal angle φv/π
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4 0.5 0.6 internal angle φv/π
−1
10
function (5.3): L2−error function (5.3): H1−error
−1
−1
ε=2
ε=2−3 −3
10
ε=2−5 −4
10
ε=2−3 −2
10
−5
ε=2
1
2
L −norm of interpolation error
10
ε=2−1 H −seminorm of interpolation error
−2
10
ε=2−7
−3
10
ε=2−7
−5
10
−9
ε=2
−9
ε=2
−6
10
−4
0
0.1
0.2
0.3
0.4
0.5 angle φv/π
0.6
0.7
0.8
0.9
1
10
0
0.1
0.2
0.3
0.4
0.5 angle φv/π
0.6
Fig. 7. The L2 -norm (left) and H 1 -seminorm (right) of the quadratic interpolation errors on an element versus the angle φv that controls the element internal angles. The element is always aligned with ∇3 u (i.e., φu = 0), and its aspect ratio is fixed as r = −1 . (i) Upper panel: for function u = y(y 2 − x2 ); (ii) middle panel: for function u = (y − x)(y 2 − 2 x2 ); (iii) bottom panel: for function (5.3) with various γ ∈ [0, arctan()].
779
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
2
2 2
function u=y(y2−ε2x2): H1−error
2
function u=y(y −ε x ): L −error −2
10
−1
ε=2
−1
ε=2
−1
H1−seminorm of interpolation error
−3
10
ε=2−3
ε=2−5
−4
10
2
L −norm of interpolation error
10
−7
ε=2
−5
10
ε=2−3
−2
−5
10
ε=2
ε=2−7 −3
10
ε=2−9
−9
ε=2 −6
10
0
0.1
0.2
0.3
0.4 0.5 0.6 internal angle φv/π
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4 0.5 0.6 internal angle φv/π
2
0
10
2 2
0.8
0.9
1
0.7
0.8
0.9
1
0.7
0.8
0.9
1
1
function u=(y−ε x)(y −ε x ): H −error
−1
ε=2
function u=(y−ε x)(y2−ε2x2): L2−error
0.7
0
10
ε=2−1
ε=2−3
−1
H1−seminorm of interpolation error
ε=2−3 −2
10
−5
ε=2 −3
10
2
L −norm of interpolation error
10
ε=2−7
−1
ε=2−5
10
ε=2−7
−2
10
−9
ε=2
−4
10
−3
10
ε=2−9
0
0.1
0.2
0.3
0.4 0.5 0.6 internal angle φv/π
0.7
0.8
0.9
1
0
0
0.1
0.2
0.3
1
10
10
function (5.3): H1−error
function (5.3): L2−error
0
ε=2−1
10
−3
ε=2 −2
10
ε=2−5 −3
10
ε=2−3
−1
10
−7
ε=2
−5
ε=2
ε=2−7 −2
10
1
2
L −norm of interpolation error
ε=2−1
H −seminorm of interpolation error
−1
10
0.4 0.5 0.6 internal angle φv/π
ε=2−9
−4
10
−3
10 ε=2−9 −5
10
0
0.1
0.2
0.3
0.4
0.5 angle φv/π
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5 angle φv/π
0.6
Fig. 8. The L2 -norm (left) and H 1 -seminorm (right) of the quadratic interpolation errors on an element versus the angle φv that controls the element internal angles. The element is always aligned with ∇3 u (i.e., φu = 0), and its aspect ratio is fixed as r = 10−1 . (i) Upper panel: for function u = y(y 2 − x2 ); (ii) middle panel: for function u = (y − x)(y 2 − 2 x2 ); (iii) bottom panel: for function (5.3) with various γ ∈ [0, arctan()].
780
(7.1)
WEIMING CAO
M (x) = Rφ
μ1
μ2
RφT ,
where φ is the angle (with the x-axis) of the orientation of ∇3 u, (7.2)
μ1 /μ2 ≈ [S3 ]2 ,
and the determinant of M (x) is taken as (7.3)
μ1 μ2 = const · (D3 /S3 )4/3
for minimizing the H 1 -seminorm of the interpolation error and as (7.4)
μ1 μ2 = const · D3 (S3 )−3/2
for minimizing the L2 -norm of the interpolation error. Some modification is clearly needed when the function is almost unidirectional (i.e., the anisotropic ratio is ∞). Finally, we would like to emphasize that the definition of orientation and anisotropic ratio for ∇3 u described in [7] and here is only one possible way to characterize the anisotropic behavior of ∇3 u. Other definitions, possibly dependent on the norms used to measure the errors, can be introduced. In addition, when problems of multicomponents are considered, the anisotropic behavior of each variable should be carefully combined to produce an ideal direction for mesh alignment and an appropriate aspect ratio. These topics need further investigation. REFERENCES [1] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., Teubner, Stuttgart, 1999. [2] I. Babuˇ ska and A. K. Aziz, On the angle condition in the finite element method, SIAM J. Numer. Anal., 13 (1976), pp. 214–226. [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] H. Borouchaki, P. L. George, and B. Mohammadi, Delaunay mesh generation governed by metric specifications. II. Applications. Finite Elem. Anal. Des., 25 (1997), pp. 85–109. [5] W. Cao, W. Huang, and R. D. Russell, A study of monitor functions for two-dimensional adaptive mesh generation, SIAM J. Sci. Comput., 20 (1999), pp. 1978–1994. [6] 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. [7] W. Cao, An interpolation error estimates in R2 based on anisotropic measures of higher order derivatives, Math. Comp., to appear. [8] L. Chen, P. Sun, and J. Xu, Optimal anisotropic meshes for minimizing interpolation errors in Lp -norm, Math. Comp., 76 (2007), pp. 179–204. [9] L. Chen and J. Xu, Optimal Delaunay triangulations, J. Comput. Math., 22 (2004), pp. 299– 308. [10] P. G. Ciarlet, The Finite Element Methods for Elliptic Problems, Classics Appl. Math. 40, SIAM, Philadelphia, 2002. [11] E. F. D’Azevedo, Optimal triangular mesh generation by coordinate transformation, SIAM J. Sci. Statist. Comput., 12 (1991), pp. 755–786. [12] E. F. D’Azevedo and R. B. Simpson, On optimal triangular meshes for minimizing the gradient error, Numer. Math., 59 (1991), pp. 321–348. [13] L. Formaggia and S. Perotto, New anisotropic a priori error estimates, Numer. Math., 89 (2001), pp. 641–667. [14] W. Huang, Measuring mesh qualities and application to variational mesh adaptation, SIAM J. Sci. Comput., 26 (2005), pp. 1643–1666.
QUADRATIC INTERPOLATION ERROR ON TRIANGLE
781
[15] W. Huang and W. Sun, Variational mesh adaptation II: Error estimates and monitor functions, J. Comput. Phys., 184 (2003), pp. 619–648. [16] E. J. Nadler, Piecewise Linear Approximation on Triangulations of a Planar Region, Ph.D. thesis, Division of Applied Mathematics, Brown University, Providence, RI, 1985. [17] S. Rippa, Long and thin triangles can be good for linear interpolation, SIAM J. Numer. Anal., 29 (1992), pp. 257–270. [18] J. R. Shewchuk, What is a good linear finite element? Interpolation, conditioning, anisotropy, and quality measure, in Proceedings of the 11th International Meshing Roundtable, Sandia National Laboratories, Albuquerque, NM, 2002, pp. 115–126. [19] R. B. Simpson, Anisotropic mesh transformations and optimal error control, Appl. Numer. Math., 14 (1994), pp. 183–198. [20] 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.