# A Fast and Efficient Numerical Algorithm for the Nonlocal Conservative Swift–Hohenberg Equation

Authors: Jingying Wang; Shuying Zhai
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, we consider a new Swift–Hohenberg equation, where the total mass of this model is conserved through a nonlocal Lagrange multiplier. Based on the operator splitting method and spectral method, a fast and efficient numerical algorithm is proposed. Three numerical examples in both two and three dimensions are provided to illustrate that the proposed algorithm is a practical, accurate, and efficient simulation tool for the nonlocal Swift–Hohenberg equation.


---

## Body

## 1. Introduction

The Swift–Hohenberg model [1] was first introduced by Swift and Hohenberg in studies of Rayleigh–Bénard convection and has become one of the paradigms of the nonlinear dynamical system leading to complex pattern formation. Recent applications of the Swift–Hohenberg equation have been extensively spanned in complex fluids and biological tissues [2, 3]. The Swift–Hohenberg model is derived by using the L2-gradient flow of the Lyapunov energy functional [4](1)Eϕ=∫Ωϕ2Δ+12ϕ+Fϕdx,where Ω is a bounded regular domain a,bn (n = 1, 2, 3) and ϕ is the density field. Fu=1/4ϕ4−η/3ϕ3−ϵ/2ϕ2 is the nonlinear smoothing potential, and η≥0 and ϵ>0 are constants with physical significance.By taking the variational derivative of (1) in L2 norm, the following Swift–Hohenberg equation is obtained:(2)ϕt=−δEϕδϕ=−Δ+12ϕ+fϕ,x,t∈Ω×0,T,where δEϕ/δϕ denotes the variational derivative and fu=F′u=ϕ3−ηϕ2−ϵϕ. Thus, the total energy Eϕ is not increasing with time(3)ddtEϕ=∫ΩδEϕδϕ∂ϕ∂tdx=−∫Ω∂ϕ∂t2dx≤0.The Swift–Hohenberg equation (2) is a fourth-order nonlinear partial differential equation, and we cannot generally obtain the analytic solution for the arbitrary initial conditions. Therefore, computer simulations play an essential role in understanding the dynamical behavior. Till now, a large number of works have been proposed to study the Swift–Hohenberg equation [5–12].Since the Swift–Hohenberg model (2) is derived by using the L2-gradient flow, the total mass of the unknown variable is not conserved with time. However, in Elder’s phase field crystal model, the phase field variable is interpreted as the density of atoms from the perspective of physical significance; thus, the total mass of the phase field variable is expected to be conserved, such as the Cahn–Hilliard equation [13, 14]. Zhang and Yang [15] very recently proposed a conservative Swift–Hohenberg equation with a nonlocal Lagrange multiplier to cancel out the variation in mass without influencing the original energy law. The conservative Swift–Hohenberg equation is as follows:(4)ϕt+Δ+12ϕ+fϕ−1Ω∫Ωϕ+fϕdx=0,x,t∈Ω×0,T,subject to the initial condition(5)ϕx,0=ϕ0x,x∈Ω¯,The periodic boundary condition is adopted here since it is usually used to remove all boundary integrals (cf. [15]). Note that when 1/Ω∫Ωϕ+fϕdx=0, equation (4) degenerates back to the standard Swift–Hohenberg model (2).The nonlocal term1/Ω∫Ωϕ+fϕdx can be understood as a nonlocal Lagrange multiplier for the mass constraint(6)ddtMϕ=ddt∫Ωϕx,tdx=0.Furthermore, system (4) still satisfies energy dissipative law (3). Taking the L2 inner product of (4) with ϕt, integrating by parts and using (6), we have(7)ddtEϕ=∫ΩϕΔ+12+fϕϕtdx=∫ΩϕtΔ+12ϕ+fϕdx−1Ω∫Ωϕ+fϕdx∫Ωϕtdx=−∫Ωϕt2dx≤0,where Eϕ is defined by (1).Compared to a large number of studies for the standard Swift–Hohenberg equation (2), there are few numerical results on the conservative Swift–Hohenberg equation (4) because the extra nonlocal term brings up substantial numerical challenges to developing fast and efficient algorithms. Zhang and Yang [15] proposed a second-order, unconditionally energy stable scheme by combining the invariant energy quadratization (IEQ) approach with the stabilization technique. The law of conservation of mass of the proposed scheme was not verified in this article [15].In this work, we aim to develop a fast explicit algorithm based on the operator splitting method [14, 16, 17]. The main idea of the method is to solve problem (4) through three parts: the linear homogeneous heat equation was first solved via the Fourier spectral method and then the nonlinear part and the nonlocal equation were solved explicitly.The rest of this paper is organized as follows. In Section2, a fast explicit operator splitting method is proposed for solving (4). In Section 3, we test the accuracy and applicability of the obtained numerical method in three examples including a 3D example. Finally, we summarize the present work and discuss the future work in Section 4.


