A Second-Order, Weakly Energy-Stable Pseudo-spectral Scheme for the Cahn–Hilliard Equation and Its Solution by the Homogeneous Linear Iteration Method

Abstract

We present a second order energy stable numerical scheme for the two and three dimensional Cahn–Hilliard equation, with Fourier pseudo-spectral approximation in space. A convex splitting treatment assures the unique solvability and unconditional energy stability of the scheme. Meanwhile, the implicit treatment of the nonlinear term makes a direct nonlinear solver impractical, due to the global nature of the pseudo-spectral spatial discretization. We propose a homogeneous linear iteration algorithm to overcome this difficulty, in which an $O(s^2)$ (where s the time step size) artificial diffusion term, a Douglas–Dupont-type regularization, is introduced. As a consequence, the numerical efficiency can be greatly improved, since the highly nonlinear system can be decomposed as an iteration of purely linear solvers, which can be implemented with the help of the FFT in a pseudo-spectral setting. Moreover, a careful nonlinear analysis shows a contraction mapping property of this linear iteration, in the discrete $\ell ^4$ norm, with discrete Sobolev inequalities applied. Moreover, a bound of numerical solution in $\ell ^\infty$ norm is also provided at a theoretical level. The efficiency of the linear iteration solver is demonstrated in our numerical experiments. Some numerical simulation results are presented, showing the energy decay rate for the Cahn–Hilliard flow with different values of $\varepsilon$.

Appendices

Proof of Lemma 2.1

For a 3-D grid function f with its discrete Fourier expansion as (5), an application of discrete Parseval equality gives

$$\begin{aligned}&\Vert f \Vert _{-1,N}^2 = \Vert ( - \Delta _N )^{-\frac{1}{2}} f \Vert _2^2 = {\mathop {\mathop {\sum }\limits _{{\ell ,m,n =-K}}}\limits _{(\ell ,m,n) \ne \mathbf{0}}^K} \lambda _{\ell , m, n}^{-1} | \hat{f}_{\ell ,m,n} |^2 , \end{aligned}$$

(99)

$$\begin{aligned}&\Vert \nabla _N f \Vert _2^2 = {\mathop {\mathop {\sum }\limits _{{\ell ,m,n =-K}}}\limits _{(\ell ,m,n) \ne \mathbf{0}}^K} \lambda _{\ell , m, n} | \hat{f}_{\ell ,m,n} |^2 , \end{aligned}$$

(100)

with $\lambda _{\ell , m, n} = 4 \pi ^2 ( \ell ^2 + m^2 + n^2 )$. This in turn yields

$$\begin{aligned} \gamma _1 \Vert f \Vert _{-1,N}^2 + \gamma _2 \Vert \nabla _N f \Vert _2^2 = {\mathop {\mathop {\sum }\limits _{{\ell ,m,n =-K}}}\limits _{(\ell ,m,n) \ne \mathbf{0}}^K} \left( \gamma _1 \lambda _{\ell , m, n}^{-1} + \gamma _2 \lambda _{\ell , m, n} \right) | \hat{f}_{\ell ,m,n} |^2 . \end{aligned}$$

(101)

A similar calculation also gives

$$\begin{aligned} \Vert f \Vert _{H_N^{\alpha _0}}^2 = {\mathop {\mathop {\sum }\limits _{{\ell ,m,n =-K}}}\limits _{(\ell ,m,n) \ne \mathbf{0}}^K} {\alpha _0} \lambda _{\ell , m, n}^{\alpha _0} | \hat{f}_{\ell ,m,n} |^2 , \quad \forall 0< \alpha _0 < 1 . \end{aligned}$$

(102)

By making comparison between (101) and (102), we conclude that (19) is a direct consequence of the following application of Young’s inequality:

$$\begin{aligned} \gamma _1 \lambda _{\ell , m, n}^{-1} + \gamma _2 \lambda _{\ell , m, n} \ge C^* \gamma _1^{\frac{1-\alpha _0}{2}} \gamma _2^{\frac{1+\alpha _0}{2}} \lambda _{\ell , m, n}^{\alpha _0} , \quad \forall -K \le \ell , m, n \le K , \end{aligned}$$

(103)

with $C^*$ only dependent on $\alpha _0$ and $\Omega$.

