Optimal estimation for the Fujino-Morley interpolation error constants

TL;DR

This paper develops a new finite element-based algorithm to accurately estimate the interpolation error constants for Fujino-Morley operators, providing concrete bounds and analyzing shape perturbations.

Contribution

It introduces a verified computational method for lower bounds of eigenvalues related to interpolation error constants, including shape variation analysis.

Findings

Provided concrete upper bounds for Fujino-Morley interpolation error constants.

Developed a new algorithm based on finite element methods with verified computation.

Analyzed how eigenvalues vary with shape perturbations of triangles.

Abstract

The quantitative estimation for the interpolation error constants of the Fujino-Morley interpolation operator is considered. To give concrete upper bounds for the constants, which is reduced to the problem of providing lower bounds for eigenvalues of bi-harmonic operators, a new algorithm based on the finite element method along with verified computation is proposed. In addition, the quantitative analysis for the variation of eigenvalues upon the perturbation of the shape of triangles is provided. Particularly, for triangles with longest edge length less than one, the optimal estimation for the constants is provided. An online demo with source codes of the constants calculation is available at http://www.xfliu.org/onlinelab/.

Tables

Table 1. Table 1: Error constants C L , 0 subscript 𝐶 𝐿 0 C_{L,0} and C L , 1 subscript 𝐶 𝐿 1 C_{L,1} for the Lagrange interpolation

$(a,b)$	$C_{L,0}$	$C_{L,1}$
$(0,1)$	0.167349	0.488767
$(1/2,\sqrt{3}/2)$	0.117134	0.318457
$(-\sqrt{2}/2,\sqrt{2}/2)$	0.245388	1.187998
$(0,0.1)$	0.108221	0.327955

Table 2. Table 2: Error constants C FM , 0 subscript 𝐶 FM 0 C_{\mbox{\tiny{FM}},0} and C FM , 1 subscript 𝐶 FM 1 C_{\mbox{\tiny{FM}},1} for the Fujino–Morley interpolation

$(a,b)$	$C_{\text{FM},0}$	$C_{\text{FM},1}$
$(0,1)$	0.090287	0.233708
$(1/2,\sqrt{3}/2)$	0.073583	0.174354
$(-\sqrt{2}/2,\sqrt{2}/2)$	0.093318	0.300773
$(0,0.1)$	0.060474	0.188918

Equations

$$({\Pi }^{\text{mbox}FM}u-u)(P)=0,P=O,A,B;$$

$${\int }_{e_i}\frac{\partial }{\partial n}({\Pi }^{\text{mbox}FM}u-u)\text{mbox}ds=0,i=1,2,3.$$

$$\Vert {\Pi }^{\text{mbox}FM}u-u{\Vert }_K\le C_0(K)|{\Pi }^{\text{mbox}FM}u-u{|}_{2,K},|{\Pi }^{\text{mbox}FM}u-u{|}_{1,K}\le C_1(K)|{\Pi }^{\text{mbox}FM}u-u{|}_{2,K}.$$

$$C_0(K)\le 0.2575h^2,C_1(K)\le C_{\text{mbox}CR}(K)\le 0.2983h.$$

$$(C_1(K)\le )C_{\text{mbox}CR}(K)\le 0.1893h.$$

$$0.07349h^2\le C_0(K)\le 0.07353h^2,0.18863h\le C_1(K)\le 0.18868h.$$

$$\Vert \Pi u-u{\Vert }_{0,K}\le d_0|u{|}_{2,K},|\Pi u-u{|}_{1,K}\le d_1|u{|}_{2,K},$$

$$V^{\text{mbox}FM}(K):=\{u\in H^2(K)|u(O)=u(A)=u(B)=0,{\int }_{e_i}\frac{\partial u}{\partial n}ds=0,i=1,2,3\}.$$

$$C_{0}(K):=\sup_{\scriptsize{\begin{array}[]{c}u\in V^{\mbox{\tiny{FM}}}(K)\\ D^{2}u\neq 0\end{array}}}\frac{\|u\|_{K}}{\|D^{2}u\|_{K}},~{}~{}~{}C_{1}(K):=\sup_{\scriptsize{\begin{array}[]{c}u\in V^{\mbox{\tiny{FM}}}(K)\\ D^{2}u\neq 0\end{array}}}\frac{\|\nabla u\|_{K}}{\|D^{2}u\|_{K}}\>.$$

$${\int }_{e_i}({\Pi }^{\text{mbox}CR}u-u)\text{mbox}ds=0,i=1,2,3.$$

$$C_{\mbox{\tiny{CR}}}(K)=\sup_{\scriptsize{\begin{array}[]{c}u\in H^{1}(K)\\ \nabla(u-\Pi^{{}^{\mbox{\tiny{CR}}}}u)\neq 0\end{array}}}\frac{\|u-\Pi^{{}^{\mbox{\tiny{CR}}}}u\|}{\|\nabla(u-\Pi^{{}^{\mbox{\tiny{CR}}}}u)\|}=\sup_{\scriptsize{\begin{array}[]{c}u\in V^{\mbox{\tiny{CR}}}(K)\\ \nabla u\neq 0\end{array}}}\frac{\|u\|}{\|\nabla u\|}\>,$$

$$V^{\text{mbox}CR}(K):=\{u\in H^1(K)|{\int }_{e_i}u\text{mbox}ds=0,i=1,2,3\}.$$

$$\underset{\text{mbox}diam(K)\le 1}{\max }{C}_{\text{mbox}CR}(K)\in [0.1890,0.1893].$$

