Discontinuous Galerkin Finite Element Methods for the Landau-de Gennes Minimization Problem of Liquid Crystals

TL;DR

This paper develops and analyzes discontinuous Galerkin finite element methods for modeling equilibrium states of liquid crystals, providing stability, error estimates, and convergence results for the nonlinear Landau-de Gennes minimization problem.

Contribution

It introduces a stable discontinuous Galerkin approach for the nonlinear Landau-de Gennes model, with rigorous error analysis and convergence proofs.

Findings

Stable discretization of the nonlinear problem

Optimal a priori error estimates in energy norm

Quadratic convergence of Newton's method

Abstract

We consider a system of second order non-linear elliptic partial differential equations that models the equilibrium configurations of a two dimensional planar bistable nematic liquid crystal device. Discontinuous Galerkin finite element methods are used to approximate the solutions of this nonlinear problem with non-homogeneous Dirichlet boundary conditions. A discrete inf-sup condition demonstrates the stability of the discontinuous Galerkin discretization of a well-posed linear problem. We then establish the existence and local uniqueness of the discrete solution of the non-linear problem. An a priori error estimate in the energy norm is derived and a best approximation property is demonstrated. Further, we prove the quadratic convergence of Newton's iterates along with complementary numerical experiments.

Tables

Table 1. Table 1: Numerical errors and orders of convergence in dG and 𝐋 2 superscript 𝐋 2 \mathbf{L}^{2} norms using piecewise linear polynomials.

$h$	${\Vert \left|\Psi -{\Psi }_{\text{dG}}^i\right|\Vert }_{\text{dG}}$	Order	${\Vert \left|\Psi -{\Psi }_{\text{dG}}^i\right|\Vert }_{{𝐋}^2}$	Order
0.3535	0.69481292E-1	1.10439646	0.37102062E-2	2.13551126
0.1767	0.29094563E-1	1.06282552	0.74601448E-3	2.08646636
0.0883	0.13383488E-1	1.04108837	0.16427963E-3	2.04993501
0.0441	0.64047595E-2	1.02761427	0.38801373E-4	2.03044486
0.0220	0.31296010E-2	1.01899219	0.94539471E-5	2.02007895
0.0141	0.19860395E-2	-	0.38377918E-5	-

Table 2. Table 2: Numerical errors and orders of convergence in dG and 𝐋 2 superscript 𝐋 2 \mathbf{L}^{2} norms using piecewise quadratic polynomials.

$h$	${\Vert \left|\Psi -{\Psi }_{\text{dG}}^i\right|\Vert }_{\text{dG}}$	Order	${\Vert \left|\Psi -{\Psi }_{\text{dG}}^i\right|\Vert }_{{𝐋}^2}$	Order
0.3535	0.20812609E-1	2.25360832	0.59563518E-3	3.22432824
0.1767	0.35369466E-2	2.17037967	0.50612758E-4	3.13307115
0.0883	0.71656832E-3	2.12009381	0.53034861E-5	3.08714692
0.0441	0.15828049E-3	2.08449035	0.60454668E-6	3.05924399
0.0220	0.36909294E-4	2.05973780	0.71945319E-7	3.04117003
0.0141	0.14720321E-4	-	0.18516670E-7	-

Table 3. Table 3: Numerical errors and orders of convergence in dG and 𝐋 2 superscript 𝐋 2 \mathbf{L}^{2} norms using piecewise cubic polynomials.

$h$	${\Vert \left|\Psi -{\Psi }_{\text{dG}}^i\right|\Vert }_{\text{dG}}$	Order	${\Vert \left|\Psi -{\Psi }_{\text{dG}}^i\right|\Vert }_{{𝐋}^2}$	Order
0.3535	0.44002858E-2	2.76328893	0.79005568E-4	3.7473338
0.1767	0.81936582E-3	2.85612226	0.76031244E-5	3.8488866
0.0883	0.12772385E-3	2.92217824	0.59938677E-6	3.9184264
0.0441	0.17487379E-4	2.95474574	0.41209120E-7	3.9524773
0.0220	0.22702291E-5	2.96925442	0.26800851E-8	3.9677981
0.0141	0.60334919E-6	-	0.45615227E-9	-

Table 4. Table 4: Initial conditions for Newton’s method.

Solution	$x=0$	$x=1$	$y=0$	$y=1$
D1	$\pi /2$	$\pi /2$	0	0
D2	$\pi /2$	$\pi /2$	$\pi$	$\pi$
R1	$\pi /2$	$\pi /2$	$\pi$	0
R2	$\pi /2$	$\pi /2$	0	$\pi$
R3	$3\pi /2$	$\pi /2$	$\pi$	$\pi$
R4	$\pi /2$	$3\pi /2$	$\pi$	$\pi$

Table 5. Table 5: Numerical energy, errors and convergence rates for D1 solution in dG and 𝐋 2 superscript 𝐋 2 \mathbf{L}^{2} norm.

$h$	Energy	${\Vert \left|{\Psi }_{\text{dG}}^n-{\Psi }_{\text{dG}}^i\right|\Vert }_{\text{dG}}$	Order	${\Vert \left|{\Psi }_{\text{dG}}^n-{\Psi }_{\text{dG}}^i\right|\Vert }_{{𝐋}^2}$	Order
0.0883	77.80650525	3.42185356	0.90770236	0.23846079E-1	1.75592834
0.0441	77.90383430	1.90966722	0.94082516	0.78385873E-2	1.83134924
0.0220	77.92229141	1.00706128	0.95848055	0.24465831E-2	1.98287313
0.0110	77.94112012	0.51823233	-	0.61895017E-3	-

Table 6. Table 6: Numerical energy, errors and convergence rates for R1 solution in dG and 𝐋 2 superscript 𝐋 2 \mathbf{L}^{2} norm.

$h$	Energy	${\Vert \left|{\Psi }_{\text{dG}}^n-{\Psi }_{\text{dG}}^i\right|\Vert }_{\text{dG}}$	Order	${\Vert \left|{\Psi }_{\text{dG}}^n-{\Psi }_{\text{dG}}^i\right|\Vert }_{{𝐋}^2}$	Order
0.0883	86.44084273	3.46656699	0.90994382	0.26417589E-1	1.79599363
0.0441	86.53931303	1.92787441	0.94166745	0.82372577E-2	1.85335905
0.0220	86.55785529	1.01547693	0.95848131	0.25178324E-2	1.99673627
0.0110	86.57670525	0.52256273	-	0.63088371E-3	-

Table 7. Table 7: Numerical errors and convergence rates for solution of ( 6.10 ) in dG and 𝐋 2 superscript 𝐋 2 \mathbf{L}^{2} norm.

Max $h$	${\Vert \left|\Psi -{\Psi }_{\text{dG}}^i\right|\Vert }_{\text{dG}}$	Order	${\Vert \left|\Psi -{\Psi }_{\text{dG}}^i\right|\Vert }_{{𝐋}^2}$	Order
0.195	0.19019778	1.02678772	0.70751443E-2	2.07685808
0.999 E-1	0.89721415 E-1	1.01083542	0.11630918E-2	1.92959391
0.333 E-1	0.29026972 E-1	0.99787440	0.13457407E-3	1.90301464
0.199 E-1	0.17390702 E-1	0.99496176	0.49154393E-4	1.86299044
0.142 E-1	0.12534273 E-1	1.00852040	0.26717823E-4	1.89491873
0.111 E-1	0.97316088 E-2	-	0.14469534E-4	-

Equations

$$\mathcal{E}({\Psi }_{ϵ})={\int }_{\Omega }({\left|\nabla {\Psi }_{ϵ}\right|}^2+{ϵ}^{-2}({\left|{\Psi }_{ϵ}\right|}^2-1{)}^2)\text{dx},$$

$$-\Delta {\Psi }_{ϵ}=2{ϵ}^{-2}(1-{\left|{\Psi }_{ϵ}\right|}^2){\Psi }_{ϵ}\text{ in }\Omega \text{ and }{\Psi }_{ϵ}=\mathbf{g}\text{ on }\partial \Omega .$$

$$\mathrm{N}({\Psi }_{ϵ};\Phi ):=A({\Psi }_{ϵ},\Phi )+B({\Psi }_{ϵ},{\Psi }_{ϵ},{\Psi }_{ϵ},\Phi )+C({\Psi }_{ϵ},\Phi )=0\text{ for all }\Phi \in \mathbf{V},$$

$$A(\Theta ,\Phi ):=a({\theta }_1,{\varphi }_1)+a({\theta }_2,{\varphi }_2),C(\Theta ,\Phi ):=c({\theta }_1,{\varphi }_1)+c({\theta }_2,{\varphi }_2),$$

$$B(\Xi ,\mathbf{\eta },\Theta ,\Phi ):=\frac{C_{ϵ}}{3}{\int }_{\Omega }\left((\Xi \cdot \mathbf{\eta })(\Theta \cdot \Phi )+2(\Xi \cdot \Theta )(\mathbf{\eta }\cdot \Phi )\right)\text{dx}=\frac{1}{3}(3b({\xi }_1,{\eta }_1,{\theta }_1,{\varphi }_1)+3b({\xi }_2,{\eta }_2,{\theta }_2,{\varphi }_2)$$

$$+2b({\xi }_2,{\eta }_1,{\theta }_2,{\varphi }_1)+2b({\xi }_1,{\eta }_2,{\theta }_1,{\varphi }_2)+b({\xi }_2,{\eta }_2,{\theta }_1,{\varphi }_1)+b({\xi }_1,{\eta }_1,{\theta }_2,{\varphi }_2)),$$

$$\text{and for }\xi ,\eta ,\theta ,\varphi \in H^1(\Omega ),a(\theta ,\varphi ):={\int }_{\Omega }\nabla \theta \cdot \nabla \varphi \text{dx},b(\xi ,\eta ,\theta ,\varphi ):=C_{ϵ}{\int }_{\Omega }\xi \eta \theta \varphi \text{dx},$$

$$\text{ and }c(\theta ,\varphi ):=-C_{ϵ}{\int }_{\Omega }\theta \varphi \text{dx}.$$

$$B({\Psi }_{ϵ},{\Psi }_{ϵ},{\Psi }_{ϵ},\Phi )=b(u,u,u,{\varphi }_1)+b(v,v,u,{\varphi }_1)+b(v,v,v,{\varphi }_2)+b(u,u,v,{\varphi }_2).$$

$${\alpha }_0{\Vert \varphi \Vert }_{L^2(\Omega )}\le {\left|\varphi \right|}_{1,\Omega }\text{ for all }\varphi \in H_0^1(\Omega ).$$

$$A(\Theta ,\Phi )\le {|||\Theta |||}_1{|||\Phi |||}_1,A(\Theta ,\Theta )\ge C_{{\alpha }_0}{|||\Theta |||}_1^2,$$

$$\text{ and }C(\Theta ,\Phi )\le C_{ϵ}{|||\Theta |||}_1{|||\Phi |||}_1,$$

$$B(\Xi ,\mathbf{\eta },\Theta ,\Phi )\lesssim {ϵ}^{-2}{|||\Xi |||}_1{|||\mathbf{\eta }|||}_1{|||\Theta |||}_1{|||\Phi |||}_1,$$

$$B(\Xi ,\mathbf{\eta },\Theta ,\Phi )\lesssim {ϵ}^{-2}{|||\Xi |||}_2{|||\mathbf{\eta }|||}_2{|||\Theta |||}_0{|||\Phi |||}_0.$$

$$b(\xi ,\eta ,\theta ,\varphi )\le C_{ϵ}{\Vert \xi \Vert }_{0,4}{\Vert \eta \Vert }_{0,4}{\Vert \theta \Vert }_{0,4}{\Vert \varphi \Vert }_{0,4}\lesssim {ϵ}^{-2}{\Vert \xi \Vert }_1{\Vert \eta \Vert }_1{\Vert \theta \Vert }_1{\Vert \varphi \Vert }_1.$$

$$b(\xi ,\eta ,\theta ,\varphi )\lesssim {ϵ}^{-2}{\Vert \xi \Vert }_{0,\infty }{\Vert \eta \Vert }_{0,\infty }{\Vert \theta \Vert }_0{\Vert \varphi \Vert }_0\lesssim {ϵ}^{-2}{\Vert \xi \Vert }_2{\Vert \eta \Vert }_2{\Vert \theta \Vert }_0{\Vert \varphi \Vert }_0\text{for all }\theta ,\varphi \in \text{V},$$

$$0<\beta :=\underset{\begin{array}{c}\Theta \in \mathbf{V}\\ {|||\Theta |||}_1=1\end{array}}{\inf }\underset{\begin{array}{c}\Phi \in \mathbf{V}\\ {|||\Phi |||}_1=1\end{array}}{\sup }⟨\mathrm{DN}({\Psi }_{ϵ})\Theta ,\Phi ⟩=\underset{\begin{array}{c}\Phi \in \mathbf{V}\\ {|||\Phi |||}_1=1\end{array}}{\inf }\underset{\begin{array}{c}\Theta \in \mathbf{V}\\ {|||\Theta |||}_1=1\end{array}}{\sup }⟨\mathrm{DN}({\Psi }_{ϵ})\Theta ,\Phi ⟩,$$

$${\mathrm{N}}_h({\Psi }_{\text{dG}};{\Phi }_{\text{dG}}):=A_{\text{dG}}({\Psi }_{\text{dG}},{\Phi }_{\text{dG}})+B_{\text{dG}}({\Psi }_{\text{dG}},{\Psi }_{\text{dG}},{\Psi }_{\text{dG}},{\Phi }_{\text{dG}})+C_{\text{dG}}({\Psi }_{\text{dG}},{\Phi }_{\text{dG}})=L_{\text{dG}}({\Phi }_{\text{dG}}),$$

