Discretization error estimates for penalty formulations of a linearized Canham-Helfrich type energy

TL;DR

This paper develops error estimates for penalty-based finite element methods solving a linearized Canham-Helfrich energy problem, crucial for modeling biomembranes with embedded particles, without needing to resolve particle boundaries.

Contribution

It provides the first discretization error estimates for penalty formulations of a fourth-order membrane energy problem with reduced regularity solutions.

Findings

Proves almost-H^{5/2} regularity of solutions.

Derives optimal error estimates for penalty finite element discretizations.

Numerical results confirm theoretical error bounds.

Abstract

This paper is concerned with minimization of a fourth-order linearized Canham-Helfrich energy subject to Dirichlet boundary conditions on curves inside the domain. Such problems arise in the modeling of the mechanical interaction of biomembranes with embedded particles. There, the curve conditions result from the imposed particle--membrane coupling. We prove almost-$H^{\frac{5}{2}}$ regularity of the solution and then consider two possible penalty formulations. For the combination of these penalty formulations with a Bogner-Fox-Schmit finite element discretization we prove discretization error estimates which are optimal in view of the solution's reduced regularity. The error estimates are based on a general estimate for linear penalty problems in Hilbert spaces. Finally, we illustrate the theoretical results by numerical computations. An important feature of the presented discretization is that it does not require to resolve the particle boundary. This is crucial in order to avoid re-meshing if the presented problem arises as subproblem in a model where particles are allowed to move or rotate.

Equations

$${\mathcal{J}}_{CHS}(\mathcal{M})={\int }_{\mathcal{M}}\frac{1}{2}\kappa H^2+{\kappa }_{G}K+\sigma \mathrm{d}{\mathcal{H}}^2\text{.}$$

$${\mathcal{J}}_{\Omega }(u)={\int }_{\Omega }\frac{1}{2}\kappa (\Delta u{)}^2+\frac{1}{2}\sigma |\nabla u{|}^2\mathrm{d}x$$

$$\Gamma =\bigcup _{i=1}^k\partial B_i=\partial {\Omega }_B\setminus \partial \Omega =\partial {\Omega }_B\cap B.$$

$$H=\{v\in H^2({\Omega }_B)|v{|}_{\partial \Omega }={\partial }_{\nu }v{|}_{\partial \Omega }=0\}$$

$$u{|}_{\partial B_i}$$

$$u{|}_{\Gamma }=f_1+\sum _{i=1}^k{\gamma }_i{\eta }_i,{\partial }_{\nu }u{|}_{\Gamma }=f_2\text{ .}$$

$$T_{\Gamma }:H^2(\Omega )$$

$$P_{\Gamma }:H^{\frac{3}{2}}(\Gamma )\times H^{\frac{1}{2}}(\Gamma )$$

$$P_{\Gamma }^1(v)=v-\sum _{i=1}^k|\partial B_i{|}^{-1}{\int }_{\partial B_i}v\mathrm{d}\sigma =v-\sum _{i=1}^k\frac{(v,{\eta }_i{)}_{L^2(\Gamma )}}{({\eta }_i,{\eta }_i{)}_{L^2(\Gamma )}}{\eta }_i,P_{\Gamma }^2(v)=v.$$

$$V_f:=\left\{v\in H_0^2(\Omega )\mid P_{\Gamma }(T_{\Gamma }v-f)=0\right\}$$

$${\mathcal{J}}_{\Omega }(u_{\epsilon })+\sum _{i=1}^2\frac{1}{{\epsilon }_i}\Vert P_{\Gamma }^i(T_{\Gamma }^i{u}_{\epsilon }-f_i){\Vert }_{L^2(\Gamma )}^2\text{.}$$

$$T_B:H^2(\Omega )$$

$$P_B:H^2(B)$$

$$P_B{T}_{B}u=P_{B}u{|}_B\text{.}$$

$${\mathcal{J}}_{\Omega }(u_{\epsilon })+\frac{1}{\epsilon }\Vert P_B(T_B{u}_{\epsilon }-u{|}_B){\Vert }_{H^s(B)}^2\text{.}$$

$$\frac{1}{2}a(u,u)-\ell (u)\text{.}$$

$$a(u,v)=\ell (v)\forall v\in U_0.$$

$$\frac{1}{2}a(u_{\epsilon }^X,u_{\epsilon }^X)-\ell (u_{\epsilon }^X)+\sum _{i=1}^m\frac{1}{2{\epsilon }_i}\Vert u_{\epsilon }^X-u{\Vert }_{b_i}^2.$$

$$a_{\epsilon }(\cdot ,\cdot )=a(\cdot ,\cdot )+\sum _{i=1}^m\frac{1}{{\epsilon }_i}{b}_i(\cdot ,\cdot ),{\ell }_{\epsilon }(\cdot )=\ell (\cdot )+\sum _{i=1}^m\frac{1}{{\epsilon }_i}{b}_i(u,\cdot )$$

$$a_{\epsilon }(u_{\epsilon }^X,v)={\ell }_{\epsilon }(v)\forall v\in X.$$

$$|a(u,v)-\ell (v)|\le \sum _{i=1}^m{c}_i\Vert v{\Vert }_{b_i}\forall v\in X.$$

$$\Vert u-u_{\epsilon }^X{\Vert }_{a_{\epsilon }}^2\le 3\underset{v\in X}{\inf }\left(\Vert u-v{\Vert }_{a_{\epsilon }}^2+\sum _{i=1}^m{\epsilon }_i{c}_i^2\right)\text{.}$$

$$\Vert e{\Vert }_{a_{\epsilon }}^2=a_{\epsilon }(e,e)=a_{\epsilon }(e,u-v)+a_{\epsilon }(e,v-u_{\epsilon }^X)$$

$$\begin{array}{rl}{a}_{\epsilon }(e,u-v)& \le \Vert e{\Vert }_{a_{\epsilon }}\Vert u-v{\Vert }_{a_{\epsilon }}\le \frac{1}{4}\Vert e{\Vert }_{a_{\epsilon }}^2+\Vert u-v{\Vert }_{a_{\epsilon }}^2.\end{array}$$

$$a_{\epsilon }(u,w)-a(u,w)$$

$$a_{\epsilon }(e,w)=a(u,w)-\ell (w)\forall w\in X.$$