$$\Omega =\{(a,b)\in {\mathbb{R}}^2|a^2+b^2\le 1,a\ge 1/2,b>0\}.$$

$$(D^{2}u,D^{2}v)=\lambda (u,v)\forall v\in V^{\text{mbox}FM}(K)$$

$$(D^{2}u,D^{2}v)=\mu (\nabla u,\nabla v)\forall v\in V^{\text{mbox}FM}(K)$$

$$C_0={\sqrt{{\lambda }_1}}^{-1},C_1={\sqrt{{\mu }_1}}^{-1}.$$

$$V^{\text{mbox}FM,h}(K):=\{u_h\in V^h|u_h(O)=u_h(A)=u_h(B)=0;$$

$${\int }_{e_i}\frac{\partial u}{\partial n}\text{mbox}ds=0\text{ for each ei of K}.\}.$$

$$(D^{2}u,D^{2}v)={\lambda }_h(u,v)\in V^{\text{mbox}FM,h}(K)$$

$$(D^{2}u,D^{2}v)={\mu }_h(\nabla u,\nabla v)\in V^{\text{mbox}FM,h}$$

$$V(h):=\{u+u_h|u\in V^{\text{mbox}FM},u_h\in V^{\text{mbox}FM,h}\},M(u,v):=(D^{2}u,D^{2}v)$$

$$N(u,v):=(u,v)\text{mbox}for{C}_0,N(u,v):=(\nabla u,\nabla v)\text{mbox}for{C}_1.$$

$$(D^2(u-P_{h}u),D^2{v}_h)=0,\forall v_h\in V^h.$$

$$\Vert u-P_{h}u\Vert \le C_{h,0}\Vert D^2(u-P_{h}u)\Vert ,\Vert \nabla (u-P_{h}u)\Vert \le C_{h,1}\Vert D^2(u-P_{h}u)\Vert .$$

$${\lambda }_1\ge \frac{{\lambda }_{h,1}}{1+{\lambda }_{h,1}{C}_{h,0}^2},{\mu }_1\ge \frac{{\mu }_{h,1}}{1+{\mu }_{h,1}{C}_{h,1}^2}.$$

$$(P_{h}u){|}_T={\Pi }^{\text{mbox}FM}(u{|}_T),\text{mbox}foreachT\in {\mathcal{T}}^h.$$

$$C_{i,h}=\underset{T\in {\mathcal{T}}^h}{\max }{C}_i(T),i=0,1.$$

$$\Vert u-P_{h}u\Vert \le 0.2575h^2\Vert D^2(u-P_{h}u)\Vert ,\Vert \nabla (u-P_{h}u)\Vert \le 0.1893h\Vert D^2(u-P_{h}u)\Vert .$$

$${\lambda }_1\ge \frac{{\lambda }_{h,1}}{1+0.2575^2{h}^4{\lambda }_{h,1}},{\mu }_1\ge \frac{{\mu }_{h,1}}{1+0.1893^2{h}^2{\mu }_{h,1}}.$$

$$Q:=\left(\begin{array}[]{cc}1&\alpha\\ 0&\beta\end{array}\right),$$

Full text

1 Introduction

The Fujino–Morley111The element discussed here is often called by “Morley element” in the existing literature. However, the same element is also proposed independently by T. Fujino in P.739 of Fujino1971 (proceedings of a conference on 1969), which is also cited by L.S.D. Morley in Morley1971 . finite element method (FEM) Fujino1971 ; Morley1968 ; Morley1971 provides a robust way to solve partial differential problems evolving bi-harmonic operators. Especially, in solving the eigenvalue problem of bi-harmonic operators, the Fujino–Morley FEM along with the Fujino–Morley interpolation operator can be utilized to find explicit lower bound for the eigenvalues; see the work of Carstensen–Gallistl Carstensen2013 and Liu Liu2015AMC . For the Fujino–Morley interpolation, there are two fundamental constants that are playing important roles in bounding the eigenvalues. In liu-you-amc-2018 , such constants are used to estimate the interpolation error constant for the quadratic Lagrange interpolation operator. Rough bounds of the two constants have been given in Carstensen2013 by using theoretical analysis. In this paper, we will propose a FEM based method to provide the optimal estimation of the constants.

Let $K$ be a triangle element with the largest edge length as $h$. The vertices of $K$ are denoted by $O$, $A$ and $B$ and the edges by $e_{1},e_{2},e_{3}$; see Fig. 1.

Over $K$, we define the Fujino–Morley interpolation operator $\Pi^{\mbox{\tiny{FM}}}$. For given $u\in H^{2}(K)$, $\Pi^{\mbox{\tiny{FM}}}u$ is a quadratic polynomial that satisfies

[TABLE]

For such an interpolation operator, the following error estimation is needed in numerical analysis, especially in solving the problem of bounding eigenvalues for bi-harmonic operators Carstensen2013 ; Liu2015AMC ; see, also, Theorem 3.1 for details of lower eigenvalue bounds evaluation:

[TABLE]

The definition of constants $C_{0}(K)$ and $C_{1}(K)$ can be found in (6). In carstensen2014guaranteed ; Carstensen2013 , to bound $C_{1}(K)$, the Crouzeix-Raviart constant $C_{\mbox{\tiny{CR}}}$ is introduced. Using the fact $C_{1}(K)\leq C_{\mbox{\tiny{CR}}}(K)$, the rough bounds for the constants are given by

[TABLE]

In Liu2015AMC , the optimal estimation of $C_{\mbox{\tiny{CR}}}(K)$ is given as

[TABLE]

