Unconditionally energy stable discontinuous Galerkin schemes for the Cahn-Hilliard equation

Hailiang Liu^† and Peimeng Yin^§ ^‡ Iowa State University, Department of Mathematics, Ames, IA 50011 [email protected] ^§ Wayne State University, Department of Mathematics, Detroit, MI 48202 [email protected]

Abstract.

In this paper, we introduce novel discontinuous Galerkin (DG) schemes for the Cahn-Hilliard equation, which arises in many applications. The method is designed by integrating the mixed DG method for the spatial discretization with the Invariant Energy Quadratization (IEQ) approach for the time discretization. Coupled with a spatial projection, the resulting IEQ-DG schemes are shown to be unconditionally energy dissipative, and can be efficiently solved without resorting to any iteration method. Both one and two dimensional numerical examples are provided to verify the theoretical results, and demonstrate the good performance of IEQ-DG in terms of efficiency, accuracy, and preservation of the desired solution properties.

Key words and phrases:

Cahn-Hilliard equation, energy dissipation, mass conservation, DG method, IEQ method.

1991 Mathematics Subject Classification:

65N12, 65N30, 35K35

1. Introduction

The Cahn-Hilliard (CH) equation, originally introduced in [6] as a model of phase separation in binary alloys, has become a fundamental equation as well as a building block in the phase field methodology for moving interface problems arising from various applications (see, e.g., [29] for the references therein).

This work is concerned with high order numerical approximations to the Cahn-Hilliard problem: find $\{u(x,t),w(x,t)\}$ for $x\in\Omega$ and $t>0$ such that

$\displaystyle u_{t}$	$\displaystyle=\nabla\cdot(M(u)\nabla w),$	(1.1)
$\displaystyle w$	$\displaystyle=-\epsilon^{2}\Delta u+F^{\prime}(u),$
$\displaystyle u(x,0)$	$\displaystyle=u_{0}(x),\quad x\in\Omega,$

where $\Omega\subseteq\mathbf{R}^{d}(d=1,2,3)$ is a bounded domain, $\epsilon$ is a positive parameter, $M(u)\geq 0$ is the mobility function, $F(u)$ is the nonlinear bulk potential, and $u_{0}(x)$ is the initial data.

The well-posedness study of the Cahn–Hilliard equation is very rich, and the existing results are mainly for two types of models: one is the constant mobility with polynomial potentials ([13]), and another is the degenerate mobility of form $M(u)=u(1-u)$ ([11]) with the potential of logarithmic types (see, e.g., [6, 12, 44, 11, 10]). Here we numerically study model (1.1) with no restrictions on the specific form of the mobility and free energy. For the degenerate CH model, we apply the established scheme to regularized systems as discussed in Section 4. While we also refer the reader to [31] for a new relaxation system to approximate the degenerate CH model.

We consider in this paper the following boundary conditions:

\text{(i)}\ u\text{ is periodic};\quad\text{or (ii)}\ \partial_{\mathbf{n}}u=M(u)\partial_{\mathbf{n}}w=0,\quad x\in\partial\Omega.

(1.2)

Here $\mathbf{n}$ stands for the unit outward normal to the boundary $\partial\Omega$ . With such boundary conditions, equation (1.1) is endowed with a gradient flow structure $u_{t}=\nabla\cdot(M(u)\frac{\delta}{\delta u}\mathcal{E})$ , dictated by the energy dissipation law

\displaystyle\frac{d}{dt}\mathcal{E}(u)=-\int_{\Omega}M(u)|\nabla w|^{2}\leq 0,

(1.3)

where the total free energy is defined by

\mathcal{E}(u)=\int_{\Omega}\left(\frac{\epsilon^{2}}{2}|\nabla u|^{2}+F(u)\right)dx.

(1.4)

The CH model is nonlinear and its analytical solutions are intractable. Also, for a gradient flow such as the CH equation steady states are of particular interest. Hence, designing accurate, efficient, and energy stable algorithms to solve it becomes essential, in particular for the accuracy of long time simulations. Keeping the energy dissipation (1.3) has been a major concern in the design of various numerical schemes, see, e.g. [5, 14, 18, 34, 32, 17, 35, 42, 7].

In this paper, we aim to develop unconditionally energy stable discontinuous Galerkin (DG) schemes for solving the CH model problem. To achieve this, we face two main challenges: (i) how to handle fourth order derivatives in the DG discretization; and (ii) how to handle the nonlinear term associated with the potential $F$ in time discretization.

For (i), several approaches have been adopted to deal with difficulties caused by the high order solution derivatives. The first one is the local DG (LDG) methods [39, 21, 33], with which the original equation is rewritten into a first order system for further DG discretization. The second one is the mixed symmetric interior penalty DG (SIPG) methods [37, 15, 16, 17], with which the penalty terms are added as interface corrections upon the global solution formulation so that the resulting scheme is stable. It is also possible to apply an ultra-weak DG discretization, such as the DG scheme in [9] for the one-dimensional biharmonic equation. The spatial DG discretization we adopt here takes advantages of both the SIPG method and the mixed DG method without interior penalty in [25, 26] for another class of fourth order PDEs.

For (ii), there are several related time discretization techniques available in the literature, including the so-called convex splitting approach [14, 38, 7] and the stabilization approach [32, 40]. The convex splitting approach, introduced in the pioneering work [14], treats the convex part of the potential implicitly and the concave part explicitly. As a result such method is unconditionally energy stable, however, it produces nonlinear schemes, thus the implementations are often complicated with potentially high computational costs. The stabilization approach, introduced in [40] for epitaxial growth models, treats the nonlinear term explicitly so that the resulting scheme becomes stable after a stabilization term is added to avoid strict time step constraints. However, a large stabilization constant is often needed to ensure energy stability.

A more feasible remedy would be to make a reasonable transform of the given potential. One technique is the Invariant Energy Quadratization (IEQ) approach introduced in [41, 45], which generalized the two types of linear energy stable schemes introduced in [19]. The IEQ method is known to provide flexibilities to treat the nonlinear terms since it only requests that the nonlinear potential be bounded from below. We note that recently a related approach, called the SAV method, has been introduced in [36, 35] with certain advantages over the IEQ. Yet efficiently solving the resulting linear system when coupled with a DG spatial discretization is subtle, see [27].

In order to design IEQ-DG schemes, our strategy is to start with the model satisfying two basic assumptions:
(i) the mobility function $M(u)$ satisfies

M(u)\geq M_{\min}>0;

(ii) there exist a constant $B>0$ , such that

F(u)>-B,

for any $u$ under consideration, and further discuss how to extend the established schemes to more general cases.

The IEQ-DG method introduced in [26] for the Swift–Hohenberg equation has several advantages, such as high order of accuracy, easy implementation, and efficiency. Unfortunately, the DG discretization in [26] when applied to the CH equation (1.1) does not seem to yield the desired energy stability. Therefore, in this paper, we exploit the direct DG (DDG) method [24] for both $u$ and $w$ , which when taking a special choice of the flux parameters can be reformulated into a mixed symmetric interior penalty DG (SIPG) scheme (see, e.g., [16]).

The IEQ approach for time discretization relies on an auxiliary variable $U=\sqrt{F(u)+B}$ , so that

U_{t}=\frac{1}{2}H(u)u_{t},\quad H(u):=F^{\prime}(u)/\sqrt{F(u)+B}.

With this transformation we update $U^{n}$ in two steps: the piecewise $L^{2}$ projection with $U_{h}^{n}=\Pi U^{n}$ , and the update step with

\frac{U^{n+1}-U_{h}^{n}}{\Delta t}=\frac{1}{2}H(u^{n}_{h})\frac{u_{h}^{n+1}-u_{h}^{n}}{\Delta t}.

The resulting IEQ-DG scheme follows from replacing the nonlinear function $F^{\prime}(u_{h}^{n+1})$ by $H(u^{n}_{h})U^{n+1}$ in the DG discretization (see the scheme formulation in (3.4)). In addition, the resulting discrete systems are linear. As a result, the methods are simple to implement and computationally efficient. With the techniques exploited in [28], the IEQ-DG schemes may be 6 extended to DG schemes of arbitrarily high order in time, yet still unconditionally energy stable.

Finally, closest to our work is [22], where the authors, building on the IEQ formulation with the LDG spatial discretization for phase field problems including the CH equation, proposed energy stable linear schemes combining with the semi-implicit spectral deferred correction to gain higher order time discretization, while the auxiliary variable is computed coupling with other unknowns. In contrast, our algorithms enable an explicit update for the auxiliary variable, hence more efficient in computation.

1.1. Our Contributions

•

We propose to solve (1.1) by simple IEQ-DG schemes, which integrate the mixed DG method for spatial discretization with the IEQ approach for time discretization, coupled with a crucial spatial projection.
•

We show that the semi-discrete DG scheme features a discrete energy dissipation law if the penalty parameter is suitably large, and present both first and second order (in time) IEQ-DG algorithms. We prove that the IEQ-DG schemes are indeed unconditionally energy stable.
•

We conduct extensive experiments to evaluate the performance of IEQ-DG. First, we present numerical results to show the high order of spatial accuracy of the proposed schemes, the energy dissipating and mass conservative properties of numerical solutions. Second, we conduct experiments on some two-dimensional pattern formation problems, all of which demonstrate the good performance of IEQ-DG.

1.2. Organization

In Section 2, we formulate a unified semi-discrete DG method for the CH equation (1.1) subject to different boundary conditions. In Section 3, we present fully discrete DG schemes and show the energy dissipation and mass conservation properties. In Section 4, we discuss extensions to the case with degenerate mobility and the logarithmic Flory-Huggins potential. In Section 5, we numerically verify the performance of IEQ-DG on different numerical examples. Finally in Section 6 some concluding remarks are given.

2. Spatial DG discretization

Let the domain $\Omega$ be a union of rectangular meshes $\mathcal{T}_{h}=\{K\}:=\bigcup_{\alpha=1}^{\mathcal{N}}K_{\alpha}$ , with $\alpha=(\alpha_{1},\cdots,\alpha_{d}),\ \mathcal{N}=(\mathcal{N}_{1},\cdots,\mathcal{N}_{d})$ and $K_{\alpha}=I_{\alpha_{1}}^{1}\times\cdots\times I_{\alpha_{d}}^{d}$ , where $I_{\alpha_{i}}^{i}=[x_{\alpha_{i}-1/2}^{i},x_{\alpha_{i}+1/2}^{i}]$ for $\alpha_{i}=1,\cdots,\mathcal{N}_{i}$ . Denote by $h_{i}=\max_{1\leq\alpha_{i}\leq N_{i}}|I_{\alpha_{i}}^{i}|$ , with $h=\max_{1\leq i\leq d}h_{i}$ . We denote the set of the interior interfaces by $\Gamma^{0}$ , the set of all boundary faces by $\Gamma^{\partial}$ , and $\Gamma_{h}=\Gamma^{0}\cup\Gamma^{\partial}$ .