$$A_{\text{dG}}(\Theta ,\Phi )$$

$$B_{\text{dG}}(\Xi ,\mathbf{\eta },\Theta ,\Phi )$$

$$=\frac{1}{3}(3b_{\text{dG}}({\xi }_1,{\eta }_1,{\theta }_1,{\varphi }_1)+3b_{\text{dG}}({\xi }_2,{\eta }_2,{\theta }_2,{\varphi }_2)+2b_{\text{dG}}({\xi }_2,{\eta }_1,{\theta }_2,{\varphi }_1)+2b_{\text{dG}}({\xi }_1,{\eta }_2,{\theta }_1,{\varphi }_2)$$

$$+b_{\text{dG}}({\xi }_2,{\eta }_2,{\theta }_1,{\varphi }_1)+b_{\text{dG}}({\xi }_1,{\eta }_1,{\theta }_2,{\varphi }_2)),$$

$$a_h(\theta ,\varphi ):=\sum _{T\in \mathcal{T}}{\int }_T\nabla \theta \cdot \nabla \varphi \text{dx},J(\theta ,\varphi ):=\sum _{E\in \mathcal{E}}{\int }_E\{\frac{\partial \theta }{\partial \eta }\}[\varphi ]\text{ds},\text{ and }{J}^{\sigma }(\theta ,\varphi ):=\sum _{E\in \mathcal{E}}{\int }_E\frac{\sigma }{h}[\theta ][\varphi ]\text{ds},$$

$$b_{\text{dG}}(\xi ,\eta ,\theta ,\varphi ):=C_{ϵ}\sum _{T\in \mathcal{T}}{\int }_T\xi \eta \theta \varphi \text{dx},c_{\text{dG}}(\theta ,\varphi ):=-C_{ϵ}\sum _{T\in \mathcal{T}}{\int }_T\theta \varphi \text{dx},$$

$$\text{ and }{l}_{\text{dG}}^i(\varphi ):=\lambda \sum _{E\in {\mathcal{E}}_D}{\int }_E\frac{\partial \varphi }{\partial \eta }{g}_i\text{ds}+\sum _{E\in {\mathcal{E}}_D}{\int }_E\frac{\sigma }{h}{g}_i\varphi \text{ds}\text{ for }1\le i\le 2.$$

$${|||\Psi -{\Psi }_{\mathrm{dG}}|||}_{\mathrm{dG}}\lesssim h.$$

$${|||\Psi -{\Psi }_{\mathrm{dG}}|||}_{\mathrm{dG}}\lesssim (1+{ϵ}^{-2})\underset{{\Theta }_{\mathrm{dG}}\in {\mathbf{V}}_h}{\min }{|||\Psi -{\Theta }_{\mathrm{dG}}|||}_{\mathrm{dG}}.$$

$${|||\Psi -{\Psi }_{\mathrm{dG}}|||}_0\lesssim h^2(1+(1+{ϵ}^{-2}{)}^2).$$

$$A_{\mathrm{dG}}({\Theta }_{\mathrm{dG}},{\Phi }_{\mathrm{dG}})\le C_A{|||{\Theta }_{\mathrm{dG}}|||}_{\mathrm{dG}}{|||{\Phi }_{\mathrm{dG}}|||}_{\mathrm{dG}},$$

$$\alpha {|||{\Phi }_{\mathrm{dG}}|||}_{\mathrm{dG}}^2\le A_{\mathrm{dG}}({\Phi }_{\mathrm{dG}},{\Phi }_{\mathrm{dG}})\text{ for all }{\Phi }_{\mathrm{dG}}\in {\mathbf{V}}_h.$$

$$B_{\mathrm{dG}}(\Xi ,\mathbf{\eta },\Theta ,\Phi )\lesssim {ϵ}^{-2}{|||\Xi |||}_{\mathrm{dG}}{|||\mathbf{\eta }|||}_{\mathrm{dG}}{|||\Theta |||}_{\mathrm{dG}}{|||\Phi |||}_{\mathrm{dG}}\text{ and, }{C}_{\mathrm{dG}}(\Theta ,\Phi )\lesssim {ϵ}^{-2}{\Vert \Theta \Vert }_{\mathrm{dG}}{\Vert \Phi \Vert }_{\mathrm{dG}},$$

$$B_{\mathrm{dG}}(\Xi ,\mathbf{\eta },\Theta ,\Phi )\lesssim {ϵ}^{-2}{|||\Xi |||}_2{|||\mathbf{\eta }|||}_2{|||\Theta |||}_0{|||\Phi |||}_0,$$

$$b_{\text{dG}}(\xi ,\eta ,\theta ,\varphi )\le C_{ϵ}{\Vert \xi \Vert }_{0,4}{\Vert \eta \Vert }_{0,4}{\Vert \theta \Vert }_{0,4}{\Vert \varphi \Vert }_{0,4}\lesssim {ϵ}^{-2}{\Vert \xi \Vert }_{\text{dG}}{\Vert \eta \Vert }_{\text{dG}}{\Vert \theta \Vert }_{\text{dG}}{\Vert \varphi \Vert }_{\text{dG}}$$

$$\text{ and, }{c}_{\text{dG}}(\theta ,\varphi )\le C_{ϵ}{\Vert \theta \Vert }_0{\Vert \varphi \Vert }_0\lesssim {ϵ}^{-2}{\Vert \theta \Vert }_{\text{dG}}{\Vert \varphi \Vert }_{\text{dG}}.$$

$$\sum _{T\in \mathcal{T}}{\int }_T(\Xi \cdot \mathbf{\eta })(\Theta \cdot \Phi )\mathrm{dx}\lesssim {|||\Xi |||}_{\mathrm{dG}}{|||\mathbf{\eta }|||}_{\mathrm{dG}}{|||\Theta |||}_{\mathrm{dG}}{|||\Phi |||}_{\mathrm{dG}}\text{ for all }\Xi ,\mathbf{\eta },\Theta ,\Phi \in {\mathbf{V}}_h.$$

$${\Vert \mathrm{v}-{\mathrm{I}}_{\mathrm{dG}}\mathrm{v}\Vert }_{H^l(T)}\le C_I{h}_T^{s-l}{\Vert \mathrm{v}\Vert }_{H^s(T)},$$

Full text

Keywords: nematic liquid crystals, energy optimization, Landau-de Gennes energy functional, discontinuous Galerkin finite element methods, error analysis, convergence rate

1 Introduction

Liquid crystals are a transitional phase of matter between the liquid and crystalline phases. They inherit versatile properties of both liquid phase (e.g. fluidity) and, crystalline phase (optical, electrical, magnetic anisotropy) which make them ubiquitous in practical applications ranging from wristwatch and computer displays, nanoparticle organizations to proteins and cell membranes. Liquid crystals present different phases as temperature varies. We consider the nematic phase formed by rod-like molecules that self-assemble into an ordered structure, such that the molecules tend to align along a preferred orientation. There are three main mathematical models for nematic liquid crystals studied in literature which are the Oseen-Frank model [35, 6, 2], Ericksen model [19] and Landau-De Gennes model [17, 39, 16].

We present some rigorous results for the discontinuous Galerkin (dG) formulation of a reduced two-dimensional Landau-de Gennes model, following the model problem studied in [37, 42]. In the reduced Landau-de Gennes theory, the state of the nematic liquid crystal is described by a tensor order parameter, $\displaystyle\textit{Q}:=s(2\mathbf{n}\otimes\mathbf{n}-\textit{I}),$ where $\textit{Q}\in\mathbf{S}_{0}:=\{\textit{Q}=(\textit{Q}_{ij})_{1\leq i,j\leq 2}\in\mathbb{R}^{2\times 2}:\textit{Q}=\textit{Q}^{T},\text{tr}\textit{Q}=0\}$ and $I$ is the $2\times 2$ identity matrix. We refer to ${}^{\prime}\mathbf{n}^{\prime}$ as the director or the locally preferred in-plane alignment direction of the nematic molecules; ${}^{\prime}s^{\prime}$ is the scalar order parameter that measures the degree of order about ${}^{\prime}\mathbf{n}^{\prime}.$ We can write ${}^{\prime}\mathbf{n}^{\prime}$ as $\mathbf{n}=(\cos\theta,\sin\theta)$, where $\theta$ is the director angle in the plane. The general Landau-de Gennes Q-tensor order parameter is a symmetric traceless $3\times 3$ matrix and we can rigorously justify our $2$-D approach in certain model situations (see [43]).

Since Q is a symmetric, traceless $2\times 2$ matrix, we can represent Q as a two dimensional quantity, $\Psi=(u,v)$ where $u=s\cos 2\theta$ and $v=s\sin 2\theta$. The stable nematic equilibria are minimizers of the Landau-de Gennes energy functional given by $\mathcal{E}(\textit{Q}):=\mathcal{E}_{B}(\textit{Q})+\mathcal{E}_{E}(\textit{Q})+\mathcal{E}_{S}(\textit{Q})-\mathcal{E}_{L}(\textit{Q})$ subject to appropriate boundary conditions. The bulk energy over the domain $\Omega$ is given by $\mathcal{E}_{B}(\textit{Q}):=\int_{\Omega}(-\frac{\alpha^{2}}{2}\text{tr}{\textit{Q}}^{2}-\frac{b^{2}}{3}\text{tr}{\textit{Q}}^{3}+\frac{c^{2}}{4}(\text{tr}{\textit{Q}}^{2})^{2})\,\textrm{dx},$ $\alpha$, $b$ and $c$ being temperature and material dependent constants. The one-constant elastic energy is $\mathcal{E}_{E}(\textit{Q}):=\int_{\Omega}\frac{L_{el}}{2}\left\lvert\nabla\textit{Q}\right\rvert^{2}\,\textrm{dx}$ with $L_{el}>0$ being an elastic constant. The surface anchoring energy is $\mathcal{E}_{S}(\textit{Q}):=\int_{\partial\Omega}W\left\lvert(\textit{Q}_{11},\textit{Q}_{12})-\mathbf{g}\right\rvert^{2}\,\textrm{ds}$ with the anchoring strength $W$ on $\partial\Omega$ and prescribed Lipschitz continuous boundary function $\mathbf{g}:\partial\Omega\rightarrow\mathbb{R}^{2}$. The electrostatic energy $\mathcal{E}_{L}(\textit{Q}):=\int_{\Omega}-C_{0}(\textit{Q}\textit{E})\cdot\textit{E}\,\textrm{dx}$ with $C_{0}$ depending on vaccum permittivity and material dependent constants and E being the electric field vector. Here ’$\cdot$’ denotes the scalar product of vectors of function entries.

In the absence of surface effects ($\mathcal{E}_{S}(\textit{Q})$) and external fields ($\mathcal{E}_{L}(\textit{Q})$), the Landau-de Gennes energy functional in the dimensionless form considered in this article is given by [37]

[TABLE]

where $\Psi_{\epsilon}=(u,v)\in\mathbf{H}^{1}(\Omega):=H^{1}(\Omega)\times H^{1}(\Omega)$ and $\epsilon$ is a small positive parameter that depends on the elastic constant, bulk energy parameters and the size of the domain. More precisely, (see [37]) after rescaling and suitable non-dimensionalization, $\epsilon=\sqrt{\frac{L_{el}}{L^{2}c^{2}}}$ is proportional to the ratio of the nematic correlation length to a characteristic domain size. The nematic correlation length (usually on the scale of nanometers) is typically related to defect core sizes. The $\epsilon\rightarrow 0$ limit describes the macroscopic limit where the domain size is much larger than the correlation length i.e. describes micron scale geometries or even larger geometries. Some typical values [37] of these physical parameters are $L_{el}=10^{-11}\,\text{N m}$, $L=8\times 10^{-5}\text{ m}$ and $c^{2}=1\times 10^{6}\,\text{N m}^{-1}$. We are concerned with the minimization of the functional $\mathcal{E}$ for the Lipschitz continuous boundary function $\mathbf{g}:\partial\Omega\rightarrow\mathbb{R}^{2}$. More precisely, the admissible space is $\mathbf{X}=\left\{\mathbf{w}\in\mathbf{H}^{1}(\Omega):\mathbf{w}=\mathbf{g}\,\,\text{on}\,\,\partial\Omega\right\}.$ The strong formulation of (1.1) seeks $\Psi_{\epsilon}\in\mathbf{H}^{1}(\Omega)$ such that

[TABLE]

The behavior of $\Psi_{\epsilon}$ as $\epsilon\rightarrow 0$ has been studied by Brezis et al in [9] with the assumption $0<\epsilon<1$ and an $\epsilon$ independent bound for ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\epsilon}|\kern-1.07639pt|\kern-1.07639pt|}_{\mathbf{H}^{2}}$ has been established. The minimizers of the energy functional (1.1) for a given smooth boundary data $\mathbf{g}$ with $\left\lvert\mathbf{g}\right\rvert=1\text{ on }\partial\Omega$ for a smooth bounded domain has been well-studied [9]. Further, the authors rigorously prove in [9] that $\Psi_{\epsilon}\rightarrow\Psi_{0}$ as $\epsilon\rightarrow 0$ in $C^{1,\alpha}(\bar{\Omega})$   $\text{ for all }\alpha<1$, where $\Psi_{0}$ is a harmonic map and is a solution of $-\Delta\Psi_{0}=\Psi_{0}\left\lvert\nabla\Psi_{0}\right\rvert^{2}\text{ on }\Omega,$ $\left\lvert\Psi_{0}\right\rvert=1$ on $\Omega,$   $\Psi_{0}=\mathbf{g}$ on $\partial\Omega$. We are interested in the dG finite element approximation of regular solutions of the boundary value problem (1.2). Let ${\rm N}:\mathbf{H}^{1}(\Omega)\rightarrow\mathbf{H}^{1}(\Omega)^{*}$ be the Ginzburg-Landau operator defined as ${\rm N}(\Psi_{\epsilon}):=-\Delta\Psi_{\epsilon}+2\epsilon^{-2}(\left\lvert\Psi_{\epsilon}\right\rvert^{2}-1)\Psi_{\epsilon}$. For $\epsilon$ small enough, the linearized operator $\displaystyle L_{\epsilon}:={\rm DN}({\Psi_{\epsilon}})$ is bijective when defined between standard spaces (e.g., $L_{\epsilon}:\mathbf{W}^{2,2}\rightarrow\mathbf{L}^{2}$) [40], although the norm of its inverse blows up as $\epsilon\rightarrow 0$.