In this paper, we provide an algorithm to evaluate the constants directly and verify that the following estimation of the above two constants holds for element $K$ of arbitrary shapes.

[TABLE]

The lower bounds of these constants means there exists element $K$ such that $C_{i}(K)$ cannot be smaller than the lower bound provided here. To have rigorous computation results, the INTLAB toolbox of interval arithmetic Ru99a is adopted for the evaluation of the constants. Particularly, the rigorous eigenvalue estimation for matrices is based on the algorithm of Behnke behnke1991calculation .

The method to be proposed in this paper for estimating constant $C_{0}$ and $C_{1}$ can also be used to estimate the constants of other interpolation operators. For example, let $\Pi$ be the Lagrange interpolation operator or the Fujino–Morley interpolation operator defined over triangle $K$, the following interpolation error estimation holds.

[TABLE]

In §6, sharp evaluation of constants $d_{0}$ and $d_{1}$ is provided; see detailed results in Table 1 and Table 2.

The rest of this paper is arranged as follows. In §2, constants $C_{0}$ and $C_{1}$ and the function space setting are introduced. The method to solve eigenvalue problems corresponding to constants is presented in §3. The theoretical analysis about the perturbation of eigenvalues on element shape variation is performed in §4. In §5, the algorithm to bound the constants for elements of arbitrary shapes is proposed and the optimal estimation is obtained. In §6, the optimal bound for $C_{0}$ and $C_{1}$ is applied in bounding error constants for other interpolation operators.

2 Preliminary

Function spaces

The standard notation for Sobolev space is used in this paper. That is, $\|\cdot\|_{\Omega}$ denotes the $L^{2}$ norm for $L^{2}$ space; $|\cdot|_{k,\Omega}$ $(k=1,2,\cdots)$ denotes the $k$th order semi-norm for functions in $H^{k}(\Omega)$. In many cases, the subscript $\Omega$ will be omitted if the domain is self-evident. The gradient operator is denoted by $\nabla$ and the second order derivative is given by $D^{2}u:=(u_{xx},u_{xy},u_{yx},u_{yy})$ for $u\in H^{2}(\Omega)$. The inner product of $L^{2}(\Omega)$ or $\left(L^{2}(\Omega)\right)^{2}$ is denoted by $(\cdot,\cdot)$.

Definition of constants

To give the definition of the interpolation constants $C_{0}$ and $C_{1}$, let us introduce the kernel space of $\Pi^{\mbox{\tiny{FM}}}$,,denoted by $V^{\mbox{\tiny{FM}}}(K)$, when the operator is applied to a triangle element $K$. That is,

[TABLE]

Over space $V^{\mbox{\tiny{FM}}}(K)$, the constants are defined by using the Rayleigh quotient.

[TABLE]

Given a reference triangle $K$ with $\text{diam}(K)=1$, let $K_{h}$ be the triangle obtained through scaling $K$ by $h$ times, that is, $\text{diam}(K_{h})=h$. It is easy to see that $C_{0}(K_{h})=h^{2}C_{0}(K)$, $C_{1}(K_{h})=hC_{1}(K)$.

Below, we also introduce the Crouzeix-Raviart interpolation constant, which will help to find the optimal estimation of the constant $C_{1}$.

Crouzeix-Raviart interpolation constant

Given $u\in H^{1}(K)$, the Crouzeix-Raviart interpolation $\Pi^{\mbox{\tiny{CR}}}u$ is a linear polynomial over $K$ such that

[TABLE]

The Crouzeix-Raviart constant associated to $\Pi^{\mbox{\tiny{CR}}}$ is defined as follows,

[TABLE]

where

[TABLE]

The Crouzeix-Raviart interpolation constant $C_{\mbox{\tiny{CR}}}$ is well investigated in Liu2015AMC :

1. a)

When vertices $O$ and $A$ of $K$ are fixed, the value of $C_{\mbox{\tiny{CR}}}(K)$ has monotonicity upon the $y-$coordinate of $B$. 2. b)

For all triangles with diameter less than $1$, the maximum value of $C_{\mbox{\tiny{CR}}}(K)$ has a rigorous bound as follows,

[TABLE] 3. c)

Numerical computation implies the maximum value is achieved when $K$ is a regular triangle.

For the relation between $C_{1}(K)$ and $C_{\mbox{\tiny{CR}}}(K)$, we have the following lemma (see, e.g., Carstensen2013 ).

Lemma 2.1

$C_{1}(K)\leq C_{\mbox{\tiny{CR}}}(K)$* for all triangle $K$.*

Proof. One can draw the conclusion by noticing that $u_{x},u_{y}\in V^{\mbox{\tiny{CR}}}(K)$ for any $u\in V^{\mbox{\tiny{FM}}}(K)$ and the inequalities $\|u_{x}\|\leq C_{\mbox{\tiny{CR}}}(K)\|\nabla u_{x}\|$, $\|u_{y}\|\leq C_{\mbox{\tiny{CR}}}(K)\|\nabla u_{y}\|$.

The task to obtain the optimal estimation of constants can be divided into two steps.

1. Step 1.

Direct evaluation for $C_{0}$ and $C_{1}$ for several sample triangles by solving the corresponding eigenvalue problem along with the Fujino-Morley FEM; see the detail in §3. 2. Step 2.

Perturbation analysis of constant $C_{0}$ and $C_{1}$ upon the change of shape of triangle element, where the monotonicity of constant $C_{0}$ and $C_{\mbox{\tiny{CR}}}$ will play an important role; see the detail in §4 and §5.

