Error-constant estimation under the maximum norm for linear Lagrange interpolation

For the linear Lagrange interpolation over a triangular domain, we propose an efficient algorithm to rigorously evaluate the interpolation error constant under the maximum norm by using the finite-element method (FEM). In solving the optimization problem corresponding to the interpolation error constant, the maximum norm in the constraint condition is the most difficult part to process. To handle this difficulty, a novel method is proposed by combining the orthogonality of the space decomposition using the Fujino–Morley FEM space and the convex-hull property of the Bernstein representation of functions in the FEM space. Numerical results for the lower and upper bounds of the interpolation error constant for triangles of various types are presented to verify the efficiency of the proposed method.


Introduction
In this paper, we consider the error estimation for the linear Lagrange interpolation over triangle elements and provide explicit values for the error constant in the error estimation under the L ∞ -norm.
Before the detailed discussion of our results, let us introduce the existing literature on the Lagrange interpolation function in a general scope.
• (1D case) Given 1-dimensional interval I = (0, 1), since H 1 (I) ⊂ C(I), we can define the Lagrange interpolation Π L u such that Π L u is a linear function satisfying (u − Π L u)(0) = (u − Π L u)(1) = 0.Then, the following results are well known as optimal estimates if u is regular enough in the sense that the right-hand sides of the inequalities are well defined: where u (2) denotes the second derivative of u, • 0,I and • ∞,I denote the L 2 -and L ∞ -norms, respectively, and | • | 1,I and | • | 2,I denote the H 1 -and H 2seminorms, respectively.The estimation presented above are optimal in the sense that there exist functions for which the equalities hold.
• (2D case) Over a triangle K with vertices p i (i = 1, 2, 3), the Lagrange interpolation function Π L u is the linear function such that (see Figure 1) (u − Π L u)(p i ) = 0, ∀i = 1, 2, 3.In case of L 2 -norm and H 1 -seminorm error estimation of Π L , one need to estimate the interpolation error constants appearing in the following inequalities: Let h be the medium edge length of K, θ the maximum angle, and αh (0 < α ≤ 1) the smallest edge length.Kikuchi and Liu [1,2] obtained a bound of C 0 and C 1 as follows: . Also, Kobayashi [3] showed that for a triangle K with edge lengths A, B, C and area S, the following holds: where the constant C 1 (K) is defined by The optimal estimation of constants C 0 (K) and C 1 (K) for a concrete K can be obtained by solving corresponding eigenvalue problems with rigorous lower eigenvalue bounds; see results of [4,5].
For L ∞ -norm error estimation under the L ∞ -norm of objective function, Waldron [6] provides the following sharp inequality: , where R is the radius of the circumscribed circle of K, d is the distance of the center c of the circumscribed circle from K, and u (2) ∞,K is defined by
A detailed discussion on the L ∞ -norm of interpolation error for a quadratic polynomial f is considered by D'Azevedo and Simpson [7].In [8], Shewchuk gives a survey of the interpolation error estimation with L ∞ -norm for both f − Π L f and ∇(f − Π L f ), along with the discussion on the relation between the interpolation error and the finite element approximation error functions.Also, the discussion on the affection of the aspect ratio of triangle element to the interpolation error can be found in Cao [9].
In this research, we consider the L ∞ -norm estimation for the Lagrange interpolation over triangle element K by using the H 2 -seminorm of the objective function, that is, Here C L (K) is the interpolation error constant to be evaluated explicitly.For example, for a unit isosceles right triangle element, we have an estimation of the optimal constant C L (K) as Estimation of C L (K) for triangle of general shapes is provided in Theorem 2.3, while sharp bounds for concrete triangles are discussed in §3.Such kind of estimation is helpful to provide explicit maximum norm error estimation for the FEM solution to boundary value problems by further applying the point-wise error estimation (see e.g., [10]), which will be considered in our succeeding work; see also classical qualitative error analysis under the maximum norm in Section 19-22 of [11]; The contribution of our paper is summarized as follows.
(1) For triangle element K of general shapes, a formula to give an upper bound of C L (K) is obtained by theoretical analysis.The bound is raw but works well for triangle element of arbitrary shapes.Particularly, our analysis tells that the value of C L (K) can be very large and tend to ∞ if the triangle element tends to degenerate to a 1D segment; see detail in §2.2.
(2) For specific triangle element K, the optimal estimation of C L (K) is obtained by solving the corresponding optimization problem over H 2 (K) under the constraint condition involving L ∞ -norm.The processing of the constraint condition with L ∞ -norm is not an easy work.We develop a novel algorithm to provide efficient and sharp estimation for the solution of the optimization problem.With a light computation, one can obtain the estimation of C L (K) with relative error less than 1%.
The rest of our paper is structured as follows.At the end of this section, we introduce the preliminary concepts and notations to be used throughout the paper.In §2, the estimation of the upper bound for C L (K) is considered using theoretical approach.The raw upper bound of the interpolation error constant is calculated for a right isosceles triangle.Also, we investigate the asymptotic behavior of the constant as the triangle tends to degenerate.In §3, using finite element method (FEM), an algorithm for the optimal estimation of the constant is proposed.Lower bounds for the constant are calculated to confirm the efficiency of the proposed algorithm.The numerical results are summarized and the conclusion is presented in §4.
Notation Let us introduce the notation for the function spaces used in this paper.In most cases, the domain Ω of functions is selected as a triangle element K.The standard notation is used for Sobolev function spaces W k,p (Ω).The associated norms and seminorms are denoted by • k,p,Ω and | • | k,p,Ω , respectively (see, e.g., Chapter 1 of [12] and Chapter 1 of [13]).Particularly, for special k and p, we use abbreviated notations as ,Ω , and L p (Ω) = W 0,p (Ω).The set of polynomials over K of up to degree k is denoted by P k (K).The second order derivative is given by D 2 u := (u xx , u xy , u yx , u yy ) for u ∈ H 2 (K).
Given a triangle K, denote each vertex by p i (i = 1, 2, 3) and the largest edge length by h K ; see Figure 2. We follow the notation introduced by Liu and Kikuchi [2] to configure a general triangle with geometric parameters.Let h, α and θ be positive constants such that Define a triangle K α,θ,h with three vertices p 1 (0, 0), p 2 (h, 0) and p 3 (αh cos θ, αh sin θ).Note that h ≤ h K .In case of h = 1, the notation K α,θ,1 is abbreviated as K α,θ .
With the above configuration of the triangle K α,θ,h , the optimal constant C L (K) in (1) can be defined as follows: By scaling of the triangle element, it is easy to confirm that C L (α, θ, h) = hC L (α, θ, 1).
In the rest of the paper, we show how to obtain explicit bounds for the error constant C L (α, θ, h).In this section, a raw upper bound of the constant is obtained through theoretical analysis.Such a bound applies to triangles of arbitrary shapes.
First, let us quote a lemma about the trace theorem, which gives estimation for the integral over edge of a triangle element.For reader's convenience, we show the proof in a concise way; refer to e.g.[14,15,16] for more detailed discussion.
Proof For any w ∈ H 1 (K), the Green theorem leads to Here, − → n is the unit outer normal direction on the boundary of K.For the term Here, H K is the height of the triangle with base as e.Thus, We can now draw the conclusion by sorting the above inequality.
Using the trace theorem, the following result provides a pointwise estimation of the interpolation error.
Lemma 2.2 Given u ∈ H 2 (K), for any point x 0 ∈ K, we have where h K is the longest edge length of K, and H K is the height of the subtriangle K = p 1 x 0 p 3 with respect to the base ẽ = p 1 x 0 (see Figure 4).
Taking the Taylor expansion of g on the segment ẽ and noting that g(p 1 ) = 0, The conclusion follows.
Liu and Kikuchi [2] considered the estimation of the constant C 1 (α, θ) for different types of triangles K = K α,θ such that where h is the medium length of K.The constant C 1 (α, θ) is used to give a bound for C L (K) as shown in the lemma below.
2.1 The case for a right isosceles triangle Using Lemma 2.3, we obtain the upper bound of the constant for right isosceles triangles: For the right isosceles triangle Then, since y = h − x, The above quantity takes its maximum value at p = h 2 , h 2 and p = (h, 0), and its maximum value is 1.Thus, for any p on p 2 p 4 , |p Hence, we obtain the error estimate for a right isosceles triangle as in (4).

