Improving the accuracy and consistency of the energy quadratization method with an energy-optimized technique

Abstract

We propose an energy-optimized invariant energy quadratization method to solve the gradient flow models in this paper, which requires only one linear energy-optimized step to correct the auxiliary variables on each time step. In addition to inheriting the benefits of the baseline and relaxed invariant energy quadratization method, our approach has several other advantages. Firstly, in the process of correcting auxiliary variables, we can directly solve linear programming problem by the energy-optimized technique, which greatly simplifies the nonlinear optimization problem in the previous relaxed invariant energy quadratization method. Secondly, we construct new linear unconditionally energy stable schemes by applying backward Euler formulas and Crank-Nicolson formula, so that the accuracy in time can reach the first- and second-order. Thirdly, comparing with relaxation technique, the modified energy obtained by energy-optimized technique is closer to the original energy, meanwhile the accuracy and consistency of the numerical solutions can be improved. Ample numerical examples have been presented to demonstrate the accuracy, efficiency and energy stability of the proposed schemes.

Invariant energy quadratization ,Gradient flow models ,energy-optimized technique ,Unconditionally energy stable ,Numerical experiments 65M12 ,35K20 ,35K35 ,35K55 ,65Z05

1 Introduction↩︎

Many physical problems such as interfacial dynamics , crystal growth , image inpainting , thin films , and polymers [9] can be modeled by gradient flow models which comply with the second law of thermodynamics. As a special kind of nonlinear partial differential equations, the gradient flow model can be derived from the mobility and the variational derivative of free energy when the total free energy is known. Generally, consider the spatial-temporal domain $$\Omega_t:=\Omega \times \left(0,T\right]$$, $$\Omega \subset \mathbb{R}^d, d=1,2,3$$ and the state variable $$\phi =\phi \left(\boldsymbol{x} , t\right) :\Omega _t\mapsto \mathbb{R}$$, the total free energy $$E(\phi)$$ consists of a quadratic term $$E_{\mathcal{L}}(\phi)$$ and a strongly nonlinear term $$E_{\mathcal{N}}(\phi)$$, namely \begin{align} \label{E0} E\left( \phi \right) =E_{\mathcal{L}}\left( \phi \right) +E_{\mathcal{N}}\left( \phi \right) =\frac{1}{2}\int_{\Omega}{\phi \mathcal{L}\phi}\, \text{d}\boldsymbol{x}+\int_{\Omega}{F\left( \phi \right)}\,\text{d}\boldsymbol{x}, \end{align}\tag{1} where $$F(\phi)$$ is a nonlinear potential function, and $$\mathcal{L}$$ is a linear self-adjoint non-negative definite operator. Then the corresponding general gradient flow system of 1 reads \label{PF} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\mathcal{G}\mu ,\\ {}&\mu =\mathcal{L}\phi +f\left( \phi \right) , \end{align}\tag{2} with the initial condition $$\phi(\boldsymbol{x},0)$$ and thermodynamically consistent boundary conditions (such as periodic and homogeneous Neumann boundaries), where $$f(\phi)=F'(\phi)$$, $$\mu =\frac{\delta E}{\delta \phi}$$ is the chemical potential, $$M$$ is a mobility constant and $$\mathcal{G}$$ is a non-negative definite operator reflecting the dissipation path. The gradient flow system $$\eqref{PF}$$ preserves the energy dissipation law \begin{align} \label{dE} \frac{dE\left( \phi \right)}{dt}=\left( \frac{\delta E}{\delta \phi},\frac{\partial \phi}{\partial t} \right) =\left( \mathcal{L}\phi +f\left( \phi \right) ,-M\mathcal{G}\mu \right) =-\left( \mu ,M\mathcal{G}\mu \right) \le 0, \end{align}\tag{3} where $$\left( \cdot ,\cdot \right)$$ denotes the $$L^2$$ product, i.e. $$\left( f,g \right) =\int_{\Omega}{fg}\,\text{d}\boldsymbol{x}, \forall f,g\in L^2\left( \Omega \right)$$.

In order to eliminate non-physical numerical solutions, the energy dissipation law must be preserved at the discrete level when designing numerical schemes for nonlinear gradient flow models. In recent years, there has been tremendous interest in designing efficient and accurate energy stable schemes for the nonlinear dissipation systems. To this end, many efficient schemes have been developed, including convex splitting method , the exponential time differencing (ETD) method, stabilized semi-implicit (SSI) approach , and the new Lagrange multiplier method [16]. Recently, a numerical approach named invariant energy quadratization (IEQ) or energy quadratization (EQ) method is proposed to solve thermodynamically consistent PDE problems, and by utilizing auxiliary variables, the original PDE problems are reformulated into equivalent PDE problems. Based on the baseline IEQ method, linear, second-order even any high-order, easy-to-implement, and unconditionally energy stable schemes are constructed for complex gradient flow models . Subsequently, the scalar auxiliary variable (SAV)-type approaches are proposed and can also be applied to solve a class of gradient flow models.

However, for the baseline IEQ method, the truncation errors are introduced in the long-time numerical simulation process, so that the numerical solutions of auxiliary variables are no longer equivalent to their original continuous definitions. We know that the baseline IEQ method preserves the modified energy dissipation law with respect to the auxiliary variables instead of the original variables, which means the original energy corresponding to the numerical solutions is not necessarily monotonically decreasing. This inconsistency also leads to a loss of accuracy when the time step is not small enough. In order to overcome the adverse effects of this inconsistency, the pioneering relaxation techniques are used to correct the numerical solutions of the auxiliary variables, so that the modified energy is closer to the original energy. In this paper, inspired by a novel and essential energy-optimal technique based on the SAV method, i.e. EOP-SAV method [33], we propose an energy-optimized approach based on the baseline IEQ method, and construct first- and second-order (in time) unconditionally energy stable schemes which can solve a large class of gradient flow problems. We effectively penalize the inconsistency between the numerical solution of auxiliary variables and their continuous definitions by using the proposed energy-optimized IEQ (EOP-IEQ) method and the main advantages of our approach are listed as follows:

• The EOP-IEQ method inherits all the advantages of the baseline IEQ method [17] and relaxed IEQ (REQ) method [30];

• Our energy-optimized technique makes the calculation process of updating auxiliary variables more simpler and efficient;

• Compared with the IEQ and REQ methods, the modified energy obtained by the energy-optimized technique is closer to the original energy.

All numerical schemes based on the EOP-IEQ method are linear unconditionally energy stable and easy-to-implement. In numerical examples of various gradient flow models, the modified energy corresponding of the numerical solutions based on the auxiliary variables almost satisfy the original energy dissipation law.

The rest of this paper is organized as follows. In Section 2, we briefly review the baseline IEQ and REQ methods for gradient flow models. In Section 3, we propose the EOP-IEQ method for general dissipative systems, and prove that the new first- and second-order (in time) numerical schemes are unconditionally energy stable. In Section 4, we verify its optimality by comparing the EOP-IEQ method with the baseline IEQ and REQ methods, and provide a large number of numerical experiments to verify the accuracy, efficiency and energy stability of the proposed schemes.

2 A brief review of the IEQ-type methods↩︎

In this section, we review the baseline IEQ method which considered by Yang in [17] and the relaxed IEQ (REQ) approach proposed in [30] by Zhao.

2.1 The baseline IEQ method↩︎

Firstly, we need to assume that the bulk free energy $$F(\phi)$$ is bounded from below which means $$F(\phi) > - C_*$$ for a positive constant $$C_*$$. Introduce an auxiliary variable \begin{align} \label{eq95q} q\left( \boldsymbol{x},t \right) := Q\left( \phi\left( \boldsymbol{x},t \right) \right) =\sqrt{F\left( \phi\left( \boldsymbol{x},t \right) \right) +C_0},\quad C_0>C_*, \end{align}\tag{4} where the positive constant $$C_0$$ makes sure $$q(\boldsymbol{x},t)$$ a well defined real function. Then we rewrite the gradient flow system 2 as the following equivalent system \label{IEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\mathcal{G}\mu ,\\ {}&\mu =\mathcal{L}\phi +\frac{q}{Q(\phi)}f\left( \phi \right) ,\\ {}&\frac{\partial q}{\partial t}=\frac{f\left( \phi \right)}{2Q(\phi)}\frac{\partial \phi}{\partial t}. \end{align}\tag{5} Using the IEQ method above, we convert the free energy 1 into the following quadratic form \begin{align} \label{E_{1}^{n+1}} \bar{E}( \phi ,q ) =\frac{1}{2}( \mathcal{L}\phi ,\phi ) +(q,q)-C_0|\Omega|, \end{align}\tag{6} It is not difficult to obtain the following modified energy dissipation law \begin{align} \label{dE_{1}^{n+1}} \frac{d\bar{E}\left( \phi,q \right)}{dt}=-M\left( \mathcal{G}\mu ,\mu \right) \le 0. \end{align}\tag{7}

Remark 1. With the IEQ method, numerical algorithms can be introduced to solve the equivalent model 5 instead of solving the original model 2 .

Now we introduce some notations that will be used throughout the paper. Let $$L^p(\Omega)$$ denote the usual Lebesgue space on $$\Omega$$ with the norm $$\|\cdot\|_{L^p}$$. $$W^{k,p}(\Omega)$$ stands for the standard Sobolev spaces equipped with the standard Sobolev norms $$\|\cdot\|_{k,p}$$. For $$p = 2$$, we write $$H^k(\Omega)$$ for $$W^{k,2}(\Omega)$$ and the corresponding norm is $$\|\cdot\|_k$$, where $$\|\cdot\|_0$$ is simply represented by $$\|\cdot\|$$.

We discretize time domain $$[0,T]$$ into uniform meshes $$0=t_0<t_1<\cdots<t_N=T$$ with $$t_n=n\tau$$, $$\tau=T/N$$, and use $$(\bullet )^{n+1}$$ to denote the numerical approximation of $$(\bullet )$$ at $$t_{n+1}$$. With above notations, we obtain the following semi-discrete $$k$$th-order IEQ/BDF$$k$$ schemes based on the backward Euler formulas.

Corollary 1 (IEQ/BDF$$k$$ schemes [17]). For initial datas, we set $$\phi ^0=\phi \left( \boldsymbol{x},0 \right)$$ and $$q^0=\sqrt{F\left( \phi ^0 \right) +C_0}$$, we can update $$(\phi^{n+1}, q^{n+1})$$ via the following step: \label{IEQ1} \begin{align} {}&\frac{\alpha _k\phi ^{n+1}-A_k\left( \phi ^n \right)}{\tau}=-M\mathcal{G}\mu ^{n+1},\\ {}&\mu ^{n+1}=\mathcal{L}\phi ^{n+1}+\frac{q^{n+1}}{Q\left( B_k\left( \phi ^n \right) \right)}f\left( B_k\left( \phi ^n \right) \right) ,\\ {}&\frac{\alpha _kq^{n+1}-A_k\left( q^n \right)}{\tau}=\frac{f\left( B_k\left( \phi ^n \right) \right)}{2Q\left( B_k\left( \phi ^n \right) \right)}\frac{\alpha _k\phi ^{n+1}-A_k\left( \phi ^n \right)}{\tau}, \end{align}\qquad{(1)} where $$Q\left( B_k\left( \phi ^n \right) \right) =\sqrt{F\left( B_k\left( \phi ^n \right) \right) +C_0}$$. The values of $$\alpha_k$$, $$A_k(\phi^n)$$ and $$B_k(\phi^n)$$ of $$k=1,2$$ are as follows:

• First-order: \begin{align} \alpha _1=1,\quad A_1\left( \phi ^n \right) =\phi ^n,\quad B_1\left( \phi ^n \right) =\phi ^n. \end{align}

• Second-order: \begin{align} \alpha _2={}&\frac{3}{2},\quad A_2\left( \phi ^n \right) =2\phi ^n-\frac{1}{2}\phi ^{n-1},\quad B_2\left( \phi ^n \right) =2\phi ^n-\phi ^{n-1}. \end{align}

Lemma 1 ([17]). The IEQ/BDF$$k$$ $$(k=1,2)$$ schemes ?? are unconditionally energy stable in the sense that

• First-order :*

$\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert q^{n+1} \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) -\lVert q^n \rVert ^2\le -\tau \left( \mathcal{G}\mu ^{n+1},\mu ^{n+1} \right)\le 0 .$