The discontinuous Galerkin finite element space can be formulated as

V_{h}=\{v\in L^{2}(\Omega):\ v|_{K}\in P^{k}(K),\ \forall K\in\mathcal{T}_{h}\},

where $P^{k}(K)$ denotes the set of polynomials of degree no more than $k$ on element $K$ . Let $K_{1}$ and $K_{2}$ be two neighboring cells. If the unit normal vector $\nu$ on element interfaces $e\in\partial K_{1}\cap\partial K_{2}$ is oriented from $K_{1}$ to $K_{2}$ , then the average $\{\cdot\}$ and the jump $[\cdot]$ operator are defined by

\{v\}=\frac{1}{2}(v|_{\partial K_{1}}+v|_{\partial K_{2}}),\quad[v]=v|_{\partial K_{2}}-v|_{\partial K_{1}},

for any function $v\in V_{h}$ , where $v|_{\partial K_{i}}\ (i=1,2)$ is the trace of $v$ on $e$ evaluated from element $K_{i}$ .

The direct DG discretization of (1.1) is to find $(u_{h}(\cdot,t),w_{h}(\cdot,t))\in V_{h}\times V_{h}$ such that for all $\phi,\ \psi\in V_{h}$ and $K\in\mathcal{T}_{h}$


$\displaystyle\int_{K}u_{ht}\phi dx=$	$\displaystyle-\int_{K}M(u_{h})\nabla w_{h}\cdot\nabla\phi dx+\int_{\partial K}M(\widehat{u_{h}})\left(\widehat{\partial_{\nu}w_{h}}\phi+(w_{h}-\widehat{w_{h}})\partial_{\nu}\phi\right)ds,$	(2.1a)
$\displaystyle\int_{K}w_{h}\psi dx=$	$\displaystyle\epsilon^{2}\left(\int_{K}\nabla u_{h}\cdot\nabla\psi dx-\int_{\partial K}\widehat{\partial_{\nu}u_{h}}\psi+(u_{h}-\widehat{u_{h}})\partial_{\nu}\psi ds\right)+\int_{K}F^{\prime}(u_{h})\psi dx,$	(2.1b)

where $\nu$ stands for the outward normal direction to $\partial K$ . On each cell interface $e\in\partial K\bigcap\Gamma^{0}$ , the numerical flux is taken as

\displaystyle\widehat{\partial_{\nu}v}=\frac{\beta_{0}[v]}{h_{e}}+\{\partial_{\nu}v\},\ \widehat{v}=\{v\},

(2.2)

for $v=w_{h},\ u_{h}$ , where $\beta_{0}>0$ is a parameter to be determined. Here $h_{e}$ is the characteristic length of interface $e$ . In case of the uniform meshes, we take $h_{e}=h_{i}$ at each interface $x_{\alpha_{i}+1/2}^{i}$ for $\alpha_{i}=0,1,\cdots,\mathcal{N}_{i}$ . The numerical fluxes on $e\in\partial K\bigcap\Gamma^{\partial}$ depend on the boundary conditions. For periodic boundary conditions, the numerical fluxes can take the same formula as those in (2.2). For homogeneous Neumann boundary conditions, the numerical fluxes on the boundary $e\in\partial K\bigcap\Gamma^{\partial}$ are defined as

\widehat{\partial_{\nu}w_{h}}=0,\ \widehat{w_{h}}=w_{h},\ \widehat{\partial_{\nu}u_{h}}=0,\ \widehat{u_{h}}=u_{h}.

(2.3)

Summation of (2.1) over all elements $K\in\mathcal{T}_{h}$ leads to a global DG formulation


$\displaystyle(u_{ht},\phi)=$	$\displaystyle-A(M(u_{h});w_{h},\phi),$	(2.4a)
$\displaystyle(w_{h},\psi)=$	$\displaystyle A(\epsilon^{2};u_{h},\psi)+\left(F^{\prime}(u_{h}),\psi\right),$	(2.4b)

where the bilinear functional is given by

	$\displaystyle A(a(x);q,v)=$	$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}a(x)\nabla q\cdot\nabla vdx+\sum_{e\in\Gamma^{0}}\int_{e}a(x)\left(\widehat{\partial_{\nu}q}[v]+[q]\{\partial_{\nu}v\}\right)ds$		(2.5)
		$\displaystyle+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\int_{e}a(x)\left(\widehat{\partial_{\nu}q}[v]+[q]\{\partial_{\nu}v\}\right)ds,$		(2.5)

where $\tau=1$ for (i) in (1.2) and $\tau=0$ for (ii) in (1.2). Here $a(x)=M(u_{h})$ in $K$ but $M(\widehat{u_{h}})$ for $x\in e$ . Note that the factor $\frac{\tau}{2}$ in (2.5) is used to indicate that for periodic boundary conditions only one end in each direction should be counted. Here each respective type of boundary conditions specified in (1.2) has been taken into account. The initial data for $u_{h}$ is taken from $V_{h}$ so that

\int_{\Omega}(u_{0}(x)-u_{h}(x,0))\phi dx=0,\quad\forall\phi\in V_{h}.

As usual we denote $u_{h}(x,0)=\Pi u_{0}(x)$ , where $\Pi$ is the piecewise $L^{2}$ projection.

We introduce the discrete energy

E(u_{h})=\frac{1}{2}A(\epsilon^{2};u_{h},u_{h})+\int_{\Omega}F(u_{h})dx,

(2.6)

and the notation

\|v\|^{2}_{DG}:=\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla v|^{2}dx+\left(\sum_{e\in\Gamma^{0}}+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\right)\int_{e}\frac{\beta_{0}}{h_{e}}[v]^{2}ds,\quad\forall v\in V_{h}.

(2.7)

We can show that if $\beta_{0}$ is suitably large, the semi-discrete DG scheme (2.4) features a discrete energy dissipation law.

Lemma 2.1.

For piecewise continuous function $a(x)>0$ , there exists $\beta^{*}>0$ such that if $\beta_{0}>\beta^{*}$ , then

A(a(x);v,v)\geq c_{0}\inf_{x\in\Omega}a(x)\|v\|_{\rm DG}^{2}\quad\forall v\in V_{h}.

(2.8)

for some $c_{0}>0$ . As a result, we have

\frac{d}{dt}E(u_{h})=-A(M(u_{h});w_{h},w_{h})\leq 0,\quad\forall t>0.

(2.9)

Proof.

(i) By Young’s inequality, we have

	$\displaystyle A(a(x);v,v)=$	$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}a(x)\|\nabla v\|^{2}dx+\left(\sum_{e\in\Gamma^{0}}+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\right)\int_{e}a(x)[v]\left(\frac{\beta_{0}}{h_{e}}[v]+2\{\partial_{\nu}v\}\right)ds$
	$\displaystyle\geq$	$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}a(x)\|\nabla v\|^{2}dx+\left(\sum_{e\in\Gamma^{0}}+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\right)\int_{e}a(x)\frac{\beta_{0}}{h_{e}}[v]^{2}ds$
		$\displaystyle-\left(\sum_{e\in\Gamma^{0}}+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\right)\left(\frac{\alpha}{h_{e}}\int_{e}a(x)[v]^{2}ds+\frac{h_{e}}{\alpha}\int_{e}a(x)\{\partial_{\nu}v\}^{2}ds\right)$
	$\displaystyle=$	$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}a(x)\|\nabla v\|^{2}dx-\frac{1}{\alpha}\left(\sum_{e\in\Gamma^{0}}+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\right)h_{e}\int_{e}a(x)\{\partial_{\nu}v\}^{2}ds$
		$\displaystyle+\left(\sum_{e\in\Gamma^{0}}+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\right)\left(\frac{\beta_{0}-\alpha}{h_{e}}\right)\int_{e}a(x)[v]^{2}ds$

for $0<\alpha<\beta_{0}$ .

Set

\beta^{*}=\sup_{v\in\tilde{V}_{h}}\frac{\left(\sum_{e\in\Gamma^{0}}+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\right)h_{e}\int_{e}a(x)\{\partial_{\nu}v\}^{2}ds}{\sum_{K\in\mathcal{T}_{h}}\int_{K}a(x)|\nabla v|^{2}dx},

(2.10)

where $\tilde{V}_{h}:=\{v\in V_{h}|\quad v\not=\text{const}\}$ . For given $a(x)$ positive, this quantity $\beta^{*}$ can be shown to be uniformly bounded. In fact, it suffices to bound each local ratio by

\frac{h_{e}\int_{\partial K}a(x)|\partial_{\nu}v|^{2}ds}{\int_{K}a(x)|\nabla v|^{2}dx}\leq\frac{\max_{x\in\partial K}a(x)}{\min_{x\in K}a(x)}\frac{h_{e}\int_{\partial K}|\partial_{\nu}v|^{2}ds}{\int_{K}|\nabla v|^{2}dx},

in which the ratio without weight $a$ is known to be bounded by a constant depending only on $k$ ; see [23, Lemma 3.1]. Thus we have

\displaystyle A(a(x);v,v)\geq\left(1-\frac{\beta^{*}}{\alpha}\right)\sum_{K\in\mathcal{T}_{h}}\int_{K}a(x)|\nabla v|^{2}dx+\left(1-\frac{\alpha}{\beta_{0}}\right)\left(\sum_{e\in\Gamma^{0}}+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\right)\int_{e}a(x)\frac{\beta_{0}}{h_{e}}[v]^{2}ds.

Set $\alpha=\sqrt{\beta_{0}\beta^{*}}$ and $c_{0}=\left(1-\sqrt{\frac{\beta^{*}}{\beta_{0}}}\right)$ , we obtain

A(a(x);v,v)\geq c_{0}\inf_{x\in\Omega}a(x)\|v\|_{\rm DG}^{2}

for $\beta_{0}>\beta^{*}$ .

(ii) Taking $\phi=w_{h},\ \psi=u_{ht}$ in (2.4), then (2.9) follows immediately. ∎

Remark 2.1.

Here $\beta^{*}$ clearly depends on the choice of $a(x)$ unless it is a constant, and the type of boundary conditions through $\tau$ . The fact that $\beta^{*}$ is increasing in $\tau$ implies that for Neumann boundary conditions it suffices to take a smaller $\beta^{*}$ .