Dependence of the constant on the shape of K
In this subsection, we consider the variation of the interpolation constant when a reference triangle, i.e., the right isosceles triangle, is transformed to a general triangle.
Remark 2.1 By using the raw bound of C L 1, π 2 ≤ 1.3712h in (4), an explicit but raw bound of C L (α, θ) is available.Later, with a sharp and rigorous estimation of C L 1, π 2 based on numerical approach, the bound can be improved as Remark 2.2 Here are two remarks on the asymptotic behavior of the constant when the triangle tends to degenerate.

Optimal estimation of the constant
In the previous section, we obtain explicit bounds for the interpolation constant for triangles of general shape.Basically, such bounds from theoretical analysis only provide raw bound for the objective constant.In this section, we propose a numerical algorithm to obtain the optimal estimation of the constant C L (K) for specific triangles.
Let us define the space V L (K) := u ∈ H 2 (K) u(p i ) = 0 (i = 1, 2, 3) .Let T h be a triangulation of the domain K and define the space ds is continuous on each interior edge e of T h .
For u h , v h ∈ V FM h (K), define the discretized H 2 -inner product and seminorm by Let us define the two quantities over the triangle K: holds.In Theorem 3.1, we describe the algorithm to bound λ by using λ h .
Given u ∈ H 2 (K), the Fujino-Morley interpolation Π FM h u is a function satisfying and at the vertices p i and edges e i of K, The Fujino-Morley interpolation has the property that (see, e.g., [4,5]) . Thus, it is easy to see that the Fujino-Morley interpolation is just the projection Below, let us introduce the theorem that provides an explicit lower bound of λ.Such a result is inspired by idea of [17] for the lower bounds of eigenvalue problems.Theorem 3.1 Suppose there exists a quantity C FM h such that, Then, the following lower found of λ(K) is available.
Proof For any u ∈ V L (K), noting that Π FM h u 2,K ≥ √ λ h Π FM h u ∞,K and applying the inequality (8), we have From the orthogonality in ( 7), we have Thus, From the definition of λ in (6), we draw the conclusion.
To apply Theorem 3.1 for bounding λ, an explicit value of C FM h is needed.Below, let us describe the way to obtain this explicit value by utilizing the raw bound of C L (α, θ).