• Second-order :*

\begin{align} {}&\frac{1}{4}\left[ \left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\left( \mathcal{L}\left( 2\phi ^{n+1}-\phi ^n \right) ,2\phi ^{n+1}-\phi ^n \right) \right] \!\\ {}&+\frac{1}{2}\left( \lVert q^{n+1} \rVert ^2+\lVert 2q^{n+1}-q^n \rVert ^2 \right)\\ {}&-\frac{1}{4}\left[ \left( \mathcal{L}\phi ^n,\phi ^n \right) +\left( \mathcal{L}\left( 2\phi ^n-\phi ^{n-1} \right) ,2\phi ^n-\phi ^{n-1} \right) \right]\\ {}&-\frac{1}{2}\left( \lVert q^n \rVert ^2+\lVert 2q^n-q^{n-1} \rVert ^2 \right)\\ \le{}& -\tau M\left( \mathcal{G}\mu ^{n+1},\mu ^{n+1} \right) \le 0. \end{align}

2.2 The REQ method↩︎

To overcome the inconsistency issue of the modified energy and original energy at the discrete level, Zhao [30] proposes a relaxation technique based on the baseline IEQ method, which called REQ method to overcome the numerical difference between $$q^{n+1}$$ and $$Q(\phi^{n+1})$$. The semi-implicit relaxed IEQ Crank-Nicolson (REQ/CN) scheme is as follows.

Corollary 2 (REQ/CN scheme [30]). For initial datas, we set $$\phi ^0=\phi \left( \boldsymbol{x},0 \right)$$ and $$q^0=\sqrt{F\left( \phi ^0 \right) +C_0}$$, then we update $$\phi^{n+1}$$, $$q^{n+1}$$ via the following two steps:

• $$\boldsymbol{Step~1} :$$ To obtain the intermediate solutions $$(\phi^{n+1}, \hat{q}^{n+1})$$ by the baseline IEQ method: \begin{align} {}&\frac{\phi ^{n+1}-\phi ^n}{\tau}=-M\mathcal{G}\mu ^{n+\frac{1}{2}},\\ {}&\mu ^{n+\frac{1}{2}}=\mathcal{L}\frac{\phi ^{n+1}+\phi ^n}{2}+\frac{\hat{q}^{n+1}+q^n}{2}\bar{H}^{n+\frac{1}{2}},\\ {}&\frac{\hat{q}^{n+1}-q^n}{\tau}=\frac{\bar{H}^{n+\frac{1}{2}}}{2}\frac{\phi ^{n+1}-\phi ^n}{\tau}, \end{align} where $$\bar{H}^{n+\frac{1}{2}}=\frac{f\left( \bar{\phi}^{n+\frac{1}{2}} \right)}{Q\left( \bar{\phi}^{n+\frac{1}{2}} \right)}$$, $$\bar{\phi}^{n+\frac{1}{2}}=\frac{3}{2}\phi ^n-\frac{1}{2}\phi ^{n-1}$$.

• $$\boldsymbol{Step~2} :$$ To update the intermediate solution $$\hat{q}^{n+1}$$ with a relaxation step, we update $$q^{n+1}$$ as: \label{REQ47CN-q} \begin{align} {}&q^{n+1}=\xi ^{n+1}\hat{q}^{n+1}+\left( 1-\xi ^{n+1} \right) Q(\phi^{n+1}),\quad\xi ^{n+1}\in V_0\cap V_1,\\ {}& V_0=\left\{ \xi \left| \xi \in \left[ 0,1 \right] \right. \right\}, \quad V_1=\left\{ \xi \left| a\xi ^2+b\xi +c \right. \le 0 \right\} , \end{align}\qquad{(2)} where the coefficients $$a$$, $$b$$ and $$c$$ are given by \begin{align} {}& a=\lVert \hat{q}^{n+1}-Q(\phi^{n+1}) \rVert ^2,\quad b=2\left( \hat{q}^{n+1}-Q\left( \phi ^{n+1} \right),Q\left( \phi ^{n+1} \right) \right), \\ {}& c=-\tau \eta \left( \mathcal{G}\mu ^{n+\frac{1}{2}},\mu ^{n+\frac{1}{2}} \right) +\lVert Q(\phi^{n+1}) \rVert ^2-\lVert \hat{q}^{n+1} \rVert ^2,\quad \eta\in[0,1]. \end{align}

Remark 2 (Optimal choice for $$\xi^{n+1}$$ [30]). Here $$\eta$$ is an artificial parameter that can be assigned. If $$a=0$$, we set $$\xi^{n+1}=0$$ which means that $$q^{n+1}=Q(\phi^{n+1})$$, and if $$a>0$$, we take $$\xi ^{n+1}=\max \left\{ 0,(-b-\sqrt{b^2-4ac})/(2a) \right\}$$.

::: {#th-E_{1}^{n+1} .lemma} Lemma 2 ([30]). The REQ/CN scheme 2 is unconditionally energy stable in the sense that $\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert q^{n+1} \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) -\lVert q^n \rVert ^2=-\tau M \left( \mathcal{G}\mu ^{n+\frac{1}{2}},\mu ^{n+\frac{1}{2}} \right) \le 0.$ :::

3 The EOP-IEQ method↩︎

The major issue is that the quxiliary variable $$q(t_{n+1})$$ is no longer equal to $$Q(\phi(t_{n+1}))$$ numerically. Thus the modified energy law [dE_{1}^{n+1}] does not necessarily satisfy the original energy law 3 at the discrete level. To better restrain the inconsistency issue between $$q(t_{n+1})$$ and $$Q(\phi(t_{n+1}))$$ (that are supposed to be equal as introduced in 4 ), Liu et.al. propose a novel modified SAV method [33], named energy-optimized SAV (EOP-SAV) approach which is unconditionally energy stable with respect to a modified energy that is closer to the original energy than the baseline SAV[26] and relaxed SAV [31] approaches. Inspired by this method, the EOP-IEQ method is considered in this paper, which has the optimal correction on modified energy compared to the IEQ and REQ methods. The main idea of the energy-optimized (EOP) technique is to correct the auxiliary variable $$q^n$$ from the perspective of original nonlinear energy $$E_{\mathcal{N}}\left( \phi^n \right)$$, so that the modified energy is closest to the original energy at the discrete level, and the numerical schemes preserve the energy dissipation law.

3.1 The EOP-IEQ/CN scheme.↩︎

Now we consider the following second-order Crank-Nicolson scheme based on the novel EOP-IEQ approach.

Corollary 3 (EOP-IEQ/CN scheme). For initial datas, we set $$\phi ^0=\phi \left( \boldsymbol{x},0 \right)$$ and $$q^0=\sqrt{F\left( \phi ^0 \right) +C_0}$$, then we update $$\phi^{n+1}$$, $$q^{n+1}$$ via the following two steps:

• $$\boldsymbol{Step~1} :$$ To obtain the intermediate solutions $$(\phi^{n+1},\hat{q}^{n+1})$$ by the baseline IEQ method: \label{EOP47CN1} \begin{align} {}&\frac{\phi ^{n+1}-\phi ^n}{\tau}=-M\mathcal{G}\mu ^{n+\frac{1}{2}},\\ {}&\mu ^{n+\frac{1}{2}}=\mathcal{L}\frac{\phi ^{n+1}+\phi ^n}{2}+\frac{\hat{q}^{n+1}+q^n}{2}\bar{H}^{n+\frac{1}{2}},\\ {}&\frac{\hat{q}^{n+1}-q^n}{\tau}=\frac{\bar{H}^{n+\frac{1}{2}}}{2}\frac{\phi ^{n+1}-\phi ^n}{\tau}, \end{align}\qquad{(3)} where $$\bar{H}^{n+\frac{1}{2}}=\frac{f\left( \bar{\phi}^{n+\frac{1}{2}} \right)}{Q\left( \bar{\phi}^{n+\frac{1}{2}} \right)}$$, $$\bar{\phi}^{n+\frac{1}{2}}=\frac{3}{2}\phi ^n-\frac{1}{2}\phi ^{n-1}$$.

• $$\boldsymbol{Step~2} :$$ To optimize the intermediate solution $$\hat{q}^{n+1}$$ with an energy-optimized step, we update the auxiliary variable $$q^{n+1}$$ as $\label{EOP47CN2} q^{n+1}=\lambda ^{n+1}Q(\phi^{n+1}),\quad \lambda ^{n+1}=\min \left\{ 1,\sqrt{\frac{E_{2}^{n+1}}{E_{1}^{n+1}}} \right\},\qquad{(4)}$ where the coefficients $$E_{1}^{n+1}$$ and $$E_{2}^{n+1}$$ are given by \begin{align} E_{1}^{n+1}=\lVert Q\left( \phi ^{n+1} \right) \rVert ^2,\quad E_{2}^{n+1}=\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) +\lVert q^n \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) . \end{align}