3. Time discretization

3.1. The IEQ reformulation

The basic idea of the IEQ methodology [43, 42] is to rewrite the energy functional into a quadratic form

E(u_{h},U)=\frac{1}{2}A(\epsilon^{2};u_{h},u_{h})+\int_{\Omega}U^{2}dx=E(u_{h})+B|\Omega|,

(3.1)

where

U=\sqrt{F(u_{h})+B}

is well-defined when $B$ is chosen so that $F(u_{h})+B>0$ . With the IEQ approach $U$ is updated by solving

U_{t}=\frac{1}{2}H(u_{h})u_{ht},

where

\displaystyle H(w)=\frac{F^{\prime}(w)}{\sqrt{F(w)+B}}.

(3.2)

For the semi-discrete DG scheme (2.4), we follow [26] where an IEQ-DG method was developed for the Swift-Hohenberg equation, to consider the following system: find $(u_{h},w_{h})\in V_{h}\times V_{h}$ and $U(x,t)$ such that


$\displaystyle U_{t}=$	$\displaystyle\frac{1}{2}H(u_{h})u_{ht},$	(3.3a)
$\displaystyle(u_{ht},\phi)=$	$\displaystyle A(M(u_{h});w_{h},\phi),$	(3.3b)
$\displaystyle(w_{h},\psi)=$	$\displaystyle A(\epsilon^{2};u_{h},\psi)+\left(H(u_{h})U,\psi\right),$	(3.3c)

for all $\phi,\psi\in V_{h}$ , subject to initial data

U(x,0)=\sqrt{F(u_{0}(x))+B},\quad u_{h}(x,0)=\Pi u_{0}(x).

Note that with the modified discrete energy (3.1) we still have the following

\frac{d}{dt}E(u_{h},U)=-A(M(u_{h});w_{h},w_{h})\leq 0.

We proceed to discretize (3.3) in time.

3.2. First order fully discrete IEQ-DG scheme

Find $(u^{n}_{h},w_{h}^{n})\in V_{h}\times V_{h}$ and $U^{n}=U^{n}(x)$ such that


$\displaystyle U^{n}_{h}=$	$\displaystyle\Pi U^{n},$	(3.4a)
$\displaystyle\frac{U^{n+1}-U_{h}^{n}}{\Delta t}=$	$\displaystyle\frac{1}{2}H(u_{h}^{n})\frac{u_{h}^{n+1}-u_{h}^{n}}{\Delta t},$	(3.4b)
$\displaystyle\left(\frac{u_{h}^{n+1}-u_{h}^{n}}{\Delta t},\phi\right)=$	$\displaystyle-A(M(u_{h}^{n});w_{h}^{n+1},\phi),$	(3.4c)
$\displaystyle(w_{h}^{n+1},\psi)=$	$\displaystyle A(\epsilon^{2};u_{h}^{n+1},\psi)+\left(H(u_{h}^{n})U^{n+1},\psi\right),$	(3.4d)

for $\phi,\psi\in V_{h}$ .

This scheme admits the following properties.

Theorem 3.1.

Suppose that $M_{\min}\leq M(u)\leq M_{\max}$ for $u$ under consideration, and set

\beta_{0}^{*}=\frac{M_{\max}}{M_{\min}}\sup_{v\in\tilde{V}_{h}}\frac{\left(\sum_{e\in\Gamma^{0}}+\frac{\tau}{2}\sum_{e\in\Gamma^{\partial}}\right)h_{e}\int_{e}\{\partial_{\nu}v\}^{2}ds}{\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla v|^{2}dx}.

(3.5)

If $\beta_{0}>\beta_{0}^{*}$ , scheme (3.4) admits a unique solution $(u_{h}^{n},w_{h}^{n})$ for any $\Delta t>0$ , and the solution $u_{h}^{n}$ satisfies the mass conservation, i.e.,

\int_{\Omega}u_{h}^{n}dx=\int_{\Omega}u_{h}^{0}dx,

(3.6)

for any $n>0$ . Moreover, for $E^{n}:=E(u_{h}^{n},U_{h}^{n})$ we have

	$\displaystyle E^{n+1}\leq E(u_{h}^{n+1},U^{n+1})=$	$\displaystyle E^{n}-\Delta tA(M(u_{h}^{n});w_{h}^{n+1},w_{h}^{n+1})$		(3.7)
		$\displaystyle-\frac{1}{2}A(\epsilon^{2};u_{h}^{n+1}-u_{h}^{n},u_{h}^{n+1}-u_{h}^{n})-\\|U^{n+1}-U_{h}^{n}\\|^{2},$		(3.7)

independent of the size of $\Delta t$ .

Proof.

Taking $\phi=1$ in (3.4c) implies (3.6). We next show the existence and uniqueness of (3.4) at each time step. Substitution of (3.4b) into (3.4c) with (3.4d) gives the following linear system


$\displaystyle(u_{h}^{n+1}/\Delta t,\phi)+A(M(u_{h}^{n});w^{n+1}_{h},\phi)=$	$\displaystyle(u_{h}^{n}/\Delta t,\phi),$	(3.8a)
$\displaystyle A(\epsilon^{2};u_{h}^{n+1},\psi)+\left(1/2H(u_{h}^{n})^{2}u_{h}^{n+1},\phi\right)-(w^{n+1}_{h},\psi)=$	$\displaystyle\left(1/2H(u_{h}^{n})^{2}u_{h}^{n},\phi\right)-(H(u_{h}^{n})U_{h}^{n},\psi).$	(3.8b)

It suffices to prove the uniqueness for this linear system. Denoting $(\tilde{u},\tilde{w})$ the difference of two possible solutions of (3.8) for fixed $(u_{h}^{n},w^{n}_{h})$ , so that


$\displaystyle(\tilde{u}/\Delta t,\phi)+A(M(u_{h}^{n});\tilde{w},\phi)=$	$\displaystyle 0,$	(3.9a)
$\displaystyle A(\epsilon^{2};\tilde{u},\psi)+\left(1/2H(u_{h}^{n})^{2}\tilde{u},\phi\right)-(\tilde{w},\psi)=$	$\displaystyle 0.$	(3.9b)

Setting $\phi=\Delta t\tilde{w},\psi=\tilde{u}$ and adding the two equations, we have

\Delta tA(M(u_{h}^{n});\tilde{w},\tilde{w})+A(\epsilon^{2};\tilde{u},\tilde{u})+\left(1/2H(u_{h}^{n})^{2}\tilde{u},\tilde{u}\right)=0.

Using (2.8) with $a=M(u_{h}^{n})$ so that $\beta_{0}^{*}>\beta^{*}$ , we have

\displaystyle 0\geq

\displaystyle\Delta tM_{\min}\|\tilde{w}\|_{DG}^{2}+\epsilon^{2}\|\tilde{u}\|_{DG}^{2}+\frac{1}{2}\int_{\Omega}H(u_{h}^{n})^{2}\tilde{u}^{2}dx,

which ensures that $\tilde{u}=const$ and $\tilde{w}=const$ . Then it follows $A(M(u_{h}^{n});\tilde{w},\phi)=A(\epsilon^{2};\tilde{u},\psi)=0$ . Thus, (3.9a) is equivalent to

(\tilde{u},\phi)=0,\quad\forall\phi\in V_{h}.

We must have $\tilde{u}=0$ . In a similar fashion, $\tilde{w}=0$ follows from (3.9b). Hence the uniqueness for (3.8) follows.

We next prove (3.7). Taking $\phi=w_{h}^{n+1}$ in (3.4c), $\psi=\frac{u_{h}^{n+1}-u_{h}^{n}}{\Delta t}$ in (3.4d) gives

\displaystyle-A(M(u_{h}^{n});w_{h}^{n+1},w_{h}^{n+1})=A(\epsilon^{2};u_{h}^{n+1},\frac{u_{h}^{n+1}-u_{h}^{n}}{\Delta t})+\left(H(u_{h}^{n})U^{n+1},\frac{u_{h}^{n+1}-u_{h}^{n}}{\Delta t}\right).

By (3.4b) and bilinearity of $A(\epsilon^{2};\cdot,\cdot)$ , the right hand side of the above equation gives

	$\displaystyle RHS=$	$\displaystyle\frac{1}{2\Delta t}\left(A(\epsilon^{2};u_{h}^{n+1},u_{h}^{n+1})-A(\epsilon^{2};u_{h}^{n},u_{h}^{n})+A(\epsilon^{2};u_{h}^{n+1}-u_{h}^{n},u_{h}^{n+1}-u_{h}^{n})\right)$
		$\displaystyle+\frac{1}{\Delta t}\left(\\|U^{n+1}\\|^{2}-\\|U_{h}^{n}\\|^{2}+\\|U^{n+1}-U_{h}^{n}\\|^{2}\right).$

Hence

	$\displaystyle E(u_{h}^{n+1},U^{n+1})=$	$\displaystyle E(u_{h}^{n},U_{h}^{n})-\Delta tA(M(u_{h}^{n});w_{h}^{n+1},w_{h}^{n+1})$		(3.10)
		$\displaystyle-\frac{1}{2}A(\epsilon^{2};u_{h}^{n+1}-u_{h}^{n},u_{h}^{n+1}-u_{h}^{n})-\\|U^{n+1}-U_{h}^{n}\\|^{2}.$		(3.10)

Implied by the fact that $\Pi$ is a contraction mapping in $L^{2}$ , we have

E(u_{h}^{n+1},U_{h}^{n+1})\leq E(u_{h}^{n+1},U^{n+1}),

(3.11)

hence (3.7) as desired.

∎

Remark 3.1.

Note that here $\beta_{0}>\beta_{0}^{*}$ serves only as a sufficient condition for stability. For $M(u)=const$ and rectangular meshes, $\beta_{0}^{*}$ can be more precisely estimated as $\beta_{0}^{*}=k^{2}$ , see [23, Lemma 3.1]. For general case, it suffices to add upon $k^{2}$ a factor of size $\max_{K}\frac{\max_{x\in\partial K}M(u_{h}^{n})}{\min_{x\in K}M(u_{h}^{n})}$ , which approaches to one as the mesh is sufficiently refined. In our numerical examples we take $\beta_{0}^{*}=3k^{2}$ .

3.3. Second order fully discrete IEQ-DG scheme