Explicit estimation of C FM h
To have an explicit value of C FM h , we first define the quantity C FM res (T ) for each element T in the triangulation T h : Here, Then, the following C FM h with an upper bound makes certain (8) holds: Remark 3.1 Let T h be a uniform triangulation of a right isosceles triangle; see a sample mesh in Figure 8.We choose an explicit upper bound of , where h is the leg length of each right triangle element.).The estimation of λ h is equivalent to finding the solution to the optimization problem where the components a ij of A are given by a ij = φ i , φ j h , {φ i } i=1,...,M are the basis functions for the Fujino-Morley space V FM h , and x ∈ R M denotes the Fujino-Morley coefficient vector of u h ∈ V FM h .To solve the optimization problem (11) is not an easy work since the L ∞ -norm of the function appears in the constraint.Here, we introduce the technique to apply Bernstein polynomials and their convex-hull property to solve the problem.Strictly speaking, a new optimization problem (12) utilizing the Bernstein polynomials will be formulated to provide a lower bound for the solution of (11).
As a preparation, let us introduce the definition of Bernstein polynomials along with the convex-hull property; refer to, e.g., [18,19] for detailed discussion.
Convex-hull property of Bernstein polynomials Given a triangle K, let (u, v, w) be barycentric coordinates for a point x in K.A Bernstein polynomial u h of degree n over a triangle K is defined by Here, J i,j,k (x) are the Bernstein basis polynomials; the coefficients d i,j,k are the control points of u h .Noticing that we can easily obtain the following convex-hull property of Bernstein polynomials: can be represented by the Bernstein basis polynomials.Let B be the N ×M matrix that transforms the Fujino-Morley coefficients x to the Bernstein coefficients d B .Note that u h is regarded as a piecewise Bernstein polynomial so that its Bernstein coefficient vector d B has the dimension N = 6 × #{elements}.From the convex-hull property of the Bernstein polynomials, the following inequality holds: Based on this inequality, we propose a new optimization by relaxing the constraint condition of (11): The solution to problem (12) provides a lower bound for (11), i.e., λ h ≥ λ h,B .Below, we propose an algorithm to solve the problem (12).Since A is positive definite, let us consider the Cholesky decomposition of A: A = R T R, where R is an M × M upper triangular matrix.Then, by letting y := Rx and B := BR −1 , problem (12) becomes The following lemma shows the solution for problem (13).
Lemma 3.1 [1] Let b T i (i = 1, . . ., N ) be the ith row of B and b T max be a row of B satisfying b max 2 = max i=1,...,N b i 2 .Then, the optimal value of problem ( 13) is given by . Hence, For any y ∈ S, from the Cauchy-Schwarz inequality, [1] Appreciation to Tamaki TANAKA and Syuuji YAMADA from Faculty of Science, Niigata University for their idea of solving this problem in an efficient way. Thus, From ( 14) and ( 15), we draw the conclusion.
Note that the diagonal elements of ).Therefore, we can solve problem (12) without performing the Cholesky decomposition of A, as shown by the following lemma.Lemma 3.2 Let D := BA −1 B T .The optimal value of ( 12) is given by , where diag(D) is the diagonal elements of D.
Theorem 3.1 gives a lower bound for λ.
, this lower bound is used to obtain an upper bound for C L (K).Below, let us summarize the procedure to obtain a lower bound for λ.