This model problem has been studied using the conforming finite element method for (1.2) for a *fixed * $\epsilon$ in [37] on a square domain with $\mathbf{g}$ consistent tangent boundary conditions. The variational formulation of a more generic Landau-de Gennes model with homogeneous Dirichlet boundary data has been studied [16] in which an abstract approach of the finite element approximation of nonsingular solution branches has been analyzed in the conforming finite element set up. However the analysis for the non-homogeneous boundary conditions has not been considered in this work. In [2] and [36], the authors have discussed a mixed finite element method for the Frank-Oseen and Ericksen-Leslie models for nematic liquid crystals, respectively.

It is well-known that Allen-Cahn equation [1] is the gradient flow equation associated with the Lyapunov energy functional $\mathcal{J}_{\epsilon}(u)=\frac{1}{2}\int_{\Omega}(\left\lvert\nabla u\right\rvert^{2}+\frac{\epsilon^{-2}}{2}(\left\lvert u\right\rvert^{2}-1)^{2})\,\textrm{dx},$ where $u:\Omega\rightarrow R$ is a scalar valued function and $\epsilon$ is a small parameter (not to be confused with the $\epsilon$ in (1.2)) known as an “interaction length” which is small compared to the characteristic dimensions on the laboratory scale. Numerical approximations of the Allen-Cahn equation have been extensively investigated in the literature. A priori error estimate for the error bounds as a function of $\epsilon$ have been analyzed and shown to be of polynomial order in $\epsilon^{-1}$ by Feng and Prohl in [22], for the conforming finite element approximation of the Allen-Cahn problem. A symmetric interior penalty discontinuous Galerkin method [21] and a posteriori error analysis [23, 7] (which has a low order polynomial dependence on $\epsilon^{-1}$) have also been studied for Allen-Cahn equation. A dG scheme has been proposed recently in [4] for the $\epsilon$- dependent stochastic Allen–Cahn equation with mild space-time noise posed on a bounded domain of $\mathbb{R}^{2}$. The Allen Cahn work is relevant to time-dependent front propagation in nematic liquid crystals, in certain reduced symmetric situations [38]. The problem considered in this article is different from the Allen Cahn equation; we have a system of two coupled nonlinear partial differential equations (PDEs) in a time independent scenario.

Our motivation comes from the planar bistable nematic device reported in [37]. This device consists of a periodic array of shallow square or rectangular wells, filled with nematic liquid crystals, subject to tangent boundary conditions on the lateral surfaces. The vertical well height is much smaller than the cross-sectional dimensions and hence, it is reasonable to assume invariance in the vertical direction and to model the profile on the bottom cross-section, taking the domain to be a square as opposed to a three-dimensional square well. This model reduction can be rigorously justified using gamma-convergence techniques [25, 43]. The tangent boundary conditions require that the nematic director, identified with the leading eigenvector $\mathbf{n}$ of the Landau-de Gennes Q-tensor order parameter, lies in the plane of the square and $\mathbf{n}$ is tangent to the square edges. Indeed, this motivates the choice of the Dirichlet conditions for Q on the square edges as described below. The tangent boundary conditions can also be phrased in terms of surface anchoring energies but this results in mixed boundary-value problems for Q which are relatively more difficult to analyse than the Dirichlet counterparts. Our results can be applied to the planar bistable nematic device; we can model a single well as a square domain with Dirichlet tangent boundary conditions on the square edges and analyze discontinuous Galerkin finite element methods (dGFEMs) to approximate regular solutions of (1.2) for a fixed $\epsilon$. This involves a semilinear system of PDEs with a cubic nonlinearity (see (1.2)) and non-homogeneous boundary conditions. As reported in [37, 42], there are six experimentally observed stable nematic equilibria, labeled as diagonal and rotated states, for this model problem. The nematic director roughly aligns along one of the square diagonals in the diagonal states whereas the director rotates by $\pi$ radians between a pair of opposite parallel edges, in a rotated state. There are two diagonal states, since there are two square diagonals, and four rotated states related to each other by a $\frac{\pi}{2}$ rotation. The dGFEMs are attractive because they are element-wise conservative, are flexible with respect to local mesh adaptivity, are easier to implement than finite volume schemes, allow for non-uniform degrees of approximations for solutions with variable regularity over the computational domain and can handle non-homogeneous boundary condition in a natural way. These methods also relax the inter-element continuity requirement in conforming FEM. An a priori error analysis of dGFEMs for general elliptic problems has been derived in [41, 28, 27]. For a comprehensive study of several dGFEMs applied to elliptic problems, see [5]. The dGFEMs are also well studied for fourth order elliptic problems [13, 24]. Recently, dGFEMs have been studied for the von Kármán equations [13] that involves a quadratic non-linearity and homogeneous boundary conditions.

To the best of our knowledge, dGFEMs have not been analysed for the nonlinear system derived from the reduced two-dimensional Landau-de Gennes energy in (1.1) and this is the primary motivation for our study. Our contributions can be summarized as follows.

We derive an elegant representation of the nonlinear operator, convergence analysis with $h$-$\epsilon$* dependency* and a priori error estimate of a dGFEM formulation for the system (1.2) with non homogeneous boundary conditions. The choice of the discretization parameter $h$ that depends on $\epsilon$ ensures that (i) the dG discretization of a linearized problem is well-posed and (ii) the corresponding discrete non-linear problem has a unique solution following applications of contraction mapping theorem. 2. 2.

We prove a best approximation result for regular solutions of the non-linear problem (1.2). 3. 3.

We prove the quadratic convergence of the Newton’s iterates to the approximate dGFEM solution. 4. 4.

Our numerical results confirm the theoretical orders of convergence and rate of convergence as a function of $h$ and $\epsilon$, in the context of the planar bistable nematic device and other representative examples.

Throughout the paper, standard notations on Sobolev spaces and their norms are employed. The standard semi-norm and norm on $H^{s}(\Omega)$ $(\text{resp.}\,W^{s,p}(\Omega))$ for $s,p$ positive real numbers, are denoted by $\left\lvert\cdot\right\rvert_{s}$ and $\left\lVert\cdot\right\rVert_{s}$ $(\text{resp. }\left\lvert\cdot\right\rvert_{s,p}\text{ and }\left\lVert\cdot\right\rVert_{s,p})$. The standard $L^{2}(\Omega)$ inner product is denoted by $(\cdot,\cdot)$. We use the notation $\mathbf{H}^{s}(\Omega)$ (resp. $\mathbf{L}^{p}(\Omega)$) to denote the product space $H^{s}(\Omega)\times H^{s}(\Omega)$ $(\text{resp. }L^{p}(\Omega)\times L^{p}(\Omega))$. The standard norms ${|\kern-1.07639pt|\kern-1.07639pt|\cdot|\kern-1.07639pt|\kern-1.07639pt|}_{s}$ ($\text{resp. }{|\kern-1.07639pt|\kern-1.07639pt|\cdot|\kern-1.07639pt|\kern-1.07639pt|}_{s,p}$) in the Sobolev spaces $\mathbf{H}^{s}(\Omega)$ ($\text{resp. }\mathbf{W}^{s,p}(\Omega)$) defined by ${|\kern-1.07639pt|\kern-1.07639pt|\Phi|\kern-1.07639pt|\kern-1.07639pt|}_{s}\!=(\left\lVert\phi_{1}\right\rVert_{s}^{2}+\left\lVert\phi_{2}\right\rVert_{s}^{2})^{\frac{1}{2}}$ for all $\Phi\!=\!(\phi_{1},\phi_{2})\in\!\mathbf{H}^{s}(\Omega)$ $(\text{resp. }{|\kern-1.07639pt|\kern-1.07639pt|\Phi|\kern-1.07639pt|\kern-1.07639pt|}_{s,p}\!\!=(\left\lVert\phi_{1}\right\rVert_{s,p}^{2}+\left\lVert\phi_{2}\right\rVert_{s,p}^{2})^{\frac{1}{2}}$ for all $\Phi\!=\!(\phi_{1},\phi_{2})\!\in\!\mathbf{W}^{s,p}(\Omega)).$ The norm on $\mathbf{L}^{2}(\Omega)$ space is defined by ${|\kern-1.07639pt|\kern-1.07639pt|\Phi|\kern-1.07639pt|\kern-1.07639pt|}_{0}\!=(\left\lVert\phi_{1}\right\rVert_{0}^{2}+\left\lVert\phi_{2}\right\rVert_{0}^{2})^{\frac{1}{2}}$ for all $\Phi\!=\!(\phi_{1},\phi_{2})\in\!\mathbf{L}^{2}(\Omega).$ Set $\textrm{V}:=H_{0}^{1}(\Omega)=\left\{\phi\in L^{2}(\Omega):\frac{\partial\phi}{\partial x},\frac{\partial\phi}{\partial y}\in L^{2}(\Omega)\,\,,\phi|_{\partial{\Omega}}=0\right\}$ and $\mathbf{V}\!=\,\!\mathbf{H}_{0}^{1}(\Omega)=H_{0}^{1}(\Omega)\times H_{0}^{1}(\Omega)$. The inequality $a\lesssim b$ abbreviates $a\leq Cb$ with the constant $C>0$ independent of mesh-size parameter ${}^{\prime}h^{\prime}$ and ${}^{\prime}\epsilon^{\prime}$. The constants that appear in various Sobolev imbedding results in the sequel are denoted using a generic notation $C_{S}$.

The paper is organised as follows. In Section 2, we present the model problem along with the weak formulation and some preliminary results. We state the dG finite element formulation of the problem in Subsection 3.1 and our main results are stated in Subsection 3.2. The existence and uniqueness of the discrete solution of the non-linear problem, error estimates, best approximation result and the convergence of Newton’s method are presented as main theorems. Section 4 contains some auxiliary results needed to prove the main results. A discrete inf-sup condition for a discrete bilinear form has been established. In Section 5, we prove the main theorems. A contraction map has been defined on the discrete space to use a fixed point argument for proving the existence and uniqueness of discrete solution. An alternative proof of the existence and uniqueness of the solution of the discrete problem using Newton-Kantorovich theorem has been given in this section. This is followed by a proof of the quadratic convergence of Newton’s method and numerical results that are consistent with the theoretical results in Section 6.

2 Preliminaries

In this section, we introduce the weak formulation for (1.2) and establish some boundedness results. The details of derivation of the weak formulation in presented in A.1.

In the weak formulation of (1.2), we seek $\Psi_{\epsilon}\in\mathbf{X}$ such that

[TABLE]

where for all $\Xi=(\xi_{1},\xi_{2}),\boldsymbol{\eta}=(\eta_{1},\eta_{2}),\Theta=(\theta_{1},\theta_{2}),\Phi=(\phi_{1},\phi_{2})\in\mathbf{H}^{1}(\Omega)$ and $C_{\epsilon}=2\epsilon^{-2},$

[TABLE]

Remark 2.1*.*

When $\Xi=\boldsymbol{\eta}=\Theta=\Psi_{\epsilon}:=(u,v)$,

[TABLE]

The operator $B(\cdot,\cdot,\cdot,\cdot)$ corresponds to the non-linear part of the system (1.2). Such representation of $B(\cdot,\cdot,\cdot,\cdot)$ yields nice properties proven in Lemma 2.3 which makes the analysis elegant.

Lemma 2.2**.**

(Poincaré inequality)**[33]** Let $\Omega$ be a bounded open Lipschitz domain in $\mathbb{R}^{2}.$ Then there exists a positive constant $\alpha_{0}=\alpha_{0}(\Omega)$ such that

[TABLE]

The boundedness and coercivity of $A(\cdot,\cdot)$ and the boundedness of $C(\cdot,\cdot)$ given below can be easily verified. For all $\Theta$, $\Phi\in\mathbf{V}$,

[TABLE]

where $C_{\alpha_{0}}$ depends on $\alpha_{0}.$ The next lemma establishes two boundedness results for $B(\cdot,\cdot,\cdot,\cdot)$.

Lemma 2.3**.**

(Boundedness of $B(\cdot,\cdot,\cdot,\cdot)$ )* For all $\Xi$, $\boldsymbol{\eta}$, $\Theta$, $\Phi\in\mathbf{V}$,*

[TABLE]

and for all $\Xi,\boldsymbol{\eta}\in\mathbf{H}^{2}(\Omega)$, $\Theta$, $\Phi\in\mathbf{V}$,

[TABLE]

Proof.

It is enough to prove the results for $b(\cdot,\cdot,\cdot,\cdot)$; then (2.6) and (2.7) follow from the definition of $B(\cdot,\cdot,\cdot,\cdot)$ and a grouping of the terms. For $\xi,\eta,\theta,\phi\in H^{1}(\Omega)$, a use of Hölder’s inequality and the Sobolev imbedding result ${H}^{1}(\Omega)\hookrightarrow L^{4}(\Omega)$ [20] leads to (2.6) as