Remark 3. In the computational simulation, the value of the intermediate variable $$\hat{q}^{n+1}$$ in $$\boldsymbol{Step~1}$$ does not need to be solved, it is only used in the following numerical analysis. According to the energy-optimized technique in $$\boldsymbol{Step~2}$$, the numerical solutions of $$q^{n+1}$$ are always non-negative. Furthermore, the elements discretized by $$q^{n+1}$$ in space are all non-negative at the fully discrete level.

::: {#th-E_{2}^{n+1} .theorem} Theorem 1. The energy-optimized technique of $$q^{n+1}$$ ?? in $$\boldsymbol{Step~2}$$ is an optimal choice to correct the modified energy dissipative law of the baseline IEQ method [17]. :::

Proof. Taking the inner products of ?? with $$\mu ^{n+\frac{1}{2}}$$, $$\frac{\phi ^{n+1}-\phi ^n}{\tau}$$ and $$\hat{q}^{n+1}+q^n$$ respectively, we obtain

\begin{align} \label{qCN} \lVert \hat{q}^{n+1} \rVert ^2=\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) +\lVert q^n \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) -\tau M\left( \mathcal{G}\mu ^{n+\frac{1}{2}},\mu ^{n+\frac{1}{2}} \right) . \end{align}\tag{8} As far as the energy dissipation law is concerned, if we correct $$q^{n+1}$$, at least the following inequality holds \begin{align} \label{qEdiss} \frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert q^{n+1} \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) -\lVert q^n \rVert ^2\le 0. \end{align}\tag{9} Then we have the following inequality \begin{align} \label{qE} 0\le \lVert q^{n+1} \rVert ^2\le \frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) +\lVert q^n \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right). \end{align}\tag{10} Let $$E_{1}^{n+1}=\lVert Q\left( \phi ^{n+1} \right) \rVert ^2$$, $$E_{2}^{n+1} = \frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) +\lVert q^n \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right)$$, we have

(1) $$E_{1}^{n+1}\le E_{2}^{n+1}$$. We update $$q^{n+1}$$ by $$q^{n+1} =Q\left( \phi ^{n+1} \right)$$ and 9 is valid obviously. It means that the discrete modified energy is totally equal to the original energy, i.e. $\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert q^{n+1} \rVert ^2=\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert Q\left( \phi ^{n+1} \right) \rVert ^2.$

(2) $$E_{1}^{n+1}> E_{2}^{n+1}$$. Combining 8 and 10 , we see that $$0\le E_{2}^{n+1}=\lVert \hat{q}^{n+1} \rVert ^2+\tau M \left( \mathcal{G}\mu ^{n+\frac{1}{2}},\mu ^{n+\frac{1}{2}} \right)$$, which means $$\lVert \hat{q}^{n+1} \rVert ^2\le E_{2}^{n+1}<\lVert Q\left( \phi ^{n+1} \right) \rVert ^2$$. In terms of the closest approximation to the original energy, $$E_{2}^{n+1}$$ is the closest value to $$\lVert Q\left( \phi ^{n+1} \right) \rVert ^2$$ and satisfies the energy dissipative law 9 . So we update $$q^{n+1}$$ by $$q^{n+1}=\sqrt{E_{2}^{n+1}/E_{1}^{n+1}}Q\left( \phi ^{n+1} \right)$$, it means that $$\lVert q^{n+1} \rVert ^2=E_{2}^{n+1}$$ and 9 is valid.

◻

Theorem 2. The EOP-IEQ/CN scheme 3 is unconditionally energy stable in the sense that \label{CN1} \begin{align} {}&\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert q^{n+1} \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) -\lVert q^n \rVert ^2\\ ={}&-\tau M\left( \mathcal{G}\mu ^{n+\frac{1}{2}},\mu ^{n+\frac{1}{2}} \right)\le 0 , \end{align}\qquad{(5)} and if $$E_{1}^{n+1}\le E_{2}^{n+1}$$, we further have the following original dissipative law $\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert Q\left( \phi ^{n+1} \right) \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) -\lVert Q\left( \phi ^n \right) \rVert ^2\le 0.$

Proof. From the equation ?? in $$\boldsymbol{Step~2}$$, we obtain $$\lVert q^{n+1} \rVert ^2\le E_{2}^{n+1}$$, so ?? is obviously valid. According to $$\boldsymbol{Step~2}$$, we also obtain $$\lVert q^n \rVert ^2\le E_{1}^{n}$$. So we have \begin{align} \label{CN2} \frac{1}{2}\left( \mathcal{L}\phi ^{n},\phi ^{n} \right) +\lVert q^{n} \rVert ^2\le \frac{1}{2}\left( \mathcal{L}\phi ^{n},\phi ^{n} \right) +\lVert Q\left( \phi ^{n} \right) \rVert ^2, \quad n\ge 0. \end{align}\tag{11} If $$E_{1}^{n+1}\le E_{2}^{n+1}$$, we obtain $$\lVert q^{n+1} \rVert ^2=E_{1}^{n+1}$$, which means \begin{align} \label{CN3} \frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert Q\left( \phi ^{n+1} \right) \rVert ^2=\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert q^{n+1} \rVert ^2. \end{align}\tag{12} Combining ?? , 11 and 12 , we obtain $\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert Q\left( \phi ^{n+1} \right) \rVert ^2\le \frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) +\lVert Q\left( \phi ^n \right) \rVert ^2.$ The proof of the theorem is now complete. ◻

Theorem 3. The EOP-IEQ/CN scheme 3 is a special case of the REQ/CN scheme 2 under the condition of the artificial parameter $$\eta=1$$ which means that the maximum relaxation is introduced.

Proof. Suppose that the intermediate solutions $$(\phi^{n+1}, \hat{q}^{n+1})$$ are calculated from the baseline IEQ/CN scheme ?? , we show that the conclusions of the two schemes are consistent with $$\eta=1$$,.

If $$\lVert \hat{q}^{n+1}-Q\left( \phi ^{n+1} \right) \rVert=0$$, we know $$\lVert \hat{q}^{n+1} \rVert =\lVert Q\left( \phi ^{n+1} \right) \rVert$$. For the the relaxation step ?? of the REQ/CN scheme 2, we show that $$\xi^{n+1}=0.$$ Since $$a=0$$, it follows that $$\xi^{n+1}=0$$. For the energy-optimized technique ?? of the EOP-IEQ/CN scheme 3, we show that $$\lambda^{n+1}=1$$. According to the definition of $$E_{2}^{n+1}$$ and 8 , we have $$E_{1}^{n+1}=\lVert Q\left( \phi ^{n+1} \right) \rVert ^2\le \lVert \hat{q}^{n+1} \rVert ^2+\tau M ( \mathcal{G}\mu ^{n+\frac{1}{2}},\mu ^{n+\frac{1}{2}} ) =E_{2}^{n+1}$$, so $$\lambda ^{n+1}=\min \left\{ 1,\sqrt{E_{2}^{n+1}/E_{1}^{n+1}} \right\} =1$$.

If $$\lVert \hat{q}^{n+1}-Q\left( \phi ^{n+1} \right) \rVert> 0$$, we know $$a>0$$, then we analyze the values of $$b$$ and $$c$$.

(1) $$b>0$$. For the relaxation step ?? , we have $$(-b-\sqrt{b^2-4ac})/(2a)<0$$ and $$\xi ^{n+1}=0$$. According to the energy-optimized technique [EOP47CN], we obtain $$\lVert Q\left( \phi ^{n+1} \right) \rVert ^2<\left( \hat{q}^{n+1},Q\left( \phi ^{n+1} \right) \right) \le \left( \lVert \hat{q}^{n+1} \rVert ^2+\lVert Q\left( \phi ^{n+1} \right) \rVert ^2 \right)/2 .$$ Combining the above equality with the definition of $$E_{2}^{n+1}$$, we obtain $$E^{n+1}_1=\lVert Q\left( \phi ^{n+1} \right) \rVert ^2<\lVert \hat{q}^{n+1} \rVert ^2\le E^{n+1}_2,$$ so $$\lambda ^{n+1}=1$$.

(2) $$b\le 0$$, $$c\le 0$$. For the relaxation step ?? , we know that $$\xi ^{n+1}=\max \left\{ 0,\left( -b-\sqrt{b^2-4ac} \right) /\left( 2a \right) \right\}=0$$, and for the energy-optimized step [EOP47CN], we obtain $$\lVert Q\left( \phi ^{n+1} \right) \rVert ^2\le \lVert \hat{q}^{n+1} \rVert ^2+\tau M ( \mathcal{G}\mu ^{n+\frac{1}{2}},\mu ^{n+\frac{1}{2}} )$$ by $$c\le 0$$, which means $$E_{1}^{n+1}\le E_{2}^{n+1}$$, so $$\lambda ^{n+1}=1$$.

(3) $$b\le 0$$, $$c>0$$. For the REQ/CN scheme ?? , we obtain $$\xi ^{n+1}=(-b-\sqrt{b^2-4ac})/(2a)$$. For the EOP-IEQ/CN scheme 3, we have $$E_{2}^{n+1}<E_{1}^{n+1}$$ by $$c>0$$, so $$\lambda ^{n+1} =\sqrt{E_{2}^{n+1}/E_{1}^{n+1}}$$.

In conclusion, when $$\xi^{n+1}=0$$, we have $$\lambda^{n+1}=1$$, so we update $$q^{n+1}$$ by $$q^{n+1}=Q(\phi^{n+1})$$; and when $$\xi^{n+1}=(-b-\sqrt{b^2-4ac})/(2a)$$, we obtain $$\lambda^{n+1}=\sqrt{E_{2}^{n+1}/E_{1}^{n+1}}$$, then we get \begin{align} q^{n+1}{}&=\frac{-b-\sqrt{b^2-4ac}}{2a}\hat{q}^{n+1}+\left( 1-\frac{-b-\sqrt{b^2-4ac}}{2a} \right) Q\left( \phi ^{n+1} \right)\\ {}& =\sqrt{\frac{E_{2}^{n+1}}{E_{1}^{n+1}}}Q\left( \phi ^{n+1} \right) . \end{align} So the REQ/CN scheme is energy-optimized with the artificial parameter $$\eta=1$$. ◻