$$\displaystyle\begin{aligned} a_{\varepsilon}(e,v-u^{X}_{\varepsilon})&=a(u,v-u^{X}_{\varepsilon})-\ell(v-u^{X}_{\varepsilon})\\ &\leq\sum_{i=1}^{m}c_{i}\|v-u^{X}_{\varepsilon}\|_{b_{i}}\\ &\leq\sum_{i=1}^{m}c_{i}\left(\|e\|_{b_{i}}+\|u-v\|_{b_{i}}\right)\\ &\leq\sum_{i=1}^{m}\Bigl{(}\frac{1}{4\varepsilon_{i}}\|e\|_{b_{i}}^{2}+\frac{1}{2\varepsilon_{i}}\|u-v\|_{b_{i}}^{2}+\frac{3}{2}\varepsilon_{i}c_{i}^{2}\Bigr{)}\\ &\leq\frac{1}{4}\|e\|_{a_{\varepsilon}}^{2}+\frac{1}{2}\|u-v\|_{a_{\varepsilon}}^{2}+\sum_{i=1}^{m}\frac{3}{2}\varepsilon_{i}c_{i}^{2}\end{aligned}$$

$$\Vert e{\Vert }_{a_{\epsilon }}^2\le \frac{1}{2}\Vert e{\Vert }_{a_{\epsilon }}^2+\frac{3}{2}\Vert u-v{\Vert }_{a_{\epsilon }}^2+\sum _{i=1}^m\frac{3}{2}{\epsilon }_i{c}_i^2$$

$${\tilde{a}}_{\epsilon }(\cdot ,\cdot )=\tilde{a}(\cdot ,\cdot )+\sum _{i=1}^m\frac{1}{{\epsilon }_i}{\tilde{b}}_i(\cdot ,\cdot ),{\tilde{\ell }}_{\epsilon }(\cdot )=\tilde{\ell }(\cdot )+\sum _{i=1}^m\frac{1}{{\epsilon }_i}{\tilde{b}}_i(u,\cdot )$$

$${\tilde{a}}_{\epsilon }({\tilde{u}}_{\epsilon }^X,v)={\tilde{\ell }}_{\epsilon }(v)\forall v\in X\text{.}$$

Full text

1. 1 Introduction

2. 2 Notation and Problem Formulations

3. 3 An Error Estimate for Linear Penalty Problems on Hilbert Spaces

4. 4 Regularity of the Hard Curve Constraint Problem

5. 5 Discretization Errors

6. 6 Numerical Example

7. 6.1 A symmetric example problem with known exact solution

8. 6.2 Discretization, quadrature, and error measurement

9. 6.3 Results for the soft curve formulation

10. 6.4 Results for the soft bulk formulation

11. 6.5 A non-symmetric example problem

12. 7 Acknowledgements

13. 8 Appendix

1. Introduction

A standard model for the behavior of biomembranes on a macroscale is the Canham–Helfrich model which describes a biological membrane mathematically as a hypersurface $\mathcal{M}\subseteq\mathbb{R}^{3}$ minimizing the Canham–Helfrich energy

[TABLE]

Here $H$ and $K$ are the mean and Gaussian curvature of $\mathcal{M}$, the coefficients $\kappa>0$ and $\kappa_{G}\geq 0$ are the corresponding bending rigidities, and $\sigma\geq 0$ is the membrane’s surface tension. Under the assumption that the membrane is “almost flat” one can justify a geometric linearization of this functional for a membrane patch, leading to the so called Monge–gauge approximation

[TABLE]