For the purpose of simplicity, the vertices of $K$ are located at $O(0,0)$, $A(1,0)$, $B(a,b)$, while $|OB|\leq 1$ and $|AB|\leq 1$. Also, due to the symmetry of the position of $B$, only the case that $a\geq 1/2$ will be considered. As a summary, we will focus on triangles with $B$ inside the following area $\Omega$ (see Fig. 2):

[TABLE]

3 Point-wise estimation of constants $C_{0}(K)$ and $C_{1}(K)$

In this section, let us describe the algorithm to estimate constants $C_{0}(K)$ and $C_{1}(K)$. The interpolation constants are determined by solving eigenvalue problems of bi-harmonic operators, and their upper bounds will be evaluated by applying the Fujino-Morley FEM. The lower bounds of the interpolation constants will be discussed in §5.3, where the conforming polynomial spaces are used.

The constants $C_{0}$ and $C_{1}$ are corresponding to the eigenvalues of the following eigenvalue problems.

Problem a)

Find $u\in V^{\mbox{\tiny{FM}}}(K)(\subset H^{2}(K))$ and $\lambda>0$ such that

[TABLE]

Problem b)

Find $u\in V^{\mbox{\tiny{FM}}}(K)(\subset H^{2}(K))$ and $\mu>0$ such that

[TABLE]

The distribution of eigenpairs of Problem a) and Problem b) has been well investigated under the theories for compact self-adjoint differential operators; see, e.g., Barbuska1991 . Let $\lambda_{1}$ and $\mu_{1}$ be the smallest eigenvalue of Problem a) and b), respectively. Then, it is easy to see that

[TABLE]

To give explicit bound of $\lambda_{1}$ and $\mu_{1}$, let introduce several finite element spaces. Let $\mathcal{T}^{h}$ be a proper triangulation of domain $K$. A general Fujino–Morley finite element space $V^{h}$ has the member function $u_{h}$ with the following properties.

1. a)

$u_{h}$ is piece-wise quadratic polynomial on each element of $\mathcal{T}^{h}$; 2. b)

$u_{h}$ is continuous on each vertex; 3. c)

$\displaystyle{\int_{e}\frac{\partial u}{\partial n}\>\mbox{d}s}$ is continuous across each interior edge $e$.

To approximate $V^{\mbox{\tiny{FM}}}$ (see definition at (5)), let us define a subspace of $V^{h}$ by introducing boundary conditions.

[TABLE]

[TABLE]

Next, let us define the approximate eigenvalue problems over $V^{\mbox{\tiny{FM}},h}(K)$. Since the function in $V^{\mbox{\tiny{FM}},h}(K)$ may not have continuity across interior edges, the differential operators $\nabla$ and $D^{2}$ are piece-wisely defined on the triangulation $\mathcal{T}^{h}$.

Problem a’)

Find $u\in V^{\mbox{\tiny{FM}},h}(K)$ and $\lambda_{h}>0$ such that

[TABLE]

Problem b’)

Find $u\in V^{\mbox{\tiny{FM}},h}(K)$ and $\mu_{h}>0$ such that

[TABLE]

Denote the smallest eigenvalue of a’) and b’) by $\lambda_{h,1}$ and $\mu_{h,1}$, respectively.

Below, we quote Theorem 2.1 of Liu2015AMC for the purpose of bounding eigenvalues, where the spaces and bilinear forms are taken as follows.

[TABLE]

[TABLE]

Particularly, $M(u,v)$ is an inner product for $V(h)$.

Theorem 3.1

Let $P_{h}:V(h)\mapsto V^{\mbox{\tiny{FM}},h}$ be the projection with respect to inner product $(D^{2}\cdot,D^{2}\cdot)$, i.e., for any $u\in V(h)$

[TABLE]

Suppose there exist quantities $C_{h,0}$ and $C_{h,1}$ such that for any $u\in V(h)$

[TABLE]

Then we have,

[TABLE]

Remark 1

The above theorem does not require $V^{h}\subset V$. Thus, we can use non-conforming finite element methods to obtain lower eigenvalue bounds.

Due to the special setting of the nonconforming space $V^{h}$, the projection $P_{h}$ here is nothing else but the Fujino–Morley interpolation defined over $\mathcal{T}^{h}$. That is,

[TABLE]

The error constants $C_{0,h}$ and $C_{1,h}$ are given by

[TABLE]

From existing estimation of $C_{0}$ and $C_{1}$, as quoted in (3) and (4), we have

[TABLE]

Therefore, the explicit lower bounds of $\lambda_{1}$ and $\mu_{1}$ are given as

[TABLE]

Remark 2

The technique of K. Kobayashi in Kobayashi-2015-CM can be applied here to provide upper bound of $C_{0}$ (lower bound of $\lambda_{1}$) directly without a prior information of $C_{0}$, i.e., the rough bound in (3) and (4).

4 Variation of constants upon perturbation of triangle

shape

In the last section, Theorem 3.1 provides a method to bound constants for a concrete triangle $K$. To bound constant $C_{0}$ and $C_{1}$ for triangles of arbitrary shapes, we need to consider the variation of constants upon the perturbation of triangle shape.

Linear perturbation of triangle $K$

Define linear mapping $Q$ by

[TABLE]

where $\alpha\approx 0$ and $\beta\approx 1$. Apply $Q$ to the point $(x,y)^{T}$ of $K$ by

[TABLE]

and $\widetilde{K}$ is the triangle with vertices $O(0,0),A(1,0)$ and $\widetilde{B}(x+\alpha y,\beta y)$. For the variation of norms of $u$ on $K$ under $Q$, we have the following lemma.

Lemma 4.2