Remark 4. Theorem 3 shows that when the artificial parameter $$\eta=1$$, the relaxation step achieves energy optimization, which means that it agrees with the results obtained by the energy-optimized step, and when $$\eta\in [0,1)$$, the energy-optimized technique is always equal to or better than the relaxation technique.

Now we establish the consistency estimate for the EOP-IEQ/CN scheme 3. The consistency errors $$\varepsilon _{1}^{n+1}$$, $$\varepsilon _{2}^{n+1}$$, $$\varepsilon _{3}^{n+1}$$ for the EOP-IEQ/CN method are determined by the following: \label{trun} \begin{align} {}&\frac{\phi _{\star}^{n+1}-\phi _{\star}^{n}}{\tau}=-M\mathcal{G}\left( \mathcal{L}\phi _{\star}^{n+\frac{1}{2}}+H\left( \phi _{\star}^{n+\frac{1}{2}} \right) \hat{q}_{\star}^{n+\frac{1}{2}} \right) +\varepsilon _{1}^{n+1},\\ {}&\frac{\hat{q}_{\star}^{n+1}-q_{\star}^{n}}{\tau}=\frac{H\left( \phi _{\star}^{n+\frac{1}{2}} \right)}{2}\left( \frac{\phi _{\star}^{n+1}-\phi _{\star}^{n}}{\tau} \right) +\varepsilon _{2}^{n+1},\\ {}&q_{\star}^{n+1}=\lambda ^{n+1}Q\left( \phi _{\star}^{n+1} \right)+\varepsilon _{3}^{n+1} ,\quad \lambda ^{n+1}=\min \left\{ 1,\sqrt{\frac{E^{n+1}_2}{E^{n+1}_1}} \right\}, \end{align}\tag{13} where $$\phi _{\star}^{n+1}=\phi \left( t_{n+1} \right)$$, $$q_{\star}^{n+1}=q\left( t_{n+1} \right)$$, $$Q_{\star}^{n+1}=Q\left( \phi(t_{n+1}) \right)$$.

Suppose that the analytical solutions ($$\phi,q$$) of the equivalent gradient flow model 5 possess the following regularity condition $\label{zheng} \left\{ \begin{array}{l} \phi \in L^{\infty}\left( 0,T;H^2\left( \Omega \right) \right) ,\quad q\in L^{\infty}\left( 0,T;L^{\infty}\left( \Omega \right) \right) ,\\ \phi _t\in L^{\infty}\left( 0,T;L^{\infty}\left( \Omega \right) \right) ,\quad q_{tt},\phi _{tt}\in L^2\left( 0,T;L^2\left( \Omega \right) \right). \end{array} \right.\tag{14}$ We first give the following lemma to establish the quantitative relation between the $$L^2$$ norm of $$H( \phi _{\star}^{n+1}) -H( \phi _{\star}^{n} )$$ and $$\phi _{\star}^{n+1}-\phi _{\star}^{n}$$ under the regularity condition 14 .

Lemma 3 ([23]). Suppose that

(1) $$F(x)$$ uniformly bounded from below: $$F(x)>-C_*$$, $$x\in(-\infty,+\infty)$$;

(2) $$F(x) \in C^2(-\infty,+\infty)$$;

(3) there exists a positive constant $$C_{\star}$$ such that $$\underset{0\le n\le N}{\max}\left( \lVert \phi^n_{\star} \rVert _{L^{\infty}} \right) \le C_{\star}$$.

We have \begin{align} \lVert H\left( \phi _{\star}^{n+1} \right) -H\left( \phi _{\star}^{n} \right) \rVert \le C_1\lVert \phi _{\star}^{n+1}-\phi _{\star}^{n} \rVert ,\quad n\ge 0, \end{align} where $$C_1$$ dependens only on $$C_*$$, $$C_0$$, $$C_{\star}$$.

Theorem 4. We consider $$\mathcal{G} = I$$ for the $$L^2$$ gradient flow for simplicity, and the the analytical solutions ($$\phi,q$$) satisfie the regularity condition 14 , then the following consistency estimate holds for the system 13 : $\lVert \varepsilon _{1}^{n+1} \rVert +\lVert \varepsilon _{2}^{n+1} \rVert +\lVert \varepsilon _{3}^{n+1} \rVert \le C\tau^2,$ where $$C$$ is independent of $$\tau$$.

Proof. Let $$R_{\phi}^{n+1}$$ be the truncation error defined by \label{err1} \begin{align} R_{\phi}^{n+1} ={}&\frac{\phi _{\star}^{n+1}-\phi _{\star}^{n}}{\tau}-\frac{\partial \phi \left( t_{n+\frac{1}{2}} \right)}{\partial t}\\ ={}&\frac{1}{2\tau}\left( \int_{t_n}^{t_{n+\frac{1}{2}}}{\left( t_n-s \right) ^2\frac{\partial ^2\phi \left( s \right)}{\partial t^2}}ds+\int_{t_{n+\frac{1}{2}}}^{t_{n+1}}{\left( t_{n+1}-s \right) ^2}\frac{\partial ^2\phi \left( s \right)}{\partial t^2}ds \right) . \end{align}\tag{15} Since \label{err2} \begin{align} {}&\frac{\phi _{\star}^{n+1}+\phi _{\star}^{n}}{2}-\phi _{\star}^{n+\frac{1}{2}}=\frac{1}{2}\left( \int_{t_n}^{t_{n+\frac{1}{2}}}{\left( s-t_n \right) \frac{\partial ^2\phi}{\partial t^2}\left( s \right)}ds+\left( t_{n+1}-s \right) \frac{\partial ^2\phi}{\partial t^2}\left( s \right) ds \right) ,\\ {}&\frac{\hat{q}_{\star}^{n+1}+\phi _{\star}^{n}}{2}-q_{\star}^{n+\frac{1}{2}}=\frac{1}{2}\left( \int_{t_n}^{t_{n+\frac{1}{2}}}{\left( s-t_n \right) \frac{\partial ^2q}{\partial t^2}\left( s \right)}ds+\left( t_{n+1}-s \right) \frac{\partial ^2q}{\partial t^2}\left( s \right) ds \right) . \end{align}\tag{16} From the Lemma 3, we obtain \begin{align} \label{err3} \lVert H\left( \phi _{\star}^{n+1} \right) -H\left( \frac{3}{2}\phi _{\star}^{n}-\frac{1}{2}\phi _{\star}^{n-1} \right) \rVert \le C_1\lVert \phi _{\star}^{n+1}-\frac{3}{2}\phi _{\star}^{n}+\frac{1}{2}\phi _{\star}^{n-1} \rVert \le C_1C_0\tau^2. \end{align}\tag{17} By the three inequalities 15 -17 and triangle inequality, we have \begin{align} \lVert \varepsilon _{1}^{n+1} \rVert \le{}& \lVert R_{\phi}^{n+1} \rVert +\varepsilon ^2\lVert \Delta \left( \frac{\phi _{\star}^{n+1}+\phi _{\star}^{n}}{2}-\phi _{\star}^{n+\frac{1}{2}} \right) \rVert \\ {}&+\lVert \frac{\hat{q}^{n+1}_{\star}+q^n_{\star}}{2}H\left( \frac{3}{2}\phi _{\star}^{n}-\frac{1}{2}\phi _{\star}^{n-1} \right) -q_{\star}^{n+\frac{1}{2}}H\left( \phi _{\star}^{n+\frac{1}{2}} \right) \rVert \\ \le{}& C\tau^2. \end{align} Similarly we can obtain that $$\lVert \varepsilon _{2}^{n+1} \rVert\le C\tau^2$$. By Theorem 3, we obtain $$q_{\star}^{n+1}=\xi ^{n+1}\hat{q}^{n+1}+\left( 1-\xi ^{n+1} \right) Q\left( \phi _{\star}^{n+1} \right) =Q\left( \phi _{\star}^{n+1} \right) +O\left( \tau ^2 \right),\forall \xi ^{n+1}\in \left[ 0,1 \right]$$. ◻

Remark 5 (Order of Accuracy [30]). The proposed EOP-IEQ/CN scheme 3 is second-order accurate in time. Notice the fact, $$\hat{q}^{n+1}=q(t_{n+1})+O(\tau^2)$$ and $$Q(\phi^{n+1})=q(t_{n+1})+O(\tau^2)$$. Thus for $$\eta=1$$, $$q^{n+1}=\xi^{n+1}\hat{q}^{n+1}+(1-\xi^{n+1})Q(\phi^{n+1})=q(t_{n+1})+O(\tau^2), \forall \xi^{n+1}\in[0,1]$$, which means the energy-optimized technique does not affect the order of accuracy for the baseline IEQ method.

3.2 The EOP-IEQ/BDF$$k$$ schemes↩︎

If we use the backward Euler formulas for the time discretization, we will have the semi-implicit EOP-IEQ/BDF$$k$$ $$(k=1,2)$$ schemes as below.

Corollary 4 ($$k$$th-order EOP-IEQ/BDF$$k$$ Scheme). For initial datas, we set $$\phi ^0=\phi \left( \boldsymbol{x},0 \right)$$ and $$q^0=\sqrt{F\left( \phi ^0 \right) +C_0}$$, then we update $$\phi^{n+1}$$, $$q^{n+1}$$ via the following two steps:

• $$\boldsymbol{Step~1} :$$ To calculate the intermediate solutions $$(\phi^{n+1}, \hat{q}^{n+1})$$ by the following semi-implicit IEQ/BDF$$k$$ schemes: \begin{align} {}&\frac{\alpha _k\phi ^{n+1}-A_k\left( \phi ^n \right)}{\tau}=-M\mathcal{G}\mu ^{n+1},\\ {}&\mu ^{n+1}=\mathcal{L}\phi ^{n+1}+\hat{q}^{n+1}H\left( B_k\left( \phi ^n \right) \right) ,\\ {}&\frac{\alpha _k\hat{q}^{n+1}-A_k\left( q^n \right)}{\tau}=\frac{H\left( B_k\left( \phi ^n \right) \right)}{2}\frac{\alpha _k\phi ^{n+1}-A_k\left( \phi ^n \right)}{\tau}, \end{align} where $$H\left( B_k\left( \phi ^n \right) \right) =\frac{f\left( B_k\left( \phi ^n \right) \right)}{Q\left( B_k\left( \phi ^n \right) \right)}$$.

• $$\boldsymbol{Step~2} :$$ To update the auxiliary variable $$q^{n+1}$$ via an energy-optimized step as:

• $$\bullet$$ First-order EOP/BDF$$1$$. $\label{EOP47BDF1} q^{n+1}=\lambda ^{n+1}Q(\phi^{n+1}),\quad \lambda ^{n+1}=\min \left\{ 1,\sqrt{\frac{E_{2}^{n+1}}{E_{1}^{n+1}}} \right\}.\qquad{(6)}$ where the coefficients $$E_{1}^{n+1}$$ and $$E_{2}^{n+1}$$ are given by \begin{align} E_{1}^{n+1}\!=\!\lVert Q( \phi ^{n+1} ) \rVert ^2,\quad E_{2}^{n+1}\!=\!\frac{1}{2}( \mathcal{L}\phi ^n,\phi ^n ) \!+\!\lVert q^n \rVert ^2 \!-\!\frac{1}{2}( \mathcal{L}\phi ^{n+1},\phi ^{n+1} ) . \end{align}

• $$\bullet$$ Second-order EOP/BDF$$2$$. \begin{align} \label{EOP47BDF2} q^{n+1}=\lambda ^{n+1}\left( Q\left( \phi ^{n+1} \right) \!-\!\frac{2}{5}q^n \right)\!+\!\frac{2}{5}q^n ,\quad \lambda ^{n+1}=\min \left\{ 1,\sqrt{\frac{\bar{E}_{2}^{n+1}}{\bar{E}_{1}^{n+1}}} \right\} , \end{align}\qquad{(7)} where the coefficients are given by \begin{align} \bar{E}_{1}^{n+1}={}&\lVert Q\left( \phi ^{n+1} \right) -\frac{2}{5}q^n \rVert ^2,\\ \bar{E}_{2}^{n+1}={}&\frac{1}{10}\left( \mathcal{L}\phi ^n,\phi ^n \right) +\frac{1}{10}\left( \mathcal{L}\left( 2\phi ^n-\phi ^{n-1} \right) ,2\phi ^n-\phi ^{n-1} \right)\\ {}&-\frac{1}{10}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) -\frac{1}{10}\left( \mathcal{L}\left( 2\phi ^{n+1}-\phi ^n \right) ,2\phi ^{n+1}-\phi ^n \right)\\ {}&+\frac{4}{25}\lVert q^n \rVert ^2+\frac{1}{5}\lVert 2q^n-q^{n-1} \rVert ^2. \end{align}