of the Canham–Helfrich energy. Here, the membrane patch is considered to be the graph $\mathcal{M}=\{(x,u(x)\mathrel{|}x\in\Omega\}$ of a function $u\colon\Omega\to\mathbb{R}$ over some reference domain $\Omega\subseteq\mathbb{R}^{2}$.

A variety of hybrid models for the coupling of embedded particles to the membrane have been considered. These hybrid models are based on a continuous surface description of the membrane while particles are described by discrete entities (see, e.g., [3, 8, 9, 12, 17]). For an overview on hybrid models we refer to [10] and the references cited therein. In the present paper we consider coupling conditions imposed on the particle boundaries and follow the notation introduced in [10]. There the reference domain $\Omega$ is split into the membrane’s domain $\Omega_{B}\subseteq\Omega$ and the particles’ domain $B=\Omega\setminus\Omega_{B}$. The model introduces membrane–particle interactions by functions which prescribe the membrane’s height profile and slope on the interface $\Gamma:=\partial\Omega_{B}\cap\partial B$. Altogether this yields an energy minimization problem subject to Dirichlet boundary conditions on $\Gamma$.

The present paper introduces and analyzes a discretization based on a penalty formulation for the corresponding boundary value problem that avoids the resolution of particle boundaries. The reason for considering penalized boundary conditions is the following: Since a biological membrane behaves like a fluid in tangential directions, particles can in principle move and rotate in plane. Any model simulating moving particles and any method computing optimal particle positions will thus have to solve multiple problems with varying particle positions. If the particle boundaries would have to be resolved this would require mesh-deformation or re-meshing which can be computationally quite expensive. An alternative is to replace the strict boundary conditions by adequate penalty terms that can be formulated without having to resolve the boundary. One such penalty approach has been introduced in [10] and is called soft curve formulation. A novel second formulation in this paper will be the soft bulk formulation. Both penalty problems will be defined later on.

Our goal is to derive discretization error estimates for those penalty problems. For second order equations such estimates are well known, see for example [2], [4] and [15]. It turns out that the fourth order problems that we are interested in can be treated sufficiently well using a simple general framework for penalty problems in Hilbert spaces. This is mostly due to the fact that the regularity of our solutions is rather limited in first place and so we can use rather simple estimates to still obtain optimal rates of convergence. Thanks to the abstract formulations we get as a byproduct a general error theory for finite element penalty problems with low regularity.

The paper is structured as follows: In Section 2 we introduce the notations and problem formulations that are used throughout this paper. Section 3 is devoted to an abstract error result for linear penalty approximations on Hilbert spaces. As a foundation for the application of this result to the problem at hand we then discuss regularity of solutions in Section 4. There we prove that a solution of the original problem lies in $H^{\frac{5}{2}-\delta}(\Omega)$ for all $\delta>0$. In Section 5 we combine the regularity with the developed abstract results to show the optimal (in view of the restricted regularity) convergence rate $O(h^{1/2-\delta})$ for a discretization with Bogner–Fox–Schmit finite elements. Finally, Section 6 illustrates our results by numerical examples that reproduce our theoretical findings.

2. Notation and Problem Formulations

We consider a membrane with $k$ embedded particles. To this end let $\Omega\subset\mathbb{R}^{2}$ be a bounded reference domain with Lipschitz boundary $\partial\Omega$ and $B_{i}\subset\Omega$, $i=1,\dots,k$ the area occupied by the $i$-th particle. For simplicity we assume that the $B_{i}$ are closed, nonempty, connected, and pairwise disjoint. Furthermore we require that each particle $B_{i}$ has a $C^{1,1}$-boundary. By $B=\bigcup_{i=1}^{k}B_{i}$ and $\Omega_{B}=\Omega\setminus B$ we denote the area occupied by the particles and the membrane, respectively. Since the $B_{i}$ are closed and disjoint the total membrane–particle interface is given by

[TABLE]

For simplicity we consider the space

[TABLE]

with homogeneous Dirichlet boundary conditions on the boundary of $\partial\Omega$. Notice that we can treat Navier and periodic boundary conditions on $\partial\Omega$ using the same techniques (cf. [10]).

We assume that the membrane–particle interaction is governed by boundary conditions on the interface. More precisely, each particle $B_{i}$ enforces a height profile given by $f^{i}_{1}\colon\partial B_{i}\to\mathbb{R}$ and a slope given by $f^{i}_{2}\colon\partial B_{i}\to\mathbb{R}$ to the membrane on its boundary $\partial B_{i}$. I. e. we consider boundary values

[TABLE]

where $\nu$ is the unit outward normal to $\partial\Omega_{B}$ and the parameter $\gamma^{i}\in\mathbb{R}$ is allowed to vary freely in order to factor out the average height on $\partial B_{i}$. This is necessary because we only want to prescribe the height profile, while the absolute or average height is not fixed. For the following we collect all such boundary data in a function $f=(f_{1},f_{2})=\sum_{i=1}^{k}f^{i}\colon\Gamma\to\mathbb{R}^{2}$ where $f^{i}=(f^{i}_{1},f^{i}_{2})\colon\partial B_{i}\to\mathbb{R}^{2}$ is extended by zero to the whole of $\Gamma$. Using this notation we consider the minimization problem:

Problem 2.1**.**

Find $u\in H$ and $\gamma\in\mathbb{R}^{k}$ minimizing $\mathcal{J}_{\Omega_{B}}(u)$ subject to

[TABLE]

Notice that the effect of the free parameter $\gamma_{i}$ is localized to $\partial B_{i}$ via the use of the indicator function $\eta_{i}:=\chi_{\partial B_{i}}\in L^{2}(\Gamma)$ which is one on $\partial B_{i}$ but vanishes on all $\partial B_{j},j\neq i$. In view of the trace theorem we from now on assume that $f=(f_{1},f_{2})\in H^{\frac{3}{2}}(\Gamma)\times H^{\frac{1}{2}}(\Gamma)$. Under this assumption it is known that Problem 2.1 admits a unique solution by application of Lax–Milgram’s theorem [10].

It is in fact possible to simplify this problem formulation. For this purpose we make use of the trace operator $T_{\Gamma}=(T_{\Gamma}^{1},T_{\Gamma}^{2})$

[TABLE]

This operator is well-defined, continuous, surjective and admits a continuous right-inverse (see e. g. [11, Theorem 1.5.1.2]). Note that we can also view $T_{\Gamma}$ as a trace operator for $H^{2}(\Omega_{B})$ due to our regularity assumptions on $\Gamma$. In order to simplify the notation for boundary values up to the average height in (2.2) we furthermore introduce the projection operator

[TABLE]

defined by $P_{\Gamma}(v_{1},v_{2})=(P_{\Gamma}^{1}(v_{1}),P_{\Gamma}^{2}(v_{2}))$ where

[TABLE]

Using this operator we can formulate the boundary conditions (2.2) as $P_{\Gamma}(T_{\Gamma}u-f)=0$ which gives rise to the space of admissible functions defined by

[TABLE]

and the corresponding minimization problem:

Problem 2.2**.**

Find $u\in V_{f}$ minimizing $\mathcal{J}_{\Omega}(u)$.

Notice that Problem 2.2 differs from Problem 2.1 not only due to the different notation for the boundary conditions on $\Gamma$, but also because it minimizes $J_{\Omega}$ for functions defined on the whole of $\Omega$ and not just on $\Omega_{B}$. However, the following result from [10] shows that a solution of Problem 2.2 also immediately yields a solution to Problem 2.1.

Proposition 2.3**.**

Let $u\in H_{0}^{2}(\Omega)$ be the solution of Problem 2.2. Define $\gamma_{i}=|\partial B_{i}|^{-1}(u-f_{1},\eta_{i})_{L^{2}(\Gamma)}$ for $i=1,\dots,k$. Then $(u|_{\Omega_{B}},\gamma)$ is the solution of Problem 2.1.

While Problem 2.2 allows for variable height of particles, the full problem considered in [10] has additional degrees of freedom. This is due to the fact that particles can move and rotate in the plane of the fluid membrane. To avoid mesh-deformation or re-meshing whenever particle positions change we will in the following drop the hard constraints at the particle boundaries in favor of a penalized approach. Replacing the hard curve constraints by penalty terms in the energy functional leads to:

Problem 2.4**.**

Find $u_{\varepsilon}\in H^{2}_{0}(\Omega)$ minimizing

[TABLE]

This formulation is more favorable than Problem 2.2 in so far as it admits a straightforward conforming finite element discretization without re-meshing in case of variable particle positions.

At this point we want to mention an alternative penalty formulation which is based on the idea that, for known shapes $B_{i}$, the solution $u|_{B}$ could be computed a-priori up to a constant per component $B_{i}$. We define the restriction operator

[TABLE]

Analogously to the curve constraint formulation we introduce the associated projection operator

[TABLE]

for $\psi_{i}:=\chi_{B_{i}}$ and some fixed $s\in[0,2]$. Using this notation the bulk constrained problem reads:

Problem 2.5**.**

Find $u\in H_{0}^{2}(\Omega)$ minimizing $\mathcal{J}_{\Omega}(u)$ subject to

[TABLE]

One quickly verifies that the solutions of Problem 2.2 and Problem 2.5 coincide. Analogously to Problem 2.4 a penalty formulation of Problem 2.5 is given by:

Problem 2.6**.**

Find $u_{\varepsilon}\in H_{0}^{2}(\Omega)$ and minimizing

[TABLE]

While the penalized formulations are more flexible in the sense that the domain $\Omega_{B}$ does not have to be resolved for discretization, one will be faced with the problem of balancing the penalty parameter and discretization errors. To this end we will first develop an abstract error estimate for penalized problems and then analyze the regularity of the solution of the hard constrained Problem 2.2. Using a suitable regularity result later allows to derive optimal $h$-dependent values of $\varepsilon_{i}$ for a discretization with mesh size $h$.

3. An Error Estimate for Linear Penalty Problems on Hilbert Spaces

In this section we derive an abstract energy error estimate for penalized discretizations of linearly constrained problems in a Hilbert space setting. Let $H$ be a Hilbert space and $U_{0}\subset H$ a closed subspace. We consider the affine closed subspace $U=U_{0}+u_{0}\subseteq H$ of $H$ for some given $u_{0}\in H$. In addition, let $a\colon H\times H\rightarrow\mathbb{R}$ be a symmetric positive semi-definite bounded bilinear form and $\ell\colon H\rightarrow\mathbb{R}$ a bounded linear form on $H$. We furthermore assume that $a(\cdot,\cdot)$ is coercive on $U_{0}$. In this setting we consider the affine constrained minimization problem:

Problem 3.1**.**

Find $u\in U$ minimizing

[TABLE]

Application of Lax–Milgram’s theorem yields the existence of a unique solution $u\in U$ of Problem 3.1 characterized by the variational equation

[TABLE]

From now on $u\in U$ denotes this unique solution.

Now let us assume that we want to approximate this problem by a penalty formulation over some closed linear subspace $X\subseteq H$ with $m\in\mathbb{N}$ penalty terms. Those penalty terms shall be given by symmetric positive semi-definite bounded bilinear forms $b_{i}\colon H\times H\rightarrow\mathbb{R}$ and penalty parameters $\varepsilon_{i}>0$. Denoting by $\|v\|_{c}:=\sqrt{c(v,v)}$ the semi-norm induced by a symmetric positive semi-definite bilinear form $c(\cdot,\cdot)$ on $H$ the penalized problem reads:

Problem 3.2**.**

Find $u^{X}_{\varepsilon}\in X$ minimizing

[TABLE]

For the sake of convenience we write

[TABLE]

for $\varepsilon=(\varepsilon_{i})_{i=1,\dots,m}\in\mathbb{R}_{+}^{m}=(0,\infty)^{m}$. We require that $a_{1}(\cdot,\cdot)$ is coercive on $X$, such that $a_{\varepsilon}(\cdot,\cdot)$ is also coercive for all $\varepsilon\in(0,1]^{m}$. Under these assumptions Lax–Milgram’s theorem implies existence of a unique solution $u^{X}_{\varepsilon}\in X$ for any $\varepsilon\in(0,1]^{m}$, characterized by the variational equation

[TABLE]

The following result states a Céa-type estimate for the error $u-u^{X}_{\varepsilon}$ resulting from penalization and discretization in $X$.

Theorem 3.3**.**

For the fixed $u\in U$ used in the definition of $\ell_{\varepsilon}$ suppose there exist constants $c_{i}>0$ for $i=1,\dots,m$ such that

[TABLE]

Then the error for the solution $u^{X}_{\varepsilon}$ of Problem 3.2 can be bounded by

[TABLE]

Proof.

Let $e=u-u^{X}_{\varepsilon}$. Then we get

[TABLE]

for all $v\in X$. Using Young’s inequality we can bound the first term by

[TABLE]

The definition (3.2) of the $u$-dependent penalized functional $\ell_{\varepsilon}$ yields

[TABLE]

Combining this identity with (3.3) we get

[TABLE]

Because of $v-u^{X}_{\varepsilon}\in X$ this implies

[TABLE]

where the second to last inequality follows from Young’s inequality.

Inserting the estimates (3.7) and (3.9) into (3.6) we obtain

[TABLE]

which proves the assertion. ∎

Next we consider the case where the evaluation of $a(\cdot,\cdot)$, $\ell(\cdot)$, and $b_{i}(\cdot,\cdot)$ is not performed exactly but approximated by some $\widetilde{a}$, $\widetilde{\ell}$ and $\widetilde{b}_{i}$, respectively. We use a notation analogous to $a_{\varepsilon}$ and $\ell_{\varepsilon}$:

[TABLE]

Instead of solving (3.3) for $u^{X}_{\varepsilon}$ directly one now computes a $\widetilde{u}^{X}_{\varepsilon}$ by solving

[TABLE]

In this setting we can prove the following Strang-type result:

Proposition 3.4**.**

Let $\left\|a\right\|$ be the continuity constant of $a$ with respect to the $H$-norm and let $a_{1}=a+\sum_{i=1}^{m}b_{i}$ and $\widetilde{a}_{1}=\widetilde{a}+\sum_{i=1}^{m}\widetilde{b}_{i}$ be coercive with respect to the constants $\alpha$ and $\widetilde{\alpha}$, respectively. Additionally to the assumptions of Theorem 3.3 suppose that for all $i\in\{1,\dots,m\}$ there exists $c_{i}\in\mathbb{R}$ such that for all $v\in X$

[TABLE]

Then there exists a constant $C=C(\alpha,\widetilde{\alpha},\left\|a\right\|)>0$ such that

[TABLE]

Proof.

See appendix ∎

4. Regularity of the Hard Curve Constraint Problem

In order to apply the results of the previous section it will be crucial to prove (3.4). For this purpose and in order to prove convergence rates by bounding the best approximation errors in Theorem 3.3 we will now derive regularity results for the solution of Problem 2.2.

Let $u\in V_{f}$ be the solution of the hard curve minimization formulation, Problem 2.2. Knowing the regularity of $u$ is central to proving discretization errors since the maximal regularity of the solution immediately reveals the optimal rates of convergence that one would expect for the corresponding discretization errors.

Our strategy is to rewrite the hard curve minimization problem as a system of elliptic partial differential equations to which we apply standard regularity theory for elliptic interface problems. To this end we define the fourth order elliptic differential operator $L:=\kappa\Delta^{2}-\sigma\Delta$ associated with the energy functional $\mathcal{J}$.

Proposition 4.1**.**

Suppose $(u_{1},u_{2})\in H^{2}(\Omega_{B})\times H^{2}(B)$ is a weak solution of the system of PDEs

[TABLE]

Then $u_{1}=u|_{\Omega_{B}}$ and $u_{2}=u|_{B}$ where $u$ is the solution of Problem 2.2. Conversely, if $u$ is the solution of Problem 2.2, then $(u|_{\Omega_{B}},u|_{B})$ solves (4.1).

Proof.

A proof for the equivalence of Problem 2.2 with a weak formulation of (4.1) is given in [10, Proposition 4.1]. ∎

Lemma 4.2**.**

Let $\Omega$ be a piecewise polygonal domain, let $\Gamma$ be smooth and suppose $(f_{1},f_{2})\in H^{\frac{7}{2}}(\Gamma)\times H^{\frac{5}{2}}(\Gamma)$. If all corners of $\Omega$ have an inner angle $\omega$ with $\omega\leq 126^{\circ}$, then $(u|_{\Omega_{B}},u|_{B})\in H^{4}(\Omega_{B})\times H^{4}(B)$.

Proof.

Let $g:=\frac{\sigma}{\kappa}\Delta u|_{\Omega_{B}}\in L^{2}(\Omega_{B})$. From Proposition 4.1 we know that $u|_{\Omega_{B}}$ is a weak solution of

[TABLE]

Since $g\in L^{2}(\Omega_{B})$ and $T_{\Gamma}^{1}u|_{\Omega_{B}}=f_{1}+\sum_{i=1}^{k}|\partial B_{i}|^{-1}(u|_{\Omega_{B}}-f_{1},\eta_{i})_{L^{2}(\Gamma)}\eta_{i}\in H^{\frac{7}{2}}(\Gamma)$ as well as $T_{\Gamma}^{2}u|_{\Omega_{B}}=f_{2}\in H^{\frac{5}{2}}(\Gamma)$ it follows from [11, Theorem 7.2.2.3] and the computations in [5] on the associated characteristic equation that $u|_{\Omega_{B}}\in H^{4}(\Omega_{B})$. The analogue argumentation on $B$ yields $u|_{B}\in H^{4}(B)$ and thus proves the assertion. ∎

Having the intrinsic regularity $u\in H^{2}(\Omega)$ and the piecewise regularity from Lemma 4.2 allows to show an improved global regularity result. The key ingredient is the following technical Lemma.

Lemma 4.3**.**

Let $v\in H^{2}(\Omega)$ such that $v|_{\Omega_{B}}\in H^{4}(\Omega_{B})$ and $v|_{B}\in H^{4}(B)$. Then $v\in H^{2+\frac{1}{2}-\delta}(\Omega)$ for all $\delta>0$.

Proof.

See appendix. ∎

Combining Lemma 4.2 and Lemma 4.3 now gives:

Corollary 4.4**.**

Let the assumptions from Lemma 4.2 hold. Then $u\in H^{\frac{5}{2}-\delta}(\Omega)$ for every $\delta>0$.

5. Discretization Errors

In this section we apply the results from the previous sections to derive discretization errors for penalized finite element approximations of the form Problem 2.4 and Problem 2.6 but with $h$-dependent penalty parameters.

We suppose from now on that $\Omega$ is a rectangular domain. In view of Lemma 4.2 this implies that $u|_{\Omega_{B}}\in H^{4}(\Omega_{B})$, $u|_{B}\in H^{4}(B)$ and $u\in H^{\frac{5}{2}-\delta}(\Omega)$ for any $\delta>0$. On the domain $\Omega$ we establish a quadrilateral grid equipped with Bogner–Fox–Schmit finite elements (see e. g. [6]). Given the set of grid nodes $N$ and the set of multi-indices $\Lambda=\{(0,0),(1,0),(0,1),(1,1)\}$, the resulting finite element space is spanned by a basis $(\psi_{p,\alpha})_{p\in N,\alpha\in\Lambda}$ of piecewise bi-cubic polynomials such that each $\psi_{p,\alpha}$ satisfies

[TABLE]

where $\partial^{\beta}$ is the partial derivative associated with the multi-index $\beta$ and $\delta_{a,b}$ is the Kronecker delta. I. e., the degrees of freedom are the values, first order partial derivatives and mixed second order partial derivatives at the vertices (see Figure 5.1).

By setting the degrees of freedom on $\partial\Omega$ to zero this yields a conforming subspace of $C^{1}(\overline{\Omega})\cap H_{0}^{2}(\Omega)$. Given a family $(S_{h})_{h\in I}$ of such finite element discretizations over $\Omega$ based on quasi-uniform grids with mesh size $h$ we will use the following well known result: There exists a constant $c>0$ such that for all $K\in[0,4]$ the approximation estimate

[TABLE]

holds true for all $h\in I$. (This is a direct consequence of classical approximation estimates, see e. g. [6, Theorem 3.1.4], and interpolation theory of function spaces, see e. g. [14, Theorem 1.1.6].)

In order to put our minimization problems into the variational framework used in Section 3 we introduce the bilinear form

[TABLE]

together with the curve penalty terms

[TABLE]

and the bulk penalty term

[TABLE]

Then the solution $u\in V_{f}\subset H_{0}^{2}(\Omega)$ of Problem 2.2 satisfies the variational equation

[TABLE]

Discretization of the penalized curve constrained Problem 2.4 in $S_{h}$ leads to the soft curve problem: Find $u_{\varepsilon}^{h}\in S_{h}$ minimizing

[TABLE]

which is equivalently characterized by the variational equation

[TABLE]

Analogously, discretization of the penalized bulk constrained Problem 2.6 in $S_{h}$ leads to the soft bulk problem: Find $u_{\varepsilon}^{h}\in S_{h}$ minimizing

[TABLE]

with the equivalent variational equation

[TABLE]

At this point we want to emphasize that the solutions of (5.4) and (5.6) are in general different. We nevertheless refer in a slight abuse of notation to both solutions as $u^{h}_{\varepsilon}$. In the following it will, however, be clear from the context whether $u^{h}_{\varepsilon}$ denotes the solution of the soft curve problem or the soft bulk problem.

In the next result we show that the soft curve formulation meets the requirement from Theorem 3.3, which we will use afterwards to get an asymptotic error estimate.

Lemma 5.1**.**

For all $v\in H_{0}^{2}(\Omega)$ the solution $u$ of Problem 2.2 satisfies

[TABLE]

Here $[w]=w|_{\Omega_{B}}-w|_{B}$ denotes the jump of the function $w$ on $\Gamma$.

Proof.

By integration by parts and as of $u\in H^{2}(\Omega)$, $u|_{\Omega_{B}}\in H^{4}(\Omega_{B})$, $u|_{B}\in H^{4}(B)$ and $Lu=0$ we are able to state for all $v\in H_{0}^{2}(\Omega)$

[TABLE]

Since $\widetilde{v}_{i}:=(\operatorname{id}_{L^{2}(\Gamma)}-P_{\Gamma}^{i})T_{\Gamma}^{i}v\in H^{\frac{5}{2}-i}(\Gamma)$ and because $\Gamma$ is sufficiently smooth there exists a function $w\in H_{0}^{2}(\Omega)$ such that $T_{\Gamma}^{i}w=\widetilde{v}_{i}$, [13, Theorem 9.4]. In particular, we have $P_{\Gamma}^{i}T_{\Gamma}^{i}w=0$ which yields $w\in V_{0}$ and

[TABLE]

by applying the variational equation (5.2) for $u$. Combined with the Cauchy–Schwarz inequality this implies

[TABLE]

which was to be shown. ∎

Now we are in the situation to show the main results, namely the discretization error estimate for the discretizations (5.3) and (5.5).

Theorem 5.2**.**

Let $u^{h}_{\varepsilon}\in S_{h}$ be the solution of the discrete soft curve problem (5.4) and assume that $\varepsilon_{1}=c_{1}h^{\lambda_{1}}$ and $\varepsilon_{2}=c_{2}h^{\lambda_{2}}$. For any $\delta>0$ define

[TABLE]

Then there exists a constant $c>0$ independent of $h$ such that

[TABLE]

In particular, $\|u-u^{h}_{\varepsilon}\|_{H^{2}(\Omega)}\in O(h^{\frac{1}{2}-\delta})$ for $\lambda_{1}\in[1-2\delta,3-2\delta]$ and $\lambda_{2}=1-2\delta$.

Proof.

In this proof we use the notation $A\lesssim B$ whenever $A\leq cB$ holds with a constant $c$ that is independent of $\varepsilon_{1}$, $\varepsilon_{2}$ and $h$.

As of Lemma 5.1 we can apply Theorem 3.3. Using this together with the coercivity of $a(\cdot,\cdot)$ on $H_{0}^{2}(\Omega)$ we conclude for all $v\in S_{h}$

[TABLE]

Note that $a(\cdot,\cdot)$ is continuous on $H^{2}(\Omega)$ and $b_{i}(\cdot,\cdot)$ is continuous on $H^{i-\frac{1}{2}+\delta}(\Omega)$. Furthermore, $u\in H^{\frac{5}{2}-\delta}(\Omega)$ according to Corollary 4.4. Choosing the interpolant $v=\overline{u}$ and applying the interpolation estimate (5.1) yields

[TABLE]

This proves the assertion. ∎

We can proceed similarly to prove convergence rates for the soft bulk formulation.

Lemma 5.3**.**

Suppose that there is a constant $c_{0}>0$ such that for every $h\in I$, $K\in[0,2]$ and $k\in[0,K]$ the inverse estimate

[TABLE]

is fulfilled. Let $u$ be the solution of Problem 2.2 and $\delta>0$. Then there exists a constant $c>0$ such that for all $h\in I$ and $v\in S_{h}$

[TABLE]

Proof.

From the proof of Lemma 5.1 we know for all $v\in S_{h}$

[TABLE]

And as in the proof of Lemma 5.1 there exists a function $w\in V_{0}\subseteq H_{0}^{2}(\Omega)$ such that $T_{\Gamma}^{i}w=T_{\Gamma}^{i}(\operatorname{id}_{B}-P_{B})T_{B}v$ for $i\in\{1,2\}$. From this we conclude

[TABLE]

as in the proof of Lemma 5.1. Furthermore note that the operators $T_{\Gamma}^{i}$ are continuous and linear over $H^{\frac{3}{2}+\delta}(B)$ and thus the corresponding operator norms ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|T_{\Gamma}^{i}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\delta}$ are bounded. It follows that