## 2. Time-Splitting Spectral Method

In this section, we present a time-splitting spectral method to simulate the asymptotic behavior of the solution of equation (4). The proposed method is based on the operator splitting method for the time and spectral method for the space. For simplicity, we only consider two-dimensional (2D) space. The 1D/3D case is defined analogously.The operator splitting method for the conservative Swift–Hohenberg equation is constructed as follows. We rewrite equation (4) in the following abstract initial value problem:(8)ut=Au+ℬu+Cu,where the operators A, ℬ, and C are defined as(9)Aϕ=−Δ+12ϕ,ℬu=−fϕ,Cu=1Ω∫Ωϕ+fϕdx.Thus, we can evolve the solution fromt to t+τ using the following first-order Lie–Trotter method [16]:(10)ϕx,y,t+τ=SAτ∘SBτ∘SCτϕx,y,t.where SAτ, SBτ, and SCτ are the exact solution operators associated with the nonlinear equation(11)ϕt=−Δ+12ϕ,the heat equation(12)ϕt=−fϕ,and the Lagrange multiplier equation(13)ϕt=1Ω∫Ωϕ+fϕdx,respectively. Higher order splitting schemes can be found in [18, 19].Before proceeding to present a fast explicit operator splitting method for (4), the exact solution operators SAτ, SBτ, and SCτ are first replaced by their numerical approximations SAτ,h, SBτ,h, and SCτ,h, respectively.Step 1 (SAτ⟶SAτ,h).
We first consider the fully discrete scheme for (11). The following spatial grid is introduced:(14)Ωhper=xi,yj=a+ih,a+jh,0≤i,j≤N−1,for the periodic condition, where h=b−a/N. Let uj,kk=uxi,yj,kτ, where τ=T/K is the time step and M is the total number of time steps.
Based on the theory of spectral method [20], the solution to equation (11) derived with periodic boundary conditions can be represented by the following band-limited Fourier series:(15)ℱN−1ϕ˜tx,y≔PNϕx,y,t=∑p=−N/2N/2∑q=−N/2N/2ϕ˜p,qtφp,qx,y,where(16)φp,q=expi2pπx−ab−a+i2qπy−ab−a,and by using the following Fourier expansion, the fast Fourier transform (FFT) may be used to calculate the Fourier coefficients {ϕ˜p,qt} from the discrete function values {ϕi,jt}:(17)ℱNϕtp,q≔ϕ˜p,qt=h2CpperCqperb−a2∑j=0N−1∑k=0N−1ϕi,jtφp,qxi,yj,and Cpper and Cqper are, respectively, defined as(18)Crper=2,r=N2,1,r<N2.
Substituting (15) into (11), we obtain the following system of ordinary differential equations for the discrete pq-th Fourier coefficients:(19)dϕ˜p,qtdt=−λp,q+12ϕ˜p,qt,where(20)λp,q=−2pπb−a2+2qπb−a2,p,q=0,±1,±2,…,±N2.
The above equation is solved analytically by the variable separation method with an initial conditionϕ^pqk, and then the solution ϕ^pqk+1 is obtained as follows:(21)ϕ˜p,qk+1=exp−τλp,q+12ϕ˜p,qk.
Hence, we have(22)ϕk+1=ℱN−1exp−τλp,q+12ℱNϕkp,q.
The errors betweenPNϕ and ϕ can be measured in Sobolev spaces. For this purpose, we denote by HpsΩ the subspace of HsΩ, which consists of functions with derivatives of order up to s−1 being periodic on Ω. The main approximation result is stated below.Theorem 1 (see [20]).
For anyϕ∈HpsΩ and s≥0, it holds that(23)PNϕ−ϕ≤CN−sϕs.This theorem indicates that the projectionPNϕ is the best approximation of ϕ in all Sobolev spaces HpsΩ. For further references about the spectral method, we refer to [20].Step 2 (SBτ⟶SBτ,h).
Now we propose a first-order difference scheme for solving equation (12) in tk,tk+1. Note that if ϕ≠0, equation (12) can be rewritten as(24)ϕtϕ=−ϕ2−ηϕ−ϵ.
Taking the integral fromtk to tk+1 on both sides of the above equation with respect to t, we have(25)lnϕk+1−lnϕk=−∫tktk+1ϕ2−ηϕ−ϵdt.
Note that∫tktk+1ϕ2−ηϕ−ϵdt=τϕk2−ηϕk−ϵ+Oτ2, and then we obtain a first-order explicit finite difference scheme for (12)(26)ϕk+1=ϕkexp−τϕk2−ηϕk−ϵ.
It is clear that ifϕ=0, the above formula is also derived.Step 3 (SCτ⟶SCτ,h).
Equation (12) involves solving a nonlocal Lagrange multiplier equation. For convenience, we denote βϕ=∫Ωϕ+fϕdx. Then, we employ the backward difference method for equation (13) in tk,tk+1 with an initial condition ϕk; this gives(27)ϕk+1−ϕkτ=1Ωβϕk+1.
Summing both sides of equation (27) over i and j, and by the property of mass conservation h2∑i,j=0N−1ϕij0=h2∑i,j=0N−1ϕijk+1, we obtain(28)βϕk+1=h2τ∑i,j=0N−1ϕij0−∑i,j=0N−1ϕijk.
A similar idea has also been mentioned in [21].
Combining (27) and (28), a first-order explicit finite difference scheme for (13) is obtained:(29)ϕk+1=ϕk+h2Ω∑i,j=0N−1ϕij0−∑i,j=0N−1ϕijk.
Then, based on equations (10), (22), (26), and (29), we obtain the following fast explicit operator splitting schemes:(30)ϕk+1,1=ℱNexp−τλp,q+12ℱN−1ϕkp,q,ϕk+1,2=ϕk+1,1exp−τϕk2−ηϕk−ϵ,ϕk+1=ϕk+1,2+h2Ω∑i,j=0N−1ϕij0−∑i,j=0N−1ϕijk,2.
The derivation of algorithm (30) shows that the temporal truncation error has the form C1τ+C2τ2. Based on this observation, we can use the following Richardson’s extrapolation technique:(31)ϕwK=2ϕ2K−ϕK,to eliminate the leading term Oτ, where ϕK and ϕ2K are the numerical solutions at the final time level calculated using time step sizes τ=T/K and τ=T/2K, respectively. Thus, we obtain a high-order numerical solution ϕwK with second-order temporal accuracy.
In addition, we define the discrete inner product(32)f,g≔h2∑i,j=0N−1fi,jgi,jand discrete ⋅∞ norm and ⋅2 norm(33)f∞=max0≤i,j≤N−1fi,j,f2=f,f.Theorem 2.
Scheme (30) satisfies the discrete mass conservation law(34)Mk=M0,k=1,2,…,K,where Mk=h2ϕi,jk.Proof.
Taking discrete inner product of the third equation of (30) with 1, summing over i and j, we have(35)ϕk+1,1=ϕk+1,2,1+h2Ω∑i,j=0N−1ϕij0−∑i,j=0N−1ϕijk,2,1=ϕ0,1.
This completes the proof.