Theorem 5. The EOP-IEQ/BDF$$k$$ $$( k = 1, 2 )$$ schemes 4 are unconditionally energy stable.

Proof. For $$k=1, 2$$, we have the following analysis process respectively.

• $$\bullet$$ First-order EOP-IEQ/BDF$$1$$ scheme. From ?? , we obtain \begin{align} \lVert q^{n+1} \rVert ^2=\left( \lambda ^{n+1} \right) ^2\lVert Q\left( \phi ^{n+1} \right) \rVert ^2\le \left( \sqrt{\frac{E_{2}^{n+1}}{E_{1}^{n+1}}} \right) ^2 E_{1}^{n+1}=E_{2}^{n+1}, \end{align} According to the definition of $$E_{2}^{n+1}$$ and above inequation, we have \begin{align} {}&\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\lVert q^{n+1} \rVert ^2-\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) -\lVert q^n \rVert ^2\\ \le{}& -\tau M \left( \mathcal{G}\mu ^{n+1},\mu ^{n+1} \right) \le 0. \end{align}

• $$\bullet$$ Second-order EOP-IEQ/BDF$$2$$ scheme. From ?? , we have $\lVert q^{n+1}-\frac{2}{5}q^n \rVert ^2=\left( \lambda ^{n+1} \right) ^2\lVert Q\left( \phi ^{n+1} \right) -\frac{2}{5}q^n \rVert ^2\le \bar{E}_{2}^{n+1}.$ By the definition of $$\bar{E}_{2}^{n+1}$$ and above inequation, we obtain \begin{align} {}&\frac{1}{4}\left[ \left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\left( \mathcal{L}\left( 2\phi ^{n+1}-\phi ^n \right) ,2\phi ^{n+1}-\phi ^n \right) \right]\\ {}&+\frac{1}{2}\left( \lVert q^{n+1} \rVert ^2+\lVert 2q^{n+1}-q^n \rVert ^2 \right)\\ {}&-\frac{1}{4}\left[ \left( \mathcal{L}\phi ^n,\phi ^n \right) +\left( \mathcal{L}\left( 2\phi ^n-\phi ^{n-1} \right) ,2\phi ^n-\phi ^{n-1} \right) \right] \\ {}&-\frac{1}{2}\left( \lVert q^n \rVert ^2+\lVert 2q^n-q^{n-1} \rVert ^2 \right)\\ \le{}& -\tau M \left( \mathcal{G}\mu ^{n+1},\mu ^{n+1} \right) . \end{align}

The proof is completed. ◻

Remark 6. It can be analyzed by a similar proof procedure of Theorem 3 that for the EOP-IEQ/BDF$$k$$ $$(k=1,2)$$ schemes, the energy-optimized technique is always equal to or superior to the relaxation technique, for any artificial parameter $$\eta=[0,1]$$.

3.3 energy-optimized techniques for the variations of the baseline IEQ method↩︎

Due to the wide applicability of the baseline IEQ method, several variants have been proposed to be suitable for different situations . Since our proposed energy-optimized technique does not restrict the form of auxiliary variables, our EOP-IEQ approach can be easily extended to gradient flow models that contain multiple auxiliary variables. Below we describe the EOP-IEQ method applied to a gradient flow with one phase field variable $$\phi$$. We consider a more general form of the free energy \begin{align} \label{ME} \mathcal{E}\left( \phi \right) =\frac{1}{2}\left( \mathcal{L}\phi ,\phi \right) +\sum_{i=1}^k{\left( F_i\left( \phi \right) ,1 \right)}, \end{align}\tag{18} where the operator $$\mathcal{L}$$ is a linear self-adjoint non-negative definite operator and $$F_i\left( \phi \right)$$, $$i=1$$, $$2$$, $$\cdots$$, $$k$$ are the bulk potentials. The general gradient flow system of the free energy 18 as follows \label{MPF} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\mathcal{G}\mu ,\\ {}&\mu =\mathcal{L}\phi +\sum_{i=1}^k{F_{i}^{'}\left( \phi \right)}, \end{align}\tag{19} which has the following energy dissipation law $\frac{dE\left( \phi \right)}{dt}=-\mathcal{G}\left( \mu ,\mu \right)\le 0 ,$ where $$M$$ is a mobility constant and the operator $$\mathcal{G}$$ is a non-negative definite operator. For the system 19 , assume that each nonlinear bulk potentials have a lower bound $$F_i(\phi)\ge C^*_i$$, $$i=1$$, $$2$$, $$\cdots$$, $$k$$. Denote $$\boldsymbol{q}=\left( q_1,q_2,\dots ,q_k \right)$$, we introduce multiple auxiliary variables $q_i( \boldsymbol{x},t ) =:Q_i(\phi( \boldsymbol{x},t )) =\sqrt{F_i( \phi( \boldsymbol{x},t ) ) +C^0_i},\quad C^0_i>C^*_i,i=1,2,\dots,k.$ Then the equivalent form of the gradient flow system 19 can be written as \label{MIEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\mathcal{G}\mu ,\\ {}&\mu =\mathcal{L}\phi +\sum_{i=1}^k{q_iH_i( \phi )},,\\ {}&\frac{\partial q_j}{\partial t}=\frac{H_j( \phi )}{2}\frac{\partial \phi}{\partial t},\quad j=1,2,\dots ,k. \end{align}\tag{20} where $$H_i( \phi ) =\frac{f_i( \phi )}{Q_i( \phi )},f_i=F_{i}^{'}\left( \phi \right),i=1,2,\dots ,k.$$ Then we can get the modified free energy as \begin{align} \bar{\mathcal{E}}\left( \phi,\boldsymbol{q} \right) =\frac{1}{2}\left( \mathcal{L}\phi ,\phi \right) +\sum_{i=1}^k{\lVert q_i \rVert ^2}-\sum_{i=1}^k{C^0_i|\Omega|}. \end{align} For the reformulated system 20 , it has the following modified energy dissipation law $\frac{d\bar{E}( \phi ,\boldsymbol{q} )}{dt}=-M\mathcal{G}( \mu ,\mu ) \le 0.$ We introduce the energy-optimized technique into the system 20 to improve the inconsistency between the modified energy and the original energy after discretization. Then a second-order EOP-IEQ scheme based on Crank-Nicolson formula can be written as follows.

Corollary 5. For initial datas, we set $$\phi ^0=\phi \left( \boldsymbol{x},0 \right)$$, $$q^0=\sqrt{F\left( \phi ^0 \right) +C_0}$$, then we update $$\phi^{n+1}$$, $$q^{n+1}$$ via the following two steps:

• $$\boldsymbol{Step~1} :$$ To obtain the numerical solution $$\phi^{n+1}$$ by the baseline IEQ method: \begin{align} {}&\frac{\phi ^{n+1}-\phi ^n}{\tau}=-M\mathcal{G}\mu ^{n+\frac{1}{2}},\\ {}&\mu ^{n+\frac{1}{2}}=\mathcal{L}\frac{\phi ^{n+1}+\phi ^n}{2}+\sum_{i=1}^k{\frac{\hat{q}_{i}^{n+1}+q_{i}^{n}}{2}H_i( \bar{\phi}^{n+\frac{1}{2}} )},\\ {}&\frac{\hat{q}_{j}^{n+1}-q_{j}^{n}}{\tau}=H_j( \bar{\phi}^{n+\frac{1}{2}}) \frac{\phi ^{n+1}-\phi ^n}{\tau},\quad j=1,2,\dots ,k. \end{align} where $H_i( \bar{\phi}^{n+\frac{1}{2}} ) =\frac{f_i( \bar{\phi}^{n+\frac{1}{2}} )}{Q_i( \bar{\phi}^{n+\frac{1}{2}} )},\quad \bar{\phi}^{n+\frac{1}{2}}=\frac{3}{2}\phi ^n-\frac{1}{2}\phi ^{n-1}.$

• $$\boldsymbol{Step~2} :$$ To update the auxiliary variables $$\boldsymbol{q}^{n+1}=\left( q_{1}^{n+1},q_{2}^{n+1},\cdots ,q_{k}^{n+1} \right)$$ with the following energy-optimized step: $q_{i}^{n+1}=\lambda ^{n+1}Q_i\left( \phi ^{n+1} \right) ,\quad \lambda ^{n+1}=\min \left\{ 1,\sqrt{\frac{E^{n+1}_2}{E^{n+1}_1}} \right\} , i=1,2,\cdots ,k,$ where $E^{n+1}_1=\sum_{i=1}^k{\lVert Q_i\left( \phi ^{n+1} \right) \rVert ^2},\quad E^{n+1}_2=\bar{E}^n-\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) .$