[TABLE]

Together with (5.7) this implies

[TABLE]

which proves the assertion. ∎

Note that (5.7) is indeed an assumption. We can not apply the standard finite element inverse estimates here as the grid is in general not matched to the domain $B$. A closer look to the proof of the standard inverse estimate reveals that the constant in (5.7) would go to infinity as $\left|E_{h}\cap B\right|/\left|E_{h}\right|\to 0$ for $I\ni h\to 0$ and some element $E_{h}$ in the grid associated to $S_{h}$. This means that the above inverse estimate over $B$ is only valid for grids that resolve $B$ sufficiently well.

Theorem 5.4**.**

Let the assumptions from Lemma 5.3 hold, $u^{h}_{\varepsilon}\in S_{h}$ the solution of the discrete soft bulk problem (5.6), and $\varepsilon=c_{0}h^{\lambda}$ for some $\delta>0$. Define

[TABLE]

Then there exists a constant $c>0$ independent of $h$ such that

[TABLE]

In particular, $\|u-u^{h}_{\varepsilon}\|_{H^{2}(\Omega)}\in\mathcal{O}(h^{\frac{1}{2}-\delta})$ for $s\leq\frac{3}{2}+\delta$ and $\lambda=4-2s$.

Proof.

As of Lemma 5.3 we can apply Theorem 3.3. Noting continuity of $b_{B}(\cdot,\cdot)$ on $H^{s}(\Omega)$ we obtain as in the proof of Theorem 5.2 that