Given $u\in H^{2}(K)$ and define $\tilde{u}=u\circ Q^{-1}\in H^{2}(\widetilde{K})$, where Q=\left(\begin{array}[]{cc}1&\alpha\\ 0&\beta\end{array}\right), and $\beta>0$. Define $\gamma=\alpha^{2}+\beta^{2}+1$. Then

1. (a)

For $L^{2}(K)$-norm, we have

[TABLE] 2. (b)

For $H^{1}(K)$-norm, we have

[TABLE] 3. (c)

For the $H^{2}(K)$-norm, we have

[TABLE]

Proof. The equality of (a) is evident. Let $K$ be the triangle with fixed vertices $O(0,0)$ and $A(1,0)$. Let us introduce the linear transform of $K$. Since $(u_{x},u_{y})^{t}=Q\cdot(\tilde{u}_{\tilde{x}},\tilde{u}_{\tilde{y}})^{t}$, we have

[TABLE]

where $\lambda_{\min}(Q^{t}Q)$ and $\lambda_{\max}(Q^{t}Q)$ denoted the minimum and maximum eigenvalues of $Q^{t}Q$, respectively. Therefore

[TABLE]

The eigenvalues of $Q^{t}Q$ are listed below:

[TABLE]

where $\gamma=\alpha^{2}+\beta^{2}+1$.

For the second order derivatives, we have

[TABLE]

where

[TABLE]

We denote the minimum and maximum eigenvalues of $T^{t}T$ by $\lambda_{\min}(T^{t}T)$ and $\lambda_{\max}(T^{t}T)$, respectively. For the $H^{2}(K)$ norm, we have

[TABLE]

where the minimum and maximum eigenvalues are

[TABLE]

$\square$

Remark 3

Note that there are repeated eigenvalue $\beta^{2}$ of positive definite matrix $T^{t}T$ with algebraic multiplicity 2. It’s easily to verify $\lambda_{\max}(T^{t}T)>\beta^{2}>\lambda_{\min}(T^{t}T)>0$ by using the fact $\gamma^{2}-4\beta^{2}=(\gamma-2\beta^{2})^{2}+4\alpha^{2}\beta^{2}>0.$

Below, we show several results for special values of $\alpha$ and $\beta$.

Lemma 4.3

Let $\alpha=0$, $\beta=1+\epsilon$. For $\epsilon>0$, we have

[TABLE]

and

[TABLE]

For $-1<\epsilon<0$, we have

[TABLE]

and

[TABLE]

Next lemma will be used in case of the perturbation of $B$ along the arc with center as $O$ and radius as $|AB|$.

Lemma 4.4

For $0<\theta<\pi$, $0<\theta+\tau<\pi$, define linear mapping $Q$ with the following $\alpha$ and $\beta$,

[TABLE]

Then, in case $\tau<0$,

[TABLE]

where $\rho$ and $\eta$ are defined by

[TABLE]

In case $\tau>0$, the above inequalities hold by exchanging the value of $\rho$ and $\eta$.

Proof

This lemma is a direct result of Lemma 4.2 by noticing the following relations.

[TABLE]

4.1 Properties of constant $C_{0}(K)$

Theorem 4.2

Let $K$ be the triangle with vertices $O(0,0)$, $A(1,0)$, $B(a,b)\in\Omega$, then $C_{0}(K)$ monotonically increases as $b$ increases.

Proof

Take $\alpha=0$, $\beta=1+\epsilon$ ($\epsilon>0$) for transformation $Q$ in (17). Then $QB$ moves $B$ along y-direction and

[TABLE]

Then we can easily draw the conclusion from the definition of the constant.

Below we consider the perturbation of $B=(\cos\theta,\sin\theta)$ along $\theta$ direction. By using Lemma 4.2 and 4.4, we can easily obtain the following result.

Theorem 4.3

Let $K$ be the triangle with vertices $O(0,0)$, $A(1,0)$, $B(\cos\theta,\sin\theta)$ ($B\in\Omega$) and $\widetilde{K}$ be the triangle with vertices $O$, $A$ and $\widetilde{B}=(\cos(\theta+\tau),\sin(\theta+\tau))$ ($\widetilde{B}\in\Omega$). We have the following two inequalities:

[TABLE]

4.2 Properties of $C_{\mbox{\tiny{CR}}}(K)$ and $C_{1}(K)$

As shown in (4), an upper bound of $C_{1}(K)$ is known already via the estimation of $C_{\mbox{\tiny{CR}}}(K)$. For the purpose of a more accurate estimation for $C_{1}(K)$, the perturbation of $C_{1}(K)$ respect to $B$ will be required in next section.

As a preparation, let us consider the perturbation analysis of $C_{\mbox{\tiny{CR}}}(K)$ upon moving $B$ along the $x-$direction.

Theorem 4.4

Let $K$ be the triangle with vertices $O(0,0)$, $A(1,0)$, $B(a,b)$ ($B\in\Omega$) and $\widetilde{K}$ be the triangle with vertices $O$, $A$ and $\widetilde{B}(a+b\epsilon,b)$ ($\widetilde{B}\in\Omega$) with the condition $|\epsilon|<1/2$. Then,

[TABLE]

Proof

Take $\alpha=\epsilon$, $\beta=1$ for transformation $Q$ in (17). Since for $\tilde{v}\in V^{\mbox{\tiny{CR}}}(\widetilde{K})$, $v=\tilde{v}\circ Q\in V^{\mbox{\tiny{CR}}}(K)$, from Lemma 4.2 , we have

[TABLE]

By using Eq. (8) and the above two equations, we have

[TABLE]