Theorem 6. The scheme 5 is unconditionally energy stable in the sense that \begin{align} {}&\frac{1}{2}\left( \mathcal{L}\phi ^{n+1},\phi ^{n+1} \right) +\sum_{i=1}^k{\lVert q_{i}^{n+1} \rVert ^2}-\frac{1}{2}\left( \mathcal{L}\phi ^n,\phi ^n \right) -\sum_{i=1}^k{\lVert q_{i}^{n} \rVert ^2}\\ \le{}& -\tau M\left( \mathcal{G}\mu ^{n+\frac{1}{2}},\mu ^{n+\frac{1}{2}} \right) \le 0. \end{align}

The proof of this theorem is not particularly difficult but will not be reproduced here.

4 Numerical implementation and results↩︎

In this section, we adopt some widely used gradient flow models including Allen-Cahn (AC) equation, Cahn-Hilliard (CH) equation, molecular beam epitaxy (MBE) model, and phase field crystal (PFC) model to verify the generality of EOP-IEQ method. Furthermore, the corresponding numerical algorithms are implemented to verify the above theoretical results. For simplicity, we consider periodic boundary condition and use the Fourier spectrum method in space in all examples, where $$N_x$$ and $$N_y$$ represent the number of Fourier modes along the $$x$$ and $$y$$ axis respectively. It is worth noting that the proposed schemes can also be applied to other thermodynamically consistent boundary conditions which satisfy energy dissipation law.

4.1 Allen-Cahn equation↩︎

To begin with, we consider the free energy [36] \begin{align} \label{E-ACH} E\left( \phi \right) =\int_{\Omega}{\frac{\alpha_0}{2}\left| \nabla \phi \right|^2}+\frac{\alpha_0}{\varepsilon ^2}F\left( \phi \right) \text{d}\boldsymbol{x,}\quad F\left( \phi \right) =\frac{\left( \phi ^2-1 \right) ^2}{4}, \end{align}\tag{21} where the interface parameters $$0<$$ $$\alpha_0$$, $$\varepsilon$$ $$\ll 1$$.

The Allen-Cahn (AC) equation [37] can be expressed in the energy variational form 2 of the free energy 21 with $$\mathcal{G} = MI$$ as follows \begin{align} \label{AC} \frac{\partial \phi}{\partial t}=-M\left( -\alpha_0 \Delta \phi +\frac{\alpha_0}{\varepsilon ^2}f\left( \phi \right) \right). \end{align}\tag{22} where $$f\left( \phi \right) =F'(\phi)=\phi ^3-\phi$$. Now we introduce the auxiliary variable $$q:=Q\left( \phi \right) =\sqrt{F\left( \phi \right) +C_0}$$, then the general model 5 is specified as \label{AC-IEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\left( -\alpha_0 \Delta \phi +\frac{\alpha_0 }{\varepsilon ^2}\frac{q}{Q\left( \phi \right)} f\left( \phi \right)\right) ,\\ {}&\frac{\partial q}{\partial t}=\frac{f\left( \phi \right)}{2Q\left( \phi \right)}\frac{\partial \phi}{\partial t}. \end{align}\tag{23}

Case A. We check the convergence rate in time of the proposed schemes for the AC equation under the $$2D$$ computational region $$\Omega = [-0.5,0.5]^2$$. The initial data is chosen to be smooth \begin{align} \label{AC-init} \phi \left( x,y,0 \right) =\tanh\frac{a+b\cos \left( 6\theta \right) -c \sqrt{x^2+y^2}}{\sqrt{2}\varepsilon}, \quad \theta =\arctan \frac{y}{x}, \end{align}\tag{24} where $$a,b,c$$ are positive constants. We set $$(N_x, N_y) =(64,64)$$, so that the spatial discrete error is negligible compared to the temporal discrete error. The parameters are $$(a,b,c) = (1.5,1.2,2\pi)$$, $$\alpha_0=0.01$$, $$\varepsilon =0.01$$, $$M=0.6$$, $$C_0=100$$. In addition, since the exact solution is unknown, we choose the EOP-IEQ/CN scheme with time step $$\tau = 1e-6$$ as the reference solution for the calculation error. In Fig. 1, we obtain the numerical errors of $$L^2$$ norm and convergence rates at $$T = 0.1$$ using various first- and second-order schemes. We observe that (i) the temporal convergence rates obtained by baseline IEQ and EOP-IEQ methods are consistent with theoretical expectations in all cases; (ii) the $$L^2$$ errors of the numerical solutions $$\phi$$ and the auxiliary variable $$q$$ in the EOP-IEQ schemes are significantly smaller than those in the baseline IEQ schemes.

Case B. We also choose the initial condition 24 to compare the accuracy of the baseline IEQ, REQ and the EOP-IEQ methods for solving the AC equation. We adjust some parameters $$(N_x, N_y) =(256,256)$$, $$\eta=0.5$$, $$M=0.51$$, $$C_0=1$$, $$\tau=1e-2$$ and $$T=10.5$$ and use the results of the EOP-IEQ/BDF2 scheme of $$\tau=1e-6$$ as the reference solution. The profiles of $$\phi$$ at various times are summarized in Fig. 2. A comparison of the modified energy (left), the energy error (middle), and the update parameters (right) for the IEQ/CN, REQ/CN, and EOP-IEQ/CN schemes shown in Fig. [F14], where the update parameters include relaxation parameters $$\xi^n$$ and energy-optimized parameters $$\lambda^n$$, where $$\lambda_1^n=|\lambda^n-1|$$. We observe that the energy-optimized parameters $$\lambda_1^n=0, \forall n\ge 1$$ in the EOP-IEQ/CN scheme, which means $$\lambda^n=1$$ and $$q^{n}=Q(\phi^n), \forall n\ge 1$$. It is essential to preserve the consistency of the modified energy with the original energy. Meanwhile, Fig. [F14] (b) shows that the errors between the corrected energy and the reference energy using the IEQ/CN and REQ/CN schemes are significantly larger. It is verified that our proposed EOP-IEQ approach ensure the modified energy closer to original energy than the REQ approach.

4.2 Cahn-Hilliard equation↩︎

Now we consider the Cahn-Hilliard (CH) equation given as \label{CH} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\Delta \mu ,\\ {}&\mu =-\alpha _0\Delta \phi +\frac{\alpha _0}{\varepsilon ^2}f\left( \phi \right) . \end{align}\tag{25} It can be expressed in the energy variational form 2 of the free energy 21 with $$\mathcal{G} = -M\Delta$$. We introduce the auxiliary variable $$q:=Q\left( \phi \right) =\sqrt{F\left( \phi \right) -\kappa\phi ^2/2+C_0}$$, then the general model in 5 is specified as \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\Delta \left( -\alpha _0\Delta \phi +\frac{\alpha _0}{\varepsilon ^2}\kappa \phi +\frac{q}{Q\left( \phi \right)}\frac{\alpha _0}{\varepsilon ^2}f_{\kappa}\left( \phi \right) \right) ,\\ {}&\frac{\partial q}{\partial t}=\frac{f_{\kappa}\left( \phi \right)}{2Q\left( \phi \right)}\frac{\partial \phi}{\partial t}, \end{align} where $$f_{\kappa}\left( \phi \right) =\phi^3-\phi-\kappa\phi$$. In addition, the CH equation satisfies the law of conservation of mass \begin{align} \label{MM} \frac{dm\left( t \right)}{dt}=0,\quad m\left( t \right) =\int_{\Omega}{\phi \left( \boldsymbol{x,}t \right)}\,\text{d}\boldsymbol{x}. \end{align}\tag{26}

Case A. The initial conditions [38] is considered as the following discrete form \begin{align} \label{CH95init} \phi \left( x,y,0 \right) =0.25+0.4rand\left( x,y \right) , \end{align}\tag{27} where $$rand(x,y)$$ represents a pseudo-random value derived from a uniform distribution within $$(-1,1)^2$$. We choose the numerical parameters $$\Omega = [-0.5,0.5]^2, (N_x,N_y)=(64,64)$$, $$\alpha_0=0.01^2$$, $$\varepsilon = 0.01$$, $$\kappa=4$$, $$M=0.001$$, $$C_0=100$$, $$T=0.1$$, and the reference solution is computed by the EOP-IEQ/BDF2 scheme for $$\tau=1e-6$$. The numerical errors and convergence rates of EOP-IEQ method are shown in Fig 4. We observed that the numerical results of the numerical solution $$\phi$$ and the auxiliary variable $$q$$ in $$L^2$$ errors obtained by using the proposed schemes are in agreement with expectations.

Case B. Considering that coarsening effect is an important phenomenon in phase field simulation, this example is used to verify the validity of using energy-optimized technique to simulates long-time dynamics. The initial data is set to be the same as 27 and we adjust some parameters to $$(N_x,N_y)=(128,128)$$, $$\tau=1e-3$$, $$M=0.1$$ and $$T=100$$. The results of the EOP-IEQ/CN scheme of $$\tau=1e-5$$ are calculated as the reference solution. Fig. 5 shows several profiles of coarsening process under the EOP-IEQ/CN scheme at $$t=1, 10, 90$$. We use the IEQ/CN and EOP-IEQ/CN schemes in Fig. 6 to investigate the modified energy, the ratio of nonlinear free energy and mass profiles during the evolution process. We observe that the modified energy obtained by the EOP-IEQ method is closer to the original energy than that obtained by the baseline IEQ method. It indicates the energy-optimized technique improves the accuracy significantly.