[TABLE]

Again choosing the interpolant $v=\overline{u}$ and applying the interpolation estimate (5.1) finally gives

[TABLE]

and thus proves the assertion. ∎

6. Numerical Example

6.1. A symmetric example problem with known exact solution

We select a problem with a known, analytically computable solution. We consider the circular domain $\hat{\Omega}=\{x\in\mathbb{R}^{2}\mathrel{|}|x|<r_{2}\}$ with one embedded particle $\hat{B}=\{x\in\mathbb{R}^{2}\mathrel{|}|x|<r_{1}\}$ for $0<r_{1}<r_{2}<1$. On the boundary of $\hat{\Omega}$ we consider homogeneous Dirichlet boundary conditions whereas on the particle boundary $\Gamma_{1}=\partial B_{1}$ we impose variable-height coupling conditions of the form (2.1), i. e.,

[TABLE]

for variable $\gamma^{1}\in\mathbb{R}$ with the given functions

[TABLE]

We selected the parameters

[TABLE]

Then the exact solution of Problem 2.1 on $\hat{\Omega}$ is the fourfold symmetric function $u\in H^{2}(\hat{\Omega})$ given by

[TABLE]

with the constants

[TABLE]

In order to discretize the circular domain $\hat{\Omega}$ with the presented Bogner–Fox–Schmit discretization, we embed $\hat{\Omega}$ into the larger rectangular domain $\Omega=[-1,1]^{2}$ and treat the condition on $\Gamma_{2}=\partial\hat{\Omega}$ like a second particle boundary condition. It is easily seen that the solution of Problem 2.1 on $\Omega$ is obtained by extending the solution given on $\hat{\Omega}$ by zero. The solution of this problem is depicted in Fig. 6.1.