In case $\epsilon>0$, we have

[TABLE]

and in case of $\epsilon<0$ and $|\epsilon|<1/2$, we have

[TABLE]

Thus we can draw the conclusion.

Remark 4

In the above proof, we use polynomial of $\epsilon$ to simplify the expression in (20). The estimation of the constant has the exact order up to linear term $\epsilon$, while the coefficients of quadratic term $\epsilon^{2}$ are overestimated for the purpose of a simple expression.

Theorem 4.5

Let $K$ be the triangle with vertices $O(0,0)$, $A(1,0)$, $B(a,b)$ ($B\in\Omega$) and $\widetilde{K}$ be the triangle with vertices $O(0,0)$, $A(1,0)$ and $\widetilde{B}(a+\alpha b,(1+\epsilon)b)$ ($\widetilde{B}\in\Omega$). Assume $|\alpha|,|\epsilon|\leq 1/2$, we have

[TABLE]

[TABLE]

Proof

Let Q=\left(\begin{array}[]{cc}1&\alpha\\ 0&1+\epsilon\end{array}\right) be the linear mapping from $K$ to $\widetilde{K}$. Usually, $\partial u/\partial n$ varies under the transform $Q$. However, for ${v}\in V^{\mbox{\tiny{FM}}}({K})$, we still have $\tilde{v}={v}\circ Q^{-1}\in V^{\mbox{\tiny{FM}}}(\widetilde{K})$. To see this, notice that the condition ${v}(O)={v}(A)={v}(B)=0$ implies $\int_{e_{i}}\partial v/\partial\tau\>\mbox{d}s=0$, where $\tau$ is the tangent vector of $e_{i}$. From the condition $\int_{e_{i}}\partial{v}/\partial n\>\mbox{d}s=0$, we have

[TABLE]

where $n$ is the unit norm vector on $e_{i}$. Hence,

[TABLE]

From Lemma 4.2, we have, with $\gamma=\alpha^{2}+(1+\epsilon)^{2}+1$,

[TABLE]

Define

[TABLE]

By taking the Taylor expansion of $\sqrt{f_{2}/f_{1}}$ at $\alpha=0,\epsilon=0$, we have, for $|\alpha|,|\epsilon|\leq 1/2$,

[TABLE]

[TABLE]

By the definition of constant $C_{1}$, we can draw the conclusion.

Remark 5

In the estimation (23) and (24), the coefficients of $|\alpha|$ and $|\epsilon|$ agree with the ones in Taylor expansion, while the coefficients of $|\alpha|^{2}$ and $|\epsilon|^{2}$ are over-estimated for the purpose of simple expression.

5 Optimal estimation of constants

From the monotonicity of $C_{0}(K)$ on $y$-coordinate of $B$ as shown in Theorem 4.2, the maximum value of $C_{0}(K)$ can only happens on the arc $r=1$, $0<\theta\leq\pi/3$. Since the value of $C_{0}(K)$ depends on the $x$- and $y$-coordinate of $B$, $C_{0}(K)$ can be regarded as a function on $x_{B}$, $y_{B}$. Fig. 3 displays the contour lines of $C_{0}(K)$ respect to $x_{B},y_{B}$ in $[0.5,1]\times(0,1]$. Numerical estimation of $C_{0}(K)$ implies that the maximum is taken when $K$ is the regular triangle.

5.1 Optimal estimation of $C_{0}$

On arc $r=1$, $0<\theta\leq\pi/3$, we perform point-wise evaluation of $C_{0}(K)$ on a subdivision of $\theta$. Then the bound of $C_{0}(K)$ on whole arc is obtained by applying perturbation theory in Theorem 4.3.

Define $\theta_{i}$ and $\tau_{i}$ by

[TABLE]

[TABLE]

For each $\theta_{i}$, we evaluate the constant $C_{0}(K)$, then using the Theorem 4.3 to give a upper bound of $C_{0}(K)$ on each sub-interval $(\theta_{i}-\tau_{i},\theta_{i}]$ with perturbation $\tau_{i}.$

Fig. 4 displays the upper bound of $C_{0}(K)$, where $x$-coordinate denote the size of $\angle BOA$ with range as $(0,\pi/3]$. The verified computing results show that

[TABLE]

5.2 Optimal estimation of $C_{1}$

The estimation of $C_{0}(K)$ is relatively easily done, which is thanks to the monotonicity of the constant $C_{0}(K)$ with respect to the $y$-coordinate of vertex $B$ of $K$. However, such property of monotonicity is not available for $C_{1}(K)$. Thus, one has to consider the case of collapsed triangles (the vertex $B$ being close to $x$-axis), which is difficult to process because the estimation of Theorem 4.5 has divergent bound for small $y$-coordinate of vertex $B(a,b)$. To avoid such difficulties, we subdivide the area into two parts: $\Omega=\Omega_{1}\cup\Omega_{2}$ (see Fig. 6). On $\Omega_{1}$, we can estimate $C_{1}$ directly; on $\Omega_{2}$, the estimation of $C_{1}$ is done through $C_{\mbox{\tiny{CR}}}(K)$ for collapsed triangles by using the relation $C_{1}(K)\leq C_{\mbox{\tiny{CR}}}(K)$.

The approximate numerical evaluation of $C_{1}(K)$ and $C_{\mbox{\tiny{CR}}}(K)$ are displayed in Fig. 5, from which we have the following information.

The maximum value of $C_{1}(K)\approx 0.18868$ happens when $B$ tends to $(1,0)$ 2. 2.