For the proof of (20), a discrete version of Sobolev embedding from $H^{\alpha _0}$ into $\ell ^4$, we have to utilize the continuous extension of f, given by (22). For simplicity of presentation, we focus our analysis in the 2-D case; for the 3-D grid function, the analysis could be carried out in a similar, yet more tedious way. And also, $\Vert \cdot \Vert$ is denoted as the standard $L^2$ norm for a continuous function.

We denote the following grid function

$$\begin{aligned} g_{i,j} = \left( f_{i,j} \right) ^2 . \end{aligned}$$

(104)

A direct calculation shows that

$$\begin{aligned} \left\| f \right\| _4 = \left( \left\| g \right\| _2 \right) ^{\frac{1}{2}} . \end{aligned}$$

(105)

Note that both norms are discrete in the above identity. Moreover, we assume the grid function g has a discrete Fourier expansion as

$$\begin{aligned} g_{i,j} = \sum _{\ell ,m=-K}^{K} (\hat{g}^N_c)_{\ell ,m} \mathrm {e}^{2\pi \mathrm{i} (\ell x_i + m y_j )} , \end{aligned}$$

(106)

and denote its continuous version as

$$\begin{aligned} G (x,y) = \sum _{\ell ,m=-K}^{K} (\hat{g}^N_c)_{\ell ,m} \mathrm {e}^{2\pi \mathrm{i} (\ell x + m y )} \in \mathcal{P}_{K} . \end{aligned}$$

(107)

With an application of the Parseval equality at both the discrete and continuous levels, we have

$$\begin{aligned} \left\| g \right\| _2^2 = \left\| G \right\| ^2 = \sum _{\ell ,m=-K}^{K} \left| (\hat{g}^N_c)_{\ell ,m} \right| ^2 . \end{aligned}$$

(108)

On the other hand, we also denote

$$\begin{aligned} H (x,y) = \left( f_S (x,y) \right) ^2 = \sum _{\ell ,m=-2K}^{2K} (\hat{h}^N)_{\ell ,m} \mathrm {e}^{2\pi \mathrm{i} (\ell x + m y )} \in \mathcal{P}_{2K} . \end{aligned}$$

(109)

The reason for $H \in \mathcal{P}_{2K}$ is because $f_N \in \mathcal{P}_{K}$. We note that $H \ne G$, since $H \in \mathcal{P}_{2K}$, while $G \in \mathcal{P}_{K}$, although H and G have the same interpolation values on at the numerical grid points $(x_i, y_j)$. In other words, g is the interpolation of H onto the numerical grid point and G is the continuous version of g in $\mathcal{P}_{K}$. As a result, collocation coefficients $\hat{g}_c^N$ for G are not equal to $\hat{h}^N$ for H, due to the aliasing error. In more detail, for $- K \le \ell , m \le K$, we have the following representations:

$$\begin{aligned} ( \hat{g}_c^N )_{\ell ,m} = \left\{ \begin{array}{l} (\hat{h}^N)_{\ell ,m} + (\hat{h}^N)_{\ell +N,m} + (\hat{h}^N)_{\ell ,m+N} + (\hat{h}^N)_{\ell +N,m+N} , \, \, \ell< 0 , m< 0 ,\\ (\hat{h}^N)_{\ell ,m} + (\hat{h}^N)_{\ell +N,m} , \, \, k< 0 , m = 0,\\ (\hat{h}^N)_{\ell ,m} + (\hat{h}^N)_{\ell +N,m} + (\hat{h}^N)_{\ell ,m-N} + (\hat{h}^N)_{\ell +N,m-N} , \, \, \ell< 0 , m> 0 ,\\ (\hat{h}^N)_{\ell ,m} + (\hat{h}^N)_{\ell -N,m} + (\hat{h}^N)_{\ell ,m-N} + (\hat{h}^N)_{\ell -N,m-N} , \, \, \ell> 0 , m> 0 ,\\ (\hat{h}^N)_{\ell ,m} + (\hat{h}^N)_{\ell -N,m} , \, \, \ell> 0 , m = 0 ,\\ (\hat{h}^N)_{\ell ,m} + (\hat{h}^N)_{\ell -N,m} + (\hat{h}^N)_{\ell ,m+N} + (\hat{h}^N)_{\ell -N,m+N} , \, \, \ell> 0 , m< 0 ,\\ (\hat{h}^N)_{\ell ,m} + (\hat{h}^N)_{\ell ,m+N} , \, \, \ell = 0 , m < 0 ,\\ (\hat{h}^N)_{\ell ,m} , \, \, \ell = 0 , m = 0 ,\\ (\hat{h}^N)_{\ell ,m} + (\hat{h}^N)_{\ell ,m-N} , \, \, \ell = 0 , m > 0 . \end{array} \right. \end{aligned}$$

(110)

With an application of Cauchy inequality, it is clear that

$$\begin{aligned} \sum _{\ell ,m=-K}^{K} \left| (\hat{g}^N_c)_{\ell ,m} \right| ^2 \le 4 \left| \sum _{\ell ,m=-2K}^{2K} (\hat{h}^N)_{\ell ,m} \right| ^2 . \end{aligned}$$

(111)

Meanwhile, an application of Parseval’s identity to the Fourier expansion (109) gives

$$\begin{aligned} \left\| H \right\| ^2 = \left| \sum _{\ell ,m=-2K}^{2K} (\hat{h}^N)_{\ell ,m} \right| ^2 . \end{aligned}$$

(112)

Its comparison with (108) indicates that

$$\begin{aligned} \left\| g \right\| _2^2 = \left\| G \right\| ^2 \le 4 \left\| H \right\| ^2 , \quad \text{ i.e. } \, \, \left\| g \right\| _2 \le 2 \left\| H \right\| , \end{aligned}$$

(113)

with the estimate (111) applied. Meanwhile, since $H (x,y) = \left( f_N (x,y) \right) ^2$, we have

$$\begin{aligned} \left\| f_N \right\| _{L^4} = \left( \left\| H \right\| _2 \right) ^{\frac{1}{2}} . \end{aligned}$$

(114)

Therefore, a combination of (105), (113) and (114) results in

$$\begin{aligned} \left\| f \right\| _4 = \left( \left\| g \right\| _2 \right) ^{\frac{1}{2}} \le \left( 2 \left\| H \right\| _{L^2} \right) ^{\frac{1}{2}} \le \sqrt{2} \left\| f_N \right\| _{L^4} . \end{aligned}$$

(115)

For the continuous function $f_N (x,y)$, we have the following estimate in Sobolev embedding (in 2-D):

(116)

Then we arrive at

$$\begin{aligned} \left\| f \right\| _4 \le \sqrt{2} \left\| f_N \right\| _{L^4} \le C_0 \Vert (-\Delta _N)^{\frac{1}{2}} f \Vert _2 . \end{aligned}$$

(117)

This finishes the proof of (20) for $d=2$.

The 3-D case could be analyzed in the same fashion, and the details are skipped for the sake of brevity. The proof of Lemma 2.1 is completed.

Proof of Lemma 2.2

For any grid function $f\in {\mathcal {G}}_N$, we recall its continuous extension, $f_S = S_N(f)\in {\mathcal {P}}_K$, as defined in (22). Since f is the point-wise grid interpolation of $f_S$, we have

$$\begin{aligned} \left\| f \right\| _\infty \le \left\| f_S \right\| _{L^\infty } . \end{aligned}$$

(118)

For the smooth function $f_S$, applying the 3-D Sobolev inequality associated to the embedding $H^2\hookrightarrow L^\infty$ and elliptic regularity, we have

$$\begin{aligned} \left\| f_S \right\| _{L^\infty } \le C \left( \left| \int _\Omega f_S \, d {\varvec{x}}\right| + \left\| \Delta f_S \right\| _{L^2} \right) . \end{aligned}$$

(119)

Subsequently, the maximum norm estimate (21) is a direct consequence of the following identities:

$$\begin{aligned} h^3 \sum _{i,j,k=0}^{N-1}f_{i,j,k} =: \overline{f} = \overline{f_S} := \int _\Omega \, f_S \, d {\varvec{x}}, \quad \Vert \Delta _N f \Vert _2 = \Vert \Delta f_S \Vert _{L^2} . \end{aligned}$$

(120)

This finishes the proof of Lemma 2.2.