## 3. Numerical Experiments

In this section, three numerical experiments are carried out to study the accuracy and efficiency of the proposed method. TheL∞-norm error (denoted as ErrL∞) and the L2-norm error (denoted as ErrL2) of the numerical solutions are computed, respectively, by the following formulas:(36)ErrL∞=ϕM+1−ϕx,y,T∞,ErrL2=ϕM+1−ϕx,y,T2.Problem 1 (convergence test).
In order to verify the temporal numerical accuracy, we consider an example on0,642×0,1 with various ϵ and η, and the initial condition is as follows:(37)ϕ0x,y=sin8π128xcos8π128y.
We useh=1/2, so the spatial discretization errors are negligible compared with the time discretization error. The computational results at T=1 for various ϵ and η are presented in Table 1. One may see that the convergence rate is close to 2 for all cases, which coincides with what we expected.Table 1
Convergence rate withh=1/2 and different M for Problem 1 at T=1.
M
=
10
20
40
80
160
ϵ = 0.05η = 0.05
Err
L
∞
8.14e − 5
1.78e − 5
4.18e − 6
1.01e − 6
2.49e − 7
Rate
L
∞
—
2.19
2.09
2.05
2.02
Err
L
2
1.30e − 3
2.95e − 4
7.04e − 5
1.72e − 5
4.25e − 6
Rate
L
2
—
2.14
2.07
2.03
2.02
ϵ = 0.25η = 2
Err
L
∞
4.61e − 3
1.10e − 3
2.66e − 4
6.52e − 5
1.62e − 5
Rate
L
∞
—
2.07
2.05
2.03
2.01
Err
L
2
3.69e − 2
9.16e − 3
2.27e − 3
5.63e − 4
1.40e − 4
Rate
L
2
—
2.01
2.01
2.01
2.01
ϵ = 0.5η = 0.5
Err
L
∞
1.25e − 4
3.30e − 5
8.43e − 6
2.13e − 6
5.34e − 7
Rate
L
∞
—
1.92
1.97
1.99
1.99
Err
L
2
1.78e − 3
4.23e − 4
1.03e − 4
2.54e − 5
6.32e − 6
Rate
L
2
—
2.08
2.04
2.02
2.01Problem 2 (effects ofϵ and η on energy stability and pattern formation).
Various related works have been proposed for the Swift–Hohenberg equation, and most of them considered the case ofη=0 except [9, 10]. In order to show the effects of ϵ and η on pattern formation, we use a benchmark simulation for the conservative Swift–Hohenberg equation, i.e., the phase transition behaviors. We set the 2D computational domain as 0,1282 and discretize it by using h=1. The initial condition reads as(38)ϕ0x,y=0.071+randx,y,where randx,y is a random number in the range of −1,1.
We first perform a numerical experiment withϵ=0.025 and η=0.1 for the stability study of the proposed method. Figure 1 shows the x−y−ϕ plots at T=2000 with different time steps, and the corresponding curves of mass M and energy E are also presented. These results suggest that the proposed method satisfies mass conservation and is unconditionally stable. The only drawback is that the method is conditionally energy stable.
Secondly, with the same parameters andτ=0.1, we present the snapshots of phase transition behaviors of the density field ϕx,y,t at different computational time steps in Figure 2.
Finally, the relationship between energy stability and the parametersϵ and η is examined. Figure 3 shows the evolution of the energies using different ϵ, η, and time steps. The results in Figure 3 suggest that the proposed method is conditionally energy stable. We can see from this figure that the energy Eϕ is mainly affected by the value of ϵ rather than η, and small time step τ is needed to keep the energy law (7) for large ϵ.Figure 1
ϕ
x
,
y
,
2000 and the curves of mass M and energy E with different time steps τ for Problem 2: (a) τ=10, (b) τ=5, (c) τ=1, (d) τ=0.5, and (e) τ=0.1.
(a)
(b)
(c)
(d)
(e)Figure 2
Snapshots ofϕ at different time steps for Problem 2: (a) ϕx,y,0, (b) ϕx,y,100, (c) ϕx,y,200, (d) ϕx,y,600, (e) ϕx,y,1000, and (f) ϕx,y,2000.
(a)
(b)
(c)
(d)
(e)
(f)Figure 3
The curves of energyE with different ϵ, η, and time steps τ for Problem 2: (a) ϵ=0.025,η=0.1, (b) ϵ=0.025,η=0.5, (c) ϵ=0.025,η=0.99, (d) ϵ=0.25,η=0.1, (e) ϵ=0.25,η=0.5, (f) ϵ=0.25,η=0.99, (g) ϵ=0.5,η=0.1, (h) ϵ=0.5,η=0.5, and (i) ϵ=0.5,η=0.99.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Problem 3 (the conservative Swift–Hohenberg equation in 3D).
We consider the evolution ofϕ in 3D space using the proposed method. The initial condition on 0,403 is(39)ϕ0x,y,z=sin8πx128cos8πy128cos8πz128.
We take the simulation parameters(40)T=1000,ϵ=0.25,η=0.05,h=1,τ=13.
In Figure4, we display the evolution of the phase field at different time periods. From the results, it can be seen that our method is applicable to 3D conservative Swift–Hohenberg equation. Moreover, the curves of mass M and energy E are also presented in Figure 5.Figure 4
Snapshots ofϕ at different time periods for Problem 3: (a) ϕx,y,10, (b) ϕx,y,100, (c) ϕx,y,300, (d) ϕx,y,400, and (e) ϕx,y,1000.
(a)
(b)
(c)
(d)
(e)Figure 5
The curves of (a) massM and (b) energy E for Problem 3.
(a)
(b)


## 4. Conclusions

In this paper, we developed a second-order, linear, and mass conservation algorithm for solving the nonlocal Swift–Hohenberg equation by combining the operator splitting technique with the Fourier spectral method. A variety of numerical experiments were presented to demonstrate the stability and the accuracy of the developed method. In particular, it is applicable to various problems with a mass conservation constraint.



---

*Source: 7012483-2020-01-06.xml*