We first obtain $u_{h}^{1},w_{h}^{1}$ and $U^{1}$ from the first order fully discrete IEQ-DG scheme (3.4). We further use the second order backward differentiation formula (BDF2) for time discretization. In other words, for $n\geq 1$ , the second order fully discrete IEQ-DG scheme is to find $(u^{n+1}_{h},w_{h}^{n+1})\in V_{h}\times V_{h}$ such that for $\forall\phi,\psi\in V_{h}$ ,


$\displaystyle U_{h}^{n}=$	$\displaystyle\Pi U^{n},$	(3.12a)
$\displaystyle\frac{3U^{n+1}-4U_{h}^{n}+U_{h}^{n-1}}{2\Delta t}=$	$\displaystyle\frac{1}{2}H(u^{n,*}_{h})\frac{3u_{h}^{n+1}-4u_{h}^{n}+u_{h}^{n-1}}{2\Delta t},$	(3.12b)
$\displaystyle\left(\frac{3u_{h}^{n+1}-4u_{h}^{n}+u_{h}^{n-1}}{2\Delta t},\phi\right)=$	$\displaystyle-A(M(u^{n,*}_{h});w_{h}^{n+1},\phi),$	(3.12c)
$\displaystyle(w_{h}^{n+1},\psi)=$	$\displaystyle A(\epsilon^{2};u_{h}^{n+1},\psi)+\left(H(u^{n,*}_{h})U^{n+1},\psi\right),$	(3.12d)

where $u^{n,*}_{h}$ is obtained using $u_{h}^{n-1}$ and $u^{n}_{h}$

\displaystyle u^{n,*}_{h}=

\displaystyle 2u_{h}^{n}-u_{h}^{n-1}.

(3.13)

Here instead of $u_{h}^{n+1}$ we use $u^{n,*}_{h}$ to avoid iteration steps in updating the numerical solution, while still maintaining second order accuracy in time.

To show the energy stability, we first present some useful identities.

Lemma 3.1.

For any symmetric bilinear functional $\mathcal{A}(\cdot,\cdot)$ , it follows

	$\displaystyle\mathcal{A}(\phi+\psi,\phi-\psi)$	$\displaystyle=\mathcal{A}(\phi,\phi)-\mathcal{A}(\psi,\psi),$
	$\displaystyle 2\mathcal{A}(\phi_{1},3\phi_{1}-2\phi_{2}-\phi_{3})$	$\displaystyle=\mathcal{A}(\phi_{1},\phi_{1})+\mathcal{A}(2\phi_{1}-\phi_{2},2\phi_{1}-\phi_{2})-\mathcal{A}(\phi_{2},\phi_{2})$
		$\displaystyle\qquad-\mathcal{A}(\phi_{3},\phi_{3})+\mathcal{A}(\phi_{1}-\phi_{3},\phi_{1}-\phi_{3}).$

Proof.

The first identity follows from a direct calculation using the symmetry of the bilinear functional $\mathcal{A}(\cdot,\cdot)$ . The second follows from a proper decomposition and using the first identity, that goes as follows:

	$\displaystyle 2\mathcal{A}(\phi_{1},3\phi_{1}-2\phi_{2}-\phi_{3})=$	$\displaystyle\mathcal{A}(\phi_{1},6\phi_{1}-4\phi_{2}-2\phi_{3})$
	$\displaystyle=$	$\displaystyle\mathcal{A}(\phi_{1},\phi_{1}+4(\phi_{1}-\phi_{2})+\phi_{1}-2\phi_{3})$
	$\displaystyle=$	$\displaystyle\mathcal{A}(\phi_{1},\phi_{1})+\mathcal{A}(2\phi_{1},2\phi_{1}-2\phi_{2})+\mathcal{A}(\phi_{1},\phi_{1}-2\phi_{3})$
	$\displaystyle=$	$\displaystyle\mathcal{A}(\phi_{1},\phi_{1})+\mathcal{A}(2\phi_{1}-\phi_{2},2\phi_{1}-\phi_{2})-\mathcal{A}(\phi_{2},\phi_{2})$
		$\displaystyle+\mathcal{A}(\phi_{1}-\phi_{3},\phi_{1}-\phi_{3})-\mathcal{A}(\phi_{3},\phi_{3}).$

∎

For the scheme (3.12), we have

Theorem 3.2.

Under the assumption in Theorem 3.1, we set $\beta_{0}^{*}$ as defined in (3.5). If $\beta_{0}>\beta_{0}^{*}$ , the second order fully discrete DG scheme (3.12) admits a unique solution $(u_{h}^{n+1},w_{h}^{n+1})$ , and the solution $u_{h}^{n}$ satisfies the mass conservation (3.6) for any $n>0$ . Moreover, for any $\Delta t>0$ it follows

\bar{E}^{n+1}-\bar{E}^{n}\leq-\Delta tA(M(u^{n,*}_{h});w_{h}^{n+1},w_{h}^{n+1})\leq 0,

(3.14)

where the modified energy is defined by

\bar{E}^{n}=\frac{E(u_{h}^{n},U_{h}^{n})+E(u_{h}^{n,*},U_{h}^{n,*})}{2},

(3.15)

with

U_{h}^{n,*}=2U_{h}^{n}-U_{h}^{n-1}.

Proof.

We first prove (3.14). Taking $\phi=2\Delta tw_{h}^{n+1}$ in (3.12c) gives

	$\displaystyle-2\Delta tA(M(u^{n,*}_{h});w_{h}^{n+1},w_{h}^{n+1})=$	$\displaystyle(w_{h}^{n+1},3u_{h}^{n+1}-4u_{h}^{n}+u_{h}^{n-1})$
	$\displaystyle\text{using}\;(\ref{FPDGFull+}d)\quad=$	$\displaystyle A(\epsilon^{2};u_{h}^{n+1},\psi)+\left(H(u^{n,*}_{h})U^{n+1},\psi\right)\quad\psi:=3u_{h}^{n+1}-4u_{h}^{n}+u_{h}^{n-1}$
	$\displaystyle\text{using}\;(\ref{FPDGFull+}b)\quad=$	$\displaystyle A(\epsilon^{2};u_{h}^{n+1},3u_{h}^{n+1}-4u_{h}^{n}+u_{h}^{n-1})+\left(3U^{n+1}-4U_{h}^{n}+U_{h}^{n-1},2U^{n+1}\right)$
	$\displaystyle=$	$\displaystyle A(\epsilon^{2};u_{h}^{n+1},3u_{h}^{n+1}-2u_{h}^{n}-u_{h}^{n,})+\left(3U^{n+1}-2U_{h}^{n}-U_{h}^{n,},2U^{n+1}\right).$

Both $A(\epsilon^{2};\cdot,\cdot)$ and $(\cdot,\cdot)$ are symmetric, by Lemma 3.1 we have

	$\displaystyle A(\epsilon^{2};u_{h}^{n+1},3u_{h}^{n+1}-2u_{h}^{n}-u_{h}^{n,*})=$	$\displaystyle\frac{1}{2}\left[A(\epsilon^{2};u_{h}^{n+1},u_{h}^{n+1})+A(\epsilon^{2};u_{h}^{n+1,},u_{h}^{n+1,})-A(\epsilon^{2};u_{h}^{n},u_{h}^{n})\right.$
		$\displaystyle\left.-A(\epsilon^{2};u_{h}^{n,},u_{h}^{n,})+A(\epsilon^{2};u_{h}^{n+1}-u_{h}^{n,},u_{h}^{n+1}-u_{h}^{n,})\right],$
	$\displaystyle\left(3U^{n+1}-2U_{h}^{n}-U_{h}^{n,*},2U^{n+1}\right)=$	$\displaystyle\\|U^{n+1}\\|^{2}+\\|2U^{n+1}-U_{h}^{n}\\|^{2}-\\|U_{h}^{n}\\|^{2}-\\|U_{h}^{n,}\\|^{2}+\\|U^{n+1}-U_{h}^{n,}\\|^{2}.$

Regrouping, we obtain

		$\displaystyle\frac{1}{2}\left[A(\epsilon^{2};u_{h}^{n+1},u_{h}^{n+1})+A(\epsilon^{2};u_{h}^{n+1,},u_{h}^{n+1,})\right]+\\|U^{n+1}\\|^{2}+\\|2U^{n+1}-U_{h}^{n}\\|^{2}$		(3.16)
		$\displaystyle=2\bar{E}^{n}-2\Delta tA(M(u^{n,}_{h});w_{h}^{n+1},w_{h}^{n+1})-\frac{1}{2}A(\epsilon^{2};u_{h}^{n+1}-u_{h}^{n,},u_{h}^{n+1}-u_{h}^{n,})-\\|U^{n+1}-U_{h}^{n,}\\|^{2}$
		$\displaystyle\leq 2\bar{E}^{n}-2\Delta tA(M(u^{n,*}_{h});w_{h}^{n+1},w_{h}^{n+1}).$

Further use of the fact that $\Pi$ is a contraction mapping in $L^{2}$ , we have

\|U_{h}^{n+1}\|^{2}\leq\|U^{n+1}\|^{2},\quad\|2U_{h}^{n+1}-U_{h}^{n}\|^{2}\leq\|2U^{n+1}-U_{h}^{n}\|^{2}.

Then the left hand side of (3.16) is bounded below by $2\bar{E}^{n+1}$ , thus (3.14) follows.

Similar to the proof of Theorem 3.1, the existence and uniqueness of the scheme (3.12) is equivalent to showing the uniqueness of $u_{h}^{n+1},w_{h}^{n+1}$ given $u_{h}^{i},w_{h}^{i},U^{i}$ with $i=n,n-1$ .

Let $(\tilde{u},\tilde{w},\tilde{U})$ be the difference of two such solutions, then


$\displaystyle\tilde{U}=$	$\displaystyle\frac{1}{2}H(u^{n,*}_{h})\tilde{u},$	(3.17a)
$\displaystyle\left(3\tilde{u},\phi\right)=$	$\displaystyle-2\Delta tA(M(u^{n,*}_{h});\tilde{w},\phi),$	(3.17b)
$\displaystyle(\tilde{w},\psi)=$	$\displaystyle A(\epsilon^{2};\tilde{u},\psi)+\left(H(u^{n,*}_{h})\tilde{U},\psi\right).$	(3.17c)

Setting $\phi=\tilde{w},\ \psi=3\tilde{u}$ , and subtracting (3.17b) from (3.17c), it follows

2\Delta tA(M(u^{n,*}_{h});\tilde{w},\tilde{w})+3A(\epsilon^{2};\tilde{u},\tilde{u})+6\|\tilde{U}\|^{2}=0,

where (3.17a) has been used to simplify the third term. By (2.8) with $a=M(u_{h}^{n})$ , we have $\beta_{0}>\beta_{0}^{*}>\beta^{*}$ , hence