The maximum value of $C_{\mbox{\tiny{CR}}}(K)\approx 0.18932$ happens at $B=(1/2,\sqrt{3}/2)$. 3. 3.

$C_{\mbox{\tiny{CR}}}(K)\approx 0.18868$ when $B$ tends to $(1,0)$. 4. 4.

The constant $C_{\mbox{\tiny{CR}}}(K)$ gives a nice bound (in the sense that the value of $C_{\mbox{\tiny{CR}}}(K)$ is less than $0.18868$) for $C_{1}(K)$ when $B$ is below $y=\sin(0.98\pi/3)$.

The above properties of $C_{1}(K)$ implies the optimal estimate value of $C_{1}$ can be done by following the strategy below:

1. Step 1.

Separate the region $\Omega$ into parts $\Omega_{1}$ and $\Omega_{2}$ as follows(see Fig. 6):

[TABLE]

2. Step 2.

In $\Omega_{1}$, by evaluating the constant $C_{1}$ for several position of $B$ and apply the perturbation result in Theorem 4.5, we will have

[TABLE] 3. Step 3.

In $\Omega_{2}$, the upper bound of $C_{1}(K)$ will be obtained through $C_{\mbox{\tiny{CR}}}(K)$,

[TABLE]

Below, we show the details to obtain (25) and (26).

Estimation of $C_{1}$ in $\Omega_{1}$

Let $K_{0}$ be the triangle with vertices at $(0,0)$, $(1,0)$, $\displaystyle B_{0}=(1/2,\sqrt{3}/2)$. Evaluation of upper bound of $C_{1}$ tells that $\displaystyle C_{1}(K_{0})<0.1744$. Notice that for all $x,y\in\Omega_{1}$, $|x-x_{B_{0}}|\leq 0.021y_{B_{0}}$ and $|y-y_{B_{0}}|\leq 0.013y_{B_{0}}$. By applying estimation (22) in Theorem 4.5 with $\alpha=0.021$ and $\epsilon=0.013$, we have

[TABLE]

Estimation of $C_{1}$ in $\Omega_{2}$

Numerical results implies $C_{1}(K)$ and $C_{\mbox{\tiny{CR}}}(K)$ have the same supremum on $\Omega_{2}$, which is reached when $B$ tends to $(1,0)$. Since $C_{1}(K)\leq C_{\mbox{\tiny{CR}}}(K)$ holds strictly, we just focus on $C_{\mbox{\tiny{CR}}}(K)$. Due to the monotonicity of $C_{\mbox{\tiny{CR}}}(K)$ on $y_{B}$, the maximum value of $C_{\mbox{\tiny{CR}}}(K)$ is taken on either of the following two boundaries.

[TABLE]

• •

On $\Gamma_{1}$, we subdivide $\Gamma_{1}$ into small closed intervals:

[TABLE]

Here, $\displaystyle{h=\frac{1}{20}\left(\cos\left(\frac{0.98\pi}{3}\right)-\frac{1}{2}\right)}$.

The estimation of $C_{\mbox{\tiny{CR}}}(K)$ on each $I_{i}$ is displayed in Fig. 7. The circles in Fig. 7 denote the point-wise evaluation of $C_{\mbox{\tiny{CR}}}(K)$ and the short bars denote the upper bound of $C_{\mbox{\tiny{CR}}}(K)$ based on Theorem 4.4. The computation results tell that

[TABLE]

• •

To estimation of $C_{\mbox{\tiny{CR}}}$ on $\Gamma_{2}$, we take the subdivision of $\theta$ as follows

[TABLE]

Define $\tau_{1}=\theta_{1}$, $\tau_{i}=\theta_{i}-\theta_{i-1}\quad(i=2,\ldots,98)$. From the estimation in Theorem 4.2 of LiuLiu2015AMC , we can estimate $C_{\mbox{\tiny{CR}}}(K)$ for $\theta$ in each interval $(0,\theta_{1}]$ and $(\theta_{i-1},\theta_{i}]$ $(i=1,2,\ldots,98)$. In Fig. 8, the point wise estimation of $C_{\mbox{\tiny{CR}}}(K)$ on $\theta_{i}$ and the upper bound of $C_{\mbox{\tiny{CR}}}(K)$ for each interval are displayed. Computation results indicate that the supremum of $C_{\mbox{\tiny{CR}}}(K)$ is reached when $B$ tends to $(1,0)$. We have the following strict bound of $C_{\mbox{\tiny{CR}}}(K)$

[TABLE]

The estimation (27), (28) and (29) validate the results in (25) and (26). Thus, we can draw the conclusion that

[TABLE]

Remark 6

The asymptotic value of $C_{1}$ when the vertex $B=(a,b)$ of $K$ tends to $(1,0)$ can be determined by theoretical analysis. Below is a sketch of determining $\lim_{\theta\to 0}C_{1}(K)$. With the analogous argument as in Theorem 5 of LIU-Kikuchi , the eigenvalue $\lambda$ corresponding to $C_{1}$ is determined by solving the following eigenvalue problem on one dimensional interval (0,1),

[TABLE]

The general solution of this ODE is given by utilizing hypergeometric function,

[TABLE]

Here, $J_{0}$ is the [math]-th order Bessel function of the first kind and $~{}_{1}F_{2},~{}_{2}F_{3}$ are hypergeometric functions. By further applying the two boundary conditions, one can obtain the following equation with $\lambda$.

[TABLE]

The above equation has infinite solutions and simple computation tells that the smallest one, denoted by $\lambda_{1}$, is given by

[TABLE]

Thus, $C_{1}(K)$ converges to $\displaystyle C_{1}=\frac{1}{\sqrt{\lambda_{1}}}=0.1886407440\cdots$ when $B$ tends to $(1,0)$.