6.2. Discretization, quadrature, and error measurement

For the finite element discretization $S_{h}\subseteq H_{0}^{2}(\Omega)$ we divide $\Omega$ into uniform squares with edges of length $h$ and equip the resulting grid with a Bogner–Fox–Schmit finite element basis. The boundary conditions over $\partial\Omega$ are enforced by setting the corresponding degrees of freedom to zero while the boundary conditions over the $\Gamma_{i}$ are replaced by penalized constraints. The resulting soft curve problem with penalized curve constraints is then given by (5.4) while the soft bulk problem with penalized bulk constraints is given by (5.6) where the area of the virtual second particle is now given by $B_{2}=\Omega\setminus\hat{\Omega}$.

We note that we drop the projection for the constraints on $\Gamma_{2}$ and $B_{2}$ because we do not allow variable height there. While our analysis is only formulated in terms of the variable-height constraints incorporated using projections, this is in fact no limitation, since all arguments directly carry over to the case without these projections.

Since the assembly of the bilinear forms $b^{i}_{\Gamma}(\cdot,\cdot)$ and $b_{B}(\cdot,\cdot)$ involve integrals over the curves $\Gamma_{i}$ or the nontrivial domains $B_{i}$, respectively, we briefly mention how these can be approximated. Integrals over full grid cells are evaluated exactly using standard tensor Gaussian quadrature rules. For quadrature over the $\Gamma_{i}$ we use a trigonometric Gauss quadrature as described in [7] which in our case is both exact for integration of finite element functions and for integration of $u$. Integrals over the non-rectangular domains $B_{i}$ are approximated using the local parameterization quadrature method introduced in [16], which is no longer exact.