Algorithm for calculating lower bound of λ(K)
a. Set up the FEM space V FM h (K) = span{φ i } M i=1 over a triangulation of the triangle domain K. b.Assemble the global matrix A = (a ij ) M ×M (a ij = φ i , φ j h ) and the transformation matrix B from Fujino-Morley coefficients to Bernstein coefficients.c.Apply Lemma 2.3 to obtain a raw bound for C FM h .d. Apply Lemma 3.1 or Lemma 3.2 to calculate λ h,B (≤ λ h ).e.The lower bound for λ is obtained through Theorem 3.1 by using λ h,B and the upper bound of C FM h .
Using uniform triangulation of a domain K, a direct estimation of the lower bound for λ without using C FM h is available.Corollary 3.1 For a uniform triangulation of K = K α,θ,h with N subdivisions for each side, the following holds: Proof Since (C L (K)) 2 = 1/λ(K) and each T ∈ T h is similar to K, we have, The conclusion is achieved by sorting the inequality.
Remark 3.2 Theoretically, for a refined uniform triangulation, the lower bound (9) using C FM h is sharper (i.e., larger) than (16).This fact can be confirmed by utilizing the following relation: For a small value of h = 1/N , we have Thus, the second equality of ( 17) holds due to C FM res (T ) < C L (T ).However, in practical computation, the raw estimate of C FM res (T ) will cause a worse bound of λ than (16).
Using Corollary 3.1, the following steps are modified from the algorithm to obtain a lower bound for λ, without using the quantity of C FM h : Revision of algorithm for calculating lower bound of λ(K) c*.Apply Lemma 3.1 or Lemma 3.2 to calculate λ h,B (≤ λ h ).d*.Solve the lower bound for λ using Corollary 3.1 along with λ h,B .
Remark 3.3 To compare the efficiencies of the two formulas (9) and ( 16), we apply them to estimate λ for a unit right isosceles K 1,π/2 .By using uniform triangulation of size h = 1/64, the estimate (9) gives λ ≥ 5.7659 and (16) gives a sharper bound as λ ≥ 5.7798.Hence, a sharper upper bound is obtained using (16) and we have the following estimation: As a comparison, the result (4) will yield a raw bound as C L (1, π/2, h) ≤ 1.3712h.
For a triangle K α,θ with two fixed vertices p 1 (0, 0), p 2 (1, 0), let us vary the vertex p 3 (x, y) and calculate the approximate value of C L (α, θ) for each position of p 3 .Note that C L can be regarded as a function with respect to the coordinate (x, y) of p 3 , which is denoted by C L (x, y).In Figure 9, we draw the contour lines of C L (x, y) where the abscissa and the ordinate denote x-and y-coordinates of p 3 , respectively.