[TABLE]

For $\xi$, $\eta\in{H}^{2}(\Omega)$, a use of the Sobolev imbedding result ${H}^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$ and the Cauchy-Schwarz inequality leads to

[TABLE]

where $"\lesssim"$ in the last two inequalities above absorbs $C_{S}$ and the constant from $C_{\epsilon}$. ∎

The existence of minimizers of (1.1) follows from the coercivity of $\mathcal{E}$ and its convexity in $\nabla\Psi_{\epsilon}$ in (1.1) and this implies the existence of a solution of the non-linear system in (2.1). The regularity result in the next lemma follows from arguments in [16, 26] and in detailed in A.2.

Lemma 2.4** (Regularity result).**

Let $\Omega$ be an open, bounded, Lipschitz and convex domain of $\mathbb{R}^{2}.$ Then for $\mathbf{g}\in\mathbf{H}^{\frac{3}{2}}(\partial\Omega),$ any solution of (1.2) belongs to $\mathbf{H}^{2}(\Omega)$.

In this article, we approximate regular solutions [32] $\Psi_{\epsilon}$ of (1.2) for a given $\epsilon$. The regularity of solution implies that the linearized operator $\langle{\rm DN}(\Psi_{\epsilon})\cdot,\cdot\rangle$ is invertible in the Banach space and is equivalent to the following inf-sup condition [18]

[TABLE]

where $\langle{\rm DN}(\Psi_{\epsilon})\Theta,\Phi\rangle:=A(\Theta,\Phi)+3B(\Psi_{\epsilon},\Psi_{\epsilon},\Theta,\Phi)+C(\Theta,\Phi)$ and the inf-sup constant $\beta$ depends on $\epsilon$. From now onwards, the subscript $\epsilon$ in $\Psi_{\epsilon}$ is suppressed in the sequel for notational brevity.

3 Discrete formulation

In this section, we derive the dGFEM formulation for (1.2) and state our main results.

3.1 The dGFEM formulation

Let $\mathcal{T}$ be a triangulation [14] of $\bar{\Omega}$ into triangles and let the discretization parameter $h$ associated with the partition $\mathcal{T}$ be defined as $h=\max_{T\in\mathcal{T}}h_{T},$ where $h_{T}=diam(T)$. Let $\mathcal{E}_{i}$( resp. $\mathcal{E}_{D}$) denote the interior (resp. boundary) edges of $\mathcal{T}$ and $\mathcal{E}:=\mathcal{E}_{i}\cup\mathcal{E}_{D}$. Also, the boundary of an element $T$ is denoted by $\partial T$ and the unit normal vector outward from $T$ is denoted by $n.$ For any interior edge $E$ shared by two triangles $T^{+}$ and $T^{-}$, let the unit normal pointing from $T^{+}$ to $T^{-}$ be $n_{+}$ [see Figure 1]. We assume the triangulation $\mathcal{T}$ be shape regular [14] in the sense that there exists $\rho>0$ such that if $h_{T}$ is the diameter of $T$, then $T$ contains a ball of radius $\rho h_{T}$ in its interior.

For a positive integer $s$, define the broken Sobolev space by $H^{s}(\mathcal{T}):=\{v\in L^{2}(\Omega):v|_{T}\in H^{s}(T)\,\text{ for all }T\in\mathcal{T}\}$ equipped with the broken norm $\left\lVert v\right\rVert_{s,\mathcal{T}}^{2}=\sum_{T\in\mathcal{T}}\left\lVert v\right\rVert_{s,T}^{2}.$ Denote $\mathbf{H}^{s}(\mathcal{T})=H^{s}(\mathcal{T})\times H^{s}(\mathcal{T})$ to be the product space with the norm ${|\kern-1.07639pt|\kern-1.07639pt|\Phi|\kern-1.07639pt|\kern-1.07639pt|}_{s,\mathcal{T}}^{2}=\left\lVert\phi_{1}\right\rVert_{s,\mathcal{T}}^{2}+\left\lVert\phi_{2}\right\rVert_{s,\mathcal{T}}^{2}\,\,\text{ for all }\Phi=(\phi_{1},\phi_{2})\in\mathbf{H}^{s}(\mathcal{T}).$ For $u\in H^{s}(\mathcal{T})$, we define the broken gradient $\nabla_{\textrm{dG}}u$ by $(\nabla_{\textrm{dG}}u)|_{T}=\nabla_{\textrm{dG}}(u|_{T})\,\,\text{ for all }T\in\mathcal{T}.$ Follow the standard convention [41] for the jump and average. Consider the finite dimensional space that consists of piecewise linear polynomials defined by $\displaystyle V_{h}=\{v\in L^{2}(\Omega):v|_{T}\in P_{1}(T)\,\,\text{for all }T\in\mathcal{T}\},$ and equip it with the mesh dependent norm defined by $\displaystyle\left\lVert v\right\rVert^{2}_{\textrm{dG}}=\left\lvert v\right\rvert^{2}_{\textrm{H}^{1}(\mathcal{T})}+J^{\sigma}(v,v)=\sum_{T\in\mathcal{T}}\int_{T}\left\lvert\nabla v\right\rvert^{2}\,\textrm{dx}+\sum_{E\in\mathcal{E}}\int_{E}\frac{\sigma}{h}[v]^{2}\,\textrm{ds},$ where $\sigma>0$ is the penalty parameter. Let $\mathbf{V}_{h}\!:=V_{h}\times V_{h}$ be equipped with the product norm defined by ${|\kern-1.07639pt|\kern-1.07639pt|\Phi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}^{2}=\left\lVert\phi_{1}\right\rVert^{2}_{\textrm{dG}}+\left\lVert\phi_{2}\right\rVert^{2}_{\textrm{dG}}\text{ for all }\Phi_{\textrm{dG}}=(\phi_{1},\phi_{2})\in\mathbf{V}_{h}.$

The discontinuous Galerkin formulation corresponding to (1.2) seeks $\Psi_{\textrm{dG}}\!\in\!\mathbf{V}_{h}$ such that for all $\Phi_{\textrm{dG}}\in\mathbf{V}_{h},$

[TABLE]

where for $\Xi=(\xi_{1},\xi_{2}),\,\boldsymbol{\eta}=(\eta_{1},\eta_{2}),\,\Theta=(\theta_{1},\theta_{2}),\,\Phi=(\phi_{1},\phi_{2})\in\mathbf{H}^{1}(\mathcal{T})$, $\displaystyle L_{\textrm{dG}}(\Phi_{\textrm{dG}})=l^{1}_{\textrm{dG}}(\phi_{1})+l^{2}_{\textrm{dG}}(\phi_{2}),$

[TABLE]

and for $\theta,\phi\in H^{1}(\mathcal{T})$, $\mathbf{g}=(g_{1},g_{2})$, $\lambda\in[-1,1],$ $\displaystyle a_{\textrm{dG}}(\theta,\phi):=a_{h}(\theta,\phi)-J(\theta,\phi)+\lambda J(\phi,\theta)+J^{\sigma}(\theta,\phi),$

[TABLE]

Remark 3.1*.*

The parameter values $\lambda=-1,0,1$ corresponds to symmetric interior penalty, incomplete interior penalty and non-symmetric interior penalty dG methods, respectively in the context of linear problems.

3.2 The main results

Our main results are given below. The proofs are presented in Sections 5 and 6. The details of the suppressed constants in ’$\lesssim$’ in the Theorems 3.2 - 3.5 will be made clear in the proofs given in Section 5.

Theorem 3.2** (Existence, uniqueness and dG norm error estimate).**

Let $\Psi$ be a regular solution of the non-linear system (2.1). For a given fixed $\epsilon>0,$ sufficiently large $\sigma$ and sufficiently small discretization parameter chosen as $h=O(\epsilon^{2+\alpha})$ for any $\alpha>0$, there exists a unique solution $\Psi_{\rm dG}$ of the discrete non-linear problem (3.1) that approximates $\Psi$ such that

[TABLE]

Theorem 3.3** (Best approximation result).**

Let $\Psi$ be a regular solution of the non-linear system (2.1). For a given fixed $\epsilon>0,$ sufficiently large $\sigma$ and sufficiently small discretization parameter chosen as $h=O(\epsilon^{2+\alpha})$ with $\alpha>0$, the unique discrete solution $\Psi_{\rm{dG}}$ of (3.1) that approximates $\Psi$ satisfies the best-approximation property

[TABLE]

Remark 3.4*.*

Theorem 3.3 shows that discrete solution $\Psi_{\rm dG}$ is "the best" approximation of $\Psi$ in $\mathbf{V}_{h}$, up to a constant. Best approximation result is mainly motivated by "the best" approximation of discrete solution of linear PDEs in Céa’s lemma [14] for conforming FEM.

We establish the $\mathbf{L}^{2}$ norm error estimate when the parameter $\lambda$ that appears in the discrete non-linear system (3.1) (in the term $A_{\rm dG}(\cdot,\cdot)$ through $a_{\textrm{dG}}(\cdot,\cdot)$) takes the value $-1$.

Theorem 3.5** ($\mathbf{L}^{2}$ norm error estimate).**

Let $\Psi$ be a regular solution of the non-linear system (2.1). For a given fixed $\epsilon>0,$ sufficiently large $\sigma$ and sufficiently small discretization parameter chosen as $h=O(\epsilon^{2+\alpha})$ for $\alpha>0$, there exists a unique solution $\Psi_{\rm dG}$ of the discrete non-linear problem (3.1) that approximates $\Psi$ such that

[TABLE]

We use Newton’s method [32] for computation of discrete solutions. It is a standard and very effective root-finding method to approximate the roots of non-linear system of PDEs.

Theorem 3.6** (Convergence of Newton’s method).**

Let $\Psi$ be a regular solution of the non-linear system (2.1) and let $\Psi_{\rm{dG}}$ solve (3.1). For a given fixed $\epsilon>0,$ sufficiently large $\sigma$ and sufficiently small discretization parameter chosen as $h=O(\epsilon^{2+\alpha})$ with $\alpha>0$, there exists $\rho_{1}>0$, independent of $h$, such that for any initial guess $\Psi^{0}_{\rm dG}$ with ${|\kern-1.07639pt|\kern-1.07639pt|\Psi^{0}_{\rm dG}-\Psi_{\rm dG}|\kern-1.07639pt|\kern-1.07639pt|}_{\rm dG}\leq\rho_{1}$, it follows ${|\kern-1.07639pt|\kern-1.07639pt|\Psi^{n}_{\rm dG}-\Psi_{\rm dG}|\kern-1.07639pt|\kern-1.07639pt|}_{\rm dG}\leq\frac{\rho_{1}}{2^{n}}$ for all $n=1,2,\dots$ and the iterates $\Psi^{n}_{\rm dG}$ of Newton’s method are well-defined and converges quadratically to $\Psi_{\rm dG};$ that is, ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\rm dG}^{n}-\Psi_{\rm dG}|\kern-1.07639pt|\kern-1.07639pt|}_{\rm dG}\leq C_{q}{|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\rm dG}^{n-1}-\Psi_{\rm dG}|\kern-1.07639pt|\kern-1.07639pt|}_{\rm dG}^{2},$ where $C_{q}$ is a constant independent of $h$.

4 Auxiliary results

This section presents some auxiliary results needed to establish the main results in Subsection 3.2. The boundedness and ellipticity results for $A_{\textrm{dG}}(\cdot,\cdot)$, and the boundedness result for $B_{\textrm{dG}}(\cdot,\cdot,\cdot,\cdot)$, $C_{\textrm{dG}}(\cdot,\cdot)$ are proved and the inf-sup conditions for a discrete linearized bilinear form and a perturbed bilinear form are established.

Lemma 4.1**.**

(Poincaré type inequality)**[34, 11]** For $\phi\in H^{1}(\mathcal{T})$, there exists a constant $C_{P}>0$ independent of $h$ and $\phi$ such that for $1\leq r<\infty$, $\left\lVert\phi\right\rVert_{L^{r}(\Omega)}\leq C_{P}\left\lVert\phi\right\rVert_{\rm dG}.$

For boundedness and coercivity results of $A_{\textrm{dG}}(\cdot,\cdot)$ and the boundedness results of $B_{\textrm{dG}}(\cdot,\cdot,\cdot,\cdot)$ and $C_{\textrm{dG}}(\cdot,\cdot)$, it is enough to prove the corresponding results for $a_{\textrm{dG}}(\cdot,\cdot)$, $b_{\textrm{dG}}(\cdot,\cdot,\cdot,\cdot)$ and $c_{\textrm{dG}}(\cdot,\cdot)$, respectively.

Lemma 4.2**.**

(Boundedness and coercivity of $A_{\rm dG}(\cdot,\cdot)$)* [41] For $\Theta_{\rm dG},\Phi_{\rm dG}\in\mathbf{V}_{h}$,*

[TABLE]

where $C_{A}$ depends on the penalty parameter $\sigma$ and the constant from trace inequality. For a sufficiently large parameter $\sigma$, there exists a positive constant $\alpha>0$ such that

[TABLE]

Lemma 4.3**.**

(Boundedness of $B_{\textrm{dG}}(\cdot,\cdot,\cdot,\cdot)$ and $C_{\textrm{dG}}(\cdot,\cdot)$)* For $\Xi,\boldsymbol{\eta},\Theta,\Phi\in\mathbf{H}^{1}(\mathcal{T})$, it holds that*

[TABLE]

and for all $\Xi,\boldsymbol{\eta}\in\mathbf{H}^{2}(\Omega)$, $\Theta$, $\Phi\in\mathbf{H}^{1}(\mathcal{T})$,

[TABLE]

where the hidden constant in $"\lesssim"$ depends on the constants from $C_{\epsilon}$, $C_{P}$ and $C_{S}$.