5.3 Lower bound of constants

To confirm the precision of obtained estimation for $C_{0}$ and $C_{1}$, we also calculate the lower bounds of constants. The task to provide lower bound is easy compared with the upper bound estimation. By evaluating $C_{0}(K)$ and $C_{1}(K)$ for a concrete element $K$ with conforming spaces, then we can have the lower bound for optimal constants.

Given a triangle $K$, the conforming finite dimensional space $W^{m}(\subset V_{0})$ over $K$ can be constructed by using polynomials.

[TABLE]

Here, $P^{m}(K)$ denotes the space of polynomial over $K$ with degree up to $m$.

For $C_{0}$, we choose $K$ as the unit regular triangle and solve Problem a) in $W^{7}$. Numerical computation tells that

[TABLE]

Similarly, for $C_{1}$, by taking the triangle with $B=(\cos 0.01,\sin 0.01)$ and solving Problem b) in $W^{6}$, we have

[TABLE]

6 Application to eigenvalue problems of Biharmonic operators

Let us apply the two fundamental constants $C_{0}$ and $C_{1}$ to estimate the error constants appearing in the Lagrange interpolation and the Fujino–Morley interpolation.

The Lagrange interpolation error constants

Given triangle element $K$ with vertices $P_{1}(0,0),$ $P_{2}(1,0)$ and $P_{3}(a,b)$. Let $\Pi_{1}$ be the Lagrange interpolation over $T$ such that, for $v\in H^{2}(K)$, $\Pi_{1}v$ is a linear polynomial and $(v-\Pi_{1}v)(P_{i})=0$, $i=1,2,3$.

Define subspace $V$ of $H^{2}(K)$ by

[TABLE]

Define constants $C_{L,0}(K),C_{L,1}(K)$ by

[TABLE]

[TABLE]

Thus, the interpolation error estimation of $\Pi_{1}$ can be given as

[TABLE]

The constants are determined by solving the following eigenvalue problems:

Problem c) Find $u\in V_{0}(K)$ and $\eta>0$ such that

[TABLE]

Problem d) Find $u\in V_{0}(K)$ and $\mu>0$ such that

[TABLE]

Denote the smallest eigenvalue of each problem by $\eta_{1}$ and $\mu_{1}$, respectively. Then $C_{L,0}=1/\sqrt{\eta_{1}}$ and $C_{L,1}=1/\sqrt{\mu_{1}}$. Choose the subspace of Fujino–Morley finite element space $V^{h}$ defined over triangulation of $K$,

[TABLE]

Let $\eta_{h,1}$ and $\mu_{h,1}$ be the smallest eigenvalues of Problem c) and d), respectively, with $V_{0}(K)$ replaced by $V_{0}^{h}(K)$.

Let us take the following setting for Theorem 2.1 of Liu2015AMC , which is similar to Theorem 3.1,

[TABLE]

With the newly obtained error constant estimation in §5, we obtain the lower bounds of $\eta_{1}$ and $\mu_{1}$ like (16).

[TABLE]

The Fujino–Morley interpolation error constants

For the Fujino–Morley interpolation $\Pi^{\mbox{\tiny{FM}}}$, let us define the following constants

[TABLE]

[TABLE]

Notice that the minimizer function for $\eta$ and $\mu$ are both orthogonal to all $P_{2}$ polynomial functions with respect to $(D^{2}\cdot,D^{2}\cdot)$. To see this, one can take the perturbation of the minimizer function with respect to any function in $P_{2}(K)$.

Let us introduce the space $U_{0}(K)$ by

[TABLE]

Then $\eta$ and $\mu$ can be characterized by Rayleigh quotients over $U_{0}(K)$,

[TABLE]

Notice that $(D^{2}\cdot,D^{2}\cdot)$ and $((I-\Pi^{FM})\cdot,(I-\Pi^{FM})\cdot)$ are both positive definite bilinear forms on $U_{0}(K)$. Thus, we can apply Theorem 2.1 of Liu2015AMC to solve corresponding eigenvalue problems.

Concrete values of interpolation error constants

In Table 1 and 2, we list the estimation of upper bounds of $C_{L,i}$ and $C_{{\mbox{\tiny{FM}}},i}$ for different shapes of triangle $K$. The underline in tables tells that the lower bound and upper bound evaluation of the constants agree with each other at the underlined digits.

Remark 7

The estimation of constants considered in this paper only concerns the largest edge length of a triangle element. By utilizing more geometric information of the element, i.e., the inner angle size and each edge length, one can have better upper bounds of interpolation constants over domain of different shapes; see the work of Liu-Kikuchi Kikuchi+Liu2007 ; LIU-Kikuchi and Kobayashi Kobayashi-2015-CM .

7 Summary

In this paper, we provide optimal estimation of two important error constants in bounding eigenvalues of bi-harmonic operators. As application of obtained estimation of the constants, the upper bounds for error constants of two interpolation operators are evaluated. Moreover, the algorithm proposed here along with the explicit constant values can be further used to give rigorous bounds for the eigenvalues of general bi-harmonic differential operators.

Acknowledgements.

The authors would like to thank for the support from the Ministry of Science and Technology, Taiwan, R.O.C. under Grant no. MOST 107-2911-M-006-506. This research is also supported by Japan Society for the Promotion of Science, Grand-in-Aid for Young Scientist (B) 26800090, Grant-in-Aid for Scientific Research (C) 18K03411 and Grant-in-Aid for Scientific Research (B) 16H03950 for the third author.