The discretization error is measured in terms of the norm $\|\Delta(u^{h}_{\varepsilon}-u)\|_{L^{2}(\Omega)}$ on $H_{0}^{2}(\Omega)$. It is important to note that the error cannot be computed accurately by simple quadrature over the grid elements, because $u$ is not smooth across the $\Gamma_{i}$. Instead we split the integration into the subdomains $B_{1}$, $B_{2}=\Omega\setminus\hat{\Omega}$, and $\Omega\setminus(B_{1}\cup B_{2})$ where we apply the quadrature rules described above.

6.3. Results for the soft curve formulation

In case of the soft curve formulation we expect convergence with order $O(h^{\frac{1}{2}})$ due to Theorem 5.2 if the penalty parameters are chosen suitably. For the numerical examples we selected penalty parameters $\varepsilon=(ch^{\lambda_{1}},ch)$ with $\lambda_{1}\in[1,3]$ according to Theorem 5.2. The constant was fixed to be $c=10^{-3}$.

Figure 6.2 shows the behavior of the $H^{2}$-discretization error over the mesh size $h$ for uniform grids and penalty parameters with $\lambda_{1}\in\{1,2,3\}$. In accordance with Theorem 5.2 we observe convergence of at least order $O(h^{1/2})$. While the observed rate seems to behave like $O(h^{1/2})$ for $\lambda_{1}\in\{2,3\}$ it seems that the rate is slightly better for $\lambda_{1}=1$. A possible explanation for this observation is that the error contribution for the term $\|\cdot\|_{b_{1}}$ in the proof of Theorem 5.2 can be improved by smaller values of $\lambda_{1}$. However, this does not lead to a better theoretical bound due to other dominating terms.

Figure 6.3 and Figure 6.4 show the $H^{1}$- and $L^{2}$-discretization errors, respectively, for the same set of example problems. In both cases we observe convergence with the order $O(h)$. While improved convergence order for weaker norms is a well-known property of many discretizations, we emphasize that this is not covered by the theory presented here and may be considered in the future with a refined analysis.

6.4. Results for the soft bulk formulation

Because we use uniform grids and do not refine with regard to the geometry of the $\Gamma_{i}$ we lose control over the inverse estimate assumption (5.7). Thus we cannot prove condition (3.4) in general. Assuming that (3.4) still holds true, we could use Theorem 5.4 in case of exact integration. However, as discussed above, we only approximate integrals using the quadrature from [16]. If this is accurate enough, we can in view of the Strang-type result in Proposition 3.4 expect convergence of order $O(h^{1/2})$ when using $\varepsilon=ch^{4-2s}$ and the penalty norm $\|\cdot\|_{H^{s}(B_{i})}$ for $s\in[0,1]$. For our numerical experiments we considered $c=10^{-3}$ and $s\in\{0,1\}$.

Figure 6.5 shows the behaviour of the discretization error over the mesh size $h$ for uniform grids and $s\in\{0,1\}$. Again the observed order is in accordance with the expected order $O(h^{1/2})$. However, especially for $s=1$ this is perturbed by strong oscillations. This is maybe due to the fact that the constant resulting from (3.4) is not uniformly bounded and can strongly vary depending on how the curve $\Gamma$ intersects mesh elements. Again the $H^{1}$- and $L^{2}$-errors depicted in Figure 6.6 and Figure 6.7, respectively, exhibit a similar behavior but the improved convergence order $O(h)$.