Proof.

We prove the boundedness results for $b_{\textrm{dG}}(\cdot,\cdot,\cdot,\cdot)$ and $c_{\textrm{dG}}(\cdot,\cdot)$. Then a use of definitions of $B_{\textrm{dG}}(\cdot,\cdot,\cdot,\cdot)$ and $C_{\textrm{dG}}(\cdot,\cdot)$, discrete Cauchy-Schwarz inequality and a grouping of the terms yields the required results. For $\xi,\,\eta,\,\theta$ and $\phi\in H^{1}(\mathcal{T})$, a use of Holder’s inequality and Cauchy-Schwarz inequality along with Lemma 4.1 leads to

[TABLE]

The proof of (4.2) follows analogously to that of the proof of (2.7) with a use of the imbedding result $H^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$ and Cauchy-Schwarz inequality. ∎

Remark 4.4*.*

A use of (4.3) leads to the following boundedness estimate

[TABLE]

The next two lemmas describe the estimates for interpolation and enrichment operators that are crucial for the error estimates.

Lemma 4.5**.**

(Interpolation estimate)**[12]** For ${\rm v}\in H^{s}(\Omega)\text{ with }s\geq 1$, there exists ${\rm{I}_{\rm dG}v}\in V_{h}$ such that for any $T\in\mathcal{T}$,

[TABLE]

for $l=0,1$ and some positive constant $C_{I}$ independent of $h$.

Lemma 4.6**.**

(Enrichment operator).* [11, 30] Let ${\rm E}_{h}:V_{h}\rightarrow V_{C}\subset H^{1}(\Omega)$ be an enrichment operator with $V_{C}$ being the Lagrange $P_{r}$ conforming finite element space associated with the triangulation $\mathcal{T}$. Then for any $\phi_{\rm dG}\in V_{h}$, the following results hold:*

[TABLE]

where $C_{en_{1}}$ and $C_{en_{2}}$ are constants independent of $h$.

For all $\Theta_{\textrm{dG}},\Phi_{\textrm{dG}}\in\mathbf{V}_{h}$, define the discrete bilinear form $\langle{\rm DN}_{h}(\Psi)\cdot,\cdot\rangle$ by

[TABLE]

Theorem 4.7**.**

Let $\Psi$ be a regular solution of the non-linear system (2.1). For a given fixed $\epsilon>0,$ a sufficiently large $\sigma$ and a sufficiently small discretization parameter chosen such that $h=O(\epsilon^{2})$, there exists a constant $\beta_{0}$ such that the following discrete inf-sup condition holds:

[TABLE]

Proof.

For $\Theta\in\mathbf{V}+\mathbf{V}_{h},$ Lemma 4.3, (2.5) and (2.7) yield $B_{\textrm{dG}}(\Psi,\Psi,\Theta,\cdot)$, $B(\Psi,\Psi,\Theta,\cdot),C_{\textrm{dG}}(\Theta,\cdot)$ and $C(\Theta,\cdot)\in\mathbf{L}^{2}(\Omega)$. Therefore, for a given $\Theta_{\textrm{dG}}\in\mathbf{V}_{h}$ with ${|\kern-1.07639pt|\kern-1.07639pt|\Theta_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}=1$, there exist $\boldsymbol{\xi}_{\epsilon}$ and $\boldsymbol{\eta}_{\epsilon}\in\mathbf{H}^{2}(\Omega)\cap\mathbf{V}$ that solve the linear systems

[TABLE]

A use of Lemmas 4.1, 4.3 and elliptic regularity leads to

[TABLE]

where the hidden constant in $"\lesssim"$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $C_{S}$ and $C_{P}$. Subtract (4.4a) from (4.4b), choose $\Phi=\boldsymbol{\eta}_{\epsilon}-\boldsymbol{\xi}_{\epsilon}$ and use (2.4) and (4.2) to obtain

[TABLE]

A use of Lemma 4.6 in (4.6) yields

[TABLE]

where the constant hidden in $"\lesssim"$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2},\alpha_{0}$ and $C_{en_{1}}$. Since $\Psi$ is regular solution of (2.1), a use of (2.8) yields that there exists $\Phi\in\mathbf{V}$ with ${|\kern-1.07639pt|\kern-1.07639pt|\Phi|\kern-1.07639pt|\kern-1.07639pt|}_{1}=1$ such that

[TABLE]

A use of (4.4b), (2.4), an introduction of intermediate terms and triangle inequality leads to

[TABLE]

Since $\displaystyle[\boldsymbol{\xi}_{\epsilon}]=0$ on $\mathcal{E}_{i}$ and $\boldsymbol{\xi}_{\epsilon}=0$ on $\mathcal{E}_{D}$, a use of second inequality in Lemma 4.6 and the triangle inequality yields

[TABLE]

Since $\Theta_{\textrm{dG}}+\textrm{I}_{\textrm{dG}}\boldsymbol{\xi}_{\epsilon}\in\mathbf{V}_{h},$ Lemma 4.2 implies that there exists $\Phi^{\epsilon}_{\textrm{dG}}\in\mathbf{V}_{h}$ with ${|\kern-1.07639pt|\kern-1.07639pt|\Phi^{\epsilon}_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}=1$ such that

[TABLE]

A use of Lemmas 4.3 and 4.6 yields the estimates for the second and third terms in (4.10) as

[TABLE]

A use of Lemmas 4.2, 4.5 and 4.6 leads to

[TABLE]

A combination of (4.11) and (4.12) in (4.10) leads to

[TABLE]

where $"\lesssim"$ includes a constant that depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2},C_{A},C_{S},C_{P},C_{en_{1}}$ and $C_{I}.$ A substitution of (4.13) in (4.9) and a use of Lemma 4.5 yields

[TABLE]

Moreover, a use of (4.13), (4.14), Lemma 4.5 and (4.7) in (4.8) yields

[TABLE]

A use of triangle inequality, (4.14), (4.15) and (4.5) leads to

[TABLE]

where the constant $C_{1}$ is independent of $h$ and $\epsilon$. Therefore, for a given $\epsilon$, the discrete inf-sup condition holds with $\beta_{0}=\frac{1}{2C_{1}}$ for $h<h_{0}:=\frac{\epsilon^{2}}{2C_{1}}$. ∎

We use the perturbed bilinear form defined as $\langle{\rm DN}_{h}(\textrm{I}_{\textrm{dG}}\Psi)\Theta_{\rm dG},\Phi_{\rm dG}\rangle:=A_{\rm dG}(\Theta_{\rm dG},\Phi_{\rm dG})+3B_{\rm dG}(\textrm{I}_{\textrm{dG}}\Psi,\textrm{I}_{\textrm{dG}}\Psi,\Theta_{\rm dG},\Phi_{\rm dG})+C_{\rm dG}(\Theta_{\rm dG},\Phi_{\rm dG})$ for all $\Theta_{\textrm{dG}},\Phi_{\textrm{dG}}\in\mathbf{V}_{h}$ in our analysis and Newton’s algorithm in our manuscript. The next lemma establishes discrete inf-sup condition for the perturbed bilinear form.

Lemma 4.8**.**

(Stability of perturbed bilinear form).* Let $\Psi$ be a regular solution of (2.1) and ${{\rm I}_{\rm dG}}\Psi$ be the interpolation of $\Psi$ from Lemma 4.5. For a given fixed $\epsilon>0,$ a sufficiently large $\sigma$ and a sufficiently small discretization parameter chosen as $h=O(\epsilon^{2})$, the perturbed bilinear form $\langle{\rm DN}_{h}({\rm I}_{\rm dG}\Psi)\cdot,\cdot\rangle$ satisfies the following discrete inf-sup condition*

[TABLE]

Proof.

For $\tilde{\Psi}=\Psi-{\rm I}_{\textrm{dG}}\Psi$,

[TABLE]

A use of (3.2) and Remark 4.4 lead to

[TABLE]

[TABLE]

A use of Theorem 4.7, (4.1) and Lemma 4.5 leads to

[TABLE]

where the constant $C_{2}$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $C_{P}$, $C_{S}$, and $C_{I}$. Hence, for a sufficiently small choice of $h<h_{2}:=\min(h_{0},h_{1})$ with $h_{1}:=\frac{\beta_{0}}{2C_{2}}\epsilon^{2}$ , the required result holds for $\frac{\beta_{0}}{2}.$ ∎

5 Proof of main results

The proofs of main results use some preliminary results that are established first. First the non-linear map $\mu^{\epsilon}_{\textrm{dG}}$ is defined on the discrete space $\mathbf{V}_{h}$ and in Theorem 5.1 it is established that $\mu^{\epsilon}_{\textrm{dG}}$ maps a closed convex set to itself and the fixed point of $\mu^{\epsilon}_{\textrm{dG}}$ is a solution of the discrete non-linear problem in (3.1) and vice-versa. For $\Phi_{\textrm{dG}}\in\mathbf{V}_{h}$, define the map, $\mu^{\epsilon}_{\textrm{dG}}:\mathbf{V}_{h}\rightarrow\mathbf{V}_{h}$ by

[TABLE]

