# Approximation of the biharmonic problem using piecewise linear finite elements

Robert Eymard, Raphaèle Herbin

► To cite this version:

Robert Eymard, Raphaèle Herbin. Approximation of the biharmonic problem using piecewise linear finite elements. Comptes Rendus. Mathématique, 2010, 348 (23-24), pp.1283-1286. 10.1016/j.crma.2010.11.002. hal-00693004

HAL Id: hal-00693004

https://hal.science/hal-00693004

Submitted on 31 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Approximation of the biharmonic problem using piecewise linear finite elements

## Approximation d'un problème biharmonique par élément fini P1

Robert Eymard<sup>a</sup>, Raphaèle Herbin<sup>b,1</sup>

### ABSTRACT

We propose an approximation of the solution of the biharmonic problem in $H_0^2(\Omega)$ which relies on the discretization of the Laplace operator using nonconforming continuous piecewise linear finite elements.

### RÉSUMÉ

Nous proposons une approximation de la solution du problème bi-harmonique dans $H_0^2(\Omega)$ basée sur la discrétisation du Laplacien par éléments finis P1 continus mais non conformes.

## 1. Introduction

A large number of discretization methods for the biharmonic operator has been proposed in the recent past. The most classical is probably the conforming finite element method. For fourth order problems, the conforming finite element space must be a finite dimensional subspace of the Sobolev space $H^2(\Omega)$. Hence elementary basis functions are sought such that the reconstructed global basis functions on $\Omega$ belong to $C^1(\bar{\Omega})$. On Cartesian meshes, such basis functions are found by generalizing the one-dimensional $P^3$ Hermite finite element to the multi-dimensional framework. This task becomes much more difficult on more general meshes and involves rather sophisticated finite elements such as the Argyris finite element on triangles in 2D, which unfortunately requires 21 degrees of freedom [7]. Hence nonconforming FEMs have also been widely studied: see e.g. [7, Section 49], [8], and references therein, and [3,4] for more recent works. Discontinuous Galerkin methods have also been recently developed and analysed [10–13]; error estimates have been derived for polynomials of degree greater or equal to two or three. Other methods which have been developed for fourth order problems include mixed methods [5] (see also references therein), [14], and compact finite difference methods [6,2,1]. All the above methods are high order methods, and therefore, rather computationally expensive and may not be so easy to implement. Recently, a cheaper low order method based on the discretization of the Laplace operator by a cell centred finite volume scheme was proposed [9]. The idea in the present paper is to show that the discretization of the Laplace operator by the piecewise linear finite element method is also possible in order to obtain a convergent low order scheme for the biharmonic operator.

<footnote><sup>a</sup> Laboratoire d'analyse et de mathématiques appliquées, UMR CNRS 8050, Université Paris-Est, 5, boulevard Descartes, Champs-sur-Marne, 77454 Marne-la-Vallée cedex 2, France</footnote>

<footnote><sup>b</sup> Laboratoire d'analyse, topologie et probabilités, UMR CNRS 6632, Université de Provence, 39, rue Joliot-Curie, 13453 Marseille, France</footnote>

<footnote>E-mail addresses: robert.eymard@univ-paris-est.fr (R. Eymard), raphaele.herbin@cmi.univ-mrs.fr (R. Herbin).</footnote>

<footnote n="1">¹ Work supported by Groupement de recherche MOMAS, PACEN/CNRS.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Definition of the domain.</figcaption></figure>

Throughout this paper, $d \in \mathbb{N} \setminus \{0\}$ denotes the space dimension, $\Omega$ is an open polygonal bounded and connected subset of $\mathbb{R}^d$, with Lipschitz-continuous boundary $\partial\Omega$, and $f \in L^2(\Omega)$, $\ell \in L^2(\Omega)$ and $g \in (L^2(\Omega))^d$. Under these assumptions, we consider the following problem: find

$$u \in H_0^2(\Omega); \quad \forall v \in H_0^2(\Omega), \\ \int_{\Omega} \Delta u(\mathbf{x}) \Delta v(\mathbf{x}) \, d\mathbf{x} = \int_{\Omega} (f(\mathbf{x})v(\mathbf{x}) + g(\mathbf{x}) \cdot \nabla v(\mathbf{x}) + \ell(\mathbf{x}) \Delta v(\mathbf{x})) \, d\mathbf{x}, \quad (1)$$

where $H_0^2(\Omega)$ is the closure in $H^2(\Omega)$ of the set $C_c^\infty(\Omega)$ of infinitely continuously differentiable functions with compact support in $\Omega$.

## 2. Definition of the scheme

