\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2017 (2017), No. 305, pp. 1--15.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2017 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2017/305\hfil Approximation of the leading singular coefficient] {Approximation of the leading singular coefficient of an elliptic fourth-order equation} \author[M. Abdelwahed, N. Chorfi, V. D. R\u{a}dulescu\hfil EJDE-2017/305\hfilneg] {Mohamed Abdelwahed, Nejmeddine Chorfi, Vicen\c{t}iu D. R\u{a}dulescu} \address{Mohamed Abdelwahed \newline Department of Mathematics, College of Sciences, King Saud University, Riyadh, Saudi Arabia} \email{mabdelwahed@ksu.edu.sa} \address{Nejmeddine Chorfi (corresponding author)\newline Department of Mathematics, College of Sciences, King Saud University, Riyadh, Saudi Arabia} \email{nchorfi@ksu.edu.sa} \address{Vicen\c{t}iu D. R\u{a}dulescu \newline AGH University of Science and Technology, Faculty of Applied Mathematics, al. Mickiewicza 30, 30-059 Krakow, Poland. \newline Department of Mathematics, University of Craiova, 200585 Craiova, Romania} \email{vicentiu.radulescu@imar.ro} \dedicatory{Communicated by Giovanni Molica Bisci} \thanks{Submitted October 11, 2017. Published December 14, 2017.} \subjclass[2010]{35J15, 78M22} \keywords{Bilaplacian equation; singularity coefficient; \hfill\break\indent dual singular method; mortar spectral element method} \begin{abstract} The solution of the biharmonic equation with an homogeneous boundary conditions is decomposed into a regular part and a singular one. The later is written as a coefficient multiplied by the first singular function associated to the bilaplacian operator. In this paper, we consider the dual singular method for finding the value of the leading singular coefficient, and we use the mortar domain decomposition technique with the spectral discretization for its approximation. The numerical analysis leads to optimal error estimates. We present some numerical results which are in perfect coherence with the analysis developed in this paper. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{remark}[theorem]{Remark} \newtheorem{example}[theorem]{Example} \allowdisplaybreaks \section{Introduction} In a polygonal domain and when the data are smooth, the solution of an elliptic differential equation is irregular. For a homogeneous problem of the bilaplacian operator, we define some singular functions contingent to the geometry of the domain. The solution is the sum of two components: a regular part and singular functions (\cite{D,G,G1,K}). The later are multiplied by appropriate coefficients called singular coefficients. To approximate the leading singularity coefficients two algorithms are deployed. We first refer to the Strang and Fix algorithm \cite{SF}, permitting to add the leading singularity function to the discrete space \cite{C1}. Secondly, the singular dual method algorithm \cite{AM,ABM}. In physics and particularly in solid mechanics (crack propagation), the leading singularity coefficient is very significant. The calculation of this coefficient was obtained by use of finite elements (Amara and Moussaoui \cite{AM,AM1}). In this paper, we propose to use the mortar spectral element method combined with the method based on the singular dual function. We decompose the domain in an union of finite number of disjoint rectangles. On each rectangle, the discrete functions are polynomials of high degree. We enforce the discrete solution to satisfy a matching condition on the interfaces. Due to the non continuity of the discrete functions, mortar spectral element technique is nonconforming. For more details on the mortar spectral element method, we relate to Bernardi et al \cite{BM1,BM2,BCM,BMP}. In this work, we prove that the order of the error estimation between the continuous leading singularity coefficient end the discrete one is optimal. This order is better than that obtained by the Strang and Fix algorithm. The paper is outlined as follows. In the section 2, we present the geometry of the domain, the continuous problem and the dual singular method which allows us to calculate the leading coefficient of the singularity. This later depends only on the solution. The approximation of the leading singularity coefficient by the mortar element spectral method and the optimality estimation of the error is described in Section 3. Finally, the results of a numerical test are given in Section 4. \section{Geometry of the domain and the continuous problem} We denote by $\Omega$ a polygonal domain in ${\mathbb{R}}^2$ such that there exists a finite number of open rectangles $\Omega_i,1 \leq i \leq I$, satisfying \begin{equation} \bar\Omega= \cup_{i=1}^{I}\bar \Omega_i,\quad \Omega_i\cap \Omega_l=\emptyset \quad \text{for } i\neq l.\label{e2.1} \end{equation} and such that the intersection of each $\bar \Omega_i,1 \leq i \leq K$, with the boundary $\partial \Omega$ is either empty or a corner or one or several entire edges of sub-domain $\Omega_i$. We choose the coordinate axes parallel to the edge of the $\Omega_i$. We denote by $\Gamma^{i,j}$, $1\leq j\leq 4$ the edges of $\Omega^i$ and by $\gamma^{i,l}$, $1\leq i\neq l\leq I$, the open segment such that $$ {\bar{\gamma}}^{i,l}= \partial \Omega^i\cap\partial \Omega^l. $$ The set $\mathcal{V}$ will be the set of all vertices of the $\Omega^i$, $1\leq i\leq I$ and $\mathcal{S}=\cup_{i=1}^I\cup_{j=1}^4{\bar{\Gamma}}^{i,j}$ the skeleton of the decomposition. We choose finite set of disjoint open segments $\gamma^k$, where $k$ belongs to a finite set $K$ such that $\mathcal{S}=\cup_{k\in K}{\bar{\gamma}}^{k}$, each $\gamma^k$, $k\in K$ is called mortar and its being a edge $\Gamma^{i(k),j(k)}$. We are interested to non-convex domains, we assume that there exists an angle equal either to $3\pi\over 2$ or to $2\pi$ (case of the crack). Handling the singular function is local process, so that there is no restriction to suppose that the non-convex corner is unique. \subsection{Notation} Let $\omega$ be the value of the non-convex angle equal either to $3\pi\over 2$ or to $2\pi$, \textbf{a} be the corresponding corner of $\Omega$ and $\Delta$ be the open domain in $\Omega$ such that $\bar\Delta$ is the union of the $\bar\Omega_i$ which contain \textbf{a}. We choose the origin of the coordinate axes at the point $\mathbf{a}$. We introduce a system of polar coordinates $(r,\theta)$ where $r$ stands for the distance from \textbf{a} and $\theta$ is such that the line $\theta=0$ contains an edge of $\partial\Omega$. For reasons which will appear later, we are lead to make the following conformity assumption: We suppose that the decomposition of the domain $\Delta$ is conforming (see figure \ref{fig1}). If $\mathbf{a}$ is a vertex of the mortar $\Gamma^{i(k),j(k)}$ which coincides with $\Gamma^l$ a side of a sub-domain $\Omega^l$, $l \ne i(k)$ then $N_{i(k)}\le N_l$, such that the restriction of a function to $\Delta$ is in $H^2(\Delta)$. \begin{figure}[ht] \centering \includegraphics[angle=270,width=0.7\hsize]{fig1} % dom.pdf \caption{Domain $\Omega$} \label{fig1} \end{figure} The following biharmonic problem with homogenous boundary conditions is the model under consideration, \begin{equation}\label{1} \begin{gathered} \Delta^2u=f \quad \text{in } \Omega\\ u=0 \quad\text{on } \partial\Omega\\ \frac{\partial u}{\partial n}=0 \quad\text{on } \partial\Omega. \end{gathered} \end{equation} This problem admits the following variational formulation: Find $u \in H_0^2(\Omega)$ such that \begin{equation}\label{f1} \forall v \in H_0^2(\Omega) \quad \int_{\Omega} \Delta u.\Delta v \,dx\,dy =\int_{\Omega} f\,v\, dx\,dy \end{equation} where $f$ is in $L^2(\Omega)$. Formulation \eqref{f1} can be written in the equivalent form: Find $u \in H_0^2(\Omega)$ such that \begin{equation}\label{f2} \forall v \in H_0^2(\Omega) \quad \sum_{i=1}^I\int_{\Omega_i} (\Delta u|_{\Omega_i})\,(\Delta v|_{\Omega_i}) \,dx\,dy =\sum_{i=1}^I\int_{\Omega_i} f|_{\Omega_i}\,v|_{\Omega_i}\, dx\,dy. \end{equation} The function $v$ is in $H_0^2(\Omega)$ if and only if $v|_{\Omega_i}\in H^2(\Omega_i)$, $1\le i \le I$, $v$ and $\partial_n v$ vanishes on $\partial \Omega$, $v|_{\Omega_i}=v|_{\Omega_l}$ and $\partial_n v|_{\Omega_i}=\partial_n v|_{\Omega_l}$ on $\gamma^{i,l}$, $1\le i\ne l \le I$. Using Lax-Milgram theorem we show that the problem is well posed. However, in a polygonal domain the global regularity of the solution depends on the angle $\omega$. For a non negative real $s$, if $f$ in $H^{s-2}(\Omega)$ then the solution $u$ of problem \eqref{1} belongs to $H^{s+2}(\Omega)$ and there exists a positive constant $C$ such that $$ \|u\|_{H^{s+2}(\Omega)} \le C \|f\|_{H^{s-2}(\Omega)}. $$ If $\omega=\frac{3\pi}{2}$, $s \le 0.544844$, and if $\omega=2\pi$, $s \le 0.5$. To enhance the regularity we decompose the solution as follows: \begin{equation}\label{decom} u=u_R+\mu\, \tau_1 \quad \text{such that } u_R\in H^{s+2}(\Omega) \end{equation} where $\mu$ is the leading singular coefficient, $\tau_1$ is the first singular function and there exists a positive constant $C$ such that $$ \|u_R\|_{H^{s+2}(\Omega)}+|\mu|\le C \|f\|_{H^{s-2}(\Omega)}. $$ $\bullet$ When $\omega=\frac{3\pi}{2}$: $s< 1.544$ and $\tau_1(r,\theta)=\chi(r,\theta) r^{1.544}\psi(\theta)$, \begin{equation}\label{f3} \begin{aligned} \psi(\theta)=&4.302\big(\cos(0.092 \theta)-\cos(1.908 \theta)\big) \\ &-1.1815\big(10.869\sin(0.092 \theta)-0.524 \sin(1.908\theta)\big). \end{aligned} \end{equation} $\bullet$ In the case where $\omega={2\pi}$: $s< 1.5$ and $\tau_1(r,\theta)=\chi(r,\theta) r^{1.5}\psi(\theta)$, \begin{equation}\label{f4} \psi(\theta)=\big(\sin(1.5\theta)-3\sin(0.5\theta) + \cos(1.5\theta)-\cos(0.5\theta)\big). \end{equation} where $\chi$ is the $C^{\infty}$ cut off function with support in $\overline{\Delta}$ which is equal to $1$ in the neighborhood of $\mathbf{a}$. To compute the leading singularity coefficient $\mu$, we use the dual singularity functions \cite{AM,F}. We consider the characteristic equation of the bilaplacian \begin{equation}\label{bilap} F(z,\omega)=\sin^2(\omega z)-z^2 \sin(\omega ^2)=0. \end{equation} The function $F$ is even with respect to $z$, then if $z_0=\xi+i\eta$ is the solution of \eqref{bilap} then it is the same for $-z_0=-\xi-i\eta$. The singular functions depend on the positive parameter $\xi=\operatorname{Re}(z)$ where $z$ is solution of the equation \eqref{bilap}. These singular functions are solutions of the following homogeneous problem: \begin{equation}\label{pb} \begin{gathered} \Delta^2 \tau_{\xi}=0 \quad \text{in } \Omega\\ \tau_{\xi}=0 \quad \text{on } \Gamma\\ \frac{\partial \tau_{\xi}}{\partial n}=0 \quad\text{on } \Gamma. \end{gathered} \end{equation} To each $\tau_{\xi}$ corresponds a function $\tau_{-\xi}$ solution of problem \eqref{pb}. \begin{remark} \label{rmk2.1} \rm We remark that the function $\tau_{\xi}$ is at least in $H^2(\Omega)$ while $\tau_{-\xi}$ is in $L^2(\Omega)$ in the neighborhood of $\mathbf{a}$. Let $\tau_1^*(r,\theta)=\chi(r,\theta)\tau_{-\xi}$ be the dual singular function and $\tau_1^*(r,\theta)=r^{1-\xi}\psi(\theta)$ in the neighborhood of $\mathbf{a}$, where $\psi$ is defined in \eqref{f3} and \eqref{f4}. \end{remark} Since $\Delta^2 (\tau_1^*)\in H^{-2}(\Omega)$, let $\varphi^*$ the solution of the variational problem: find $\varphi^*\in H^2_0(\Omega)$ such that \begin{equation}\label{dualp} \int_{\Omega}\Delta \varphi^* \Delta v \, dx =\langle \Delta^2(\tau_1^*),v\rangle ,\quad \forall v \in H^2_0(\Omega). \end{equation} Then the singularity coefficient is \begin{equation}\label{sing} \mu=c\int_{\Omega} f (\tau_1^*-\varphi^*)dx\,dy. \end{equation} where $c$ is given by \cite{D,F} \begin{equation}\label{defc} c=8\xi(\omega)(\xi(\omega)+1)\int_0^\omega \exp(-(z+i)\theta)\psi(\theta)\, d\theta \end{equation} where $\xi(\omega)=\sup\big\{\operatorname{Re}(z): z \text{ solution of \eqref{bilap}},\; z\ne \pm1 \Big\}$, $\psi$ is defined in \eqref{f3} and \eqref{f4}. \section{Discrete problem} The discretization is based on the Galerkin method. The goal is to construct the discrete space which is not included in the continuous one because the method is not conform. Since our problem is of order $4$, we have to deal with two boundary conditions. Then we need two mortar functions; one for the trace and the other for the normal derivative. Let $\delta=(N_i)_{1\le i \le I}$ the discretization parameter where $N_i$ are the degrees of the approximation polynomials in each sub-domain $\Omega_i$, $1\le i \le I$ ($\mathbb{P}_N(\Omega)$ is the space of polynomial functions of degree less or equal to $N$). We introduce the space $\mathcal M_\delta$ called the space of mortar functions made of couples $(\varphi_0,\varphi_1)$ such that their restriction to each $\Gamma^{i(k),j(k)}$ is a polynomial of degree less than $N_i(k)$ and such that the following properties hold: at each vertex $e$ of a sub-domain, each couple $(\varphi_0,\varphi_1)$ defines a unique value $\varphi^e$, a unique derivative $\varphi_x^e$ with respect to $x$, a unique derivative $\varphi_y^e$ with respect to $y$ and a unique mixed $\varphi_{xy}^e$ with respect to $x$ and $y$. We define the discrete space $X_\delta$ as the space of functions $v_\delta$ such that: \begin{itemize} \item[(i)] for each $i$, $1\le i \le I$, the restriction of $v_{\delta}$ to $\Omega_i$ belongs to $\mathbb P_{N_i}(\Omega_i)$; \item[(ii)] $v_\delta$ and its normal derivative vanish on $\partial \Omega$; \item[(iii)] there exists a mortar couple $(\varphi_0,\varphi_1)$ in $\mathcal M_\delta$ such that for $1\le i \le I$, \begin{gather*} \forall e \in \mathcal{V}, \quad v_\delta|_{\Omega_i}(e)=\varphi^e, \quad \frac{\partial}{\partial x}( v_\delta|_{\Omega_i})(e)=\varphi_x^e,\\ \frac{\partial}{\partial y}( v_\delta|_{\Omega_i})(e)=\varphi_y^e , \quad \frac{\partial^2}{\partial x \partial y}( v_\delta|_{\Omega_i})(e)=\varphi_{xy}^e, \end{gather*} \begin{gather}\label{mat1} \forall \psi \in \mathbb P_{N_i-4}(\Gamma^{ij}), \quad \int_{\Gamma^{ij}}(v_{\delta}|_{\Omega_i}-\varphi_0)(\eta)\psi(\eta)\, d \eta=0,\\ \label{mat2} \int_{\Gamma^{ij}}(\frac{\partial v_{\delta}}{\partial n}|_{\Omega_i} -\varphi_1)(\eta)\psi(\eta)\, d \eta=0, \quad 1\le j \le 4. \end{gather} \end{itemize} The choice of $\mathbb P_{N_i-4}(\Gamma^{ij})$ is justified by the fact that four conditions are enforced on the vertices of $\Gamma^{ij}$, one on the function and one on its normal derivative at each vertex. In the case of the problem of order four and to take into account the boundary conditions, it is more appropriate to use a quadrature formula which uses the function values on the extremities as well as the values of its normal derivative. The following lemma defines this quadrature formula (see \cite{BM3,SS} for the proof) \begin{lemma} \label{lem3.1} Let $N\ge 2$ an integer, there exists a unique set of points $\xi_j$, $1\le j \le N-1$, a unique set of positif reals $\rho_j$, $1\le j \le N-1$, $\rho_+$, $\rho_-$ such that for all polynomials $\varphi$ in $\mathbb P_{2N-1}(]-1,1[)$ \begin{equation}\label{lem41} \int_{-1}^1 \varphi(x)\, dx = \sum_{j=1}^{N-1}\varphi(\xi_j)\rho_j +\varphi(-1)\rho_- + \varphi(1)\rho_+ \end{equation} \end{lemma} \begin{remark} \label{rmk3.2} \rm The nodes $\xi_j$; $1\le j \le N-1$, are the zeros of the derivative of the Legendre polynomial $L_N$. We refer to \cite{BM3} for the calculus of $\xi_j$ and $\rho_j$, $1\le j \le N-1$. \end{remark} Given two functions $u$, $v$ continuous on $\overline \Omega =[-1,1]\times[-1,1]$ and vanishes on its boundary, we define the following discrete scalar product $$ (u,v)_N=\sum_{j=1}^{N-1}\sum_{l=1}^{N-1}u(\xi_j,\xi_l) v(\xi_j,\xi_l) \rho_j \rho_l. $$ If $T^i$ is the bijection from $]-1,1[^2$ in $\Omega_i$, we define $$ (u,v)_{N_i}=\frac{|\Omega_i|}{4}\sum_{j=1}^{N_i-1} \sum_{l=1}^{N_i-1}(u\circ T^i)(\xi_j,\xi_l) (v\circ T^i)(\xi_j,\xi_l) \rho_j \rho_l. $$ Hence, for each value of $\delta$, the discrete problem is written: Find $u_{\delta} $ in $ X_{\delta}$ such that for all $v_{\delta} \in X_{\delta}$, \begin{equation}\label{f5} \sum_{i=1}^{I}(\Delta u_{\delta}|_{\Omega_i},\Delta v_{\delta}|_{\Omega_i})_{N_i} =\sum_{i=1}^{I}(f,v_{\delta}|_{\Omega_i})_{N_i}. \end{equation} See \cite{Z2} for the numerical analysis and the implementation of problem \eqref{f5} using the mortar spectral element method. We present in this work a method based on dual singular function which will allow us to approximate the leading singularity coefficient of the bilaplacian operator. This method has high precision compared to the Strang and Fix algorithm \cite{SF}. Let $X_{\delta}^* = X_{\delta}+\mathbb R \tau_1$ be the augmented discrete space, which is a Banach space by the following discrete norm, for all $u_{\delta}^*=u_{\delta}+\mu\tau_1 \in X_{\delta}^*$, $$ \|u_{\delta}^*\|_{1*} =\sum_{i=1}^I \Big(\|u_{\delta}/_{\Omega_i} \|_{H^2(\Omega_i)}^2 +|\mu|^2 \| \tau_1/_{\Omega_i}\|_{H^2(\Omega_i)}^2 \Big)^{1/2}. $$ We ask the following two discrete problems: (1) find the function $u_{\delta}^*=u_{\delta}+\mu\tau_1 \in X_{\delta}^*$ such that for all $v_{\delta}^*=v_{\delta}+\xi\tau_1 \in X_{\delta}^*$ we have \begin{equation}\label{f6} a_{\delta}^*(u_{\delta}^*,v_{\delta}^*) =\sum_{i=1}^{I}\int_{\Omega_i}f|_{\Omega_i} v_{\delta}^*|_{\Omega_i}\, dx\,dy; \end{equation} (2) find $\varphi_{\delta}^*$ in $X_{\delta}^*$ such that for all $\psi_{\delta}^*\in X_{\delta}^*$, \begin{equation}\label{f7} a_{\delta}^*(\varphi_{\delta}^*,\psi_{\delta}^*) =\sum_{i=1}^{I}\int_{\Omega_i}\Delta^2 \tau_1^*|_{\Omega_i} \psi_{\delta}^*|_{\Omega_i}\, dx\,dy. \end{equation} The bilinear form $a_{\delta}^*$ is defined by \begin{equation} \begin{aligned} a_{\delta}^*(u_{\delta}^*,v_{\delta}^*) &=\sum_{i=1}^{I}\Big((\Delta u_\delta|_{\Omega_i}, \Delta v_\delta|_{\Omega_i})_{N_i}+\mu \int_{\Omega_i} \Delta \tau_1|_{\Omega_i} \Delta v_\delta|_{\Omega_i}\, dx\,dy \\ &\quad+ \xi \int_{\Omega_i} \Delta \tau_1|_{\Omega_i} \Delta u_\delta|_{\Omega_i}\, dx\,dy + \mu \xi \int_{\Omega_i} (\Delta \tau_1|_{\Omega_i})^2 \, dx\,dy\Big). \end{aligned} \end{equation} We refer to \cite{ACR} for the numerical analysis of this problem and to \cite{ACOR} for its implementation. The following proposition gives us the expression of the discrete leading singularity coefficient. \begin{proposition} \label{prop3.3} Let $u$, $\varphi^*$, $u_\delta^*$ and $\varphi_{\delta}^*$ be respectively the solutions of \eqref{1}, \eqref{dualp}, \eqref{f6} and \eqref{f7}, we have (i) \begin{equation}\label{f7b} (\frac{1}{c})\mu_\delta=\int_{\Omega} f \tau_1^*\, dx\,dy + \int_{\Omega} u_\delta^* \Delta^2 \tau_1^*\, dx\,dy = \int_{\Omega} f (\tau_1^*-\varphi_\delta^*)\, dx\,dy \end{equation} (ii) \begin{equation}\label{f8} \begin{aligned} &(\frac{1}{c})(\mu-\mu_\delta) \\ &=\sum_{i=1}^I\int_{\Omega_i} \Delta(u-u_\delta^*)|_{\Omega_i} \Delta(\varphi^*-\varphi_\delta^*)|_{\Omega_i}\, dx\,dy\\ &\quad +\sum_{1\le i \ne l \le I}\int_{\Gamma^{il}} \Big[\frac{\partial (\Delta u)}{\partial n_i} (\varphi_{\delta}^*|_{\Omega_i}-\varphi_{\delta}^*|_{\Omega_l}) - \frac{\partial (\Delta \varphi^*)}{\partial n_i} (u_{\delta}^*|_{\Omega_i} -u_{\delta}^*|_{\Omega_l}) \Big] \end{aligned} \end{equation} where $c$ is defined in \eqref{defc} \end{proposition} \begin{proof} We consider $D$ the intersection of the domain $\Omega$ and the ball of center $\mathbf{a}$ and radius $R$. We consider that the cut-off function $\chi$ is equal $1$ in $D$ then we choose $R$ such that $\Delta^2\tau_1=\Delta^2\tau_1^*=0$ then from \eqref{1} and \eqref{decom} we have $$ \int_{\Omega} f \tau_1^* \, dx\,dy = \int_{D} -\Delta^2 u_{\delta} \tau_1^* \, dx\,dy + \int_{\Omega\backslash D}-\Delta^2 u_{\delta}^* \tau_1^* \, dx\,dy. $$ By double integration by parts we conclude \begin{align*} &\int_{\Omega} f \tau_1^* \, dx\,dy + \int_{D} u_\delta^* \Delta^2 \tau_1^* \, dx\,dy \\ &=\int_0^\omega \big(\partial_r(\Delta(u_\delta^*-u_\delta))\tau_1^* -(u_\delta^*-u_\delta)\partial_r (\Delta \tau_1^*)\big)(R,\theta) r\, d\theta \end{align*} Since $u_\delta^*-u_\delta=\mu_\delta \tau_1$, we have $$ c\Big(\int_{\Omega}f \tau_1^*\, dx\,dy + \int_\Omega u_\delta^* \Delta^2 \tau_1^*\, dx\,dy\Big)=\mu_\delta. $$ To obtain the second equality we replace $v$ by $u$ in problem \eqref{dualp}. To show \eqref{f8} we use \eqref{sing} and \eqref{f7}. We obtain \begin{equation} (\frac{1}{c})(\mu-\mu_\delta) =\int_{\Omega}f (\varphi_\delta^*-\varphi^*)\, dx\,dy = \sum_{i=1}^I \int_{\Omega_i}\Delta^2 u | _{\Omega_i}(\varphi^* -\varphi_{\delta}^*)|_{\Omega_i}\, dx\,dy. \end{equation} By double integration by parts we have \begin{equation} \begin{aligned} &(\frac{1}{c})(\mu-\mu_\delta) \\ &=\sum_{i=1}^I \int_{\Omega_i} \Delta u \Delta(\varphi^*-\varphi_\delta^*)\, dx\,dy + \sum_{1\le i 0$. The error between $\mu$ and $\mu_{\delta}$ satisfies the following estimate, for $\varepsilon>0$, $$ |\mu-\mu_\delta|\le C N^{-2} \Big(\sum_{1\le i \le I} N_i^{-\sigma_i}\Big) \| f \|_{H^{s-2}(\Omega)}, $$ $N=\inf_{1\le i \le I} N_i$ and $$ \sigma_i=\begin{cases} s-2 &\text{if $\Omega_i$ does not contain any vertices of $\Omega$},\\ \inf (s-2, 2 \eta_1 (\frac{\pi}{2}) - \varepsilon) &\text{if $\overline{\Omega_i}$ contains one vertex of $\Omega $ other than } \mathbf{a},\\ \inf (s-2, 2 \eta_1 (\omega) - \varepsilon) &\text{if } \overline{\Omega_i} \text{ contains } \mathbf{a}, \end{cases} $$ where $\eta_1(\omega)$ is the second real solution of equation \eqref{bilap} in the band $0 < \operatorname{Re}(z)