Case C. In the following examples, we will use the second-order EOP-IEQ/CN scheme [EOP47CN] for numerical simulation if not specified. We study the evolution of isolated ellipses and circles [32]. The initial condition is $\phi \left( x,y,0 \right) =\tanh\left( \frac{dist\left( \left( x,y \right) ,\varGamma \right)}{\sqrt{2}\varepsilon} \right) ,$ where $$dist\left( \left( x,y \right) ,\varGamma \right)$$ represents the signed distance from point $$\left( x,y \right)$$ to the ellipse and circular interfaces $$\varGamma$$. The center of the ellipse is located at $$(-0.1, -0.1 )$$, its major and minor axes are $$\sqrt{2}/5$$, $$\sqrt{2}/10$$, the central position of the circle is $$(0.25,0.25)$$, and the radius is set to $$0.1$$. Other parameters are $$(N_x,N_y)=(128,128)$$, $$\alpha_0=0.01^2$$, $$\tau=1e-4$$, $$\varepsilon = 0.01$$, $$\kappa=4$$, $$M=1$$, $$C_0=100$$ and $$T=1$$. We use the results of the EOP-IEQ/CN scheme of $$\tau=1e-6$$ as the reference solution. Fig. 7 shows the morphological evolution under the IEQ/CN scheme with or without energy-optimized technique, and it can also be seen that the small circle is gradually absorbed by the ellipse until it completely disappears. We observe that the numerical results of the EOP-IEQ/CN scheme are closer to the reference solution from the enlarged figures. At the same time, it can be observed that the modified energy using EOP-IEQ/CN scheme is closer to the reference energy from the energy curve in Fig. 8 (a), (b), while the baseline IEQ method has significant deviation. In addition, we observe that the baseline IEQ method can maintain the conservation of mass regardless of the energy-optimized technique in Fig. 8 (c).

4.3 Phase field crystal model↩︎