2\Delta tM_{\min}\|\tilde{w}\|_{DG}^{2}+3\epsilon^{2}\|\tilde{u}\|^{2}_{DG}+6\|\tilde{U}\|^{2}\leq 0,

which ensures that $\tilde{u}=const$ , $\tilde{w}=const$ and $\tilde{U}=0$ . Thus, using (3.17) again, we must have $\tilde{u}=\tilde{w}=0$ . Thus leads to the existence and uniqueness of the scheme (3.12).

Taking $\phi=1$ in (3.12c), it follows

\int_{\Omega}u_{h}^{n+1}dx=\frac{1}{3}\int_{\Omega}4u_{h}^{n}-u_{h}^{n-1}dx.

(3.18)

From Theorem 3.1, we have

\int_{\Omega}u_{h}^{1}dx=\int_{\Omega}u_{h}^{0}dx,

(3.19)

which when combined with (3.18) gives the mass conservation (3.6). ∎

3.4. Algorithms

The details related to the schemes implementation are summarized in the following algorithms.

3.4.1. Algorithm for the first order fully discrete IEQ-DG scheme (3.4)

•
Step 1 (Initialization) From the given initial data $u_{0}(x)$
1. (1)
  
  Generate $u_{h}^{0}=\Pi u_{0}(x)\in V_{h}$ ;
2. (2)
  
  Generate $U^{0}=\sqrt{F(u_{0}(x))+B}$ , where $B$ is a priori chosen so that $\inf(F(v)+B)>0$ .
•
Step 2 (Evolution)
1. (1)
  
  Project $U^{n}$ back into $V_{h}$ , namely compute $U^{n}_{h}=\Pi U^{n}$ ;
2. (2)
  
  Solve the linear system (3.8) for $u_{h}^{n+1},w_{h}^{n+1}$ ;
3. (3)
  
  Update $U^{n+1}$ using (3.4b), then return to (1) in Step 2.

3.4.2. Algorithm for the second order fully discrete IEQ-DG scheme (3.12)

•
Step 1 (Initialization) From the given initial data $u_{0}(x)$
1. (1)
  
  Generate $u_{h}^{0}=\Pi u_{0}(x)\in V_{h}$ ;
2. (2)
  
  Generate $U^{0}=\sqrt{F(u_{0}(x))+B}$ , where $B$ is a priori chosen so that $\inf(F(v)+B)>0$ ; and
3. (3)
  
  Solve for $u_{h}^{1},w_{h}^{1}$ and $U^{1}$ for $\forall x\in S$ through Algorithm 3.4.1 for the first order fully discrete IEQ-DG scheme (3.4).

•

Step 2 (Evolution)

(1)

Project $U^{n}$ back into $V_{h}$ , namely compute $U^{n}_{h}=\Pi U^{n}$ ;

(2)

Solve the linear system for $u_{h}^{n+1},w_{h}^{n+1}$ ,

	$\displaystyle\left(\frac{3u_{h}^{n+1}}{2\Delta t},\phi\right)+A(M(u^{n,*}_{h});w_{h}^{n+1},\phi)=$	$\displaystyle\left(\frac{4u_{h}^{n}-u_{h}^{n-1}}{2\Delta t},\phi\right),$
	$\displaystyle A(\epsilon^{2};u_{h}^{n+1},\psi)+\frac{1}{2}\left(\left(H(u^{n,*}_{h})\right)^{2}u_{h}^{n+1},\psi\right)-(w_{h}^{n+1},\psi)=$	$\displaystyle RHS,$

where $RHS=-\left(H(u^{n,*}_{h})\frac{4U_{h}^{n}-U_{h}^{n-1}}{3}-\frac{1}{2}\left(H(u^{n,*}_{h})\right)^{2}\frac{4u_{h}^{n}-u_{h}^{n-1}}{3},\psi\right)$

(3)

Update $U^{n+1}$ through (3.12b), i.e.,

U^{n+1}=\frac{1}{2}H(u^{n,*}_{h})u_{h}^{n+1}+\left(\frac{4U_{h}^{n}-U_{h}^{n-1}}{3}-\frac{1}{2}H(u^{n,*}_{h})\frac{4u_{h}^{n}-u_{h}^{n-1}}{3}\right),

then return to (1) in Step 2.

Remark 3.2.

Higher order (in time) IEQ discretization is possible. We omit the details here due to space limitation. Interested readers are referred to [20] for some arbitrarily high-order linear schemes for gradient flow models.

4. Extensions

It is known that solving the Cahn-Hilliard equation with degenerate mobility and/or logarithmic potential is more difficult since it requires a point-wise control of the numerical solution. We discuss how our schemes can be applied by a proper modification.

4.1. Mobility

Though the mobility is often taken as a constant for simplicity, a thermodynamically reasonable choice is actually the degenerate mobility $M(u)=u(1-u)$ (see e.g., [11]). There is hope that solutions which initially take values in the interval $[0,1]$ will do so for all positive time (which is not true for fourth-order parabolic equations without degeneracy). We remark that only values in the interval $[0,1]$ are physically meaningful. Such degeneracy leads to numerical difficulties.

Here, we follow [11, 3] by considering the modified mobility