The map $\mu^{\epsilon}_{\textrm{dG}}$ is well-defined and this follows from the inf-sup conditions of $\langle{\rm DN}_{h}({{\rm I}_{\rm dG}\Psi})\cdot,\cdot\rangle$ in Lemma 4.8. We use the abbreviation $\mu_{\textrm{dG}}:=\mu^{\epsilon}_{\textrm{dG}}$ in the rest of the article for notational convenience. Let $\mathbb{B}_{R}(\textrm{I}_{\textrm{dG}}\Psi):=\{\Phi_{\textrm{dG}}\in\mathbf{V}_{h}:{|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{\textrm{dG}}\Psi-\Phi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq R\}.$ The existence and uniqueness result of the discrete solution $\Psi_{\textrm{dG}}$ in Theorem 3.2 is an application of Brouwer’s fixed point theorem. Then, the energy norm error estimate is presented. As an alternative approach, the existence and uniqueness of discrete solution is established using Newton-Kantorovich theorem and is followed by the best approximation result.

Theorem 5.1** (Mapping of ball to ball).**

Let $\Psi$ be a regular solution of the non-linear system (2.1). For a given fixed $\epsilon>0,$ a sufficiently large $\sigma$ and a sufficiently small discretization parameter chosen as $h=O(\epsilon^{2+\alpha})$ with $\alpha>0$, there exists a positive constant $R(h)$ such that $\mu_{\rm dG}$ maps the ball $\mathbb{B}_{R(h)}({\rm I}_{\rm dG}\Psi)$ to itself;

[TABLE]

Proof.

The solution of (2.1) that belongs to $\Psi\in\mathbf{H}^{2}(\Omega)\cap\mathbf{V}$ satisfies the dG formulation

[TABLE]

A use of the definition and linearity of $\langle{\rm DN}_{h}({{\rm I}_{\rm dG}\Psi})\cdot,\cdot\rangle$, (5.1) and (5.2) leads to

[TABLE]

A use of Lemmas 4.2, 4.3 and 4.5 yields

[TABLE]

A rearrangement of the terms in $B_{\rm dG}(\cdot,\cdot,\cdot,\cdot)$ and a use Lemmas 4.3 and 4.5 lead to

[TABLE]

Set $\tilde{\mathbf{e}}=\Theta_{\textrm{dG}}-\textrm{I}_{\textrm{dG}}\Psi$ and use definition of $B_{\textrm{dG}}(\cdot,\cdot,\cdot,\cdot)$, a regrouping of terms and Remark 4.4 to obtain

[TABLE]

The hidden constants in $"\lesssim"$ in the estimates of $T_{1},T_{2}$, $T_{3}$ and $T_{4}$ depend on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $C_{A},$ $C_{P}$, $C_{S}$, and $C_{I}$. A substitution of the estimates for $T_{1}$, $T_{2}$, $T_{3}$ and $T_{4}$ in (5) with a use of ${|\kern-1.07639pt|\kern-1.07639pt|{\rm I}_{\rm dG}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{\rm dG}\lesssim{|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$ leads to

[TABLE]

A use of Lemma 4.8 yields that there exists a $\Phi_{\textrm{dG}}\in\mathbf{V}_{h}$ with ${|\kern-1.07639pt|\kern-1.07639pt|\Phi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}=1$ such that

[TABLE]

A use of (5.4) and ${|\kern-1.07639pt|\kern-1.07639pt|\tilde{\mathbf{e}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq R(h)$ in (5.5) leads to

[TABLE]

where $C_{3}$ is a positive constant that depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $\beta_{0}$, $C_{A}$, $C_{S}$, $C_{P}$ and $C_{I}$. Assume $h\leq\epsilon^{2+\alpha}$ with $\alpha>0$ so that $h\epsilon^{-2}\leq h^{\frac{\alpha}{2+\alpha}}$. Choose $R(h)=2C_{3}h.$ For $h<h_{4}:=\min({h_{2},h_{3}})$ with $h_{3}^{\frac{\alpha}{2+\alpha}}<\frac{1}{2(1+4C^{2}_{3})^{2}}<\frac{1}{2}$, we obtain

[TABLE]

Since $h<h_{3}<\frac{1}{2^{\frac{2+\alpha}{\alpha}}}<1$ and $\frac{8C_{3}^{3}}{(1+4C_{3}^{2})^{2}}<1,$ ${|\kern-1.07639pt|\kern-1.07639pt|\textrm{I}_{\textrm{dG}}\Psi-\mu_{\rm dG}(\Theta_{\textrm{dG}})|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq 2C_{3}h=R(h).$ This completes the proof of Theorem 5.1 . ∎

Remark 5.2*.*

We derive error estimates with $h$-$\epsilon$ dependency, for a given fixed $\epsilon$. This provides a sufficient condition on the choice of the discretization parameter $h$ for a given fixed $\epsilon$ that ensures convergence. A large value of $\alpha$ would require a very small value of $h$. Equally, a very small choice of $\alpha$ would require $h\to 0$ from the estimates above. In this respect, the choices $\alpha\ll 1$ and $\alpha\gg 1$ lead to computationally expensive scenarios and we only focus on bounded $O(1)$ values of $\alpha$ in this manuscript.

Lemma 5.3** (Contraction result).**

For a given fixed $\epsilon>0,$ a sufficiently large $\sigma$ and a sufficiently small discretization parameter chosen as $h=O(\epsilon^{2+\alpha})$ with $\alpha>0$, the following contraction result holds: $\text{for }\Theta_{1},\Theta_{2}\in\mathbb{B}_{R(h)}({\rm I}_{\rm dG}\Psi)$,

[TABLE]

with some positive constant $C$ independent of $h$ and $\epsilon$ .

Proof.

Let $\Theta_{1}$ and $\Theta_{2}\in\mathbb{B}_{R(h)}(\textrm{I}_{\textrm{dG}}\Psi)$. For $\Phi_{\textrm{dG}}\in\mathbf{V}_{h}$, a use of (5.1), the definition and linearity of $\langle{\rm DN}_{h}({{\rm I}_{\rm dG}\Psi})\cdot,\cdot\rangle$, an elementary manipulation and grouping of term lead to

[TABLE]

Set $\tilde{\mathbf{e}}_{1}=\textrm{I}_{\textrm{dG}}\Psi-\Theta_{1}$, $\tilde{\mathbf{e}}_{2}=\textrm{I}_{\textrm{dG}}\Psi-\Theta_{2}$ and $\mathbf{e}=\Theta_{1}-\Theta_{2}$. A use of Remark 4.4, ${|\kern-1.07639pt|\kern-1.07639pt|\tilde{\mathbf{e}}_{1}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq 2C_{3}h$ and ${|\kern-1.07639pt|\kern-1.07639pt|\tilde{\mathbf{e}}_{2}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq 2C_{3}h$ yields

[TABLE]

A use of (5) and Lemma 4.8 yields that there exists a $\Phi_{\textrm{dG}}\in\mathbf{V}_{h}$ with ${|\kern-1.07639pt|\kern-1.07639pt|\Phi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}=1$ such that

[TABLE]

The hidden constant in $"\lesssim"$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $\beta_{0}$, $C_{P}$ and $C_{3}$. A use of $h\leq\epsilon^{2+\alpha}$ with $\alpha>0$ in (5.8) implies $h\epsilon^{-2}\leq h^{\frac{\alpha}{2+\alpha}}$ and this completes the proof of Lemma 5.3. ∎

Remark 5.4*.*

Note that we require only $h=O(\epsilon^{2})$ to prove the discrete inf-sup conditions in Theorem 4.7 and Lemma 4.8, whereas we need $h=O(\epsilon^{2+\alpha})$ with $\alpha>0$ to prove Theorem 5.1 and Lemma 5.3.

Now we present the proof of results stated in Subsection 3.2.

Proof of Theorem 3.2.

Let $\Psi$ be a regular solution of the non-linear system (2.1) and let $\Psi_{\rm{dG}}$ solve (3.1). A use of Theorem 5.1 yields that $\mu_{\rm dG}$ maps a non-empty convex closed subset $\mathbb{B}_{R(h)}(\textrm{I}_{\textrm{dG}}\Psi)$ of a finite dimensional vector space $\mathbf{V}_{h}$ to itself. Also, $\mu_{\rm dG}$ is continuous. Therefore an application of the Brouwer fixed point theorem [33] yields that $\mu_{\rm dG}$ has at least one fixed point, say $\Psi_{\textrm{dG}}$ in this ball $\mathbb{B}_{R(h)}(\textrm{I}_{\textrm{dG}}\Psi)$ (for details see A.3). That is,

[TABLE]

The contraction result in Lemma 5.3 establishes the uniqueness of the solution of (3.1) for a sufficiently small $h$. The proof of error estimate is straightforward using Lemma 4.5 and (5.9). ∎

An alternative approach using Newton-Kantorovich theorem also provides an explicit formula for the radius $R(h)$ of the ball $\mathbb{B}_{R(h)}(\textrm{I}_{\textrm{dG}}\Psi)$ and proves the existence and uniqueness of discrete solution.

The Newton scheme below is motivated by $\Psi_{\rm dG}^{n}$ and $\Psi_{\rm dG}^{n-1}$ in place of $\Psi_{\rm dG}$ and $\textrm{I}_{\textrm{dG}}\Psi$, respectively, a substitution of $\mu_{\textrm{dG}}(\Psi_{\textrm{dG}}^{n})=\Psi_{\textrm{dG}}^{n}$ in (5.1) and a use of definition of $\langle{\rm DN}_{h}({{\rm I}_{\rm dG}\Psi})\cdot,\cdot\rangle$ in Lemma 4.8. These substitutions yield

[TABLE]

Theorem 5.5**.**

[29, 44, 32]** Suppose that the mapping ${\rm N}_{h}:D\subset\mathbf{V}_{h}\rightarrow\mathbf{V}_{h}$ is Fréchet differentiable on a open convex set $D$, and the derivative ${\rm DN}_{h}(\cdot)$ is Lipschitz continuous on $D$ with Lipschitz constant $L$. For a fixed starting point $x^{0}_{\rm dG}\in D,$ the inverse ${\rm DN}_{h}(x^{0}_{\rm dG})^{-1}$ exists as a continuous operator on $\mathbf{V}_{h}.$ The real numbers $a$ and $b$ are chosen such that

[TABLE]

and $h*:=abL\leq\frac{1}{2}.$ Also, the first approximation $x_{\rm dG}^{1}:={x^{0}_{\rm dG}}-{\rm DN}_{h}(x^{0}_{\rm dG})^{-1}{\rm N}_{h}(x^{0}_{\rm dG})$ has a property that the closed ball $\bar{U}(x_{\rm dG}^{1};r)$ lies within the domain of definition $D,$ where $r=\frac{1-\sqrt{1-2h^{*}}}{aL}-b.$ Then the following are true.

Existence and uniqueness. There exists a solution $x_{\rm dG}\in\bar{U}(x_{\rm dG}^{1};r)$ and the solution is unique on $\bar{U}(x^{0}_{\rm dG};r*)\cap D,$ that is on a suitable neighborhood of the initial point ${x^{0}_{\rm dG}}$ with $r*=\frac{1+\sqrt{1-2h^{*}}}{aL}.$ 2. 2.

Convergence of Newton’s method. The Newton’s scheme with initial iterate $x^{0}_{\rm dG}$ leads to a sequence $x^{n}_{\rm dG}$ in $\bar{U}(x^{0}_{\rm dG};r*)$, which converges to $x_{\rm dG}$ with error bound

[TABLE]

Theorem 5.6** (Existence and uniqueness of discrete solution).**

Let $\Psi$ be a regular solution of the continuous non-linear system ${\rm N}(\Psi;\Phi)=0$ for all $\Phi\in\mathbf{V}.$ For a given fixed $\epsilon>0,$ a sufficiently large $\sigma$ and a sufficiently small discretization parameter chosen as $h=O(\epsilon^{2+\alpha})$ with $\alpha>0$, the following results hold true:

there exists a solution $\Psi_{\rm dG}\in\mathbf{V}_{h}$ to ${\rm N}_{h}(\Psi_{\rm dG};\Phi_{\rm dG})=0$ for all $\Phi_{\rm dG}\in\mathbf{V}_{h}$ such that ${|\kern-1.07639pt|\kern-1.07639pt|\Psi-\Psi_{\rm dG}|\kern-1.07639pt|\kern-1.07639pt|}_{\rm dG}\leq\rho,$ where $\rho=C_{\rho}(h+b+r)$ with $b,r$ defined in Theorem 5.5. The constant $C_{\rho}$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$ and $C_{I}$, 2. 2.

there is at most one solution $\Psi_{\rm dG}$ to ${\rm N}_{h}(\Psi_{\rm dG};\Phi_{\rm dG})=0$ for all $\Phi_{\rm dG}\in\mathbf{V}_{h}$ in $\bar{U}({\rm I}_{\rm dG}\Psi;r*)\cap D$ with $r*=\frac{1+\sqrt{1-C_{K}h\epsilon^{-2}(1+h\epsilon^{-2})}}{2\beta_{0}^{-1}C_{L}\epsilon^{-2}},$ where the constant $C_{K}$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $C_{A}$, $C_{S}$, $C_{P}$, $C_{I}$, $\beta_{0}$ and $C_{L}$, a constant from Lipschitz continuity of ${\rm DN}_{h}(\cdot),$ 3. 3.

the sequence $\Psi^{n}_{\rm dG}$ of iterates converges to $\Psi_{\rm dG}$ and the error bound is given by

[TABLE]

Proof.

A use of the definition of ${\rm DN}_{h}(\cdot)$, Lemma 4.3 and a simple manipulation leads to the fact that ${\rm DN}_{h}(\cdot)$ is Lipschitz continuous on $D$ with Lipschitz constant $L=C_{L}\epsilon^{-2},$ where $C_{L}$ is a constant independent of $\epsilon.$ For a choice of $x^{0}_{\textrm{dG}}={\rm I}_{\rm dG}\Psi$, Lemma 4.8 yields

[TABLE]

with $a=2\beta_{0}^{-1}.$ Given $\Phi_{\textrm{dG}}\in\mathbf{V}_{\textrm{dG}}$ with ${|\kern-1.07639pt|\kern-1.07639pt|\Phi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}=1,$ (5.2) leads to

[TABLE]

A use of the estimates of $T_{1},T_{2},T_{3}$ from (5) yields

[TABLE]

where the constant $C_{N}$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $C_{A}$, $C_{S}$, $C_{P}$ and $C_{I}$. A use of (5.14) and (5.15) yields

[TABLE]

where $b=2C_{N}\beta_{0}^{-1}h(1+h\epsilon^{-2}).$ A sufficiently small choice of $h\leq\epsilon^{2+\alpha}$ with $\alpha>0,$ leads to $h^{*}=abL=4C_{N}C_{L}\beta_{0}^{-2}h\epsilon^{-2}(1+h\epsilon^{-2})\leq 4C_{N}C_{L}\beta_{0}^{-2}h^{\frac{\alpha}{2+\alpha}}(1+h^{\frac{\alpha}{2+\alpha}})<C_{K}h^{\frac{\alpha}{2+\alpha}}$ with $C_{K}=8C_{N}C_{L}\beta_{0}^{-2}$. For a choice of $h<h_{6}:=\min(h_{4},h_{5})$ with $h_{5}^{\frac{\alpha}{2+\alpha}}<\frac{1}{2C_{K}},$ $h^{*}<\frac{1}{2}.$ Therefore, an application of Theorem 5.5 yields the existence of the discrete solution $\Psi_{\rm dG}$ with $r=\frac{1-\sqrt{1-C_{K}h\epsilon^{-2}(1+h\epsilon^{-2})}}{2\beta_{0}^{-1}C_{L}\epsilon^{-2}}-2C_{N}\beta_{0}^{-1}h(1+h\epsilon^{-2})$ and $r*=\frac{1+\sqrt{1-C_{K}h\epsilon^{-2}(1+h\epsilon^{-2})}}{2\beta_{0}^{-1}C_{L}\epsilon^{-2}}$. A use of the second part of Theorem 5.5 yields the error bound ${|\kern-1.07639pt|\kern-1.07639pt|{\rm I}_{\rm dG}\Psi-\Psi_{\rm dG}|\kern-1.07639pt|\kern-1.07639pt|}_{\rm dG}\leq\frac{1-\sqrt{1-C_{K}h\epsilon^{-2}(1+h\epsilon^{-2})}}{2\beta_{0}^{-1}C_{L}\epsilon^{-2}}=O(h)$ which justifies the choice of $R(h)$ in Theorem 5.1. A use of triangle inequality, Lemma 4.5 and the second estimate in (5.11) for the first Newton correction leads to

[TABLE]

A substitution of $h^{*},a$ and $L$ in (5.12) yields the error bound in the Newton’s convergence given by (5.13). ∎

Proof of Theorem 3.3.

Let $\Psi^{*}_{\textrm{dG}}$ be the best approximation of $\Psi$ in $\mathbf{V}_{h}$. Then

[TABLE]

Set $\mathbf{e}_{\textrm{dG}}=\Psi^{*}_{\rm dG}-\Psi_{\rm dG}\in\mathbf{V}_{h}.$ A use of Theorem 4.7 yields that there exists a $\Phi_{\rm dG}\in\mathbf{V}_{h}$ with ${|\kern-1.07639pt|\kern-1.07639pt|\Phi_{\rm dG}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}=1$ such that

[TABLE]

Set $\mathbf{e}_{1}=\Psi-\Psi_{\rm dG}$. A use of Taylor series expansion of ${\rm N}_{h}(\cdot;\cdot)$ around $\Psi$ leads to

[TABLE]

A use of (3.1) and (5.2) lead to

[TABLE]

Rewrite $\mathbf{e}_{1}=(\Psi-\Psi^{*}_{\rm dG})+\mathbf{e}_{\textrm{dG}}$. A use of the linearity of $\langle{\rm DN}_{h}(\Psi)\cdot,\cdot\rangle$, (5.16) and the above equality leads to

[TABLE]

Since $\langle{\rm DN}_{h}(\Psi)\mathbf{e}_{1},\Phi_{\rm dG}\rangle=A_{\rm dG}(\mathbf{e}_{1},\Phi_{\rm dG})+3B_{\rm dG}(\Psi,\Psi,\mathbf{e}_{1},\Phi_{\rm dG})+C_{\rm dG}(\mathbf{e}_{1},\Phi_{\rm dG}),$ $\langle{\rm D}^{2}{\rm N}_{h}(\Psi)(\mathbf{e}_{1})\mathbf{e}_{1},\Phi_{\rm dG}\rangle=6B_{\textrm{dG}}(\mathbf{e}_{1},\mathbf{e}_{1},\Psi,\Phi_{\rm dG})$ and $\langle{\rm D}^{3}{\rm N}_{h}(\Psi)(\mathbf{e}_{1})(\mathbf{e}_{1})\mathbf{e}_{1},\Phi_{\rm dG}\rangle=6B_{\textrm{dG}}(\mathbf{e}_{1},\mathbf{e}_{1},\mathbf{e}_{1},\Phi_{\rm dG}),$ a use of Lemmas 4.2 and 4.3 leads to

[TABLE]

The constant in $"\lesssim"$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $C_{A}$, $C_{S}$ and $C_{P}$. A use of (5.16) and (5.18) in (5.17) yields

[TABLE]

The triangle inequality and (5.19) leads to

[TABLE]

For a sufficiently small choice of the discretization parameter $h=O(\epsilon^{2+\alpha})$ with $\alpha>0$, use ${|\kern-1.07639pt|\kern-1.07639pt|\mathbf{e}_{1}|\kern-1.07639pt|\kern-1.07639pt|}_{\rm dG}\leq C_{e}h$ in Theorem 3.2 and (5.20) to obtain

[TABLE]

where the constant $C_{4}$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $C_{A}$, $C_{S}$, $C_{P}$ and $\beta_{0}$. Since $h\epsilon^{-2}\leq h^{\frac{\alpha}{2+\alpha}}$, a sufficiently small choice of $h<h_{8}:=\min(h_{6},h_{7})$ with $h_{7}^{\frac{\alpha}{2+\alpha}}=\frac{C_{4}}{2(C_{e}+C_{e}^{2})}$ completes the proof of Theorem 3.3 for $C_{B}=\frac{2}{C_{4}}.$ ∎

The next two lemmas are required to prove Theorem 3.5.

Lemma 5.7**.**

For $\lambda=-1$, any $\Psi\in\mathbf{H}^{2}(\Omega)$, $\boldsymbol{\chi}\in\mathbf{H}^{2}(\Omega)\cap\mathbf{V}$ and the interpolation ${\rm I}_{\rm dG}\Psi\in\mathbf{V}_{h}$ satisfy

[TABLE]

Proof.

A use of definition of $A_{\rm dG}(\cdot,\cdot)$, $[\boldsymbol{\chi}]=0$ on $\mathcal{E}_{i},$ $\boldsymbol{\chi}=0$ on $\mathcal{E}_{D}$ and integration by parts yields

[TABLE]

We have $\sum_{T\in\mathcal{T}}\int_{\partial T}\frac{\partial\boldsymbol{\chi}}{\partial\eta}\cdot({\rm I}_{\rm dG}\Psi-\Psi)\,\textrm{ds}=\sum_{E\in\mathcal{E}_{i}}\int_{E}[\frac{\partial\boldsymbol{\chi}}{\partial\eta}]\cdot\{{\rm I}_{\rm dG}\Psi-\Psi\}\,\textrm{ds}+\sum_{E\in\mathcal{E}}\int_{E}\{\frac{\partial\boldsymbol{\chi}}{\partial\eta}\}\cdot[{\rm I}_{\rm dG}\Psi-\Psi]\,\textrm{ds}.$ Since $\boldsymbol{\chi}\in\mathbf{H}^{2}(\Omega),$ $[\nabla{\boldsymbol{\chi}}]=0$ for all $E\in\mathcal{E}_{i}.$ A use of Cauchy-Schwartz inequality and Lemma 4.5 leads to,

[TABLE]

where the constant suppressed in $"\lesssim"$ depends on $C_{I}.$ This concludes the proof. ∎

For given $G\in\mathbf{L}^{2}(\Omega)$ and $\Psi$ that solves (1.2), consider the linear dual problem:

[TABLE]

The weak formulation that corresponds to (5.21) seeks $\boldsymbol{\chi}\in\mathbf{V}$ such that

[TABLE]

A use of (2.8) establishes the well-posedness of (5.22).

Lemma 5.8**.**

A solution $\boldsymbol{\chi}$ of (5.22) belongs to $\mathbf{H}^{2}(\Omega)\cap\mathbf{V}$ and satisfies ${|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\lesssim(1+\epsilon^{-2}){|\kern-1.07639pt|\kern-1.07639pt|G|\kern-1.07639pt|\kern-1.07639pt|}_{0},$ where the hidden constant in $"\lesssim"$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2},C_{S}$ and $\beta.$

Proof.

A use of (2.8) and (5.22) yields

[TABLE]

A use of (2.5) and (2.7) implies that ${|\kern-1.07639pt|\kern-1.07639pt|-3B(\Psi,\Psi,\cdot,\boldsymbol{\chi})-C(\cdot,\boldsymbol{\chi})+G|\kern-1.07639pt|\kern-1.07639pt|}_{0}\lesssim\epsilon^{-2}{|\kern-1.07639pt|\kern-1.07639pt|\boldsymbol{\chi}|\kern-1.07639pt|\kern-1.07639pt|}_{1}+{|\kern-1.07639pt|\kern-1.07639pt|G|\kern-1.07639pt|\kern-1.07639pt|}_{0}$. Hence, the elliptic regularity [20] with a boot-strapping argument and (5.23) completes the proof. ∎

Proof of Theorem 3.5.

Set $G=\textrm{I}_{\textrm{dG}}\Psi-\Psi_{\textrm{dG}}$ in (5.21). Multiply (5.21) by $\Phi_{\rm dG}=\textrm{I}_{\textrm{dG}}\Psi-\Psi_{\textrm{dG}}$ and integrate by parts to obtain

[TABLE]

Here, $\langle{\rm DN}_{h}(\Psi)\textrm{I}_{\textrm{dG}}\Psi-\Psi_{\textrm{dG}},\boldsymbol{\chi}\rangle=A_{\rm dG}(\textrm{I}_{\textrm{dG}}\Psi-\Psi_{\textrm{dG}},\boldsymbol{\chi})+B_{\rm dG}(\Psi,\Psi,\textrm{I}_{\textrm{dG}}\Psi-\Psi_{\textrm{dG}},\boldsymbol{\chi})+C_{\rm dG}(\textrm{I}_{\textrm{dG}}\Psi-\Psi_{\textrm{dG}},\boldsymbol{\chi}).$ Note that the terms that involve $[\boldsymbol{\chi}]$ in the definition of $A_{\rm dG}(\cdot,\cdot)$ are zero. However, they are retained for the ease of further manipulations. A use of (3.1) and the fact that $\Psi\in\mathbf{H}^{2}(\Omega)$ satisfies the discrete formulation (3.1) leads to

[TABLE]

Lemmas 4.2, 4.3, 4.5 and 5.7 and Theorem 3.2 yield

[TABLE]

Lemmas 4.3 and 4.5 lead to

[TABLE]

Set $\mathbf{e}_{3}=\Psi_{\textrm{dG}}-\Psi$ and estimate $T_{8}$ as in $T_{4}$ of Theorem 5.1 and use Theorem 3.2 to obtain

[TABLE]

A combination of the estimates for $T_{5}$, $T_{6}$, $T_{7}$ and $T_{8}$ in (5.25) and a use of Lemma 5.8 yields

[TABLE]

A use of (5.26) and triangle inequality yields

[TABLE]

where the constants suppressed in $"\lesssim"$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2}$, $\beta_{0}$, $C_{A}$, $C_{S}$, $C_{P}$ and $C_{I}$. This completes the proof of Theorem 3.5. ∎

6 Numerical Implementation

In Subsection 6.1, we prove Theorem 3.6. The second subsection discusses results of numerical experiments that justify the theoretical estimates.

6.1 Convergence of Newton’s method

Now we establish that Newton iterates in (5) converges quadratically to the discrete solution. The proof follows by modification of the approach used in [13]. While the linearized operator in [13] has a system of bilinear and trilinear forms; in this case we have system of bilinear and quadrilinear forms that leads to modification in the bounds. Moreover, the choice of the radius of the ball $\rho_{1},$ in which the initial guess $\Psi_{\rm dG}^{0}$ needs to be chosen so as the Newton’s method converges depends on the non-linearity and needs to be chosen carefully. The effect of the parameter ${}^{\prime}\epsilon^{\prime}$ is also considered in this proof.

Proof of Theorem 3.6.

Following the proof of Lemma 4.8, for a sufficiently small choice of the discretization parameter $h=O(\epsilon^{2}),$ there exists a sufficiently small constant $\delta>0$ independent of $h$ such that for each $Z_{\textrm{dG}}\in\mathbf{V}_{h}$ that satisfies ${|\kern-1.07639pt|\kern-1.07639pt|Z_{\textrm{dG}}-\textrm{I}_{\textrm{dG}}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\delta$, the bilinear form

[TABLE]

satisfies discrete inf-sup condition given by

[TABLE]

For a sufficiently small choice of the discretization parameter $h=O(\epsilon^{2+\alpha})$ with $\alpha>0$, (5.9) leads to

[TABLE]

Assume that the initial guess $\Psi_{\textrm{dG}}^{0}$ satisfies ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{0}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\rho_{1}\leq\frac{\delta}{2}.$ A use of this and (6.3) along with a triangle inequality leads to

[TABLE]

Therefore, $Z_{\textrm{dG}}=\Psi_{\textrm{dG}}^{0}$ in (6.1) leads to the discrete inf-sup condition for $\langle DN_{h}(\Psi_{\textrm{dG}}^{0})\cdot,\cdot\rangle$, and this implies that there exists a unique $\Psi_{\textrm{dG}}^{1}\in\mathbf{V}_{h}$ ( which is the first Newton iterate in (5) ) satisfying the well-posed system

[TABLE]

The discrete inf-sup condition (6.2) implies the existence of a $\Phi_{\textrm{dG}}\in\mathbf{V}_{h}$ with ${|\kern-1.07639pt|\kern-1.07639pt|\Phi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}=1$ such that

[TABLE]

A use of (6.1), (6.5) and (3.1) in (6.6) yields

[TABLE]

The right hand side of (6.1) is estimated analogous to $T_{4}$ in Theorem 5.1. This followed by a use of the triangle inequality and (6.3) leads to

[TABLE]

where the constant $C_{5}$ depends on ${|\kern-1.07639pt|\kern-1.07639pt|\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{2},$ $\beta_{0}$, $C_{3},$ $C_{P}$ and $C_{S}.$ A choice of $\rho_{1}<\min(\frac{\delta}{2},\frac{-1+\sqrt{1+2\epsilon^{2}C_{5}^{-1}}}{2})$ yields ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{1}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\frac{1}{2}{|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{0}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\frac{\rho_{1}}{2}$. A use of this estimate, (6.3) and triangle inequality yields ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{1}-\textrm{I}_{\textrm{dG}}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq{|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{1}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}+{|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}-\textrm{I}_{\textrm{dG}}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\delta.$ Also, ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{1}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq C_{6}{|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{0}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}^{2},$ where $C_{6}$ depends on $C_{5}$, $\rho_{1}$ and $\epsilon$.

Therefore, ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{j}-\textrm{I}_{\textrm{dG}}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\delta$ and ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{j}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\frac{\rho_{1}}{2^{j}}$ satisfies for $j=1.$ Assume that this holds for some $j\in\mathbb{N}.$ Then, $Z_{\textrm{dG}}=\Psi_{\textrm{dG}}^{j}$ in (6.1) leads to the discrete inf-sup condition for $\langle{\rm DN}_{h}(\Psi_{\textrm{dG}}^{j})\cdot,\cdot\rangle$ which implies that there exists a unique $\Psi_{\textrm{dG}}^{j+1}\in\mathbf{V}_{h}$ satisfying the well-posed system

[TABLE]

Then, a use of (6.2) and following the proof for $j=1$, we obtain

[TABLE]

Since $\rho_{1}<\min(\frac{\delta}{2},\frac{-1+\sqrt{1+2\epsilon^{2}C_{5}^{-1}}}{2})$, ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{j+1}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\frac{1}{2}{|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{j}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\frac{\rho_{1}}{2^{j+1}}$ and ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{j+1}-\textrm{I}_{\textrm{dG}}\Psi|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq\delta$ with a quadratic convergence rate given by ${|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{j+1}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}\leq C_{6}{|\kern-1.07639pt|\kern-1.07639pt|\Psi_{\textrm{dG}}^{j}-\Psi_{\textrm{dG}}|\kern-1.07639pt|\kern-1.07639pt|}_{\textrm{dG}}^{2}.$ This completes the proof using mathematical induction. ∎

6.2 Numerical experiments

In this subsection, the computational error and convergence rate of discrete solutions for dGFEM are illustrated for some benchmark problems. We study the convergence of discrete solutions, for various values of $\epsilon$. These numerical experiments have been implemented using FEniCS [3] library and the results support the theoretical findings (for the details of implementation procedure see A.4). Let $e_{i}$ and $h_{i}$ be the error and the mesh parameter at the $i$-th level, respectively. The $i$-th level experimental order of convergence is defined by $\displaystyle\alpha_{i}:=log(e_{n}/e_{i})/log(h_{n}/h_{i})$ for $i=1,\ldots,n-1$ and $n$ corresponds to the final iteration considered in numerical experiments.

Example 6.2.1*.*

For the problem (1.2), set $\Omega=(0,1)\times(0,1)$ and the parameter value $\epsilon=0.2$. Compute the right hand side for the manufactured exact solution $u=x(1-x)y(1-y)$ and $v=x(1-x)y(1-y).$

We discretise the domain $\Omega=(0,1)\times(0,1)$ into triangles and in the uniform refinement process, each triangle $T$ is divided into four similar triangles. Tables 1, 2 and 3 present the numerical errors and orders of convergence in energy and $\mathbf{L}^{2}$ norms computed using piecewise polynomials of degree $1,$ $2$ and $3$, respectively.

Remark 6.1*.*

This numerical example verifies that the theoretical convergence rates obtained in energy norm (resp. $\mathbf{L}^{2}$ norm) are $1$, $2$ and $3$ (resp. $2,$ $3$ and $4$) for piecewise $P_{1}$, $P_{2}$ and $P_{3}$ polynomial approximations. Improved rate of convergences suggest an improvement in the theoretical estimates if the exact solution is smooth. The convergence analysis for $h$-$p$ dGFEM taking into account the effect of $h$-$\epsilon$-$p$ dependency is a problem of future interest.

Example 6.2.2*.*

Consider the problem (1.2) in $\Omega=(0,1)\times(0,1).$ We will approximate the system (1.2) using the Dirichlet boundary condition [37] given by

[TABLE]

where the parameter $d=3\epsilon$ and the trapezoidal shape function $\textit{T}_{d}:[0,1]\rightarrow{\mathbb{R}}$ is defined to be

[TABLE]

This example is motivated by the benchmark problem studied in [37]. There are two classes of stable experimentally observable configurations [see Figure (2)]: the diagonal states for which the director is aligned along the cell/square diagonals and the rotated states, for which the directors rotate by $\pi$-radians across the width of cell.

We use Newton’s method to approximate the six different solutions $\Psi=(u,v)$ corresponding to the D1, D2, R1, R2, R3 and R4 states. We compute six different initial conditions by solving the Laplace equation [37] with Dirichlet boundary conditions as specified in Table 4. For example, in the D1 case, the system $\displaystyle\Delta\theta=0\text{ in }\Omega,\theta(x,0)=\theta(x,1)=0\text{ and }\theta(0,y)=\theta(1,y)=\frac{\pi}{2}$ is solved using dGFEM. Then the corresponding initial condition for Newton’s method (for the Landau-de Gennes mdodel) is defined to be

[TABLE]

where $s=1$ at the interior nodes and $s=\left\lvert\mathbf{g}\right\rvert$ at the boundary nodes. We obtain the expected theoretical convergence rates using piecewise linear polynomial in dG and $\mathbf{L}^{2}$ norm as $O(h)$ and $O(h^{2})$, respectively. In Tables 5 and 6, we record the computational errors and convergence rates of solutions for the diagonal state D1 and rotated state R1 respectively for $\epsilon=0.02$. Similar errors and optimal convergence rates are obtained for D2, R2, R3 and R4 solutions. In Figure 2, we plot the converged director plots and scalar order parameter for the six states, D1, D2, R1, R2, R3 and R4, respectively. In Figures 3 - 5, we plot the level curves of the corresponding converged D1 and D2 diagonal, R1, R2, R3 and R4 rotated solutions.

Figure 6 plots the dG error vs the discretization parameter for various values of $\epsilon$. Note that errors are sensitive to the choice of discretization parameter as $\epsilon$ decreases. It is difficult to verify the exact dependence of $h$ and $\epsilon$ from these numerical results, except that as $\epsilon\rightarrow 0$, smaller mesh-sizes are needed for convergence.

Remark 6.2*.*

In [37], the authors study the model problem of nematic-filled square wells in the conforming FEM set up and present numerical errors and convergence rates for conforming FEM, in $\mathbf{H}^{1}$ and $\mathbf{L}^{2}$ norm to be of $O(h)$ and $O(h^{2}),$ respectively, for the six different stable states. However, they do not study the convergence trends as a function of $\epsilon$ and do not establish the theoretical order of convergence.

Example 6.2.3*.*

Consider an annular domain filled with the compound $5CB$, a standard liquid crystal material [15]. We consider the weak formulation of the Euler-Lagrange equations for a reduced two-dimensional Landau-de Gennes energy as in [15]

[TABLE]

where $\mathcal{C}=\frac{\tilde{C}}{\left\lvert A\right\rvert}$ is a constant that depends on the material parameters $A$ and $\tilde{C}$ of the liquid crystal. The parameter values are $A=0.172\times 10^{6}N/m^{2}$ and $\tilde{C}=1.73\times 10^{6}N/m^{2}$ for the compound $5CB$. The domain has an outer radius $L_{0}$ and the inner radius $0.5L_{0}$, where $L_{0}$ is the characteristic length scale. The ratio ($\tilde{\rho}$) of the inner to the outer radius has been chosen to be $0.5$ to capture the radial solution analyzed by Bethuel et al. [8]. In [8], the authors study the Oseen-Frank radial solution in this annulus within the one-constant framework with energy given by $\mathcal{E}_{OF}(\mathbf{n})=K\int_{\Omega}\left\lvert\nabla\mathbf{n}\right\rvert^{2}\,\textrm{dx}$, where $K>0$ is an elastic constant. The radial solution is simply given by $\mathbf{n}=\left(\frac{x}{\sqrt{x^{2}+y^{2}}},\frac{y}{\sqrt{x^{2}+y^{2}}}\right)$. In the reduced Landau-de Gennes framework, the corresponding manufactured solution, $\Psi=(u,v)$, is defined by $u=\frac{2x^{2}}{(x^{2}+y^{2})}-1\text{ and }v=\frac{2xy}{(x^{2}+y^{2})}$, or equivalently $\Psi=\left(2\mathbf{n}\otimes\mathbf{n}-I\right)$, $\mathbf{n}$ is the Oseen-Frank radial solution studied in [15] and $I$ is the $2\times 2$ identity matrix. The data $\mathbf{f}$ and $\mathbf{g}$ are calculated by substituting this manufactured solution into the left-hand side of (6.10).

We use dGFEM formulation based on piecewise linear polynomials over triangles, to study the convergence of the discrete solution to this manufactured solution. The mesh is generated using mshr, the mesh generation component in FEniCS. The curved boundary is approximated using polygons. The errors and order of convergences in both energy and $\mathbf{L}^{2}$ norms are tabulated in Table 7 and the results are consistent with our theoretical estimates. In Figure 7a, we plot the initial triangulation for the annular domain and in Figures 7b - 7c, we plot the level curves of the corresponding converged radial solution.

Remark 6.3*.*

In [15], (6.10) is implemented using conforming FEM and the numerical orders of convergences in $\mathbf{L}^{\infty},$ $\mathbf{L}^{2}$ and $\mathbf{H}^{1}$ norms are observed to be of orders of $O(h^{2}),$ $O(h^{2})$ and $O(h)$, respectively. However, the theoretical error estimates are not derived in [15].

7 Conclusion

The paper focuses on a rigorous dG formulation of the model problem studied in [37, 42]. It is not clear if the dG formulation gives us any additional physical insight for this model problem, at least for micron-sized wells. However, dG methods combined with a posteriori error analysis generally work well for problems with less regularity and with non-homogeneous boundary conditions or for singular liquid crystal potentials as proposed by Katriel et.al [31] and studied in [6], or when $\epsilon\rightarrow 0$ (sharp interface limit). Hence, we expect our work to be foundational for subsequent cutting-edge numerical studies of interfacial phenomena, higher-dimensional defects, quenching phenomena for LCs which are necessarily characterized by singular behavior. A particularly interesting application would be the dG formulation of the continuum theory of ferronematics [10]. Ferronematics in confinement naturally exhibit domain walls and interior point defects and it would be interesting to see how dG results compare to numerical results from confoming FEM, particularly near structural defects e.g. nematic vortices and boundary layers.

Acknowledgements

R.M. gratefully acknowledges support from institute Ph.D. fellowship and the Visiting Postgraduate Scholar scheme run by the University of Bath. A.M. gratefully acknowledges hospitality from IIT Bombay. A.M. is supported by an EPSRC Career Acceleration Fellowship EP/J001686/1 and EP/J001686/2, University of Bath Internationalisation schemes and an OCIAM Visiting Fellowship, the Keble Advanced Studies Centre. N.N. gratefully acknowledges support from the Faculty of Science at the University of Bath and the support by DST SERB MATRICS grant MTR/2017/000 199. The authors gratefully acknowledge Prof. M. Vanninathan for the insightful discussions and also the HPC facility of IIT Bombay for the computational resources. The authors thank the referees for the constructive comments that has improved the quality of the paper significantly.

Appendix A Appendix

A.1 Derivation of weak formulation

The Landau-de Gennes energy functional considered in this article is given by

[TABLE]

with Dirichlet boundary condition $\Psi=\mathbf{g}$ on $\partial\Omega.$ Consider the real valued function $i[\tau]:=\mathcal{E}(\Psi+\tau\Phi)\,\,(\tau\in\mathbb{R})$ with $\Phi\in\mathbf{V}.$ Since $\Psi$ is a minimizer of $\mathcal{E}(\cdot)$ and $\Psi+\tau\Phi=\Psi=\mathbf{g}$ on $\partial\Omega,$ $\mathcal{E}(\cdot)$ has a minimum at $\tau=0.$ Therefore, $i^{\prime}(0)=\lim_{\tau\rightarrow 0}\dfrac{\mathcal{E}(\Psi+\tau\Phi)-\mathcal{E}(\Psi)}{\tau}=0.$ A manipulation using the definition of $\mathcal{E}(\cdot)$ leads to

[TABLE]

which yields the weak formulation (2.1) in operator form. The operator form in (2.1) has a representation of the nonlinear part in terms of $B(\cdot,\cdot,\cdot,\cdot)$ and this is crucial for an elegant analysis.

A.2 Regularity result

Theorem A.1**.**

[26]** [33] Let $\Omega$ be a convex, bounded and open subset of $\mathbb{R}^{2}$ with smooth boundary $\partial\Omega.$ Then, for $F\in\mathbf{L}^{2}(\Omega)$ and $\mathbf{g}\in\mathbf{H}^{\frac{3}{2}}(\partial\Omega)$, there exists a $\Psi\in\mathbf{H}^{2}(\Omega)$ such that $-\Delta\Psi=F\text{ in }\Omega\text{ and }\Psi=\mathbf{g}\text{ on }\partial\Omega.$

Proof of Lemma 2.4.

Rewrite the system (1.2) as: $-\Delta\Psi=F_{1}(\Psi)\text{ in }\Omega\text{ and }\Psi=\mathbf{g}\text{ on }\partial\Omega,$ where $F_{1}(\Psi)=-\dfrac{2}{\epsilon^{2}}(\left\lvert\Psi\right\rvert^{2}-1)\Psi$. Expand the expression for $F_{1}(\Psi)$, where $\Psi=(u,v)\in\mathbf{H}^{1}(\Omega)$, use the Sobolev embedding result $H^{1}(\Omega)\hookrightarrow L^{p}(\Omega)$ for all $p\geq 1$, and the Hölder’s inequality to prove that $F_{1}(\Psi)\in\mathbf{L}^{2}(\Omega).$ Now a use of Theorem A.1 and a bootstrapping argument [20] imply that the solution $\Psi$ of (1.2) belongs to $\mathbf{H}^{2}(\Omega)$. ∎

A.3 Details of proof of Theorem 3.2

Theorem A.2** (Brouwer fixed point theorem).**

[33]** Let $Y$ be a finite dimensional Hilbert space and $f:K\rightarrow K$ be a continuous map from a non-empty, compact and convex subset $K$ of $Y$ which maps into $K.$ Then $f$ has a fixed point on $K$.

A use of Lemma 5.3 yields that $\mu_{\textrm{dG}}$ is Lipschitz continuous on the ball $\mathbb{B}_{R(h)}({\rm I}_{\rm dG}\Psi)$. Set $K=\mathbb{B}_{R(h)}({\rm I}_{\rm dG}\Psi)$ and $Y=\mathbf{V}_{h}$. Theorem 5.1 yields $\mu_{\textrm{dG}}$ maps $K$ to itself. Therefore, an application of Theorem A.2 yields that $\mu_{\textrm{dG}}$ has a fixed point, say $\Psi_{\rm dG}$. Now we prove any fixed point of $\mu_{\textrm{dG}}$ is a solution of (3.1). Since $\Psi_{\rm dG}$ is a fixed point of the map $\mu_{\textrm{dG}}$, we have

[TABLE]

A use of (5.1) on the left hand side and definition of $\langle DN_{h}({{\rm I}_{\rm dG}\Psi})\cdot,\cdot\rangle$ on the right hand side leads to

[TABLE]

This implies $\Psi_{\rm dG}$ solves (3.1). The uniqueness of the solution $\Psi_{\rm dG}$ follows from the contraction result in Lemma 5.3.

A.4 Implementation procedure

We compute the approximate solutions $\Psi_{\rm dG}$ using Newton’s method. To compute the initial guess $\Psi^{0}_{\rm dG}$ in Example 6.2.2, first solve the Oseen-Frank system

[TABLE]

with the corresponding Dirichlet boundary condition, say $g_{D}$. The details of boundary functions corresponding to each diagonal and rotated solution can be found in Table 4. The dG formulation corresponding to (A.1) seeks $\theta_{\textrm{dG}}\in V_{h}$ such that for all $\phi_{\textrm{dG}}\in V_{h},$

[TABLE]

**Algorithm

**

In Examples 6.2.1 and 6.2.3, where the exact solutions are known, the initial guesses for the Newton’s method are chosen as the solutions of the corresponding linear systems.