6.5. A non-symmetric example problem

While the first example allowed to compute the exact solution due to its symmetry properties, we will now assess the discretization for a non-symmetric example with several particles but without known exact solution.

We consider the domain $\Omega=[-1,1]^{2}$ with homogeneous Dirichlet boundary conditions on $\partial\Omega$ and we again select the parameters

[TABLE]

This time we consider four elliptical particles embedded non-symmetrically into the membrane domain $\Omega$. The ellipses’ major and minor axes have length $0.4$ and $0.2$, respectively. Their positions and orientations are depicted in Figure 6.8.

Regarding the particle–membrane coupling conditions we assume that all particles have a constant height profile $f^{i}_{1}=0$ and a constant slope $f^{i}_{2}$ varying for different particles. More precisely, we use

[TABLE]

The solution resulting from that setup is depicted in Fig. 6.8.

Since we do not know the exact solution for this problem we estimate the error by comparing our discrete solutions to a reference solution computed on a finer grid. We computed the solutions for the different problem formulations on uniform grids with mesh size $h=2^{-k}$ for $k=3,\dots,6$ while the reference solution was computed on a uniform grid with mesh size $h=2^{-8}$ using the soft curve formulation with $\varepsilon=10^{-3}(h^{3},h)$.

The (approximate) discretization errors for soft curve and soft bulk formulation are shown in Figure 6.9 and Figure 6.10, respectively. As before we used $\varepsilon=(ch^{\lambda},ch)$ with $\lambda\in\{1,2,3\}$ for the soft curve formulation and $\varepsilon=ch^{4-2s}$ with $s\in\{0,1\}$ for the soft bulk formulation. In both cases we selected $c=10^{-3}$.

For the soft bulk formulation the projection $P_{B}$ may lead to a densely populated matrix for those degrees of freedom located inside of the particles. To avoid the resulting computational effort for $s=1$ we used the $H^{1}(B_{i})$-norm in the form

[TABLE]

Then the penalty term that incorporates the $H^{1}(B_{i})$–projections takes the form

[TABLE]

which results in a sparse matrix again.

For both formulations, soft curve and soft bulk, we observe a rate that is slightly better than the expected order of convergence, namely $O(h^{1/2})$. Note that we restricted the mesh size $h=2^{-k}$ to powers of two in this example to simplify comparison with the reference solution. Furthermore, in the case of soft bulk constraints with $s=0$ we were not able to carry out the computations for $h=2^{-6}$ due to hardware restrictions. This explains the much smaller number of data points in the plots.

As before we also computed the $H^{1}$- and $L^{2}$-errors for the same set of example problems. For the soft curve formulation we again observe that the order of the $H^{1}$- and $L^{2}$-error depicted in Figure 6.11 and Figure 6.12, respectively, is essentially squared in comparison with the $H^{2}$-error leading to a convergence order which is approximately $O(h)$. The situation is different for the soft bulk formulation. Here, the $H^{1}$-error depicted in Figure 6.13 is of order $O(h^{3/2})$ and the $L^{2}$-error depicted in Figure 6.14 is of order $O(h^{5/2})$. While the improved $H^{1}$- and $L^{2}$-order is not covered by the presented theory anyway, we also cannot explain the surprising difference in the observed order for the symmetric and non-symmetric example problems.

7. Acknowledgements

This research has been funded by Freie Universität Berlin and was carried out in the associated project AP1 ”Particles in lipid bilayers” of the CRC 1114 ”Scaling Cascades in Complex Systems”.

8. Appendix

Proof of Lemma 4.3.

Let $u_{1}:=u|_{\Omega_{B}}$ and $u_{2}:=u|_{B}$ as well as $\Omega_{1}:=\Omega_{B}$ and $\Omega_{2}:=B$. Define $\varphi_{\alpha}(x,y):=\frac{|D^{\alpha}u(x)-D^{\alpha}u(y)|^{2}}{\|x-y\|^{2s+2}}$. For $s\in(0,\frac{1}{2})$ the Sobolev–Slobodeckij norm of $u$ is given by

[TABLE]

Splitting the domain of integration yields

[TABLE]

Using the layercake principle we write

[TABLE]

By the assumptions from Section 2 it holds that $\Omega_{1}$ and $\Omega_{2}$ are Lipschitz domains. By virtue of the Sobolev embedding theorems we conclude $u_{i}\in C^{2}(\overline{\Omega_{i}})$, see for example [1, Theorem 5.4]. Consequently, the constant $C:=\max_{|\alpha|=2}\sup_{\Omega_{1}\times\Omega_{2}}|D^{\alpha}u_{1}(x)-D^{\alpha}u_{2}(y)|^{2}\in\mathbb{R}_{>0}$ exists and thus the implication

[TABLE]

holds for all $(x,y)\in\Omega_{1}\times\Omega_{2}$. Let $\Gamma:=\partial\Omega_{1}\cap\partial\Omega_{2}$. From

[TABLE]

we infer

[TABLE]

This expression is finite for $s<\frac{1}{2}$. This implies $u\in H^{2+\frac{1}{2}-\delta}(\Omega)$ for all $\delta>0$ as was to be shown. ∎

Proof of Proposition 3.4.

For brevity we define $\widetilde{u}:=\widetilde{u}^{X}_{\varepsilon}$ and assume without loss of generality that $\varepsilon=(\varepsilon_{i})_{i=1,\dots,m}\leq 1$ holds component-wise. Let $v\in X$ and define $\widetilde{w}:=v-\widetilde{u}$. We have by (3.11)

[TABLE]

Using the continuity of $a$, the coercivity of $\widetilde{a}_{\varepsilon}$ and $\varepsilon\leq 1$ this leads to the inequality

[TABLE]

On the other hand we have the equality

[TABLE]

which by application of (3.10), Cauchy–Schwarz inequality, coercivity of $a_{\varepsilon}$, definition of $\widetilde{a}_{\varepsilon}$ and $\widetilde{\ell}_{\varepsilon}$ and $\varepsilon\leq 1$ yields

[TABLE]

Then from the triangle inequality and (8.2) inserted into (8.1) we get

[TABLE]

which proves the statement by taking the infimum over all $v\in X$ and replacing $\widetilde{w}=v-\widetilde{u}$ by the supremum over all $w\in X$. ∎