\widetilde{M}(u)=\left\{\begin{array}[]{rl}M(\sigma)&u\leq\sigma,\\ M(u)&\sigma<u<1-\sigma,\\ M(1-\sigma)&u\geq 1-\sigma,\\ \end{array}\right.

(4.1)

It is obvious that for given $\sigma$ ,

\widetilde{M}(u)\geq M_{\min}>0,

and it is well-defined for $u\in(-\infty,\infty)$ . Numerically, we apply our scheme using this modified mobility with a small $\sigma$ .

4.2. Flory-Huggins potential

A practical choice for the potential is the Logarithmic Flory-Huggins function [4, 6, 8]

F(v)=\frac{\theta}{2}\left(v\ln v+(1-v)\ln(1-v)\right)+\frac{\theta_{c}}{2}v(1-v),\quad\ v\in[0,1],

(4.2)

where $\theta,\theta_{c}>0$ are physical parameters. This function is non-convex with double wells for $\theta_{c}>2\theta$ , and it only has a single well and admits only a single phase for $\theta_{c}\leq 2\theta$ [37].

The domain of the logarithmic potential (4.2) is $(0,1)$ , which requires the numerical solution be strictly inside $(0,1)$ . For some numerical schemes, such solution bounds can be established (see, e.g., [10, 11, 30, 8, 7]).

For high order DG schemes it is rather difficult to preserve the numerical solution within $(0,1)$ . We choose to regularize the logarithmic Flory-Huggins potential (4.2) by extending its domain from $(0,1)$ to $(-\infty,\infty)$ . Such regularization technique is commonly used to remove the numerical overflow; see, e.g., [8, 1, 11, 2, 42]. Specifically, it can be replaced by the twice continuously differentiable function

\widetilde{F}(v)=\left\{\begin{array}[]{rl}\frac{\theta}{2}\left(v\ln v+(1-v)\ln\sigma+\frac{(1-v)^{2}}{2\sigma}-\frac{\sigma}{2}\right)+\frac{\theta_{c}}{2}v(1-v),&v\geq 1-\sigma,\\ \frac{\theta}{2}\left(v\ln v+(1-v)\ln(1-v)\right)+\frac{\theta_{c}}{2}v(1-v),&\sigma<v<1-\sigma,\\ \frac{\theta}{2}\left((1-v)\ln(1-v)+v\ln\sigma+\frac{v^{2}}{2\sigma}-\frac{\sigma}{2}\right)+\frac{\theta_{c}}{2}v(1-v),&v\leq\sigma,\\ \end{array}\right.

and thus $\widetilde{F}(v)$ is well defined for $v\in(-\infty,\infty)$ . It was argued in [11] that the solution with regularized $\widetilde{M}(u)$ and $\widetilde{F}(u)$ converges to the solution to the original problem as $\sigma\to 0$ . This treatment has been applied in numerical simulations, for example in [3]. In this paper, we apply our IEQ-DG schemes to problems formulated with the modified mobility and the regularized potential.

5. Numerical examples

In this section, we will carry out several numerical tests in both 1D and 2D to demonstrate both temporal and spatial accuracy of the IEQ-DG schemes (3.4) and (3.12), the mass conservation and energy dissipation properties. For the spatial accuracy, we will choose $\Delta t$ sufficiently small such that the spatial discretization error is dominant. Likewise, for the temporal accuracy, we will set spatial meshes sufficiently refined such that temporal discretization error is dominant. In the following numerical examples, the parameter $\beta_{0}=k^{2}+0.5k$ for problems with constant mobility and $\beta_{0}=3k^{2}+0.5k$ for other cases. The parameter $B=1$ as default unless specified.

Example 5.1.

(1D spatial accuracy test) Consider the Cahn-Hilliard equation (1.1) with $M=1$ and double-well potential $F(u)=\frac{1}{4}(u^{2}-1)^{2}$ in $\Omega=[0,2\pi]$ with periodic boundary conditions. Here, we follow Example 5.2 in [33] by adding a source term

s(x,t)=-e^{-t}\sin x\left(3e^{-2t}\cos 2x+3e^{-2t}\cos^{2}x+1\right)

(5.1)

to the Cahn-Hilliard equation (1.1), so that the exact solution is

u(x,t)=e^{-t}\sin x.

(5.2)

We use the fully discrete IEQ-DG scheme (3.12) with a term $\left(s(x,t^{n+1}),\phi\right)$ added to the right hand side of (3.12c), and we test the DG scheme based on $P^{k}$ polynomials, with $k=1,2,3$ . Both errors and orders of accuracy at $T=1$ are reported in Table 1. These results show that $(k+1)$ th order of accuracy in both $L^{2}$ and $L^{\infty}$ norms are obtained.

Table 1. 1D

L^{2},\;L^{\infty}

errors and orders of accuracy at

T=1

$k$	$\Delta t$		N=10	N=20		N=40		N=80
$k$	$\Delta t$		error	error	order	error	order	error	order
1	1e-3	$\\|u-u_{h}\\|_{L^{2}}$	3.09646e-02	8.07876e-03	1.94	2.03575e-03	1.99	5.10124e-04	2.00
1	1e-3	$\\|u-u_{h}\\|_{L^{\infty}}$	1.68270e-02	4.58886e-03	1.87	1.16103e-03	1.98	2.91198e-04	2.00
2	1e-4	$\\|u-u_{h}\\|_{L^{2}}$	3.56585e-04	4.17179e-05	3.10	5.12149e-06	3.03	6.35139e-07	3.01
2	1e-4	$\\|u-u_{h}\\|_{L^{\infty}}$	4.34261e-04	5.50274e-05	2.98	6.89646e-06	3.00	8.63616e-07	3.00
3	1e-5	$\\|u-u_{h}\\|_{L^{2}}$	2.57355e-05	1.66983e-06	3.95	1.05343e-07	3.99	6.62827e-09	3.99
3	1e-5	$\\|u-u_{h}\\|_{L^{\infty}}$	1.86828e-05	1.27898e-06	3.87	8.12040e-08	3.98	5.31052e-09	3.93

Note that the test error in Table 1 is expected to be of order $O(\Delta t^{2}+h^{k+1})$ . In order to observe the desired order in space, time step needs to be smaller so that $\Delta t\leq O(h^{k+1}/2)$ when $k$ gets larger. This comment applies to other cases as well.

Example 5.2.

(2D spatial accuracy test with constant mobility and double-well potential) For the Cahn-Hilliard equation (1.1) with $M(u)=1$ and double-well potential $F(u)=\frac{1}{4}(u^{2}-1)^{2}$ in $\Omega$ with appropriate boundary conditions, we add a source term

s(x,y,t)=-\frac{w(x,y,t)}{4}+\frac{\epsilon^{2}w(x,y,t)}{4}-\frac{3w(x,y,t)v(x,y,t)}{2}+\frac{3w(x,y,t)^{3}}{2}-\frac{w(x,y,t)}{2}

to the right hand side of (1.1), where

	$\displaystyle w(x,y,t)=$	$\displaystyle 0.1e^{-t/4}\sin(x/2)\sin(y/2),$
	$\displaystyle v(x,y,t)=$	$\displaystyle\left(0.1e^{-t/4}\cos(x/2)\sin(y/2)\right)^{2}+\left(0.1e^{-t/4}\sin(x/2)\cos(y/2)\right)^{2},$

so that the exact solution is

u(x,y,t)=w(x,y,t).

Here the parameter $\epsilon=0.1$ . We test this example by DG scheme (3.4) with a term $\left(s(x,y,t^{n+1}),\phi\right)$ added to the right hand side of (3.4c), and the DG scheme is based on polynomials of degree $k$ with $k=1,2,3$ on rectangular meshes.

Test case 1. (Periodic BC) In this test case, we take $\Omega=[0,4\pi]^{2}$ and consider periodic boundary conditions. Both errors and orders of accuracy at $T=0.01$ are reported in Table 2. These results show that $(k+1)$ th order of accuracy in both $L^{2}$ and $L^{\infty}$ are obtained.

Table 2. 2D

L^{2},\;L^{\infty}

errors at

T=0.01

with mesh

N\times N

$k$	$\Delta t$		N=8	N=16		N=32		N=64
$k$	$\Delta t$		error	error	order	error	order	error	order
1	1e-3	$\\|u-u_{h}\\|_{L^{2}}$	3.16822e-02	8.03463e-03	1.98	2.02336e-03	1.99	5.04024e-04	2.01
1	1e-3	$\\|u-u_{h}\\|_{L^{\infty}}$	1.38669e-02	3.74776e-03	1.89	9.59555e-04	1.97	2.40239e-04	2.00
2	1e-4	$\\|u-u_{h}\\|_{L^{2}}$	4.52729e-03	5.75115e-04	2.98	7.33589e-05	2.97	9.21578e-06	2.99
2	1e-4	$\\|u-u_{h}\\|_{L^{\infty}}$	2.32640e-03	2.95229e-04	2.98	4.06866e-05	2.86	5.26926e-06	2.95
3	1e-5	$\\|u-u_{h}\\|_{L^{2}}$	4.46670e-04	2.97916e-05	3.91	1.89117e-06	3.98	1.18585e-07	4.00
3	1e-5	$\\|u-u_{h}\\|_{L^{\infty}}$	3.20555e-04	1.80104e-05	4.15	1.02204e-06	4.14	6.16224e-08	4.05

Test case 2. (Neumann BC) Considering $\Omega=[-\pi,3\pi]^{2}$ with homogenous Neumann boundary conditions (1.2(ii)), both errors and orders of accuracy at $T=0.01$ are reported in Table 3. These results also show $(k+1)$ th order of accuracy in both $L^{2}$ and $L^{\infty}$ .

Table 3. 2D

L^{2},\;L^{\infty}

errors at

T=0.01

with mesh

N\times N

$k$	$\Delta t$		N=8	N=16		N=32		N=64
$k$	$\Delta t$		error	error	order	error	order	error	order
1	1e-3	$\\|u-u_{h}\\|_{L^{2}}$	3.16822e-02	8.03463e-03	1.98	2.02336e-03	1.99	5.04024e-04	2.01
1	1e-3	$\\|u-u_{h}\\|_{L^{\infty}}$	1.38669e-02	3.74776e-03	1.89	9.59555e-04	1.97	2.40239e-04	2.00
2	1e-4	$\\|u-u_{h}\\|_{L^{2}}$	4.52729e-03	5.75115e-04	2.98	7.33591e-05	2.97	9.18427e-06	3.00
2	1e-4	$\\|u-u_{h}\\|_{L^{\infty}}$	2.32640e-03	2.95229e-04	2.98	4.06885e-05	2.86	5.08342e-06	3.00
3	1e-5	$\\|u-u_{h}\\|_{L^{2}}$	4.46670e-04	2.97916e-05	3.91	1.89102e-06	3.98	1.18133e-07	4.00
3	1e-5	$\\|u-u_{h}\\|_{L^{\infty}}$	3.20555e-04	1.80104e-05	4.15	1.02406e-06	4.14	6.40520e-08	4.00

Example 5.3.

(2D spatial accuracy test with constant mobility and logarithmic potential) We consider the Cahn-Hilliard equation (1.1) with constant mobility $M(u)=1$ , the logarithmic Flory-Huggins potential (4.2) with $\theta=\theta_{c}=2$ , the parameters $\epsilon=1$ and $B=10$ . We add an appropriate source term $s(x,y,t)$ to the right hand side of (1.1) such that the exact solution is

u(x,y,t)=\frac{1}{10}e^{-t/4}\sin(x/4)\sin(y/4)+\frac{1}{2}.

We test this example by DG scheme (3.12) with a term $\left(s(x,y,t^{n+1}),\phi\right)$ added to the right hand side of (3.12c), and the DG scheme is also based on polynomials of degree $k$ with $k=1,2,3$ on rectangular meshes.

Test case 1. (Periodic BC) In this test case, we take $\Omega=[0,8\pi]^{2}$ and consider periodic boundary conditions. Both errors and orders of accuracy at $T=0.01$ are reported in Table 4. These results show that $(k+1)$ th order of accuracy in both $L^{2}$ and $L^{\infty}$ are obtained.

Table 4. 2D

L^{2},\;L^{\infty}

errors at

T=0.01

with mesh

N\times N

$k$	$\Delta t$		N=8	N=16		N=32		N=64
$k$	$\Delta t$		error	error	order	error	order	error	order
1	1e-3	$\\|u-u_{h}\\|_{L^{2}}$	6.34010e-02	1.62047e-02	1.97	4.04183e-03	2.00	1.00777e-03	2.00
1	1e-3	$\\|u-u_{h}\\|_{L^{\infty}}$	1.38744e-02	3.74858e-03	1.89	9.55245e-04	1.97	2.39967e-04	1.99
2	1e-4	$\\|u-u_{h}\\|_{L^{2}}$	9.39224e-03	1.18059e-03	2.99	1.46853e-04	3.01	1.83323e-05	3.00
2	1e-4	$\\|u-u_{h}\\|_{L^{\infty}}$	2.45698e-03	3.14143e-04	2.97	3.74571e-05	3.07	4.54860e-06	3.04
3	5e-6	$\\|u-u_{h}\\|_{L^{2}}$	1.09183e-03	6.72768e-05	4.02	4.09870e-06	4.04	2.54225e-07	4.01
3	5e-6	$\\|u-u_{h}\\|_{L^{\infty}}$	2.30167e-04	1.58541e-05	3.86	1.02039e-06	3.96	6.42180e-08	3.99

Test case 2. (Neumann BC) In this test case, we take $\Omega=[-2\pi,2\pi]^{2}$ and consider Neumann boundary conditions. Both errors and orders of accuracy at $T=0.01$ are reported in Table 5. These results show that $(k+1)$ th order of accuracy in both $L^{2}$ and $L^{\infty}$ are obtained.

Table 5. 2D

L^{2},\;L^{\infty}

errors at

T=0.01

with mesh

N\times N

$k$	$\Delta t$		N=8	N=16		N=32		N=64
$k$	$\Delta t$		error	error	order	error	order	error	order
1	1e-3	$\\|u-u_{h}\\|_{L^{2}}$	1.27997e-01	3.55296e-02	1.85	9.55174e-03	1.90	2.13203e-03	2.16
1	1e-3	$\\|u-u_{h}\\|_{L^{\infty}}$	5.54685e-02	1.49970e-02	1.89	3.92808e-03	1.93	9.67492e-04	2.02
2	1e-4	$\\|u-u_{h}\\|_{L^{2}}$	1.87014e-02	2.35480e-03	2.99	2.94393e-04	3.00	3.69614e-05	2.99
2	1e-4	$\\|u-u_{h}\\|_{L^{\infty}}$	1.05130e-02	1.30587e-03	3.01	1.62919e-04	3.00	2.04032e-05	3.00
3	5e-6	$\\|u-u_{h}\\|_{L^{2}}$	2.23974e-03	1.24902e-04	4.16	7.57937e-06	4.04	4.99051e-07	3.92
3	5e-6	$\\|u-u_{h}\\|_{L^{\infty}}$	1.47731e-03	8.22721e-05	4.17	4.36207e-06	4.24	3.29965e-07	3.72

Example 5.4.

(2D spatial accuracy test with degenerate mobility and logarithmic potential) We consider the Cahn-Hilliard equation (1.1) with degenerate mobility $M(u)=u(1-u)$ , the logarithmic Flory-Huggins potential (4.2) with $\theta=\theta_{c}=2$ , the parameters $\epsilon=1$ and $B=10$ . We add an appropriate source term $s(x,y,t)$ to the right hand side of (1.1) such that the exact solution is

u(x,y,t)=\frac{2}{5}e^{-t/4}\sin(x/2)\sin(y/2)+\frac{1}{2}.

We test this example by DG scheme (3.4) with a term $\left(s(x,y,t^{n+1}),\phi\right)$ added to the right hand side of (3.4c), and the DG scheme is also based on polynomials of degree $k$ with $k=1,2,3$ on rectangular meshes.

Test case 1. (Periodic BC) In this test case, we take $\Omega=[0,4\pi]^{2}$ and consider periodic boundary conditions. Both errors and orders of accuracy at $T=0.01$ are reported in Table 6. These results show that $(k+1)$ th order of accuracy in both $L^{2}$ and $L^{\infty}$ are obtained.

Table 6. 2D

L^{2},\;L^{\infty}

errors at

T=0.01

with mesh

N\times N

$k$	$\Delta t$		N=8	N=16		N=32		N=64
$k$	$\Delta t$		error	error	order	error	order	error	order
1	1e-3	$\\|u-u_{h}\\|_{L^{2}}$	1.31235e-01	3.29574e-02	1.99	8.27934e-03	1.99	2.08160e-03	1.99
1	1e-3	$\\|u-u_{h}\\|_{L^{\infty}}$	5.56010e-02	1.49372e-02	1.90	3.81584e-03	1.97	9.59510e-04	1.99
2	1e-4	$\\|u-u_{h}\\|_{L^{2}}$	2.05688e-02	2.51806e-03	3.03	3.05650e-04	3.04	3.79714e-05	3.01
2	1e-4	$\\|u-u_{h}\\|_{L^{\infty}}$	1.13806e-02	1.32194e-03	3.11	1.48147e-04	3.16	1.77820e-05	3.06
3	5e-6	$\\|u-u_{h}\\|_{L^{2}}$	2.82305e-03	1.48385e-04	4.25	8.56909e-06	4.11	5.53886e-07	3.95
3	5e-6	$\\|u-u_{h}\\|_{L^{\infty}}$	1.58906e-03	9.24779e-05	4.10	4.63277e-06	4.32	3.35743e-07	3.79

Test case 2. (Neumann BC) In this test case, we take $\Omega=[-\pi,3\pi]^{2}$ and consider Neumann boundary conditions. Both errors and orders of accuracy at $T=0.01$ are reported in Table 7. These results show that $(k+1)$ th order of accuracy in both $L^{2}$ and $L^{\infty}$ is obtained.

Table 7. 2D

L^{2},\;L^{\infty}

errors at

T=0.01

with mesh

N\times N

$k$	$\Delta t$		N=8	N=16		N=32		N=64
$k$	$\Delta t$		error	error	order	error	order	error	order
1	1e-3	$\\|u-u_{h}\\|_{L^{2}}$	1.31235e-01	3.29574e-02	1.99	8.27934e-03	1.99	2.08160e-03	1.99
1	1e-3	$\\|u-u_{h}\\|_{L^{\infty}}$	5.56010e-02	1.49372e-02	1.90	3.81584e-03	1.97	9.59510e-04	1.99
2	1e-4	$\\|u-u_{h}\\|_{L^{2}}$	2.05688e-02	2.51806e-03	3.03	3.05650e-04	3.04	3.79715e-05	3.01
2	1e-4	$\\|u-u_{h}\\|_{L^{\infty}}$	1.13806e-02	1.32194e-03	3.11	1.48147e-04	3.16	1.77820e-05	3.06
3	5e-6	$\\|u-u_{h}\\|_{L^{2}}$	2.82305e-03	1.48385e-04	4.25	8.56909e-06	4.11	5.59243e-07	3.94
3	5e-6	$\\|u-u_{h}\\|_{L^{\infty}}$	1.58906e-03	9.24779e-05	4.10	4.63278e-06	4.32	3.42344e-07	3.76

Example 5.5.

(Temporal Accuracy Test) Following the test case 2 in Example 5.3, we produce numerical solutions at $T=1$ using DG schemes (3.4) and (3.12) based on $P^{2}$ polynomails with time steps $\Delta t=2^{-m}$ with $2\leq m\leq 5$ and appropriate meshes. The $L^{2},L^{\infty}$ errors and orders of convergence are shown in Table 8, and these results confirm that DG schemes (3.4) and (3.12) are first order and second order in time, respectively.

Table 8.

L^{2},L^{\infty}

errors and EOC at

T=1

with time step

\Delta t

Scheme	Mesh		$\Delta t=2^{-2}$	$\Delta t=2^{-3}$		$\Delta t=2^{-4}$		$\Delta t=2^{-5}$
Scheme	Mesh		error	error	order	error	order	error	order
(3.4)	$32^{2}$	$\\|u-u_{h}\\|_{L^{2}}$	4.21032e-03	2.06620e-03	1.03	1.02380e-03	1.01	5.09705e-04	1.01
(3.4)	$32^{2}$	$\\|u-u_{h}\\|_{L^{\infty}}$	7.48743e-04	3.67246e-04	1.03	1.81917e-04	1.01	9.05192e-05	1.01
(3.4)	$64^{2}$	$\\|u-u_{h}\\|_{L^{2}}$	4.21016e-03	2.06606e-03	1.03	1.02364e-03	1.01	5.09522e-04	1.01
(3.4)	$64^{2}$	$\\|u-u_{h}\\|_{L^{\infty}}$	7.48925e-04	3.67359e-04	1.03	1.82003e-04	1.01	9.05934e-05	1.01
(3.12)	$64^{2}$	$\\|u-u_{h}\\|_{L^{2}}$	1.32995e-03	3.18993e-04	2.06	7.69932e-05	2.05	1.88763e-05	2.03
(3.12)	$64^{2}$	$\\|u-u_{h}\\|_{L^{\infty}}$	2.24427e-04	5.35331e-05	2.07	1.27796e-05	2.07	3.12208e-06	2.03
(3.12)	$128^{2}$	$\\|u-u_{h}\\|_{L^{2}}$	1.34130e-03	3.18230e-04	2.08	7.69137e-05	2.05	1.88416e-05	2.03
(3.12)	$128^{2}$	$\\|u-u_{h}\\|_{L^{\infty}}$	2.27726e-04	5.31968e-05	2.10	1.27386e-05	2.06	3.09279e-06	2.04

Example 5.6.

Following [37], we consider the Cahn-Hilliard equation (1.1) with constant mobility $M(u)=1$ , the logarithmic Flory-Huggins potential

F(u)=600\left(u\ln u+(1-u)\ln(1-u)\right)+1800u(1-u),

and the parameters $\epsilon=1$ and $B=10^{2}$ . The equation is subject to the initial condition

u_{0}(x,y)=\left\{\begin{array}[]{rl}0.71,&(x,y)\in\Omega_{1},\\ 0.69,&(x,y)\in\Omega_{2},\\ \end{array}\right.

where the square domain

\Omega=[-0.5,0.5]\times[-0.5,0.5],\quad\Omega_{1}=[-0.2,0.2]\times[-0.2,0.2],\quad\Omega_{2}=\Omega\backslash\Omega_{1}.

The boundary conditions are taken as Neumann BCs, (ii) in (1.2).

Test case 1. We first solve this problem by the first order fully discrete IEQ-DG scheme (3.4) based on $P^{1}$ and $P^{2}$ polynomials with time step $\Delta t=10^{-7}$ and meshes $40\times 40$ and $80\times 80$ , respectively. The contours at $T=8\times 10^{-5}$ are shown in Figure 1, and the corresponding energy and mass evolutions are shown in Figure 2. From Figure 1, we find that the solution structure is well resolved even on coarser mesh and lower order $P^{1}$ polynomials, and the scheme (3.4) using $P^{2}$ polynomials gives a better resolution than that using $P^{1}$ polynomials on coarser meshes $40\times 40$ , but there is no noticeable difference with solution on refined meshes $80\times 80$ or higher order polynomial $P^{2}$ as shown in Figure 1(b)-(d). The pattern structure is well consistent with that obtained in [37]. Figure 2(a) shows that the numerical solution of the scheme (3.4) satisfies the energy dissipation law, Figure 2(b) and 2(c) show that the numerical solution conserves the total mass $\int_{\Omega}u_{h}^{n}dx=0.6932$ under an appropriate tolerance.

Refer to caption — Figure 1. The contours of numerical solution for the scheme (3.4).

Test case 2. We solve this problem again by the second order fully discrete IEQ-DG scheme (3.12) based on $P^{1}$ and $P^{2}$ polynomials. In Figure 3, we show the contours at $T=8\times 10^{-5}$ obtained based on $P^{1}$ polynomials with mesh $40\times 40$ and time steps $\Delta t=10^{-7},8\times 10^{-8},5\times 10^{-8},2\times 10^{-8}$ , respectively. From Figure 3, we find the pattern structure is comparable to that in Figure 1(b)-(d) even with time step $\Delta t=10^{-7}$ and lower order $P^{1}$ polynomials.

Figure 4 shows that the numerical solution of the scheme (3.12) satisfies the energy dissipation law (3.14), but we do find that the modified energy (3.15) better approximate the original energy with a smaller time step $\Delta t$ , a smaller mesh size $h$ or polynomials of a higher degree. Figure 5 implies the numerical solutions with different time steps $\Delta t$ conserve the total mass $\int_{\Omega}udx=0.6932$ under an appropriate tolerance.

Example 5.7.

Following [37], we further consider the Cahn-Hilliard equation (1.1) with degenerate mobility $M(u)=u(1-u)$ , the logarithmic Flory-Huggins potential

F(u)=3000\left(u\ln u+(1-u)\ln(1-u)\right)+9000u(1-u),

and the parameters $\epsilon=1$ and $B=10^{3}$ . The initial condition is

u_{0}(x,y)=0.63+0.05\text{rand}(x,y),

where $\text{rand}(x,y)$ is the random perturbation function in $[-1,1]$ and has zero mean. For the boundary conditions, we take Neumann BCs (ii) in (1.2).

We solve this problem by the scheme (3.12) based on $P^{2}$ polynomials with meshes $64\times 64$ and time step $\Delta t=10^{-8}$ . The evolution of the concentration field is shown in Figure 6. The corresponding energy and mass evolutions are shown in Figure 7. Figure 6 clearly shows the two phases of the concentration evolution. The first phase is governed by spinodal decomposition and phase separation, which is roughly corresponding to the first three figures of Figure 6, this period is basically terminated as soon as the local concentration is driven to either value of the two binodal points. The second phase is governed by grain coarsening, approximately from $t=8\times 10^{-6}$ onwards the generated patterns cluster and grains tend to coarsen, which is a very slow process. Figure 6 shows statistically similar patterns in the numerical solution as those in [37]. Figure 7 further confirms the numerical solution of the scheme (3.12) satisfies the energy dissipation law and conserves the total mass $\int_{\Omega}udx=0.63$ .

6. Conclusion

In this paper, we integrate the mixed DG method with the IEQ method to design both first and second order fully discrete DG schemes that inherit the energy dissipation law and mass conservation of the continuous equation irrespectively of the mesh and time steps. The spatial discretization is based on the mixed DG method, and the temporal discretization is based on the IEQ approach introduced in [41] for treating nonlinear potentials. Coupled with a spatial projection, the resulting IEQ-DG algorithms are easy to implement without resorting to any iteration method, and proven to be unconditionally energy stable and mass conservative. We have presented numerical examples to verify our theoretical results, and demonstrated the good performance of the scheme in terms of efficiency, accuracy, and preservation of solution properties such as energy dissipation and mass conservation.

Acknowledgments

This work was partially supported by the National Science Foundation under Grant DMS1812666.

References

[1] J. W. Barrett and J. F. Blowey. An error bound for the finite element approximation of the Cahn-Hilliard equation with logarithmic free energy. Numer. Math., 72:1–20, 1995.
[2] J. W. Barrett and J. F. Blowey. Finite element approximation of the Cahn-Hilliard equation with concentration dependent mobility. Math. Comp., 68(22):487–517, 1999.
[3] J. W. Barrett, J. F. Blowey and H. Garcke. Finite element approximation of the Cahn–Hilliard equation with degenerate mobility. SIAM J. Numer. Anal., 37(1):286–318, 1999.
[4] J. F. Blowey and C. M. Elliott. The Cahn-Hilliard gradient theory for phase separation with non-smooth free energy Part I: Mathematical analysis. Euro. J. Appl. Math., 2:233–279, 1991.
[5] J. F. Blowey and C. M. Elliott. The Cahn-Hilliard gradient theory for phase separation with non-smooth free energy Part II: Numerical analysis. Euro. J. Appl. Math., 3:147–179, 1992.
[6] J. W. Cahn and J. E. Hilliard. Free energy of a nonuniform system. I. interfacial free energy. J. Chem. Phys., 28:258–267, 1958.
[7] W. Chen, C. Wang, X. Wang, S. M. Wise. Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential. J. Comput. Phys.: X, 3:100031, 2019.
[8] M. I. M. Copetti and C. M. Elliott. Numerical analysis of the Cahn-Hilliard equation with a logarithmic free energy. Numer. Math., 63(1):39–65, 1992.
[9] Y. Cheng and C.-S. Shu. A discontinuous Galerkin finite element method for time dependent partial differential equations with higher order derivatives. Math. Comp., 77(262):699–730, 2008.
[10] A. Debussche and L. Dettori. On the Cahn-Hilliard equation with a logarithmic free energy. Nonlinear Anal., 24:1491–1514, 1995.
[11] C. M. Elliott and H. Garcke. On the Cahn-Hilliard equation with degenerate mobility. SIAM J. Math. Anal., 27(2):404–423, 1996.
[12] C. M. Elliott and S. Luckhaus. A generalized diffusion equation for phase separation of a multi-component mixture with interfacial energy. SFB 256 Preprint 195. University of Bonn, 1991.
[13] C. M. Elliott and S. Zheng. On the Cahn-Hilliard equation. Arch. Rat. Mech. Anal., 96:339–357, 1986.
[14] D. J. Eyre. Unconditionally gradient stable time marching the Cahn-Hilliard equation. In Computational and mathematical models of microstructural evolution (San Francisco, CA, 1998), volume 529 of Mater. Res. Soc. Sympos. Proc., pages 39–46. MRS, 1998.
[15] X. Feng and O. A. Karakashian. Fully discrete dynamic mesh discontinuous Galerkin methods for the Cahn-Hilliard equation of phase transition. Math. Comp., 76:1093–1117, 2007.
[16] X. Feng, Y. Li and Y. Xing. Analysis of mixed interior penalty discontinuous Galerkin methods for the Cahn-Hilliard equation and the Hele-Shaw flow. SIAM J. Numer. Anal., 54(2):825-–847, 2016.
[17] A. S.-Filibelioǧlu, B. Karasözen and M. Uzunca. Energy stable interior penalty discontinuous Galerkin finite element method for Cahn-Hilliard equation. Int. J. Nonlinear Sci. Numer. Simul., 18(5):303-–314, 2017.
[18] D. Furihata. A stable and conservative finite difference scheme for the Cahn-Hilliard equation. Numer. Math., 87(4):675–699, 2001.
[19] F. Guillen-Gonzalez and G. Tierra. On linear schemes for a Cahn-Hilliard diffuse interface model. J. Comput. Phys., 234:140–171, 2013.
[20] Y. Gong, J. Zhao and Q. Wang. Arbitrarily high-order linear energy stable schemes for gradient flow models. J. Comput. Phys., 419:109610, 2020.
[21] R. Guo and Y. Xu. Efficient solvers of discontinuous Galerkin discretization for the Cahn-Hilliard equations. J. Sci. Comput., 58(2):380–408, 2014.
[22] R. Guo and Y. Xu. Semi-implicit spectral deferred correction method based on the invariant energy quadratization approach for phase field problems. Commun. Comput. Phys., 26(1):87–113, 2019.
[23] H. Liu. Optimal error estimates of the direct discontinuous Galerkin method for convection-diffusion equations. Math. Comp., 84:2263–2295, 2015.
[24] H. Liu and J. Yan. The direct discontinuous Galerkin (DDG) method for diffusion with interface corrections. Commun. Comput. Phys., 8(3):541–564, 2010
[25] H. Liu and P. Yin. A mixed discontinuous Galerkin method without interior penalty for time-dependent fourth order problems. J. Sci. Comput., 77(1):467-501, 2018.
[26] H. Liu and P. Yin. Unconditionally energy stable DG schemes for the Swift-Hohenberg equation. J. Sci. Comput., 81(2):789-819, 2019.
[27] H. Liu and P. Yin. On the SAV-DG method for a class of fourth order gradient flows. arXiv preprint, arXiv:2008.11877, 2020.
[28] H. Liu and P. Yin. Energy stable Runge-Kutta discontinuous Galerkin schemes for fourth order gradient flows. arXiv preprint, arXiv:2101.00152, 2020.
[29] G. B. McFadden. Phase field models of solidification. Contemporary Mathematics, 295:107–145, 2002.
[30] A. Miranville and S. Zelik. Robust exponential attractors for Cahn-Hilliard type equations with singular potentials. Math. Meth. Appl. Sci., 27:545–582, 2004.
[31] B. Perthame, A. Poulain. Relaxation of the Cahn-Hilliard equation with degenerate double-well potential ad degenerate mobility. arXiv preprint, arXiv:1908.11294.
[32] J. Shen and X. Yang. Numerical approximations of Allen-Cahn and Cahn-Hilliard equations. Discrete Contin. Dyn. Syst. A, 28:1669–1691, 2010.
[33] H. Song and C.-W. Shu. Unconditional energy stability analysis of a second order implicit-explicit local discontinuous Galerkin method for the Cahn-Hilliard equation. J. Sci. Comput., 73:1178–1203, 2017.
[34] Z.-Z. Sun. A second-order accurate linearized difference scheme for the two-dimensional Cahn-Hilliard equation. Math. Comp., 64:1463–1471, 1995.
[35] J. Shen, J. Xu, and J. Yang. A new class of efficient and robust energy stable schemes for gradient flows. SIAM Rev., 61: 474–506, 2019.
[36] J. Shen, J. Xu and X. Yang. The scalar auxiliary variable (SAV) approach for gradient flows. J. Comput. Phys., 352:407–417, 2018.
[37] G. N. Wells. E. Kuhl and K. Garikipati. A discontinuous Galerkin method for the Cahn-Hilliard equation. J. Comput. Phys., 218:860–877, 2006.
[38] X. Wu, G. J. v. Zwieten and K. G. v. d. Zee. Stabilized second-order convex splitting schemes for Cahn-Hilliard models with application to diffuse-interface tumor-growth models. Int. J. Numer. Meth. Biomed. Engng., 30:180–203, 2014.
[39] Y. Xia, Y. Xu and C.-W. Shu. Local discontinuous Galerkin methods for the Cahn-Hilliard type equations. J. Comput. Phys., 227:472–491, 2007.
[40] C. Xu and T. Tang. Stability analysis of large time-stepping methods for epitaxial growth models. SIAM. J. Num. Anal., 44:1759–1779, 2006.
[41] X. Yang. Linear, first and second order and unconditionally energy stable numerical schemes for the phase field model of homopolymer blends. J. Comput. Phys., 302: 509–523, 2016.
[42] X. Yang, J. Zhao. On linear and unconditionally energy stable algorithms for variable mobility Cahn-Hilliard type equation with logarithmic Flory-Huggins potential. Commun. Comput. Phys., 25(3):703–728, 2019.
[43] X. Yang, J. Zhao, Q. Wang and J. Shen. Numerical approximations for a three components Cahn-Hilliard phase-field model based on the invariant energy quadratization method. Math. Models Methods Appl. Sci., 27(11):1993–2030, 2017.
[44] J. Yin. On the existence of nonnegative continuous solutions of the Cahn-Hilliard equation. J. Differential Equations, 97:310–327, 1992.
[45] J. Zhao, Q. Wang and X. Yang. Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach. Int. J. Numer. Methods Eng., 110(3):279–300, 2017.

	$\displaystyle A(\epsilon^{2};u_{h}^{n+1},3u_{h}^{n+1}-2u_{h}^{n}-u_{h}^{n,*})=$	$\displaystyle\frac{1}{2}\left[A(\epsilon^{2};u_{h}^{n+1},u_{h}^{n+1})+A(\epsilon^{2};u_{h}^{n+1,},u_{h}^{n+1,})-A(\epsilon^{2};u_{h}^{n},u_{h}^{n})\right.$
		$\displaystyle\left.-A(\epsilon^{2};u_{h}^{n,},u_{h}^{n,})+A(\epsilon^{2};u_{h}^{n+1}-u_{h}^{n,},u_{h}^{n+1}-u_{h}^{n,})\right],$
	$\displaystyle\left(3U^{n+1}-2U_{h}^{n}-U_{h}^{n,*},2U^{n+1}\right)=$	$\displaystyle\\|U^{n+1}\\|^{2}+\\|2U^{n+1}-U_{h}^{n}\\|^{2}-\\|U_{h}^{n}\\|^{2}-\\|U_{h}^{n,}\\|^{2}+\\|U^{n+1}-U_{h}^{n,}\\|^{2}.$