Lower bound of the constant
To confirm the precision of the obtained estimation for the Lagrange interpolation constant, the lower bounds of the constants are calculated.Let u h be the function obtained by numerical computation solving the minimization problem.To obtain the lower bound, an appropriate polynomial f over K of higher degree d is selected by solving the minimization problem below: where p i denote the nodes of the triangulation of K. From the definition of λ(K) in ( 6) and the relation C L (K) = 1/ λ(K), we have a lower bound of C L (K) as follows: Remark 3.4 For the unit right isosceles triangle K 1,π/2 , the upper bound for the constant is obtained by solving the optimization problem with mesh size 1/64.Meanwhile, the lower bound of the constant is obtained by using a polynomial of degree 9.The two-side bounds reads:

Numerical results and conclusion
In this section, we perform numerical computation to obtain the estimation of the interpolation error constant C L (K) for triangles of various shapes.First, let us confirm the shape of the function u h that solves the minimization problem for λ h,B in case K being the unit isosceles right triangle.The contour lines of u h are displayed in Figure 10.The numerical computation tells that the maximum value of u h happens on the midpoint of the hypotenuse of K.Note that the maximum value of u h is around 0.95 while the maximum of its Bernstein coefficients is above 1.
Let us also compare the lower bounds of λ obtained through Theorem 3.1 and Corollary 3.1 for various triangles.Table 1 tells that the values obtained using Corollary 3.1 gives a sharper estimate of λ.
Table 2 summarizes the results for the lower and upper bounds of the constant for different types of triangle K 1,θ with the mesh size as h = 1/32 and h = 1/64.Figure 11 demonstrates the convergency of the upper and lower bounds of the interpolation error constant as the mesh is refined.It implies that the convergency order of upper bounds depends on the shape of the triangles.The theoretical analysis on the efficiency and the convergency of the algorithm in solving the optimization problem is beyond the scope of this paper and will be systematically investigated in our succeeding research.Conclusion In this research, we provide explicit estimates for the L ∞ -norm error constant C L of the linear Lagrange interpolation function over triangular elements.The formula in Theorem 2.1 provides a bound of C L that holds for triangle of arbitrary shapes.Theorem 3.1 in §3 proposes a numerical approach to obtain optimal bounds for the constant C L over a concrete triangle.The optimization problem corresponding to C L is novelly solved by utilizing the convex-hull property of Bernstein polynomials.In the near future, the convergency of the numerical approach to solve the optimization problems involving the maximum norm will be systematically considered.

Figure 1 A
Figure 1 A linear Lagrange interpolation function Π L u defined on a triangle K

Figure 3 A
Figure 3 A triangle K with base e and height H K

Figure 4 A
Figure 4 A subtriangle K in a triangle K

Figure 5 A
Figure 5 A right isosceles triangle K 1, π 2 ,h

Figure 6 A
Figure 6 A triangle K α,θ with angle θ close to π

Figure 7 A
Figure 7 A right triangle K α, π 2 with one leg length close to 0.

Figure 8 A
Figure 8 A uniform triangulation of a right isosceles triangle

Table 1
The lower bounds for λ through Theorem 3.1 and Corollary 3.1.

Table 2
The lower and upper bounds of C L (1, θ) for triangles of different shapes