Next, we investigate the phase field crystal (PFC) model , considering the free energy \begin{align} \label{E-PFC} E\left( \phi \right) =\int_{\Omega}{\left( \frac{1}{4}\phi ^4+\frac{\beta-\varepsilon }{2}\phi ^2-\beta\left| \nabla \phi \right|^2+\frac{1}{2}\left( \Delta \phi \right) ^2 \right) \,\text{d}\boldsymbol{x}}, \end{align}\tag{28} where $$\beta$$ the positive parameter and $$\beta>\varepsilon$$. Then the PFC equation can be expressed in the energy variational form 2 of the free energy 28 with $$\mathcal{G} = -M\Delta$$ as follows \begin{align} {}&\phi _t=-M\Delta \mu, \\ {}&\mu =\phi ^3+\left( \beta-\varepsilon \right) \phi +2\beta\Delta \phi +\Delta ^2\phi . \end{align} For the baseline IEQ method, we introduce the auxiliary variable $$q:=Q\left( \phi \right) =\phi^2$$, then the general system 5 is specified as \label{PFC-IEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=M\Delta \mu ,\\ {}&\mu =q\phi +\left( \beta-\varepsilon \right) \phi +2\beta\Delta \phi +\Delta ^2\phi ,\\ {}&\frac{\partial q}{\partial t}=2\phi \frac{\partial \phi}{\partial t}. \end{align}\tag{29} In addition, the PFC equation also satisfies the law of conservation of mass, i.e. 26 is valid.

We simulate the complex dynamics of the growth of a polycrystal in a two-dimensional supercooled liquid by using the EOP-IEQ method to further highlight its power for simulating long-time dynamics, involving the motion of the liquid crystal interface and the grain boundaries separating the crystals. To define the initial configuration, we set up three microcrystals with different orientations as shown in Fig. 9 (a) in the computing domain $$[0,800]^2$$. Other parameters are $$(N_x,N_y)=(512,512)$$, $$\varepsilon =0.25$$, $$\beta=1$$, $$M=1$$, $$\tau=1e-2$$ and $$T=1500$$. We use the results of the EOP-IEQ/CN scheme of $$\tau=1e-4$$ as the reference solution. Fig. 9 shows microstructure evolution dynamics for triangular crystal phases driven by the PFC model using the EOP-IEQ/CN scheme. The results also show that the different arrangements of crystals lead to defects and dislocations, similar to those reported in the literature [27]. The comparison of the numerical results of discrete modified energies is shown in Fig. 10 (a). The EOP-IEQ/CN scheme is more accurate than the baseline IEQ/CN scheme in Fig. 10 (b). Most importantly, the ratio of nonlinear free energy $$\lVert Q\left( \phi ^n \right) \rVert ^2/\lVert q^n \rVert ^2=1$$ by using the energy-optimized technique in Fig. 10 (b), which shows that the EOP-IEQ method preserves the consistency of $$Q( \phi ^n)$$ and $$q^n$$. Finally, Fig. 10 (c) indicates that the IEQ/CN scheme can ensure the conservation of mass regardless of the energy-optimized technique.

4.4 Molecular beam epitaxy model↩︎

Finally, we focus on the molecular beam epitaxy (MBE) model without slope section. Considering the Ehrlich-Schwoebel energy \begin{align} \label{E-MBE} E\left( \phi \right) =\int_{\Omega}{\left( \frac{\varepsilon ^2}{2}\left( \Delta \phi \right) ^2-\frac{1}{2}\ln \left( 1+\left| \nabla \phi \right|^2 \right) \right) \,\text{d}\boldsymbol{x}}. \end{align}\tag{30} Then the MBE equation can be expressed in the energy variational form 2 of the free energy 30 with $$\mathcal{G} = MI$$, which reads as $\frac{\partial \phi}{\partial t}=-M\left( \varepsilon ^2\Delta ^2\phi +\nabla \cdot \left( \frac{1}{1+\left| \nabla \phi \right|^2} \nabla \phi \right) \right) .$ By introducing the auxiliary variable $$q :=Q\left( \phi \right) =\sqrt{\ln \left( 1+\left| \nabla \phi \right|^2 \right) +C_0}$$, the MBE model can be equivalently written in quadratic form as follows \label{MBE-IEQ} \begin{align} {}&\frac{\partial \phi}{\partial t}=-M\left( \varepsilon ^2\Delta ^2\phi +\nabla \cdot \left( q\boldsymbol{H}(\phi) \right) \right) ,\\ {}&\frac{\partial q}{\partial t}=\boldsymbol{H}(\phi)\cdot \nabla \frac{\partial \phi}{\partial t}, \end{align}\tag{31} where $\boldsymbol{H}(\phi)=\frac{\nabla \phi}{\left( 1+\left| \nabla \phi \right|^2 \right) Q(\phi) }.$ Then the modified energy is $\bar{E}\left( \phi ,q \right) =\frac{\varepsilon ^2}{2}\lVert \Delta \phi \rVert ^2-\frac{1}{2}\lVert q \rVert ^2+\frac{1}{2}C_0\left| \Omega \right|.$

The energy-optimized technique for the MBE model without slope section 31 is somewhat different from that of the general scheme [EOP47CN], so we reconsider the EOP-IEQ method for the equivalent system 31 based on the Crank-Nicolson formula.

Corollary 6 (EOP-IEQ/CN scheme for the MBE model without slope section). For initial datas, we set $$\phi ^0=\phi \left( \boldsymbol{x},0 \right)$$, $$q^0=\sqrt{\ln \left( 1+\left| \nabla \phi^0 \right|^2 \right) +C_0}$$, then we update $$\phi^{n+1}$$, $$q^{n+1}$$ via the following two steps:

• $$\boldsymbol{Step~1}$$ : To obtain the numerical solution $$\phi^{n+1}$$ by the baseline IEQ method: \begin{align} {}&\frac{\phi ^{n+1}-\phi ^n}{\tau}=-M\left( \varepsilon ^2\Delta ^2\frac{\phi ^{n+1}+\phi ^n}{2}+\nabla \cdot \left( \frac{\hat{q}^{n+1}+q^n}{2}\boldsymbol{\bar{H}}^{n+\frac{1}{2}} \right) \right) ,\\ {}&\frac{\hat{q}^{n+1}-q^n}{\tau}=\boldsymbol{\bar{H}}^{n+\frac{1}{2}}\cdot \nabla \frac{\phi ^{n+1}-\phi ^n}{\tau}, \end{align} where $$\boldsymbol{\bar{H}}^{n+\frac{1}{2}}=\boldsymbol{H}\left( \frac{3}{2}\phi ^n-\frac{1}{2}\phi ^{n-1} \right) .$$

• $$\boldsymbol{Step~2}$$ : To update the auxiliary variable $$q^{n+1}$$ with an energy-optimized step: \begin{align} q^{n+1}=\lambda ^{n+1}Q(\phi^{n+1}),\quad \lambda ^{n+1}=\left\{ \begin{array}{l} 1,E^{n+1}_2\ge 0,\\ \min \left\{ 1,\sqrt{\frac{-E^{n+1}_2}{E^{n+1}_1}} \right\} ,E^{n+1}_2<0, \end{array} \right. \end{align} where \begin{align} E^{n+1}_1=\lVert Q( \phi ^{n+1} ) \rVert ^2,\quad E_{2}^{n+1}=\varepsilon ^2\lVert \Delta \phi ^n \rVert ^2-\lVert q^n \rVert ^2-\varepsilon ^2\lVert \Delta \phi ^{n+1} \rVert ^2. \end{align}

Theorem 7. The scheme 6 is unconditionally energy stable in the sense that \begin{align} \frac{\varepsilon ^2}{2}\lVert \Delta \phi ^{n+1} \rVert ^2-\frac{1}{2}\lVert q^{n+1} \rVert ^2 -\frac{\varepsilon ^2}{2}\lVert \Delta \phi ^n \rVert ^2+\frac{1}{2}\lVert q^n \rVert ^2 =-\frac{1}{M\tau}\lVert \phi ^{n+1}-\phi ^n \rVert ^2\le 0. \end{align}

The proof of the Theorem 7 is not particularly difficult but will be omitted. We define the roughness measurement function $$W(t)$$ [43] as follows $W\left( t \right) =\frac{1}{\left| \Omega \right|^{\frac{1}{2}}}\lVert \phi \left( \boldsymbol{x},t \right) -\bar{\phi}\left( t \right) \rVert , \quad \bar{\phi}\left( t \right) =\frac{1}{|\Omega|}\int_{\Omega}{\phi \left( \boldsymbol{x},t \right)}\,\text{d}\boldsymbol{x}.$

In the last example, we also use the EOP-IEQ method and the baseline IEQ method to solve the MBE equation. Numerical simulations are carried out in the rectangular region $$[-\pi, \pi] ^2$$. We start from the classical benchmark problem $\phi(x,y,0) = 0.1(\sin{3x}\sin{3y}+\cos{5x}\cos{5y}).$ The parameters are $$(Nx,Ny) = (128,128)$$, $$M=0.5$$, $$\varepsilon = 0.1$$, $$C_0=1$$, $$\tau = 1e-2$$ and $$T=50$$. We use the results of the EOP-IEQ/CN scheme of $$\tau=1e-4$$ as the reference solution. Several profiles of simulation process under the EOP-IEQ/CN scheme at $$t=1, 10, 45$$ in Fig. 11. We summarize the normalized numerical results of the modified energy evolution in Fig. 12 (a). It can be observed that the EOP-IEQ/CN scheme provides significantly more accurate results than the baseline IEQ/CN scheme in Fig. 12 (b). Furthermore, Fig. 12 shows that the modified energy and roughness decay very quickly at the beginning, and after a long period of time, the roughness begins to increase. In general, the results are consistent with those in [41].

5 Conclusion↩︎

In this paper, we propose an EOP-IEQ method from the perspective of original energy, and construct a new class of first- and second-order (in time) unconditionally energy stable schemes, which inherit the advantages of the baseline IEQ and REQ schemes. The novel energy-optimized technique avoids the nonlinear optimization problems caused by the correction of auxiliary variables in REQ method, so our process of updating auxiliary variables is simpler and more efficient. Ample numerical examples of gradient flow models, including AC equation, CH equation, PFC equation and MBE equation, verify the accuracy, efficiency and energy stability of the proposed schemes. Furthermore, our proposed numerical schemes are not only linearly solvable, but also have no significant deviation between the numerical solutions and the reference solutions. On the one hand, we observe that the modified energy using energy-optimized technique is significantly more closer to the original energy than the relaxation technique and baseline calculated step. On the other hand, in most cases of long-term numerical simulation, the ratio of nonlinear free energy $$q^n=Q( \phi ^n ), n\ge 0$$ by using the EOP-IEQ method, indicating that the numerical solutions obtained by the EOP-IEQ method preserve the original energy dissipation law. Overall, the main conclusion of this paper is that the energy-optimized technique is applicable to all available IEQ schemes in the literature with the natures of ease to use and improved accuracy.

CRediT authorship contribution statement↩︎

Xiaoqing Meng: Methodology, Software, Validation, Formal analysis, Investigation, Writing-Original Draft, Visualization. Aijie Cheng: Validation, Formal analysis, Resources, Writing-Review & Editing, Supervision. Zhengguang Liu: Conceptualization, Methodology, Validation, Formal analysis, Writing-Review & Editing, Supervision.

Declaration of competing interest↩︎

The authors declare that there are no conflicts of interest.

Data availability↩︎

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Acknowledgment↩︎

This work was supported by Natural Science Outstanding Youth Fund of Shandong Province (Grant No: ZR2023YQ007).

References↩︎

[1]
D. M. Anderson, G. B. McFadden, A. A. Wheeler, Diffuse-interface methods in fluid mechanics, Annual Review of Fluid Mechanics30(1997) 139–165.
[2]
D. Hou, M. Azaiez, C. Xu, A variant of scalar auxiliary variable approaches for gradient flows, Journal of Computational Physics395(2019) 307–332.
[3]
R. Kobayashi, Modeling and numerical simulations of dendritic crystal growth, Physica D: Nonlinear Phenomena63(1993) 410–423.
[4]
K. R. Elder, M. Katakowski, M. Haataja, M. Grant, Modeling elasticity in crystal growth, Physical Review Letters88(2002) 245701.
[5]
X. Bai, J. Sun, J. Shen, W. Yao, Z. Guo, A Ginzburg-Landau-$$H^{-1}$$ model and its SAV algorithm for image inpainting, Journal of Scientific Computing96(2023) 40.
[6]
X. Zhang, D. Zhai, T. Li, Y. Zhou, Y. Lin, Image inpainting based on deep learning: A review, Information Fusion90(2023) 74–94.
[7]
L. Giacomelli, F. Otto, Variatonal formulation for the lubrication approximation of the Hele-Shaw flow, Calculus of Variations and Partial Differential Equations13(2001) 377–403.
[8]
C. Wang, X. Wang, S. M. Wise, Unconditionally stable schemes for equations of thin film epitaxy, Discrete and Continuous Dynamical Systems28(2010) 405–423.
[9]
J. G. E. M. Fraaije, Dynamic density functional theory for microphase separation kinetics of block copolymer melts, Journal of Chemical Physics99(1993) 9202–9212.
[10]
C. M. Elliott, A. M. Stuart, The global dynamics of discrete semilinear parabolic equations, SIAM Journal on Numerical Analysis30(1993) 1622–1663.
[11]
A. Baskaran, J. S. Lowengrub, C. Wang, S. M. Wise, Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation, SIAM Journal on Numerical Analysis51(2013) 2851–2873.
[12]
M. Hochbruck, A. Ostermann, Explicit exponential Runge-Kutta methods for semilinear parabolic problems, SIAM Journal on Numerical Analysis43(2005) 1069–1090.
[13]
M. Hochbruck, A. Ostermann, Exponential integrators, Acta Numerica19(2010) 209–286.
[14]
X. Feng, T. Tang, J. Yang, Stabilized Crank-Nicolson/Adams-Bashforth schemes for phase field models, East Asian Journal on Applied Mathematics3(2013) 59–80.
[15]
J. Shen, X. Yang, Numerical approximations of Allen-Cahn and Cahn-Hilliard equations, Discrete and Continuous Dynamical Systems28(2010) 1669–1691.
[16]
Q. Cheng, C. Liu, J. Shen, A new Lagrange multiplier approach for gradient flows, Computer Methods in Applied Mechanics and Engineering367(2020) 113070.
[17]
X. Yang, Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends, Journal of Computational Physics327(2016) 294–316.
[18]
X. Yang, J. Zhao, X. He, Linear, second order and unconditionally energy stable schemes for the viscous Cahn-Hilliard equation with hyperbolic relaxation using the invariant energy quadratization method, Journal of Computational and Applied Mathematics343(2018) 80–97.
[19]
X. Yang, J. Zhao, Q. Wang, J. Shen, Numerical approximations for a three components Cahn-Hilliard phase-field model based on the invariant energy quadratization method, Mathematical Models and Methods in Applied Sciences27(2017) 1993–2030.
[20]
J. Zhao, Q. Wang, X. Yang, Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach, International Journal for Numerical Methods in Engineering110(2017) 279–300.
[21]
J. Zhao, X. Yang, Y. Gong, Q. Wang, A novel linear second order unconditionally energy stable scheme for a hydrodynamic Q-tensor model of liquid crystals, Computer Methods in Applied Mechanics and Engineering318(2017) 803–825.
[22]
J. Zhao, X. Yang, Y. Gong, X. Zhao, J. Li, Q. Wang, A general strategy for numerical approximations of non-equilibrium models-Part I: Thermodynamical systems, International Journal of Numerical Analysis and Modeling15(2018) 884–918.
[23]
X. Yang, G. Zhang, Convergence analysis for the Invariant Energy Quadratization (IEQ) Schemes for solving the Cahn-Hilliard and Allen-Cahn equations with general nonlinear potential, Journal of Scientific Computing82(2020) 1–28.
[24]
Y. Gong, J. Zhao, Energy-stable Runge-Kutta schemes for gradient flow models using the energy quadratization approach, Applied Mathematics Letters94(2019) 224–231.
[25]
Y. Gong, J. Zhao, Q. Wang, Arbitrarily high-order unconditionally energy stable schemes for thermodynamically consistent gradient flow models, SIAM Journal on Scientific Computing42(2020) B135–B156.
[26]
J. Shen, J. Xu, J. Yang, A new class of efficient and robust energy stable schemes for gradient flows, SIAM Review61(2019) 474–506.
[27]
Z. Liu, X. Li, The exponential scalar auxiliary variable (E-SAV) approach for phase field models and its explicit computing, SIAM Journal on Scientific Computing42(2020) B630–B655.
[28]
F. Huang, J. Shen, Z. Yang, A highly efficient and accurate new scalar auxiliary variable approach for gradient flows, SIAM Journal on Scientific Computing42(2020) A2514–A2536.
[29]
M. Jiang, Z. Zhang, J. Zhao, Improving the accuracy and consistency of the scalar auxiliary variable (SAV) method with relaxation, Journal of Computational Physics456(2022) 110954.
[30]
J. Zhao, A revisit of the energy quadratization method with a relaxation technique, Applied Mathematics Letters120(2021) 107331.
[31]
Y. Zhang, X. Li, Efficient and accurate exponential SAV algorithms with relaxation for dissipative system, Communications in Nonlinear Science and Numerical Simulation127(2023) 107530.
[32]
Q. Huang, C. Yuan, G. Zhang, L. Zhang, A computationally optimal relaxed scalar auxiliary variable approach for solving gradient flow systems, Computers & Mathematics with Applications156(2024) 64–73.
[33]
Z. Liu, Y. Zhang, X. Li, A novel energy-optimal scalar auxiliary variable (EOP-SAV) approach for gradient flows, 2023. http://arxiv.org/abs/2304.11288.
[34]
Y. Gong, J. Zhao, Energy-stable Runge-Kutta schemes for gradient flow models using the energy quadratization approach, Applied Mathematics Letters94(2019) 224–231.
[35]
X. Yang, L. Ju, Efficient linear schemes with unconditional energy stability for the phase field elastic bending energy model, Computer Methods in Applied Mechanics and Engineering315(2017) 691–712.
[36]
F. Huang, J. Shen, Z. Yang, A highly efficient and accurate new scalar auxiliary variable approach for gradient flows, SIAM Journal on Scientific Computing42(2020) A2514–A2536.
[37]
S. M. Allen, J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metallurgica27(1979) 1085–1095.
[38]
Q. Cheng, C. Liu, J. Shen, Generalized SAV approaches for gradient systems, Journal of Computational and Applied Mathematics394(2021) 113532.
[39]
Z. Zhang, Y. Ma, Z. Qiao, An adaptive time-stepping strategy for solving the phase field crystal model, Journal of Computational Physics249(2013) 204–215.
[40]
X. Yang, D. Han, Linearly first- and second-order, unconditionally energy stable schemes for the phase field crystal model, Journal of Computational Physics330(2017) 1116–1134.
[41]
X. Yang, J. Zhao, Q. Wang, Numerical approximations for the molecular beam epitaxial growth model based on the invariant energy quadratization method, Journal of Computational Physics333(2017) 104–127.
[42]
W. Chen, Z. Chen, J. Cheng, Y. Wang, New epitaxial thin-film models and numerical approximation, Mathematical Methods in the Applied Sciences40(2015) 3855–3881.
[43]
B. Li, J. Liu, Thin film epitaxy with or without slope selection, European Journal of Applied Mathematics14(2003) 713–743.