Let $\Omega$ be a polyhedral domain in $\mathbb{R}^d$, with $d \in \mathbb{N}^*$. We consider a classical simplicial mesh $\mathcal{T}$ of $\Omega$ [7, Section 6]; we denote by $h_T$ the supremum of the diameters of all $S \in \mathcal{T}$, and by $\theta_T$ the infimum of the angles of all triplets of vertices of the simplex $S$, for any simplex $S \in \mathcal{T}$.

Let $\mathcal{V}$ be the finite set of the vertices of the mesh, and $\mathcal{V}_S$ the set of the vertices of a given simplex $S \in \mathcal{T}$. Let $\mathcal{V}_{\text{ext}}$ be the set of all vertices of the mesh lying on the boundary and $\mathcal{V}_{\text{int}} = \mathcal{V} \setminus \mathcal{V}_{\text{ext}}$ the set of the interior vertices. For all $S \in \mathcal{T}$ and all $\mathbf{z} \in \mathcal{V}_S$, we denote by $K_{S,z}$ the subset of $S$ of all points whose barycentric coordinate related to $\mathbf{z}$ is greater than that related to any $\mathbf{z}' \in \mathcal{V}_S$ with $\mathbf{z}' \neq \mathbf{z}$ (see Fig. 1). We then denote for all $\mathbf{z} \in \mathcal{V}$ by $K_\mathbf{z}$ the union of all $K_{S,z}$, for all $S \in \mathcal{T}$ such that $\mathbf{z} \in \mathcal{V}_S$, which leads to $\bigcup_{\mathbf{z}' \in \mathcal{V}} K_{\mathbf{z}'} = \Omega$. We denote by $V_\mathbf{z}$ the set of all $\mathbf{y} \in \mathcal{V}$ such that there exists $S \in \mathcal{T}$ with $\{\mathbf{y}, \mathbf{z}\} \subset \mathcal{V}_S$. For any $\mathbf{z} \in \mathcal{V}$, let $\varphi_\mathbf{z}$ be the piecewise affine basis function of the P1 finite element, such that $\varphi_\mathbf{z}(\mathbf{z}) = 1$ and $\varphi_\mathbf{z}(\mathbf{z}') = 0$ for all $\mathbf{z}' \in \mathcal{V} \setminus \{\mathbf{z}\}$. We then denote by $V_T$ the vector space spanned by all functions $\varphi_\mathbf{z}$, $\mathbf{z} \in \mathcal{V}$ (then any $u \in V_T$ is such that $u = \sum_{\mathbf{z} \in \mathcal{V}} u(\mathbf{z})\varphi_\mathbf{z}$).

Next, the linear form $\Delta_\mathbf{z}: V_T \to \mathbb{R}$ for $\mathbf{z} \in \mathcal{V}$ is defined by

$$\forall u, v \in V_T, \quad \int_{\Omega} \nabla u(\mathbf{x}) \cdot \nabla v(\mathbf{x}) \, d\mathbf{x} = - \sum_{\mathbf{z} \in \mathcal{V}} |K_\mathbf{z}| v(\mathbf{z}) \Delta_\mathbf{z} u,$$

and the discrete Laplace operator $\Delta_T: V_T \to L^2(\Omega)$ is defined by

$$\forall u \in V_T, \quad \Delta_T u(\mathbf{x}) = \sum_{\mathbf{z} \in \mathcal{V}} \Delta_\mathbf{z} u \mathbf{1}_{K_\mathbf{z}}(\mathbf{x}), \quad \text{for a.e. } \mathbf{x} \in \Omega.$$

Note that, if we denote

$$T_{zy} = - \int_{\Omega} \nabla \varphi_{\mathbf{z}}(\mathbf{x}) \cdot \nabla \varphi_{\mathbf{y}}(\mathbf{x}) \, d\mathbf{x},$$

we get, using the property $\sum_{\mathbf{y} \in \mathcal{V}} T_{zy} = 0$ since $\sum_{\mathbf{y} \in \mathcal{V}} \varphi_{\mathbf{y}}(\mathbf{x}) = 1$ for all $\mathbf{x} \in \Omega$,

$$\Delta_\mathbf{z} u = \frac{1}{|K_\mathbf{z}|} \sum_{\mathbf{y} \in \mathcal{V}_\mathbf{z}} T_{zy}(u(\mathbf{y}) - u(\mathbf{z})), \quad \forall \mathbf{z} \in \mathcal{V}.$$
We then have

$$ \forall u \in V_T, \forall v \in V_T, \quad \int_{\Omega} \nabla u(\mathbf{x}) \cdot \nabla v(\mathbf{x}) \, d\mathbf{x} = - \int_{\Omega} P_T v(\mathbf{x}) \Delta_T u(\mathbf{x}) \, d\mathbf{x}, \qquad (2) $$

where $P_T v$ denotes the piecewise constant function defined by the values $v_z$ in $K_z$, for all $z \in \mathcal{V}$. The discrete approximation space is then defined by

$$ V_{T,0} = \{u \in V_T, u(\mathbf{z}) = 0 \text{ for all } \mathbf{z} \in \mathcal{V}_{\text{ext}}\}, $$

and the approximation scheme consists in finding

$$ u \in V_{T,0}; \forall v \in V_{T,0}, \quad \int_{\Omega} \Delta_T u(\mathbf{x}) \Delta_T v(\mathbf{x}) \, d\mathbf{x} = \int_{\Omega} (f(\mathbf{x})v(\mathbf{x}) + g(\mathbf{x}) \cdot \nabla v(\mathbf{x}) + \ell(\mathbf{x}) \Delta_T v(\mathbf{x})) \, d\mathbf{x}. \qquad (3) $$

### 3. Convergence analysis

The convergence analysis relies on the following sequence of results, the proof of which is given in a forthcoming paper:

**Lemma 3.1 (Existence, uniqueness and estimate on the solution of (3)).** There exists $C > 0$, only depending on $\Omega$ and $\theta$, such that for any $u \in V_{T,0}$ satisfying (3), then

$$ \begin{aligned} \|u\|_{L^2(\Omega)} &\le C (\|f\|_{L^2(\Omega)} + \|g\|_{L^2(\Omega)^d} + \|\ell\|_{L^2(\Omega)}), \\ \|\nabla u\|_{L^2(\Omega)^d} &\le C (\|f\|_{L^2(\Omega)} + \|g\|_{L^2(\Omega)^d} + \|\ell\|_{L^2(\Omega)}), \\ \|\Delta_T u\|_{L^2(\Omega)} &\le C (\|f\|_{L^2(\Omega)} + \|g\|_{L^2(\Omega)^d} + \|\ell\|_{L^2(\Omega)}). \end{aligned} $$

As a consequence, there exists one and only one $u \in V_{T,0}$ such that (3) holds.

The key point for the proof of the preceding result is the relation (2) with $v = u$, together with the inequality

$$ \|u - P_T u\|_{L^2(\Omega)} \le h_T \|\nabla u\|_{L^2(\Omega)^d}. $$

**Lemma 3.2 (Compactness of a sequence of approximate solutions).** Let $(T_m)_{m \in \mathbb{N}}$ be a sequence of classical simplicial discretizations of $\Omega$ [7, Section 6], such that $h_{T_m}$ tends to 0 as $m \to \infty$. Assume that there exists $\theta > 0$ with $\theta < \theta_{T_m}$ for all $m \in \mathbb{N}$. Let $(u_m)_{m \in \mathbb{N}}$ be a sequence of functions such that $u_m \in V_{T_m,0}$ for all $m \in \mathbb{N}$. Assume that the sequence $(\Delta_{T_m} u_m)_{m \in \mathbb{N}}$ is bounded in $L^2(\Omega)$; then there exists a subsequence of $(T_m)_{m \in \mathbb{N}}$, again denoted $(\tilde{T}_m)_{m \in \mathbb{N}}$, and $u \in H_0^2(\Omega)$, such that the corresponding subsequence $(u_m)_{m \in \mathbb{N}}$ satisfies:

(i) $u_m \to u \in L^2(\Omega)$,

(ii) $\nabla u_m \to \nabla u \in L^2(\Omega)^d$,

(iii) $\Delta_{T_m} u_m \to \Delta u$ weakly in $L^2(\Omega)$,

as $m \to \infty$.

The difficult point in the proof of the above lemma is the $H_0^2(\Omega)$ regularity of $u$. This point is proved using that, prolonging $u$ by 0 in $\mathbb{R}^d \setminus \Omega$, then $u \in H^1(\mathbb{R}^d)$ together with $\Delta u \in L^2(\mathbb{R}^d)$.

**Lemma 3.3 (Interpolation of regular functions with compact support).** Let $\theta > 0$ with $\theta < \theta_T$. Let $\varphi \in C_c^2(\Omega)$ and let $a = d(\text{support}(\varphi), \partial \Omega)$. Then there exists $\tilde{\mathcal{I}}_T \varphi \in V_{T,0}$ and $C > 0$ only depending on $\theta$ such that

(i) $\|\tilde{\mathcal{I}}_T \varphi - \varphi\|_{L^2(\Omega)} \le Ch_T \frac{|\varphi|_2}{a^2}$,

(ii) $\|\nabla \tilde{\mathcal{I}}_T \varphi - \nabla \varphi\|_{L^2(\Omega)^d} \le Ch_T \frac{|\varphi|_2}{a^2}$,

(iii) $\|\Delta_T \tilde{\mathcal{I}}_T \varphi - \bar{\Delta}_T \varphi\|_{L^2(\Omega)} \le Ch_T \frac{|\varphi|_2}{a^2}$, where $|\varphi|_2 = \max_{i,j=1,d} \|\partial_{ij}^2 \varphi\|_{L^\infty(\Omega)}$, and $\bar{\Delta}_T \varphi$ is the piecewise constant function equal in $K_z$ to $\frac{1}{|K_z|} \int_{K_z} \Delta \varphi(\mathbf{x}) \, d\mathbf{x}$, for all $z \in \mathcal{V}$.

In the preceding lemma, the function $\tilde{\mathcal{I}}_T \varphi$ is obtained by multiplying the approximate solution in $V_{T,0}$ of the Dirichlet problem with right-hand side $-\Delta \varphi$ by a regularization function whose second order derivatives are bounded by $C/a^2$.
**Theorem 3.4 (Convergence of the scheme).** Let $u \in H_0^2(\Omega)$ be the solution of Problem (1); let $\mathcal{T}$ be a conformal simplicial mesh of $\Omega$ and $u_{\mathcal{T}} \in V_{\mathcal{T},0}$ be the solution of (3). Then, as $h_{\mathcal{T}}$ tends to 0 with $0 < \theta \le \theta_{\mathcal{T}}$:

(i) $u_{\mathcal{T}}$ converges in $L^2(\Omega)$ to $u$,

(ii) $\nabla u_{\mathcal{T}}$ converges in $L^2(\Omega)^d$ to $\nabla u$,

(iii) $\Delta_{\mathcal{T}} u_{\mathcal{T}}$ converges in $L^2(\Omega)$ to $\Delta u$.

The proof of the convergence theorem is obtained by first extracting a sequence of approximate solutions thanks to Lemmas 3.1 and 3.2, then setting in (3) $v = \tilde{I}_{\mathcal{T}}\varphi$ given by Lemma 3.3 and letting $h_{\mathcal{T}}$ tend to 0.

**Remark 3.5 (Order of convergence).** An error estimate of order 1 may be proven if the exact solution has a compact support in $\Omega$. Using some triangular meshes, numerical results show an order 2 error on $u$, an order 1 error on $\nabla u$ and an order 1/2 on $\Delta u$.

**References**

[1] M. Ben-Artzi, I. Chorev, J.-P. Croisille, D. Fishelov, A compact difference scheme for the biharmonic equation in planar irregular domains, SIAM J. Numer. Anal. 47 (4) (2009) 3087-3108.

[2] M. Ben-Artzi, J.-P. Croisille, D. Fishelov, A fast direct solver for the biharmonic problem in a rectangular grid, SIAM J. Sci. Comput. 31 (1) (2008) 303-333.

[3] C. Bi, L. Li, Mortar finite volume method with Adini element for biharmonic problem, J. Comput. Math. 22 (3) (2004) 475-488.

[4] S.C. Brenner, L.-Y. Sung, $C^0$ interior penalty methods for fourth order elliptic boundary value problems on polygonal domains, J. Sci. Comput. 22/23 (2005) 83-118.

[5] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991.

[6] G. Chen, Z. Li, P. Lin, A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow, Adv. Comput. Math. 29 (2) (2008) 113-133.

[7] P.G. Ciarlet, The finite element method, in: P.G. Ciarlet, J.-L. Lions (Eds.), Handbook of Numerical Analysis, vol. III, Part I, North-Holland, Amsterdam, 1991.

[8] P. Destuynder, M. Salaun, Mathematical Analysis of Thin Plate Models, Mathématiques & Applications [Mathematics & Applications], vol. 24, Springer-Verlag, Berlin, 1996.

[9] R. Eymard, T. Gallouët, Herbin, Finite volume schemes for the biharmonic problem on general meshes, 2010, submitted for publication.

[10] E.H. Georgoulis, P. Houston, Discontinuous Galerkin methods for the biharmonic problem, IMA J. Numer. Anal. 29 (3) (2009) 573-594.

[11] T. Gudi, N. Nataraj, A.K. Pani, Mixed discontinuous Galerkin finite element method for the biharmonic equation, J. Sci. Comput. 37 (2) (2008) 139-161.

[12] I. Mozolevski, E. Süli, P.R. Bösing, hp-version a priori error analysis of interior penalty discontinuous Galerkin finite element approximations to the biharmonic equation, J. Sci. Comput. 30 (3) (2007) 465-491.

[13] E. Süli, I. Mozolevski, hp-version interior penalty DGFEMs for the biharmonic equation, Comput. Methods Appl. Mech. Engrg. 196 (13-16) (2007) 1851-1863.

[14] T. Wang, A mixed finite volume element method based on rectangular mesh for biharmonic equations, J. Comput. Appl. Math. 172 (1) (2004) 117-130.