# A linear adaptive second-order backward differentiation formulation scheme
for the phase field crystal equation∗

Dianming Hou1 Zhonghua Qiao2

###### Abstract.

In this paper, we present and analyze a linear fully discrete second order
scheme with variable time steps for the phase field crystal equation. More
precisely, we construct a linear adaptive time stepping scheme based on the
second order backward differentiation formulation (BDF2) and use the Fourier
spectral method for the spatial discretization. The scalar auxiliary variable
approach is employed to deal with the nonlinear term, in which we only adopt a
first order method to approximate the auxiliary variable. This treatment is
extremely important in the derivation of the unconditional energy stability of
the proposed adaptive BDF2 scheme. However, we find for the first time that
this strategy will not affect the second order accuracy of the unknown phase
function $\phi^{n}$ by setting the positive constant $C_{0}$ large enough such
that $C_{0}\geq 1/\tau.$ The energy stability of the adaptive BDF2 scheme is
established with a mild constraint on the adjacent time step radio
$\gamma_{n+1}:=\tau_{n+1}/\tau_{n}\leq 4.8645$. Furthermore, a rigorous error
estimate of the second order accuracy of $\phi^{n}$ is derived for the
proposed scheme on the nonuniform mesh by using the uniform $H^{2}$ bound of
the numerical solutions. Finally, some numerical experiments are carried out
to validate the theoretical results and demonstrate the efficiency of the
fully discrete adaptive BDF2 scheme.

###### Key words and phrases:

Phase field crystal equation, linear adaptive BDF2 scheme, scalar auxiliary
variable approach, unconditional energy stability, convergence analysis

###### 1991 Mathematics Subject Classification:

35K55, 65M12, 65M15, 65F30

∗The work of D. Hou is supported by NSFC grant 12001248, the NSF of the
Jiangsu Province grant BK20201020, the NSF of Universities in Jiangsu Province
of China grant 20KJB110013 and the Hong Kong Polytechnic University grant
1-W00D. Z. Qiao’s work is partially supported by the Hong Kong Research Grant
Council RFS grant RFS2021-5S03 and GRF grant 15302919, the Hong Kong
Polytechnic University internal grant 4-ZZLS, and the CAS AMSS-PolyU Joint
Laboratory of Applied Mathematics.
1School of Mathematics and Statistics, Jiangsu Normal University, 221116
Xuzhou, China. Current address: Department of Applied Mathematics, The Hong
Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong.
(dmhou@stu.xmu.edu.cn).
2Corresponding author. Department of Applied Mathematics, The Hong Kong
Polytechnic University, Hung Hom, Kowloon, Hong Kong. (zqiao@polyu.edu.hk).

##  1\. Introduction

The phase field crystal (PFC) model was firstly proposed in [5, 6] and has
been frequently used for modelling the crystal growth at the atomic sale in
space but on diffusive scale in time. It can simulate nucleated crystallites
at arbitrary locations and orientations with the elastic and plastic
deformations. The PFC model can account for incorporating the periodic nature
of the crystal lattices through the following Swift-Hohenberg type free
energy, defined by

$E(\phi)=\int_{\Omega}\Big{(}\frac{1}{2}\phi(\Delta+\beta)^{2}\phi+\frac{1}{4}\phi^{4}-\frac{\varepsilon}{2}\phi^{2}\Big{)}d{\tenbi
x},$ (1.1)

which is minimized by periodic density fields. Here
$\Omega\subset\mathbb{R}^{d},d=2$ or $3$, $\phi$ is the atom density field,
$\beta$ and $\varepsilon$ are two positive constants such that
$0<\varepsilon<\beta^{2}$. The operator $(\beta+\Delta)^{2}\phi$ is defined by
$(\Delta+\beta)^{2}\phi=\Delta^{2}\phi+2\beta\Delta\phi+\beta^{2}\phi.$ Thus
the conserved phase field crystal equation reads

$\begin{array}[]{r@{}l}\begin{cases}\displaystyle\frac{\partial\phi}{\partial
t}=\Delta\mu,&\mbox{in}\quad\Omega\times(0,T),\\\
\displaystyle\mu=(\Delta+\beta)^{2}\phi+\phi^{3}-\varepsilon\phi,&\mbox{in}\quad\Omega\times(0,T),\\\
\displaystyle\phi({\tenbi x},0)=\phi_{0}({\tenbi
x}),&\mbox{in}\quad\Omega.\end{cases}\end{array}$ (1.2)

Due to the gradient structure of (1.2), the above PFC system satisfies the
following energy dissipation law

$\frac{d}{dt}E(\phi)=-\|\nabla\mu\|^{2}\leq 0.$ (1.3)

This indicates the free energy $E(\phi)$ always decreases in time. In
numerical simulations, it is crucial to design numerical schemes preserving
the above energy dissipation law at the discrete level, see e.g., [1, 7, 8,
10, 4, 30, 34, 26, 27, 18, 17, 22] and references therein. These schemes are
often called energy stable schemes in practice. A direct implicit/explicit
treatment for the nonlinear term in (1.2) usually yields a severe stability
restriction on the time step sizes to preserve the discrete form of the energy
law (1.3). In the last two decades, many efforts have been made to relax this
restriction. For example, Wise, Wang and Lowengrub developed a first order
convex splitting scheme for the PFC equation in [29]. The unconditional energy
stability and the corresponding error estimate of their proposed scheme were
rigorously established. The second order convex splitting method for the PFC
model based on the Crank-Nicolson formula was also studied in their subsequent
work [16]. They have proved that the Crank-Nicolson convex splitting scheme
was unconditionally energy stable in the sense of the discrete free energy
bounded by the initial one, which indicates the uniform $L^{\infty}$ bound of
the numerical solutions. Then the corresponding error estimate of the proposed
scheme could be easily derived with the help of the uniform bound of the
numerical solutions, seeing [3]. In [11], Guo and Xu presented first order and
second order unconditionally energy stable schemes based on a local
discontinuous Galerkin method with the convex splitting approach. An
unconditionally energy stable BDF2 convex splitting scheme has been developed
for the PFC model by Li, Mei and You in [19]. They rigorously established the
mass conservation, the unconditional energy stability and the convergence of
the proposed numerical scheme. Though all the above convex splitting schemes
are uniquely solvable due to the gradient of a strictly convex function, their
resulted nonlinear systems must be solved at each time step, which is
relatively time-consuming. Thus it is highly desired to develop linear high
order energy stable schemes for the PFC model (1.2). In [31], Yang et al.
proposed two linear second order unconditionally energy stable schemes for the
PFC equation using the invariant energy quadratization (IEQ) approach.
Recently, Li and Shen have adopted the scalar auxiliary variable (SAV)
approach to construct a linear Crank-Nicolson scheme for the PFC model in
[20], where the unconditional energy stability and the corresponding error
estimate of the fully discrete scheme were rigorously proved.

Another important feature of the PFC model is that its evolution process
usually takes rather long time before it reaches the steady state. Moreover,
in this process, the solution of the PFC equation usually undergoes both fast
and slow changes at different time. Therefore, it is reasonable to employ some
time adaptive strategies in the simulation, where unconditionally energy
stable numerical schemes with variable time steps are very needed. Combining
the second order Crank-Nicolson convex splitting scheme, an adaptive time-
stepping strategy based on the variation of the free energy was adopted for
the PFC model in [33], where both the steady-state solution and the dynamical
changes of the solution have been resolved accurately and efficiently. As
known, the BDF2 schemes usually achieve stronger stability than the one based
on the Crank-Nicolson formula, especially for the equations with strong
stiffness [12, 14, 32]. Recently, an adaptive fully implicit BDF2 time-
stepping method was studied for the PFC model in [21]. Under a restriction on
the adjacent time-step ratio $\gamma_{n}\leq 3.561$, the energy stability and
the corresponding $L^{2}$ error estimate of the proposed adaptive scheme were
rigorously derived by using the discrete orthogonal convolution kernels and
convolution inequalities. However, these adaptive time-stepping schemes in
[21, 33] are also nonlinear and the nonlinear iteration is unavoidable at each
time step. Very recently, combining the nonuniform BDF2 formula and the SAV
approach, we successfully derived a linear second order unconditionally energy
stable time stepping scheme for gradient flows with the adjacent time-step
ratio $\gamma_{n}\leq 4.8465$ in [13]. Though it avoids solving the nonlinear
system for the phase function $\phi^{n}$ by explicit treatment for the
nonlinear term, there still exits a nonlinear algebraic equation to be solved
for the auxiliary variable at each time step.

The goal of this paper is to develop and analyze a linear second order
adaptive time-stepping scheme based on the BDF2 formula for the PFC model
(1.2) using the SAV approach proposed in [26, 27]. The linear BDF2 scheme
based on the SAV approach has been proved unconditionally energy stable on the
uniform temporal mesh for the PFC model in [31]. However, it seems to be
rather difficult to prove similar results on nonuniform temporal meshes. In
this paper, we construct a linear unconditionally energy stable BDF2 scheme
for the PFC model on the nonuniform temporal mesh for the first time by using
a first order approximation on the auxiliary variable. The first order
treatment of the auxiliary variable will not affect the second order accuracy
of the unknown atom density field $\phi^{n}$ by setting the positive constant
$C_{0}$ large enough such that $C_{0}\geq 1/\tau$. Moreover, the constructed
adaptive BDF2 scheme avoids solving the nonlinear algebraic equation for the
auxiliary variable as needed in our previous work [13]. More precisely, there
are only two sixth-order equations with constant coefficients to be solved at
each time step, which can be solved efficiently by exiting fast Poisson
solvers. Furthermore, under a mild restriction on the adjacent time-step ratio
$\gamma_{n}\leq 4.8465$, we rigorously establish the unconditional energy
stability and the corresponding error estimate of the proposed fully discrete
BDF2 scheme on the nonuniform temporal mesh.

The rest of the paper is organized as follows. In section 2, we review the SAV
method proposed in [26, 27], and construct the linear BDF2 scheme on the
nonuniform temporal mesh. The unconditional energy stability for the proposed
scheme is rigorously proved in the sense of a modified discrete energy. In
section 3, a rigorous error estimate of the fully discrete adaptive BDF2
scheme is derived for the PFC model, which shows the second order accuracy of
the numerical solution $\phi^{n}$ in time with larger enough constant $C_{0}$
such that $C_{0}\geq 1/\tau.$ Several numerical examples are given in section
4 to validate the theoretical results and the efficiency of the proposed
method. Finally, the paper ends with some concluding remarks in section 5.

##  2\. Scalar auxiliary variable approach

We first introduce some notations which will be used throughout the paper.
Denote $\|\cdot\|_{m}$ the standard $H^{m}$ norm. The norm and inner product
of $L^{2}(\Omega):=H^{0}(\Omega)$ are denoted by $\|\cdot\|$ and
$(\cdot,\cdot)$, respectively. For the sake of simplicity, we consider the
periodic boundary condition for (1.2).

The key idea of the SAV approach [12, 26] is to transform the nonlinear term
to a simple quadratic form by introducing a scalar auxiliary variable $r(t)$,
defined by

$r(t)=\sqrt{E_{1}(\phi)+C_{0}}:=\sqrt{\int_{\Omega}F(\phi)d{\tenbi x}+C_{0}},$
(2.1)

where $F(\phi)=\frac{1}{4}\phi^{4}-\frac{S+\varepsilon}{2}\phi^{2}$ with
$S\geq 0$, and $C_{0}$ is a positive constant such that $E_{1}(\phi)+C_{0}>0.$
In this paper, we further assume that $C_{0}$ is large enough such that

$\frac{C_{0}}{2}\leq E_{1}(\phi)+C_{0}\leq 2C_{0}.$ (2.2)

Then, an equivalent system of the original problem (1.2) with scalar auxiliary
variable is given as

$\begin{array}[]{r@{}l}\begin{cases}\displaystyle\frac{\partial\phi}{\partial
t}=\displaystyle\Delta\mu,\\\\[9.0pt]
\displaystyle\mu\displaystyle=(\Delta+\beta)^{2}\phi+S\phi+\frac{r(t)}{\sqrt{E_{1}(\phi)+C_{0}}}F^{\prime}(\phi),\\\\[15.0pt]
\displaystyle\frac{dr}{dt}\displaystyle=\frac{1}{2\sqrt{E_{1}(\phi)+C_{0}}}\int_{\Omega}F^{\prime}(\phi)\frac{\partial\phi}{\partial
t}d{\tenbi x}.\end{cases}\end{array}$ (2.3)

(2.3) satisfies the following energy dissipation law

$\displaystyle\frac{d\overline{E}(\phi,r)}{dt}=-\|\nabla\mu\|^{2}\leq 0,$
(2.4)

where
$\overline{E}(\phi,r):=\frac{1}{2}\|(\Delta+\beta)\phi\|^{2}+\frac{S}{2}\|\phi\|^{2}+r^{2}$.
Moreover, we have

$\overline{E}(\phi,r)\leq{\color[rgb]{0,0,0}\overline{E}}(\phi_{0}({\tenbi
x}),r(0)),$ (2.5)

which shows the $H^{2}$ bound of the solution of (2.3). From the definition of
the auxiliary variable $r(t)$ in (2.1), it follows that
$\overline{E}(\phi,r)=E(\phi)+C_{0}$, which indicates that the energy
dissipation law (2.4) is the same as (1.3).

We are now ready to construct efficient time stepping schemes for (2.3) to
numerically approximate the solution $\phi$ on the nonuniform temporal mesh.
Let $\tau_{n}:=t_{n}-t_{n-1},n=1,2\dots,M$ be the time step size,
$\gamma_{n+1}={\tau_{n+1}}/{\tau_{n}}$ be the adjacent time–step ratio, and
$\tau=\max\\{\tau_{n},n=1,2,\dots,M\\}$ be the maximum time step size of the
temporal mesh.

###  2.1. A second-order BDF scheme on nonuniform temporal mesh

We firstly recall the first order scheme for solving (2.3) reported in [20]

$\displaystyle\displaystyle{}\frac{\phi^{n+1}-\phi^{n}}{\tau_{n+1}}$
$\displaystyle\displaystyle=\Delta\mu^{n+1},$ (2.6a)
$\displaystyle\displaystyle\mu^{n+1}$
$\displaystyle\displaystyle=(\Delta+\beta)^{2}\phi^{n+1}+S\phi^{n+1}+\frac{r^{n+1}}{\sqrt{E_{1}^{n}+C_{0}}}F^{\prime}(\phi^{n}),$
(2.6b) $\displaystyle\displaystyle\frac{r^{n+1}-r^{n}}{\tau_{n+1}}$
$\displaystyle\displaystyle=\frac{1}{2\sqrt{E_{1}^{n}+C_{0}}}\int_{\Omega}F^{\prime}(\phi^{n})\frac{\phi^{n+1}-\phi^{n}}{\tau_{n+1}}d{\tenbi
x},$ (2.6c)

where $\phi^{n}$ is an approximation to $\phi(t_{n})$ and
$E^{n}_{1}:=E_{1}(\phi^{n})$. As stated in [20], the scheme (2.6) is
unconditionally energy stable in the sense that the following discrete energy
law holds

$\overline{E}(\phi^{n+1},r^{n+1})-\overline{E}(\phi^{n},r^{n})\leq 0.$ (2.7)

To construct a second-order BDF scheme for (2.3), we firstly introduce a
second order approximation $F^{j+\sigma}_{2}\phi$ to $\phi^{\prime}(t)$ at
$t=t_{j+\sigma}$ as follows, seeing also [13, 15]. Here,
$t_{j+\sigma}:=t_{n}+\sigma\tau_{n+1}$ with $\frac{1}{2}\leq\sigma\leq 1.$

$\displaystyle
F^{j+\sigma}_{2}\phi:=\frac{1}{\tau_{j+1}}\Big{[}\frac{1+2\sigma\gamma_{j+1}}{1+\gamma_{j+1}}\phi^{j+1}-(1+(2\sigma-1)\gamma_{j+1})\phi^{j}+\frac{(2\sigma-1)\gamma^{2}_{j+1}}{1+\gamma_{j+1}}\phi^{j-1}\Big{]},$
(2.8)

where $\phi^{j}$ is an approximation to $\phi(t_{j})$. Note that when
$\sigma=1$, it becomes the classical second order BDF with variable time
steps:

$\displaystyle
F^{j+1}_{2}\phi=\frac{1}{\tau_{j+1}}\Big{[}\frac{1+2\gamma_{j+1}}{1+\gamma_{j+1}}\phi^{j+1}-(1+\gamma_{j+1})\phi^{j}+\frac{\gamma^{2}_{j+1}}{1+\gamma_{j+1}}\phi^{j-1}\Big{]},$
(2.9)

and Crank-Nicolson formula

$\displaystyle F^{j+1/2}_{2}\phi=\frac{\phi^{j+1}-\phi^{j}}{\tau_{j+1}}$

for $\sigma=\frac{1}{2}$.

Now, we are ready to construct the second-order BDF scheme on nonuniform
temporal mesh for the PFC system (2.3),

$\displaystyle\displaystyle F^{n+\sigma}_{2}\phi=\Delta\mu^{n+\sigma},$
(2.10a)
$\displaystyle\displaystyle\mu^{n+\sigma}=(\Delta+\beta)^{2}\phi^{n+\sigma}+S\phi^{n+\sigma}+\frac{r^{n+1}}{\sqrt{E_{1}^{n}+C_{0}}}F^{\prime}(\phi^{*,n+\sigma}),$
(2.10b)
$\displaystyle\displaystyle\frac{r^{n+1}-r^{n}}{\tau_{n+1}}=\frac{1}{2\sqrt{E_{1}^{n}+C_{0}}}\int_{\Omega}F^{\prime}(\phi^{*,n+\sigma})\frac{\phi^{n+1}-\phi^{n}}{\tau_{n+1}}d{\tenbi
x},$ (2.10c)

where $\phi^{n+\sigma}:=\sigma\phi^{n+1}+(1-\sigma)\phi^{n}$ and
$\phi^{*,n+\sigma}:=\phi^{n}+\sigma\gamma_{n+1}(\phi^{n}-\phi^{n-1})$. Though
only the first order finite difference approximation is used in (2.10c), we
will see in section 3, the proposed scheme (2.10) is actually of second order
accuracy for $\phi$ with $C_{0}\geq 1/\tau$.

Before rigorously proving the energy stability of the scheme (2.10), we need
the following essential inequality, seeing also [13, 15],

$\begin{array}[]{r@{}l}\displaystyle\delta_{t}\phi^{n+1}\cdot
F^{n+\sigma}_{2}\phi\geq\displaystyle\Big{[}g(\gamma_{n+2})+\frac{1}{2}G(\gamma_{n+1},\gamma_{n+2})\Big{]}\frac{|\delta_{t}\phi^{n+1}|^{2}}{\tau_{n+1}}-g(\gamma_{n+1})\frac{|\delta_{t}\phi^{n}|^{2}}{\tau_{n}},\end{array}$
(2.11)

where $\delta_{t}\phi^{n+1}=\phi^{n+1}-\phi^{n}$. And in (2.11), for any
$0<s,z\leq\gamma_{**}(\sigma)$,

$g(\gamma_{n+1}):=\frac{(2\sigma-1)\gamma_{n+1}^{\frac{3}{2}}}{2(1+\gamma_{n+1})}>0$

and

$G(s,z):=\frac{2+4\sigma
s-(2\sigma-1)s^{\frac{3}{2}}}{1+s}-\frac{(2\sigma-1)z^{\frac{3}{2}}}{1+z}\geq
0.$

Here $\gamma_{**}(\sigma)$ is the positive root of $G(z,z)=0$. Moreover, the
root function $\gamma_{**}(\sigma)$ is decreasing for
$\sigma\in[\frac{1}{2},1]$ with $\gamma_{**}(1)\approx 4.8645$ and
$\gamma_{**}(\sigma)\rightarrow+\infty,$ as $\sigma\rightarrow\frac{1}{2}.$
Thus $\gamma_{**}(\sigma)\geq 4.8645$ for any $\sigma\in[\frac{1}{2},1].$
Furthermore, for $4\leq\gamma_{*}(\sigma)<\gamma_{**}(\sigma)$, it holds

$G(s,z)\geq\min\\{G(0,\gamma_{*}(\sigma)),G(\gamma_{*}(\sigma),\gamma_{*}(\sigma))\\}\geq
G(\gamma_{*}(\sigma),\gamma_{*}(\sigma))>0$

for any $0<s,z\leq\gamma_{*}(\sigma)<\gamma_{**}(\sigma)$.

To prove the following theorem on the unconditional energy stability of the
scheme (2.10), we further introduce the inverse Laplace operator
$\psi=(-\Delta)^{-1}u$ for $u\in L^{2}(\Omega)$ and $\int_{\Omega}ud{\tenbi
x}=0$, as defined in [28]:

$\begin{array}[]{r@{}l}\begin{cases}-\Delta\psi=u,\\\
\displaystyle\int_{\Omega}\psi d{\tenbi x}=0.\end{cases}\end{array}$

######  Theorem 2.1.

For $1/2\leq\sigma\leq 1$ and
$0<\gamma_{n+1}\leq\gamma_{**}(\sigma),n=1,\cdots,M-1$, it holds for the
scheme (2.10) that

$\displaystyle\widetilde{E}^{n+1}-\widetilde{E}^{n}\leq 0,$ (2.12)

where the discrete modified energy $\widetilde{E}^{n}$ is defined by:

$\displaystyle\widetilde{E}^{n}:=\displaystyle\overline{E}(\phi^{n},r^{n})+\frac{(2\sigma-1)\gamma_{n+1}^{\frac{3}{2}}}{2+2\gamma_{n+1}}\frac{\|\nabla^{-1}(\phi^{n}-\phi^{n-1})\|^{2}}{\tau_{n}}.$

###### Proof.

Taking the inner products of (2.10a) and (2.10b) with
$(-\Delta)^{-1}\delta_{t}\phi^{n+1}$, $\delta_{t}\phi^{n+1}$ respectively, and
multiplying (2.10c) with $2r^{n+1}$, we obtain

$\begin{array}[]{r@{}l}&\displaystyle\Big{(}\nabla^{-1}F^{n+\sigma}_{2}\phi,\nabla^{-1}\delta_{t}\phi^{n+1})\Big{)}=-(\mu^{n+\sigma},\delta_{t}\phi^{n+1}),\\\\[8.0pt]
&\displaystyle\Big{(}\mu^{n+\sigma},\delta_{t}\phi^{n+1}\Big{)}=\big{(}(\Delta+\beta)\phi^{n+\sigma},(\Delta+\beta)\delta_{t}\phi^{n+1}\big{)}+S\big{(}\phi^{n+\sigma},\delta_{t}\phi^{n+1}\big{)}\\\
&\hskip
79.6678pt\displaystyle+\frac{r^{n+1}}{\sqrt{E_{1}^{n}+C_{0}}}\Big{(}F^{\prime}(\phi^{*,n+1}),\delta_{t}\phi^{n+1}\Big{)},\\\\[9.0pt]
&\displaystyle\frac{2r^{n+1}(r^{n+1}-r^{n})}{\tau_{n+1}}=\frac{r^{n+1}}{\sqrt{E_{1}^{n}+C_{0}}}\int_{\Omega}F^{\prime}(\phi^{*,n+1})\frac{\delta_{t}\phi^{n+1}}{\tau_{n+1}}d{\tenbi
x}.\end{array}$

From the above equations, it follows that

$\begin{array}[]{r@{}l}&\Big{(}\nabla^{-1}F^{n+\sigma}_{2}\phi,\nabla^{-1}\delta_{t}\phi^{n+1}\Big{)}+\big{(}(\Delta+\beta)\phi^{n+\sigma},(\Delta+\beta)\delta_{t}\phi^{n+1}\big{)}\\\
\\\
&+S\big{(}\phi^{n+\sigma},\delta_{t}\phi^{n+1}\big{)}+2r^{n+1}(r^{n+1}-r^{n})=0.\end{array}$
(2.13)

Combining the inequality (2.11) and

$\begin{array}[]{r@{}l}&2a^{k+\sigma}(a^{k+1}-a^{k})=|a^{k+1}|^{2}-|a^{k}|^{2}+(2\sigma-1)|a^{k}-a^{k-1}|^{2},\\\\[9.0pt]
&2a^{k+1}(a^{k+1}-a^{k})=|a^{k+1}|^{2}-|a^{k}|^{2}+|a^{k+1}-a^{k}|^{2},\end{array}$
(2.14)

we can easily derive the modified discrete energy law (2.12) from
(LABEL:eqt1). Then we complete the proof. ∎

###  2.2. Fully discrete scheme with Fourier spectral method

In this work, we only present the results in the two dimensional space. The
results in the three dimensional space can be obtained in a similar way.

Let $\Omega:=(0,L)^{2}$ and we use the uniform Fourier mode number of $2N+1$
for each direction. The trial and test function space of the Fourier spectral
method is defined by

$X_{N}:=\mbox{span}\\{e^{i2(kx+ly)\pi/L}:-N\leq k,l\leq N\\},$

where $i=\sqrt{-1}$. By taking the Fourier spectral discretization in space,
we can have the following fully discrete scheme of the semi-discrete approach
(2.10): finding $\phi^{n+1}_{N}\in X_{N}$ and $r^{n+1}$ such that

$\displaystyle\displaystyle(F^{n+\sigma}_{2}\phi_{N},p_{N})=$
$\displaystyle\displaystyle-(\nabla\mu^{n+\sigma}_{N},\nabla
p_{N}),~{}~{}~{}\forall p_{N}\in X_{N},$ (2.15a)
$\displaystyle\displaystyle(\mu^{n+\sigma}_{N},q_{N})=$
$\displaystyle\displaystyle\big{(}(\Delta+\beta)\phi^{n+\sigma}_{N},(\Delta+\beta)q_{N})+S(\phi^{n+\sigma}_{N},q_{N})$
$\displaystyle\displaystyle+\frac{r^{n+1}}{\sqrt{E_{1,N}^{n}+C_{0}}}(F^{\prime}(\phi^{*,n+\sigma}_{N}),q_{N}),~{}~{}~{}\forall
q_{N}\in X_{N}$ (2.15b)
$\displaystyle\displaystyle\frac{r^{n+1}-r^{n}}{\tau_{n+1}}=$
$\displaystyle\displaystyle\frac{1}{2\sqrt{E_{1,N}^{n}+C_{0}}}\int_{\Omega}F^{\prime}(\phi^{*,n+\sigma}_{N})\frac{\phi^{n+1}_{N}-\phi^{n}_{N}}{\tau_{n+1}}d{\tenbi
x},$ (2.15c)

where $E^{n}_{1,N}:=E_{1}(\phi^{n}_{N})$,
$\phi^{n+\sigma}_{N}:=\sigma\phi^{n+1}_{N}+(1-\sigma)\phi^{n}_{N}$ and
$\phi^{*n+\sigma}_{N}:=\phi^{n}_{N}+\sigma\gamma_{n+1}(\phi^{n}_{N}-\phi^{n-1}_{N})$.

Since the proof of the unconditional energy stability for (2.15) is
essentially similar to that of the semi-discrete scheme (2.10), we only state
the following theorem and omit the proof.

######  Theorem 2.2.

For the fully discrete scheme (2.15), if $1/2\leq\sigma\leq 1$ and
$0<\gamma_{n+1}\leq\gamma_{**}(\sigma),n=1,2,\cdots,M-1$, we have

$\displaystyle\widetilde{E}^{n+1}_{N}-\widetilde{E}^{n}_{N}\leq 0,$ (2.16)

where $\widetilde{E}^{n}_{N}$ is defined by:

$\displaystyle\widetilde{E}^{n}_{N}:=\overline{E}(\phi^{n}_{N},{\color[rgb]{0,0,0}r^{n}})+\displaystyle\frac{(2\sigma-1)\gamma_{n+1}^{\frac{3}{2}}}{2+2\gamma_{n+1}}\frac{\|\nabla^{-1}(\phi^{n}_{N}-\phi^{n-1}_{N})\|^{2}}{\tau_{n}}.$

##  3\. Error analysis

In this section, we establish a rigorous error estimate for the fully discrete
scheme (2.15) with a mild smoothness requirement on the exact solution. For
simplicity, we only focus our attention on the investigation of the numerical
solutions for the special case with $\sigma=1$ in what follows, and the
obtained results are also right for the general parameter $\sigma\in[1/2,1].$
We denote $\gamma_{*}:=\gamma_{*}(1)$. From (2.16), we can deduce that for any
$\phi_{0}({\tenbi x})\in H^{2}(\Omega)$ and $0<\gamma_{n}\leq\gamma_{*}$,

$\|(\Delta+\beta)\phi^{n}_{N}\|+\|\phi^{n}_{N}\|+|r^{n}|\leq\mathcal{M},~{}~{}~{}n=1,2,\cdots,M$
(3.1)

for the fully discrete second order scheme (2.15). Here $\mathcal{M}$ is a
positive constant only depending on $\Omega$ and the initial condition
$\phi_{0}({\tenbi x})$. This indicates the $H^{2}$ bound of the numerical
solution and the uniform bound of the auxiliary variable $r^{n}$ for all
$n=1,2,\cdots,M.$

We recall the $L^{2}$ orthogonal projection $\Pi_{N}$:
$L^{2}(\Omega)\rightarrow X_{N}$, defined by

$(\Pi_{N}u-u,v)=0,~{}~{}\forall v\in X_{N}.$ (3.2)

It is obvious that $\Pi_{N}u$ is the truncated Fourier series, namely

$\Pi_{N}u=\sum_{k,l=-N}^{N}\widehat{u}_{k,l}e^{i2(kx+ly)\pi/L},$

with

$\widehat{u}_{k,l}=(u,e^{i2(kx+ly)\pi/L})=\frac{1}{|\Omega|}\int_{\Omega}u(x,y){\color[rgb]{0,0,0}e^{-i2(kx+ly)\pi/L}}dxdy.$

For this $L^{2}$ orthogonal projection $\Pi_{N}$, we have the following
optimal estimate, seeing also [24],

$\|\Pi_{N}u-u\|_{\mu}\leq CN^{\mu-m}\|u\|_{m},$ (3.3)

for any $u\in H^{m}_{per}(\Omega)$ and $0\leq\mu\leq m$. Here
$H^{m}_{per}(\Omega)$ is the subspace of $H^{m}(\Omega)$, which consists of
functions with derivatives of order up to $m-1$ being $L$-periodic. We denote
$\overline{e}^{n}_{u}:=u^{n}_{N}-\Pi_{N}u(t_{n})$ with
$\overline{e}^{0}_{u}=0$, $\widehat{e}^{n}_{u}:=\Pi_{N}u(t_{n})-u(t_{n})$,
$e^{n}_{r}:=r^{n}-r(t_{n})$ and
$\phi(t_{*,n+1}):=\phi(t_{n})+\gamma_{n+1}(\phi(t_{n})-\phi(t_{n-1})).$ Then
we have

$e^{n}_{\phi}:=\phi^{n}_{N}-\phi(t_{n})=\overline{e}^{n}_{\phi}+\widehat{e}^{n}_{\phi},~{}e^{n}_{\mu}:=\mu^{n}_{N}-\mu(t_{n})=\overline{e}^{n}_{\mu}+\widehat{e}^{n}_{\mu}.$

For simplicity, we also denote
$\delta_{t}\overline{e}^{n}_{u}:=\overline{e}^{n}_{u}-\overline{e}^{n-1}_{u}$.
The error estimates of $\phi^{n}_{N}$ and $r^{n}$ are derived in the following
theorem.

######  Theorem 3.1.

Assume that $\phi^{0}:=\phi({\tenbi x},0)\in
H^{2}(\Omega),0<\gamma_{n}\leq\gamma_{*},$ $\tau_{1}\leq C\tau^{\frac{4}{3}}$,
and $F(\cdot)\in C^{3}(\mathbb{R})$. Provided that

$\begin{array}[]{r@{}l}&\displaystyle\phi\in
L^{\infty}(0,T;H^{m}_{per}),\phi_{t}\in L^{\infty}(0,T;H^{m}_{per})\cap
L^{2}(0,T;H^{1}),\\\\[7.0pt] &\displaystyle{\phi_{tt}\in
L^{2}(0,T;H^{1}),\phi_{ttt}\in L^{2}(0,T;H^{-1})}\end{array}$ (3.4)

for $m\geq 2$, it holds

$\begin{array}[]{r@{}l}&\displaystyle\frac{\gamma^{\frac{3}{2}}_{n+2}}{2(1+\gamma_{n+2})}\frac{\|\nabla^{-1}(e^{n+1}_{\phi}-e^{n}_{\phi})\|^{2}}{\tau_{n+1}}+\frac{1}{2}\|(\Delta+\beta)e^{n+1}_{\phi}\|^{2}+\frac{S}{2}\|e^{n+1}_{\phi}\|^{2}+|e^{n+1}_{r}|^{2}\\\\[9.0pt]
\leq&\displaystyle
C\exp(T)\Big{[}\frac{\tau^{2}}{C_{0}}\int_{0}^{T}\big{(}\|\phi_{t}(s)\|^{2}_{H^{1}}+\|\phi_{tt}(s)\|^{2}\big{)}ds+\tau^{4}\int_{0}^{T}\big{(}\|\phi_{tt}(s)\|^{2}_{H^{1}}\\\\[9.0pt]
&\displaystyle+\|\phi_{ttt}(s)\|^{2}_{H^{-1}}\big{)}ds+T\Big{(}N^{-2m}\|\phi_{t}\|^{2}_{L^{\infty}(0,T;H^{m})}+N^{4-2m}\|\phi\|_{L^{\infty}(0,T;H^{m})}\Big{)}\Big{]}.\end{array}$
(3.5)

Moreover, we have

$\begin{array}[]{r@{}l}&\displaystyle\frac{\gamma^{\frac{3}{2}}_{n+2}}{2(1+\gamma_{n+2})}\frac{\|\nabla^{-1}(e^{n+1}_{\phi}-e^{n}_{\phi})\|^{2}}{\tau_{n+1}}+\frac{1}{2}\|(\Delta+\beta)e^{n+1}_{\phi}\|^{2}+\frac{S}{2}\|e^{n+1}_{\phi}\|^{2}\\\\[11.0pt]
\leq&\displaystyle
C\exp(T)\Big{[}\frac{\tau^{2}}{C^{2}_{0}}\int_{0}^{T}\big{(}\|\phi_{t}(s)\|^{2}_{H^{1}}+\|\phi_{tt}(s)\|^{2}\big{)}ds+\tau^{4}\int_{0}^{T}\big{(}\|\phi_{tt}(s)\|^{2}_{H^{1}}\\\\[11.0pt]
&\displaystyle+\|\phi_{ttt}(s)\|^{2}_{H^{-1}}\big{)}ds+T\Big{(}N^{-2m}\|\phi_{t}\|^{2}_{L^{\infty}(0,T;H^{m})}+N^{4-2m}\|\phi\|_{L^{\infty}(0,T;H^{m})}\Big{)}\Big{]},\end{array}$
(3.6)

and if we set the positive constant $C_{0}$ large enough such that
$C_{0}\geq\frac{1}{\tau}$, it holds

$\begin{array}[]{r@{}l}&\displaystyle\frac{\gamma^{\frac{3}{2}}_{n+2}}{2(1+\gamma_{n+2})}\frac{\|\nabla^{-1}(e^{n+1}_{\phi}-e^{n}_{\phi})\|^{2}}{\tau_{n+1}}+\frac{1}{2}\|(\Delta+\beta)e^{n+1}_{\phi}\|^{2}+\frac{S}{2}\|e^{n+1}_{\phi}\|^{2}\\\\[11.0pt]
\leq&\displaystyle
C\exp(T)\Big{[}\tau^{4}\int_{0}^{T}\big{(}\|\phi_{t}(s)\|^{2}_{H^{1}}+\|\phi_{tt}(s)\|^{2}_{H^{1}}+\|\phi_{ttt}(s)\|^{2}_{H^{-1}}\big{)}ds\\\\[11.0pt]
&\displaystyle+T\Big{(}N^{-2m}\|\phi_{t}\|^{2}_{L^{\infty}(0,T;H^{m})}+N^{4-2m}\|\phi\|_{L^{\infty}(0,T;H^{m})}\Big{)}\Big{]}.\end{array}$
(3.7)

###### Proof.

From (2.5), (3.1) and the Sobolev embedding theorem, $H^{2}\subset
L^{\infty}$, it follows that $\|\phi(t)\|_{L^{\infty}},$
$\|\phi^{n}\|_{L^{\infty}}\leq C,$ where the positive constant $C$ only
depends on $\phi^{0},\Omega$ and $T$. Combining the definition of $F(\cdot)$,
we have

$\|F(\cdot)\|_{L^{\infty}},\|F^{\prime}(\cdot)\|_{L^{\infty}},\|F^{\prime\prime}(\cdot)\|_{L^{\infty}},\|F^{\prime\prime\prime}(\cdot)\|_{L^{\infty}}\leq
C,$ (3.8)

for $\phi$ and $\phi^{n},n=0,1,2,\cdots,M$. In what follows, we use
$\mathcal{A}\lesssim\mathcal{B}$ to represent that $\mathcal{A}\leq
C\mathcal{B}$, where the general positive constant C is independent of
$C_{0}$. From the definition of $r(t)$ in (2.1), we can derive

$r_{tt}=-\frac{\big{(}\int_{\Omega}F^{\prime}(\phi)\phi_{t}d{\tenbi
x}\big{)}^{2}}{4\sqrt{(E_{1}(\phi)+C_{0})^{3}}}+\frac{\int_{\Omega}\big{[}F^{\prime\prime}(\phi)\phi^{2}_{t}+F^{\prime}(\phi)\phi_{tt}\big{]}d{\tenbi
x}}{2\sqrt{E_{1}(\phi)+C_{0}}}.$

Moreover, together with (2.2), (LABEL:sm_cd) and (3.8), and using Hölder
inequality, we can deduce that

$\int_{t_{n}}^{t_{n+1}}|r_{tt}|^{2}dt\lesssim\frac{1}{C_{0}}\int_{t_{n}}^{t_{n+1}}(\|\phi_{t}\|^{2}_{L^{4}}+\|\phi_{tt}\|^{2})dt\lesssim\frac{1}{C_{0}}\int_{t_{n}}^{t_{n+1}}(\|\phi_{t}\|^{2}_{H^{1}}+\|\phi_{tt}\|^{2})dt,$
(3.9)

where we have used the Sobolev embedding theorem, $H^{1}\subset L^{4}.$ From
(2.3), (2.10) and the definition of the $L^{2}$ orthogonal projection
$\Pi_{n}\phi(t_{n})$ in (3.2), we derive the following error equation for
$n\geq 1$:

$\displaystyle\displaystyle(F^{n+1}_{2}\overline{e}_{\phi},p_{N})+(\nabla\overline{e}^{n+1}_{\mu},\nabla
p_{N})=(T_{1}^{n},p_{N}),$ (3.10a)
$\displaystyle\displaystyle(\overline{e}^{n+1}_{\mu},q_{N})=\big{(}(\Delta+\beta)\overline{e}^{n+1}_{\phi},(\Delta+\beta)q_{N}\big{)}+S(\overline{e}^{n+1}_{\phi},q_{N})$
$\displaystyle\hskip
54.06006pt+\frac{e^{n+1}_{r}}{\sqrt{E^{n}_{1,N}+C_{0}}}\big{(}F^{\prime}(\phi^{*,n+1}_{N}),q_{N}\big{)}+\displaystyle(J^{n}_{1}+T^{n}_{2},q_{N}),$
(3.10b) $\displaystyle\displaystyle
e^{n+1}_{r}-e^{n}_{r}=\frac{1}{2\sqrt{E_{1}^{n}+C_{0}}}\int_{\Omega}F^{\prime}(\phi^{*,n+1}_{N})(\overline{e}^{n+1}_{\phi}-\overline{e}^{n}_{\phi})d{\tenbi
x}$ $\displaystyle\hskip
54.06006pt+\int_{\Omega}J^{n}_{2}\cdot\Big{(}\Pi_{N}\big{(}\phi(t_{n+1})-\phi(t_{n})\big{)}\Big{)}d{\tenbi
x}+J^{n}_{3}-v^{n}_{1}+v^{n}_{2},$ (3.10c)

where

$\begin{array}[]{r@{}l}J^{n}_{1}:=&\displaystyle
r(t_{n+1})\Big{[}\frac{F^{\prime}(\phi^{*,n+1}_{N})}{\sqrt{E^{n}_{1,N}+C_{0}}}-\frac{F^{\prime}(\phi(t_{*,n+1}))}{\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}}\Big{]},\\\\[11.0pt]
J^{n}_{2}:=&\displaystyle\frac{1}{2}\Big{[}\frac{F^{\prime}(\phi^{*,n+1}_{N})}{\sqrt{E^{n}_{1,N}+C_{0}}}-\frac{F^{\prime}(\phi(t_{n+1}))}{\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}}\Big{]},\\\
J^{n}_{3}:=&\displaystyle\int_{\Omega}\frac{F^{\prime}(\phi(t_{n+1}))}{2\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}}(\widehat{e}^{n+1}_{\phi}-\widehat{e}^{n}_{\phi})d{\tenbi
x},\end{array}$

and the truncation errors are defined by

$\displaystyle T^{n}_{1}$
$\displaystyle\displaystyle=\phi_{t}(t_{n+1})-\partial_{t}(\Pi_{2,n}\phi(t))|_{t=t_{n+1}},~{}T^{n}_{2}\displaystyle=F^{\prime}(\phi(t_{*,n+1}))-F^{\prime}(\phi(t_{n+1}))$
(3.11a) $\displaystyle v^{n}_{1}$
$\displaystyle\displaystyle=r(t_{n+1})-r(t_{n})-\tau_{n+1}r_{t}(t_{n+1})=\int_{t_{n}}^{t_{n+1}}(t_{n}-s)r_{tt}(s)ds,$
(3.11b) $\displaystyle v^{n}_{2}$
$\displaystyle\displaystyle=\Big{(}\frac{F^{\prime}(\phi(t_{n+1}))}{2\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}},\int_{t_{n}}^{t_{n+1}}(t_{n}-s)\phi_{tt}(s)ds\Big{)}.$
(3.11c)

Here, $\Pi_{2,n}\phi(t)$ denotes the quadratic interpolation polynomial
associated with $(t_{n-1},\phi(t_{n-1})),$ $(t_{n},\phi(t_{n}))$ and
$(t_{n+1},\phi(t_{n+1}))$. Taking
$p_{N}=(-\Delta)^{-1}(\delta_{t}\overline{e}^{n+1}_{\phi})$ and
$q_{N}=\delta_{t}\overline{e}^{n+1}_{\phi}$ into (3.10a) and (3.10b),
respectively, and multiplying (3.10c) with $2e^{n+1}_{r}$, we sum them up to
derive

$\begin{array}[]{r@{}l}&\displaystyle\big{(}\nabla^{-1}F^{n+1}_{2}\overline{e}_{\phi},\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\big{)}+\big{(}(\Delta+\beta)\overline{e}^{n+1}_{\phi},(\Delta+\beta)\delta_{t}\overline{e}^{n+1}_{\phi})+S(\overline{e}^{n+1}_{\phi},\delta_{t}\overline{e}^{n+1}_{\phi})\\\\[8.0pt]
&+2e^{n+1}_{r}(e^{n+1}_{r}-e^{n}_{r})\\\\[8.0pt]
=&\displaystyle-(J^{n}_{1}+T^{n}_{2},\delta_{t}\overline{e}^{n+1}_{\phi})+\big{(}\nabla^{-1}T^{n}_{1},\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\big{)}\\\
&\displaystyle+2e^{n+1}_{r}\int_{\Omega}J^{n}_{2}\cdot\Pi_{N}(\phi(t_{n+1})-\phi(t_{n}))d{\tenbi
x}+2e^{n+1}_{r}J^{n}_{3}+2e^{n+1}_{r}(-v^{n}_{1}+v^{n}_{2}).\end{array}$

Using (2.11), the identities (LABEL:Id), and Young’s inequality, we derive

$\begin{array}[]{r@{}l}&\displaystyle\frac{1}{2}\Big{[}\Big{(}\frac{\gamma^{\frac{3}{2}}_{n+2}}{1+\gamma_{n+2}}+G(\gamma_{*},\gamma_{*})\Big{)}\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\|^{2}}{\tau_{n+1}}-\frac{\gamma^{\frac{3}{2}}_{n+1}}{1+\gamma_{n+1}}\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n}_{\phi}\|^{2}}{\tau_{n}}\Big{]}\\\\[9.0pt]
&\displaystyle+\frac{1}{2}\big{[}\|(\Delta+\beta)\overline{e}^{n+1}_{\phi}\|^{2}-\|(\Delta+\beta)\overline{e}^{n}_{\phi}\|^{2}\big{]}+\frac{S}{2}\big{[}\|\overline{e}^{n+1}_{\phi}\|^{2}-\|\overline{e}^{n}_{\phi}\|^{2}\big{]}+|e^{n+1}_{r}|^{2}-|e^{n}_{r}|^{2}\\\\[11.0pt]
\leq&\displaystyle C(\gamma_{*})\tau_{n+1}\big{[}\|\nabla
J^{n}_{1}\|^{2}+\|\nabla^{-1}T^{n}_{1}\|^{2}+\|\nabla
T^{n}_{2}\|^{2}\big{]}+G(\gamma_{*},\gamma_{*})\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\|^{2}}{2\tau_{n+1}}\\\\[11.0pt]
&\displaystyle+C\tau_{n+1}\|\phi_{t}\|_{L^{\infty}(0,T;L^{2})}\big{[}|e^{n+1}_{r}|^{2}+\|J^{n}_{2}\|^{2}\big{]}+\frac{C\tau_{n+1}}{\sqrt{C_{0}}}\big{[}|e^{n+1}_{r}|^{2}\\\
&\displaystyle+N^{-2m}\|\phi_{t}\|^{2}_{L^{\infty}(0,T;H^{m})}\big{]}+\tau_{n+1}|e^{n+1}_{r}|^{2}+\frac{2}{\tau_{n+1}}(|v^{n}_{1}|^{2}+|v^{n}_{2}|^{2}),\end{array}$

where we have used

$\begin{array}[]{r@{}l}&\displaystyle 2e^{n+1}_{r}\cdot J^{n}_{3}\\\\[7.0pt]
\leq&\displaystyle\frac{\|F^{\prime}(\phi(t_{n+1}))\|_{L^{\infty}}}{\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}}\int_{\Omega}e^{n+1}_{r}(\widehat{e}^{n+1}_{\phi}-\widehat{e}^{n}_{\phi})d{\tenbi
x}\\\\[9.0pt]
\leq&\displaystyle\frac{C}{\sqrt{C_{0}}}\Big{[}\tau_{n+1}|e^{n+1}_{r}|^{2}+\frac{1}{\tau_{n+1}}\int_{\Omega}|\widehat{e}^{n+1}_{\phi}-\widehat{e}^{n}_{\phi}|^{2}d{\tenbi
x}\Big{]}\\\\[9.0pt]
=&\displaystyle\frac{C}{\sqrt{C_{0}}}\Big{[}\tau_{n+1}|e^{n+1}_{r}|^{2}+\frac{1}{\tau_{n+1}}\int_{\Omega}\big{|}\Pi_{N}\big{(}\phi(t_{n+1})-\phi(t_{n})\big{)}-\big{(}\phi(t_{n+1}-\phi(t_{n}))\big{)}\big{|}^{2}d{\tenbi
x}\Big{]}\\\\[9.0pt]
\leq&\displaystyle\frac{C}{\sqrt{C_{0}}}\Big{[}\tau_{n+1}|e^{n+1}_{r}|^{2}+\tau_{n+1}N^{-2m}\|\phi_{t}\|^{2}_{L^{\infty}(0,T;H^{m})}\Big{]}.\end{array}$

Combining with (LABEL:sm_cd), we have following error inequality

$\begin{array}[]{r@{}l}&\displaystyle\frac{1}{2}\Big{[}\frac{\gamma^{\frac{3}{2}}_{n+2}}{1+\gamma_{n+2}}\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\|^{2}}{\tau_{n+1}}-\frac{\gamma^{\frac{3}{2}}_{n+1}}{1+\gamma_{n+1}}\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n}_{\phi}\|^{2}}{\tau_{n}}\Big{]}+\frac{1}{2}\big{[}\|(\Delta+\beta)\overline{e}^{n+1}_{\phi}\|^{2}\\\\[9.0pt]
&\displaystyle-\|(\Delta+\beta)\overline{e}^{n}_{\phi}\|^{2}\big{]}+\frac{S}{2}\big{[}\|\overline{e}^{n+1}_{\phi}\|^{2}-\|\overline{e}^{n}_{\phi}\|^{2}\big{]}+|e^{n+1}_{r}|^{2}-|e^{n}_{r}|^{2}\\\\[11.0pt]
\lesssim&\displaystyle\tau_{n+1}\Big{[}\|\nabla
J^{n}_{1}\|^{2}+\|J^{n}_{2}\|^{2}+|e^{n+1}_{r}|^{2}+\|\nabla^{-1}T^{n}_{1}\|^{2}+\|\nabla
T^{n}_{2}\|^{2}\\\
&\displaystyle+N^{-2m}\|\phi_{t}\|^{2}_{L^{\infty}(0,T;H^{m})}\Big{]}+\frac{2(|v^{n}_{1}|^{2}+|v^{n}_{2}|^{2})}{\tau_{n+1}}.\end{array}$
(3.12)

Moreover, we have

$\displaystyle\nabla J^{n}_{1}=-r(t_{n+1})\Big{[}\frac{\nabla
I^{n}_{1}}{\sqrt{E^{n}_{1,N}+C_{0}}}+\nabla F^{\prime}(\phi(t_{*,n+1}))\cdot
I^{n}_{2}\Big{]},$ (3.13)

where

$I^{n}_{1}:=F^{\prime}(\phi^{*,n+1}_{N})-F^{\prime}(\phi(t_{*,n+1})),\ \
I^{n}_{2}:=\frac{1}{\sqrt{E^{n}_{1,N}+C_{0}}}-\frac{1}{\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}}.$
(3.14)

From (3.8) and $\phi(t_{*,n+1})\in H^{2}$, it follows that

$\begin{array}[]{r@{}l}\|\nabla
F^{\prime}(\phi(t_{*,n+1}))\|=&\displaystyle\|F^{{}^{\prime\prime}}(\phi(t_{*,n+1}))\nabla\phi(t_{*,n+1})\|\\\\[7.0pt]
\leq&\displaystyle\|F^{{}^{\prime\prime}}(\phi(t_{*,n+1}))\|_{L^{\infty}}\|\nabla\phi(t_{*,n+1})\|\leq
C.\end{array}$

Therefore, together with (2.1) and (2.2), we deduce that

$\|\nabla J^{n}_{1}\|^{2}\lesssim[\|\nabla
I^{n}_{1}\|^{2}+C_{0}|I^{n}_{2}|^{2}].$ (3.15)

For the term $J^{n}_{2}$, using triangle inequality and (2.2), we have

$\begin{array}[]{r@{}l}\|J^{n}_{2}\|^{2}=&\displaystyle\frac{1}{2}\Big{\|}\frac{I^{n}_{1}}{\sqrt{E^{n}_{1,N}+C_{0}}}+\frac{T^{n}_{2}}{\sqrt{E^{n}_{1,N}+C_{0}}}+F^{\prime}(\phi(t_{n+1}))\cdot
I^{n}_{2}\Big{\|}^{2}\\\\[11.0pt]
\lesssim&\displaystyle\frac{\|I^{n}_{1}\|^{2}}{C_{0}}+\frac{\|T^{n}_{2}\|^{2}}{C_{0}}+|I^{n}_{2}|^{2}.\\\\[9.0pt]
\end{array}$ (3.16)

Then, we substitute the estimates (3.15) and (3.16) into the error inequality
(LABEL:equ2) to get

$\begin{array}[]{r@{}l}&\displaystyle\frac{\gamma^{\frac{3}{2}}_{n+2}}{1+\gamma_{n+2}}\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\|^{2}}{2\tau_{n+1}}-\frac{\gamma^{\frac{3}{2}}_{n+1}}{1+\gamma_{n+1}}\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n}_{\phi}\|^{2}}{2\tau_{n}}+\frac{1}{2}\big{[}\|(\Delta+\beta)\overline{e}^{n+1}_{\phi}\|^{2}\\\\[9.0pt]
&\displaystyle-\|(\Delta+\beta)\overline{e}^{n}_{\phi}\|^{2}\big{]}+\frac{S}{2}\big{[}\|\overline{e}^{n+1}_{\phi}\|^{2}-\|\overline{e}^{n}_{\phi}\|^{2}\big{]}+|e^{n+1}_{r}|^{2}-|e^{n}_{r}|^{2}\\\\[9.0pt]
\lesssim&\displaystyle\tau_{n+1}\Big{[}|e^{n+1}_{r}|^{2}+\frac{1}{C_{0}}\|I^{n}_{1}\|^{2}+\|\nabla
I^{n}_{1}\|^{2}+C_{0}|I^{n}_{2}|^{2}+\|\nabla^{-1}T^{n}_{1}\|^{2}+\frac{1}{C_{0}}\|T^{n}_{2}\|^{2}\\\
&\displaystyle+\|\nabla
T^{n}_{2}\|^{2}+N^{-2m}\|\phi_{t}\|^{2}_{L^{\infty}(0,T;H^{m})}\Big{]}+\frac{2}{\tau_{n+1}}(|v^{n}_{1}|^{2}+|v^{n}_{2}|^{2}).\end{array}$
(3.17)

Moreover, from (3.14), (3.8), (3.3) and (2.2),

$\begin{array}[]{r@{}l}\|I^{n}_{1}\|^{2}=&\displaystyle\|F^{\prime}(\phi^{*,n+1}_{N})-F^{\prime}(\phi(t_{*,n+1}))\|^{2}\lesssim\|e^{n}_{\phi}\|^{2}+\|e^{n-1}_{\phi}\|^{2},\\\\[8.0pt]
\lesssim&\|\overline{e}^{n}_{\phi}+\widehat{e}^{n}_{\phi}\|^{2}+\|\overline{e}^{n-1}_{\phi}+\widehat{e}^{n-1}_{\phi}\|^{2}\\\\[8.0pt]
\lesssim&\displaystyle\|\overline{e}^{n}_{\phi}\|^{2}+\|\overline{e}^{n-1}_{\phi}\|^{2}+N^{-2m}\|\phi\|_{L^{\infty}(0,T;H^{m})},\\\\[8.0pt]
|I^{n}_{2}|^{2}=&\displaystyle\Big{|}\frac{E_{1}(\phi(t_{n+1}))-E^{n}_{1,N}}{\mathcal{G}(C_{0})}\Big{|}^{2}\lesssim\displaystyle\frac{1}{C^{3}_{0}}\big{|}E_{1}(\phi(t_{n+1}))-E^{n}_{1,N}\big{|}^{2}\\\\[8.0pt]
\lesssim&\displaystyle\frac{1}{C^{3}_{0}}\Big{(}\big{|}E_{1}(\phi(t_{n+1}))-E_{1}(\phi(t_{n}))\big{|}^{2}+\big{|}E_{1}(\phi(t_{n}))-E^{n}_{1,N}\big{|}^{2}\Big{)}\\\\[8.0pt]
\lesssim&\displaystyle\frac{1}{C^{3}_{0}}\Big{(}\tau_{n+1}\int_{t_{n}}^{t_{n+1}}\|\phi_{t}(s)\|^{2}ds+\|e^{n}_{\phi}\|^{2}\Big{)}\\\\[8.0pt]
\lesssim&\displaystyle\frac{1}{C^{3}_{0}}\Big{(}\tau_{n+1}\int_{t_{n}}^{t_{n+1}}\|\phi_{t}(s)\|^{2}ds+\|\overline{e}^{n}_{\phi}\|^{2}+N^{-2m}\|\phi\|_{L^{\infty}(0,T;H^{m})}\Big{)},\end{array}$
(3.18)

where the $\mathcal{G}(C_{0})$ is defined by

$\mathcal{G}(C_{0}):=\sqrt{E^{n}_{1,N}+C_{0}}\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}\Big{(}\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}+\sqrt{E^{n}_{1,N}+C_{0}}\Big{)}.$

From the definition of $I^{n}_{1}$ in (3.14), we have

$\|\nabla
I^{n}_{1}\|^{2}=\displaystyle\|I^{n}_{1,1}+I^{n}_{1,2}\|^{2}\leq\|I^{n}_{1,1}\|^{2}+\|I^{n}_{1,2}\|^{2},$
(3.19)

where $I^{n}_{1,1}$ and $I^{n}_{1,2}$ are defined by

$\begin{array}[]{r@{}l}&\displaystyle
I^{n}_{1,1}:=F^{\prime\prime}(\phi^{*,n+1}_{N})\nabla\big{(}\phi^{*,n+1}_{N}-\phi(t_{*,n+1})\big{)},\\\\[7.0pt]
&\displaystyle
I^{n}_{1,2}:=\big{(}F^{\prime\prime}(\phi^{*,n+1}_{N})-F^{\prime\prime}(\phi(t_{*,n+1}))\big{)}\nabla\phi(t_{*,n+1}).\end{array}$

Furthermore, using (3.8), the Hölder’s inequality and the Sobolev embedding
theorem, $H^{1}\subset L^{6}$, we have

$\begin{array}[]{r@{}l}\|I^{n}_{1,1}\|^{2}\lesssim&\displaystyle\|\nabla
e^{n}_{\phi}\|^{2}+\|\nabla e^{n-1}_{\phi}\|^{2}\\\\[7.0pt]
\|I^{n}_{1,2}\|^{2}\lesssim&\|(\phi^{*,n+1}_{N}-\phi(t_{*,n+1}))\nabla\phi(t_{*,n+1})\|^{2}\\\\[7.0pt]
\lesssim&\|\phi^{*,n+1}_{N}-\phi(t_{*,n+1})\|^{2}_{L^{6}}\|\nabla\phi(t_{*,n+1})\|^{2}_{L^{3}}\\\\[7.0pt]
\lesssim&\|\phi^{*,n+1}_{N}-\phi(t_{*,n+1})\|^{2}_{H^{1}}\|\phi(t_{*,n+1})\|^{2}_{H^{2}}\\\\[7.0pt]
\lesssim&\|e^{n}_{\phi}\|^{2}_{H^{1}}+\|e^{n-1}_{\phi}\|^{2}_{H^{1}}.\end{array}$
(3.20)

Combining (3.19) and (3.20) gives

$\begin{array}[]{r@{}l}\|\nabla I^{n}_{1}\|^{2}\lesssim&\displaystyle\|\nabla
e^{n}_{\phi}\|^{2}+\|\nabla
e^{n-1}_{\phi}\|^{2}+\|e^{n}_{\phi}\|^{2}+\|e^{n-1}_{\phi}\|^{2}\\\\[7.0pt]
\lesssim&\displaystyle\|\Delta\overline{e}^{n}_{\phi}\|^{2}+\|\Delta\overline{e}^{n-1}\|^{2}+\|\overline{e}^{n}_{\phi}\|^{2}+\|\overline{e}^{n-1}_{\phi}\|^{2}+N^{4-2m}\|\phi\|^{2}_{L^{\infty}(0,T;H^{m})},\\\\[7.0pt]
\lesssim&\displaystyle\|(\Delta+\beta)\overline{e}^{n}_{\phi}\|^{2}+\|(\Delta+\beta)\overline{e}^{n-1}_{\phi}\|^{2}+\|\overline{e}^{n}_{\phi}\|^{2}+\|\overline{e}^{n-1}_{\phi}\|^{2}+N^{4-2m}\|\phi\|^{2}_{L^{\infty}(0,T;H^{m})}\end{array}$
(3.21)

in which we have used

$\begin{array}[]{r@{}l}&\displaystyle\|\nabla e^{n+1}_{\phi}\|^{2}=(-\Delta
e^{n+1}_{\phi},e^{n+1}_{\phi})\leq\frac{\|\Delta
e^{n+1}_{\phi}\|^{2}+\|e^{n+1}_{\phi}\|^{2}}{2},\\\\[7.0pt]
&\displaystyle\|\Delta
e^{n+1}_{\phi}\|^{2}=\|\Delta\overline{e}^{n+1}_{\phi}+\Delta\widehat{e}^{n+1}_{\phi}\|^{2}\leq\|\Delta\overline{e}^{n+1}_{\phi}\|^{2}+CN^{4-2m}\|\phi\|^{2}_{L^{\infty}(0,T;H^{m})}.\end{array}$

For the truncation errors of $T^{n}_{1},T^{n}_{2},v^{n}_{1}$, and $v^{n}_{2}$,
the following estimates hold, seeing also [2, 25],

$\begin{array}[]{r@{}l}\|\nabla^{-1}T^{n}_{1}\|^{2}\lesssim&\displaystyle(\tau_{n}+\tau_{n+1})^{3}\int_{t_{n-1}}^{t_{n+1}}\|\phi_{ttt}(s)\|^{2}_{H^{-1}}ds,\\\
\|T^{n}_{2}\|^{2}\lesssim&\displaystyle\|\phi(t_{n+1})-\phi(t_{*,n+1})\|^{2}\lesssim(\tau_{n}+\tau_{n+1})^{3}\int_{t_{n-1}}^{t_{n+1}}\|\phi_{tt}(s)\|^{2}ds,\\\
\|\nabla
T^{n}_{2}\|^{2}\lesssim&\displaystyle\|\phi(t_{n+1})-\phi(t_{*,n+1})\|^{2}+\|\nabla(\phi(t_{n+1})-\phi(t_{*,n+1}))\|^{2}\\\\[9.0pt]
\lesssim&\displaystyle(\tau_{n}+\tau_{n+1})^{3}\int_{t_{n-1}}^{t_{n+1}}\|\phi_{tt}(s)\|^{2}_{H^{1}}ds,\\\
|v^{n}_{1}|^{2}\leq&\displaystyle\frac{1}{3}\tau^{3}_{n+1}\int_{t_{n}}^{t_{n+1}}|r_{tt}(s)|^{2}ds\lesssim\frac{\tau^{3}_{n+1}}{C_{0}}\int_{t_{n}}^{t_{n+1}}(\|\phi_{t}\|^{2}_{H^{1}}+\|\phi_{tt}\|^{2})dt,\\\\[11.0pt]
|v^{n}_{2}|^{2}\lesssim&\displaystyle\tau_{n+1}^{3}\Big{\|}\frac{F^{\prime}(\phi(t_{n+1}))}{2\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}}\Big{\|}^{2}\int_{t_{n}}^{t_{n+1}}\|\phi_{tt}(s)\|^{2}ds\\\\[11.0pt]
\lesssim&\displaystyle\frac{\tau^{3}_{n+1}}{C_{0}}\int_{t_{n}}^{t_{n+1}}\|\phi_{tt}(s)\|^{2}ds,\end{array}$
(3.22)

where the estimate for $r_{tt}$ in (3.9) is used. Then we take the above
estimates (3.18)-(3.22) into (LABEL:equ5), and sum up it from 1 to $n$ to get

$\begin{array}[]{r@{}l}&\displaystyle\frac{\gamma^{\frac{3}{2}}_{n+2}}{2(1+\gamma_{n+2})}\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\|^{2}}{\tau_{n+1}}+\frac{1}{2}\|(\Delta+\beta)\overline{e}^{n+1}_{\phi}\|^{2}+\frac{S}{2}\|\overline{e}^{n+1}_{\phi}\|^{2}+|e^{n+1}_{r}|^{2}\\\\[11.0pt]
\lesssim&\displaystyle\frac{\gamma^{\frac{3}{2}}_{2}}{2(1+\gamma_{2})}\frac{\|\nabla^{-1}\overline{e}^{1}_{\phi}\|^{2}}{\tau_{1}}+\frac{1}{2}\|(\Delta+\beta)\overline{e}^{1}_{\phi}\|^{2}+\frac{S}{2}\|\overline{e}^{1}_{\phi}\|^{2}+|e^{1}_{r}|^{2}\\\\[8.0pt]
&\displaystyle+\sum_{k=1}^{n}\tau_{k+1}\big{[}|e^{n+1}_{r}|^{2}+\|(\Delta+\beta)\overline{e}^{n}_{\phi}\|^{2}+\|(\Delta+\beta)\overline{e}^{n-1}_{\phi}\|^{2}+\|\overline{e}^{n}_{\phi}\|^{2}\\\\[11.0pt]
&\displaystyle+\|\overline{e}^{n-1}_{\phi}\|^{2}\big{]}+t_{n+1}\big{[}N^{-2m}\|\phi_{t}\|^{2}_{L^{\infty}(0,T;H^{m})}+N^{4-2m}\|\phi\|_{L^{\infty}(0,T;H^{m})}\big{]}\\\\[8.0pt]
&\displaystyle+\frac{\tau^{2}}{C_{0}}\int_{t_{1}}^{t_{n+1}}\big{(}\|\phi_{t}(s)\|^{2}_{H^{1}}+\|\phi_{tt}(s)\|^{2}\big{)}ds+\tau^{4}\int_{0}^{t_{n+1}}\big{(}\|\phi_{tt}(s)\|^{2}_{H^{1}}+\|\phi_{ttt}(s)\|^{2}_{H^{-1}}\big{)}ds.\end{array}$
(3.23)

We use the first order scheme (2.6) for the initial time step to get the
following error equation for $n=0$:

$\begin{array}[]{r@{}l}&\displaystyle(\frac{\overline{e}^{1}_{\phi}}{\tau_{1}},p_{N})+(\nabla\overline{e}^{1}_{\mu},\nabla
p_{N})=-(\tau_{1}^{-1}\int_{0}^{t_{1}}s\phi_{tt}(s)ds,p_{N}),\\\
&\displaystyle(\overline{e}^{1}_{\mu},q_{N})=\big{(}(\Delta+\beta)\overline{e}^{1}_{\phi},(\Delta+\beta)q_{N}\big{)}+S(\overline{e}^{1}_{\phi},q_{N})+\frac{e^{1}_{r}}{\sqrt{E^{0}_{1,N}+C_{0}}}(F^{\prime}(\phi^{0}_{N}),q_{N})\\\\[11.0pt]
&\hskip
31.2982pt\displaystyle+\Big{(}\frac{r(t_{1})F^{\prime}(\phi^{0}_{N})}{\sqrt{E^{0}_{1,N}+C_{0}}}-\frac{r(t_{1})F^{\prime}(\phi_{0}({\tenbi
x}))}{\sqrt{E_{1}(\phi(t_{1}))+C_{0}}},q_{N}\Big{)}-\Big{(}F^{\prime}(\phi_{0}({\tenbi
x}))-F^{\prime}(\phi(t_{1})),q_{N}\Big{)},\\\\[13.0pt] &\displaystyle
e^{1}_{r}=\frac{\int_{\Omega}F^{\prime}(\phi^{0}_{N})\overline{e}^{1}_{\phi}d{\tenbi
x}}{2\sqrt{E^{0}_{1,N}+C_{0}}}+\frac{1}{2}\int_{\Omega}\Big{[}\frac{F^{\prime}(\phi^{0}_{N})}{\sqrt{E^{0}_{1,N}+C_{0}}}-\frac{F^{\prime}(\phi(t_{1}))}{\sqrt{{\color[rgb]{0,0,0}E_{1}(\phi(t_{1}))+C_{0}}}}\Big{]}\Pi_{N}\big{(}\phi(t_{1})-\phi(t_{0})\big{)}d{\tenbi
x}\\\\[13.0pt] &\hskip
28.45274pt\displaystyle+\int_{\Omega}\frac{F^{\prime}(\phi(t_{1}))}{\sqrt{{\color[rgb]{0,0,0}E_{1}(\phi(t_{1}))+C_{0}}}}(\widehat{e}^{1}_{\phi}-\widehat{e}^{0}_{\phi})d{\tenbi
x}-v^{0}_{1}+v^{0}_{2},\end{array}$

where $\phi^{0}_{N}:=\Pi_{N}\phi_{0}({\tenbi x}).$ Similar to the argument for
the case $n\geq 1$, the following estimate can be easily derived with the
assumption $\tau_{1}\leq C\tau^{4/3}$

$\begin{array}[]{r@{}l}\displaystyle\frac{\|\nabla^{-1}\overline{e}^{1}_{\phi}\|^{2}}{2\tau_{1}}+\frac{1}{2}\|(\Delta+\beta)\overline{e}^{1}_{\phi}\|^{2}+\frac{S}{2}\|\overline{e}^{1}_{\phi}\|^{2}+|e^{1}_{r}|^{2}\lesssim&\displaystyle\tau_{1}^{3}+\tau_{1}N^{4-2m}\|\phi_{0}({\tenbi
x})\|_{H^{m}}\\\\[8.0pt]
\lesssim&\displaystyle\tau^{4}+\tau_{1}N^{4-2m}\|\phi\|_{L^{\infty}(0,T;H^{m})}.\end{array}$
(3.24)

Combining with (LABEL:equ9),(3.24) and (3.3), together with the discrete
Gronwall’s lemma and the triangle inequality, the desired estimate
(LABEL:equa3) is obtained.

Next, we will prove the estimate (LABEL:equa4), which can show that the
numerical solution $\phi^{n}$ is of second-order accuracy in time with
$C_{0}\geq{1}/{\tau}$. The error equation for $n\geq 1$ reads

$\begin{array}[]{r@{}l}&(F^{n+1}_{2}\overline{e}_{\phi},p_{N})+(\nabla(\Delta+\beta)\overline{e}^{n+1}_{\phi},\nabla(\Delta+\beta)p_{N})+S(\nabla\overline{e}^{n+1}_{\phi},\nabla
p_{N})\\\\[9.0pt]
=&\displaystyle\big{(}1-\frac{r^{n+1}}{\sqrt{E^{n}_{1,N}+C_{0}}}\big{)}(\nabla
F^{\prime}(\phi^{*,n+1}_{N}),\nabla
p_{N})-\big{(}\nabla(I^{n}_{1}+T^{n}_{2}),\nabla
p_{N}\big{)}-(T^{n}_{1},p_{N}).\end{array}$

Taking $p_{N}=(-\Delta)^{-1}(\delta_{t}\overline{e}^{n+1}_{\phi})$, and using
the Young’s inequality, we have

$\begin{array}[]{r@{}l}&\displaystyle\frac{1}{2}\Big{[}\big{(}\frac{\gamma^{\frac{3}{2}}_{n+2}}{1+\gamma_{n+2}}+G(\gamma_{*},\gamma_{*})\big{)}\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\|^{2}}{\tau_{n+1}}-\frac{\gamma^{\frac{3}{2}}_{n+1}}{1+\gamma_{n+1}}\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n}_{\phi}\|^{2}}{\tau_{n}}\Big{]}\\\\[9.0pt]
&\displaystyle+\frac{1}{2}\big{[}\|(\Delta+\beta)\overline{e}^{n+1}_{\phi}\|^{2}-\|(\Delta+\beta)\overline{e}^{n}_{\phi}\|^{2}\big{]}+\frac{S}{2}\big{[}\|\overline{e}^{n+1}_{\phi}\|^{2}-\|\overline{e}^{n}_{\phi}\|^{2}\big{]}\\\\[11.0pt]
\leq&\displaystyle
C(\gamma_{*})\tau_{n+1}\Big{[}\Big{|}1-\frac{r^{n+1}}{\sqrt{E^{n}_{1,N}+C_{0}}}\Big{|}^{2}+\|\nabla
I^{n}_{1}\|^{2}+\|\nabla^{-1}T^{n}_{1}\|^{2}+\|\nabla
T^{n}_{2}\|^{2}\Big{]}\\\
&\displaystyle+G(\gamma_{*},\gamma_{*})\frac{\|\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\|^{2}}{2\tau_{n+1}}.\end{array}$
(3.25)

Combined the estimate (3.18) for $I^{n}_{2}$ with (LABEL:equa3), it gives

$\begin{array}[]{r@{}l}&\displaystyle\Big{|}1-\frac{r^{n+1}}{\sqrt{E^{n}_{1,N}+C_{0}}}\Big{|}^{2}\\\
=&\displaystyle\Big{|}\frac{r^{n+1}}{\sqrt{E^{n}_{1,N}+C_{0}}}-\frac{r(t_{n+1})}{\sqrt{E_{1}(\phi(t_{n+1}))+C_{0}}}\Big{|}^{2}=\Big{|}\frac{e^{n+1}_{r}}{\sqrt{E^{n}_{1,N}+C_{0}}}+r(t_{n+1})\cdot
I^{n}_{2}\Big{|}^{2}\\\\[11.0pt]
\lesssim&\displaystyle\frac{|e^{n+1}_{r}|^{2}}{C_{0}}+\frac{1}{C^{2}_{0}}\Big{(}\tau_{n+1}\int_{t_{n}}^{t_{n+1}}\|\phi_{t}(s)\|^{2}ds+\|\overline{e}^{n}_{\phi}\|^{2}+N^{-2m}\|\phi\|_{L^{\infty}(0,T;H^{m})}\Big{)}\\\
\lesssim&\displaystyle\frac{\tau^{2}}{C^{2}_{0}}\int_{0}^{T}\big{(}\|\phi_{t}(s)\|^{2}_{H^{1}}+\|\phi_{tt}(s)\|^{2}\big{)}ds+\frac{\tau^{4}}{C_{0}}\int_{0}^{T}\big{(}\|\phi_{tt}(s)\|^{2}_{H^{1}}+\|\phi_{ttt}(s)\|^{2}_{H^{-1}}\big{)}ds\\\
&\displaystyle+\frac{1}{C^{2}_{0}}\Big{(}\tau_{n+1}\int_{t_{n}}^{t_{n+1}}\|\phi_{t}(s)\|^{2}ds+\|\overline{e}^{n}_{\phi}\|^{2}+N^{-2m}\|\phi\|_{L^{\infty}(0,T;H^{m})}\Big{)}.\end{array}$
(3.26)

Consequently, taking (LABEL:equa5), the estimates of $\|I^{n}_{1}\|^{2}$,
$\|T^{n}_{1}\|^{2}$ and $\|T^{n}_{2}\|^{2}$ into (LABEL:esti_phi), summing it
up from 1 to n, and applying the discrete Gronwall’s lemma, we can deduce the
desired estimate (LABEL:equa4) and (LABEL:equa4_1). Then the proof is
completed. ∎

######  Remark 3.1.

Our error norms contain the terms $\|\nabla^{-1}e^{1}_{\phi}\|^{2}/\tau_{1}$
for $n=0$ in (3.24) and
$\|\nabla^{-1}\delta_{t}\overline{e}^{n+1}_{\phi}\|/\tau_{n+1}$ for $n\geq 1$
in (LABEL:equa4_1). Thus, in order to guarantee the second order accuracy of
the fully discrete scheme (2.15) with this type of error norm in
(LABEL:equa4_1), it is necessary to employ the first order scheme (2.6) with
$\tau_{1}\leq\tau^{4/3}$ to start the time marching. It is worth noting that
the constant $C\exp(T)$ in the error estimates (3.5)-(3.7) would be rather
large, which may be not optimal estimates particularly for the long-time
simulation case. This exponential constant is commonly appeared in the error
analysis for the time-stepping schemes using the discrete Gronwall’s
inequality as an analysis tool. To avoid using the discrete Gronwall’s
inequality, another possible way to establish the corresponding error
estimates is to adopt the mathematical induction method. However, it is not
easy to give a detailed uniform bound for the general positive constant $C$
such as $C\leq 1$ in the derivation of the error estimate, which is a key
issue in successfully deriving the error estimates by the mathematical
induction approach. It is also an open, interesting problem now and needs more
further effort.

##  4\. Numerical results

This section is devoted to numerically validating the theoretical results in
terms of stability and accuracy. For simplicity, we set
$\beta=1,S=\varepsilon,$ and the errors are computed by $L^{\infty}$-norm
throughout the numerical tests.

###  4.1. Test of the convergence order

In this subsection, we investigate the accuracy of second order in time and
spectral convergence in space for the fully discrete scheme (2.15). Here, we
consider the PFC model with $\varepsilon=0.025$ and the initial data, given by

$\phi_{0}(x,y)=\sin\Big{(}\frac{\pi x}{16}\Big{)}\cos\Big{(}\frac{\pi
y}{16}\Big{)}.$

The computational domain is set to be $(0,32)\times(0,32)$. At $T=1$, a
reference solution of the phase function $\phi^{n}$ is computed by the fully
discrete scheme (2.15) with a fixed small uniform time step size $\Delta
t=1e-5$ and a $256\times 256$ space mesh. We first check the temporal accuracy
on the uniform temporal mesh. In Figure 1 (a), the temporal errors of the
numerical solutions of the phase function $\phi^{n}$ with different time step
sizes are plotted in log-log scale. It is shown that the scheme (2.15) with
$C_{0}=1/\tau$ achieves the expected second order accuracy in time. We also
display the second order approximation of the computed
$r^{n+1}/\sqrt{E^{n}_{1,N}+C_{0}}$ to $1$ in Figure 1 (b). Then we study the
accuracy of the spatial discretization. In Figure 1 (c), we plot the
$L^{\infty}$-errors in semi-log scale with respect to the Fourier mode number
$N$. The error curve is almost a straight line in the semi-log scale, which
indicates that the numerical solutions are exponentially convergent with
respect to $N$. This is consistent with our theoretical results.

Next, we investigate the temporal error behaviors of the fully discrete scheme
(2.15) on the nonuniform temporal mesh to validate our theoretical error
estimates in Theorem 3.1. The nonuniform temporal mesh
$\\{\widehat{t}_{n}\\}_{n=0}^{N}$ used here is produced by $40\%$ random
perturbation of the uniform mesh $\\{t_{n}=n\tau\\}_{n=0}^{N}$, i.e., the
random time grids are given by

$\widehat{t}_{n}:=t_{n}+0.4\tau\cdot rand(\cdot),\quad n=0,1,\cdots,N,$

where $rand(\cdot)$ is a random data from -1 to 1. The convergence order is
computed by

$\mbox{Order}=\frac{\log_{10}(\|e^{M}_{\phi}\|_{\infty}/\|e^{2M}_{\phi}\|_{\infty})}{\log_{10}(\tau(M)/\tau(2M))},$

where $e^{M}_{\phi}$ and $\tau(M)$ denote the $L^{\infty}$-norm error of
$\phi$ at the finial time T and the maximum time-step size for a total of M
subintervals on $[0,T]$, respectively. In Table 1, we present convergence
rates in time for the phase function $\phi^{n}$ in (2.15) at $T=1$ with
$\sigma=$ 1/2, 2/3 and 1. It is observed the proposed scheme (2.15) achieves
about second order accuracy for the phase function $\phi^{n}$ for all tested
$\sigma$ on the nonuniform mesh, even for the cases with the adjacent time
step radio $\gamma_{n}>4.8645$. Thus it indicates that $\gamma_{**}=4.8645$
may be not the optimal upper bound on the adjacent time step ratios, and
deserves further investigation. In addition, compared to the numerical
convergence results for the case of the uniform temporal meshes in Figure 1,
the convergence order of the proposed scheme shown in Table 1 achieves obvious
oscillations due to the randomness of the random nonuniform meshes used.

(a) Convergence rates in time

(b) $\displaystyle\max_{n}\big{|}r^{n+1}/\sqrt{E^{n}_{1,N}+C_{0}}-1\big{|}$

(c) Convergence rates in space

Figure 1. Convergence rates at $T=1$ for the scheme (2.15) with $C_{0}=1/\tau$.  Table 1.  Convergence rates in time at $T=1$ for the phase function $\phi^{n}$ in (2.15) on nonuniform temporal meshes with $C_{0}=1/\tau$. temporal mesh | $\sigma=\frac{1}{2}$ | $\sigma=\frac{2}{3}$ | $\sigma=1$
---|---|---|---
$M$ | $\tau$ | $\max\\{\gamma_{n}\\}$ | error | order | error | order | error | order
$20$ | 7.91e-2 | 3.64 | 8.30e-5 | - | 1.01e-4 | - | 1.45e-4 | -
$40$ | 4.12e-2 | 4.79 | 1.42e-6 | 2.71 | 1.93e-5 | 2.55 | 2.95e-5 | 2.44
$80$ | 2.21e-2 | 4.16 | 4.49e-6 | 1.84 | 4.72e-6 | 2.26 | 6.11e-6 | 2.52
$160$ | 1.11e-2 | 5.41 | 7.53e-7 | 2.58 | 8.57e-7 | 2.47 | 1.28e-6 | 2.25
$320$ | 5.55e-3 | 6.40 | 8.55e-8 | 3.16 | 1.05e-7 | 3.00 | 1.52e-7 | 3.10
$640$ | 2.70e-3 | 5.24 | 2.18e-8 | 1.89 | 2.84e-8 | 1.81 | 4.09e-8 | 1.81
$1280$ | 1.39e-3 | 7.06 | 5.71e-9 | 2.03 | 7.12e-9 | 2.10 | 1.00e-8 | 2.12

###  4.2. Time adaptive strategy

In this part, the efficiency of the fully discrete scheme (2.15) with a time
adaptive strategy is investigated through a long time phase transition
process. As we know, the phase transition process usually goes through
different stages, changing quickly at the beginning and then rather slowly
afterwards, until it reaches a steady state after a long time evolution. Thus
it would be very efficient to apply some time adaptive strategy in the
simulation based on our proposed nonuniform BDF2 scheme (2.15).

We consider the phase transition process governed by the PFC equation (1.2)
with $\varepsilon=0.1$ and the initial data $\phi^{0}_{i,j}=0.08+\eta_{i,j}$,
where $\eta_{i,j}$ is a random number between $-0.08$ and $0.08$. Here, the
computational domain is set to be $\Omega=(0,256)^{2}$. The time adaptive
strategy proposed in [23] is adopted in the simulation with our fully discrete
scheme (2.15):

$\displaystyle\tau_{n+1}=\min\Big{(}\max\big{(}\tau_{min},\frac{\tau_{max}}{\sqrt{1+\gamma|E^{{}^{\prime}}(t)|^{2}}}\big{)},4.8645\cdot\tau_{n}\Big{)},$
(4.1)

where $\tau_{min},\tau_{max}$ are predetermined minimum and maximum time step
sizes, and $\gamma$ is a positive constant. This type of the time adaptive
strategy is based on the energy variation. It will automatically select small
time steps when the energy decays rapidly and choose larger time steps
otherwise. In this simulation, we set $\tau_{min}=0.01$, $\tau_{max}=5$ and
$\gamma=10^{5}$. Furthermore, $512\times 512$ Fourier modes are used for the
spatial discretization.  In order to capture the characteristics of the phase
transition process with rapid energy changes where the small time step size
$\tau_{min}$ is used, we set $C_{0}=1/\tau_{min}$ to improve the local
accuracy of our proposed numerical scheme at these times of rapid energy
changes. It has also been checked that the chosen $C_{0}$ is large enough such
that $E_{1}(\cdot)+C_{0}>0.$

In Figure 2, we plot the evolutions of the phase function $\phi^{n}$ up to
$T=5000$ with four different types of temporal meshes, viz., the fixed large
time step mesh with $\tau=5$, two different ones from the time adaptive
strategy (4.1), and the uniform small time step mesh with $\tau=0.01$. In
Figure 3 (a), it shows that there is no obvious difference at $t=20$, and the
corresponding energies consist with each other at this stage. After the energy
goes through a large variation stage, the fixed large time step with $\tau=5$
yields inaccurate numerical solutions, while both two tested adaptive time
strategies with $\tau_{max}=1$ and $5$ give correct phase transition patterns
consistent with the results obtained by the small time step case $\tau=0.01$.
It is observed in Figure 3 (a) and (b) that both the computed modified energy
$\widetilde{E}^{n}_{N}$ and the computed original energy $E(\phi^{n}_{N})$ are
dissipative in time for the tested four different types of temporal meshes.
Moreover, as shown in Figure 3 (b), the energy dissipation of the proposed
scheme (2.15) with the adaptive strategy (4.1) consists very well with the one
using the fixed small time step $\tau=0.01$. The corresponding time step sizes
are also plotted in Figure 3 (c), which shows the efficiency of our proposed
scheme (2.15) with the adaptive strategy (4.1) during this long time
simulation. Furthermore, using a smaller $\tau_{max}$ in adaptive time
strategy (4.1) certainly improves the accuracy more or less, but may lead to
considerable increase in the computational cost, as shown in Figure 3 (c).

(a) solutions with fixed large time step size $\tau=5$.

(b) solutions with adaptive time steps with $\tau_{min}=0.01,\tau_{max}=5$ and
$\gamma=10^{5}.$

(c) solutions with adaptive time steps with $\tau_{min}=0.01,\tau_{max}=1$ and
$\gamma=10^{5}.$

(d) solutions with fixed small time step size $\tau=0.01$.

Figure 2. Solution snapshots of the phase transition process for the PFC
equation at around $t=20,50,500,1500,$ and $5000$, respectively.

(a) the computed modified energy $\widetilde{E}^{n}_{N}-C_{0}$

(b) the computed energy $E(\phi^{n})$

(c) adaptive time step sizes with $\tau_{max}=1$ and $\tau_{max}=5$

Figure 3. Evolution in time of the modified energy, the original energy, and
the adaptive time step sizes.

###  4.3. Crystal growth in a supercooled liquid

In this subsection, we simulate the growth of a polycrystal in a supercooled
liquid. Initially, three crystallites with different orientations are given by

$\phi_{0}(x_{l},y_{l})=0.285+0.446\Big{(}\cos\Big{(}\frac{0.66}{\sqrt{3}}y_{l}\Big{)}\cos(0.66x_{l})-0.5\cos\Big{(}\frac{1.32}{\sqrt{3}}y_{l}\Big{)}\Big{)},~{}~{}l=1,2,3,$

where $x_{l}$ and $y_{l}$ define a local system of cartesian coordinates that
is oriented with the crystallite lattice. Similar numerical examples can be
found in [6, 9, 16, 20, 31]. The computed domain is $(0,800)^{2}$,
$\varepsilon=0.25$, and the initial three crystallites are located in the
blocks with length of 40 in the computed domain, seeing Figure 4 for the case
of $t=0$. In order to obtain three different orientations of the crystallites,
we define the local coordinates $(x_{l},y_{l})$ by using the following affine
transformation of the global coordinates $(x,y)$,

$\begin{array}[]{r@{}l}\begin{cases}x_{l}=x\cos(\theta)-y\sin(\theta),\\\
y_{l}=x\sin(\theta)+y\cos(\theta),\end{cases}\end{array}$

where $\theta=-\pi/4,0,\pi/4$ for $l=1,2,3$, respectively. This produces a
rotation given by the angle $\theta$.

The simulation is performed by the second order scheme (2.15) with $\sigma=1$
and $1024\times 1024$ Fourier modes. Here, we use two types of temporal
meshes: the fixed small time step with $\tau=0.01$ and the adaptive strategy
(4.1) with $\tau_{min}=0.01,\tau_{max}=1$ and $\gamma=10.$ In Figure 4, we
display the evolution of the crystal growth up to $T=1500$ for these two types
of temporal meshes, respectively. It is observed that the snapshots of the
crystal growth consist with each other at same time stages for these two
temporal meshes. The growth of the crystalline phase and the motion of well-
defined crystal–liquid interfaces are observed. Moreover, it shows that the
different alignment of the crystallites causes defects and dislocations. This
phenomenon has been also observed in [6, 9, 16, 20, 31]. In Figure 5, the
evolution of the modified energy, the original energy, and the corresponding
adaptive time step sizes are plotted to validate the stability and the
efficiency of our proposed scheme (2.15) with the adaptive strategy (4.1).

######  Remark 4.1.

As shown in Figure 3 and Figure 5, the evolutions of the computed original
energy with the tested time adaptive strategies are in agreement with the one
computed by the small time step $\tau=0.01$ better than the case of the
computed modified energy. This may be due to the computed auxiliary variable
$r^{n}$ and the extra term
$\frac{(2\sigma-1)\gamma_{n+1}^{3/2}}{2+2\gamma_{n+1}}\frac{\|\nabla^{-1}(\phi^{n}_{N}-\phi^{n-1}_{N})\|^{2}}{\tau_{n}}$
contained in the modified energy in (2.16).

(a) Solutions with fixed time step size $\tau=0.01$

(b) Solutions with adaptive time steps with $\tau_{min}=0.01,\tau_{max}=1$ and
$\gamma=10.$

Figure 4. Solution snapshots of crystal growth for the PFC equation at around
$t=0,100,200,400,800$ and $1500$, respectively.

(a) the computed modified energy $\widetilde{E}^{n}_{N}-C_{0}$

(b) the computed energy $E(\phi^{n})$

(c) adaptive time step sizes with $\gamma=10$

Figure 5. Evolution in time of the energy and the adaptive time step size for
crystal growth.

##  5\. Concluding remarks

We have developed and analyzed a linear second order nonuniform BDF scheme
based on the SAV approach for the PFC model. The energy stability of the
proposed scheme is established for the nonuniform temporal mesh with a mild
restriction on the adjacent time step radio $\gamma_{n}\leq 4.8645.$ Although
we only use a first order method to approximate the auxiliary variable $r(t)$,
the second order accuracy of the phase function $\phi^{n}$ can be guaranteed
by setting a large enough positive constant $C_{0}$ such that $C_{0}\geq
1/\tau$. Moreover, a rigorous error estimate of our proposed fully discrete
nonuniform BDF2 scheme is established with some proper assumptions on the
regularity of the solution. Finally, a series of numerical tests have been
carried out to validate the theoretical results and the efficiency of the
proposed scheme combined with the time adaptive strategy.

## References

  * [1] A. Baskaran, J. S. Lowengrub, C. Wang, and S. M. Wise.  Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation.  SIAM J. Numer. Anal., 51(5):2851–2873, 2013.
  * [2] W. Chen, X. Wang, Y. Yan, and Z. Zhang.  A second order BDF numerical scheme with variable steps for the Cahn–Hilliard equation.  SIAM J. Numer. Anal., 57(1):495–525, 2019.
  * [3] L. Dong, W. Feng, C. Wang, S. M. Wise, and Z. Zhang.  Convergence analysis and numerical implementation of a second order numerical scheme for the three-dimensional phase field crystal equation.  Comput. Math. Appl., 75(6):1912–1928, 2018.
  * [4] Q. Du, L. Ju, X. Li, and Z. Qiao.  Maximum bound principles for a class of semilinear parabolic equations and exponential time differencing schemes.  SIAM Rev., 63(2):317–359, 2021.
  * [5] K. R. Elder and M. Grant.  Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals.  Phys. Rev. E, 70(5):051605, 2004.
  * [6] K. R. Elder, M. Katakowski, M. Haataja, and M. Grant.  Modeling elasticity in crystal growth.  Phys. Rev. Lett., 88(24):245701, 2002.
  * [7] C. M. Elliott and A. M. Stuart.  The global dynamics of discrete semilinear parabolic equations.  SIAM J. Numer. Anal., 30(6):1622–1663, 1993.
  * [8] D. J. Eyre.  Unconditionally Gradient Stable Time Marching the Cahn-Hilliard Equation.  Mater. Res. Soc. Sympos. Proc., 529:39–46, 1998.
  * [9] H. Gomez and X. Nogueira.  An unconditionally energy-stable method for the phase field crystal equation.  Comput. Methods Appl. Mech. Eng., 249-252:52–61, 2012.
  * [10] Z. Guan, J. S. Lowengrub, C. Wang, and S. M. Wise.  Second order convex splitting schemes for periodic nonlocal Cahn-Hilliard and Allen-Cahn equations.  J. Comput. Phys., 277:48–71, 2014.
  * [11] R. Guo and Y. Xu.  Local discontinuous galerkin method and high order semi-implicit scheme for the phase field crystal equation.  SIAM J. Sci. Comput., 38(1):A105–A127, 2016.
  * [12] D. Hou, M. Azaiez, and C. Xu.  A variant of scalar auxiliary variable approaches for gradient flows.  J. Comput. Phys., 395:307–332, 2019.
  * [13] D. Hou and Z. Qiao.  An implicit–explicit second–order BDF numerical scheme with variable steps for gradient flows.  J. Sci. Comput., 94(2):39, 2023.
  * [14] D. Hou and C. Xu.  Robust and stable schemes for time fractional molecular beam epitaxial growth model using SAV approach.  J. Comput. Phys, 445(15):110628, 2021.
  * [15] D. Hou and C. Xu.  A second order energy dissipative schemes for time fractional $L^{2}$ gradient flows using SAV approach.  J. Sci. Comput., 90:25, 2022.
  * [16] Z. Hu, S. M. Wise, C. Wang, and J. S. Lowengrub.  Stable and efficient finite-difference nonlinear-multigrid schemes for the phase field crystal equation.  J. Comput. Phys., 228(15):5323–5339, 2009.
  * [17] L. Ju, X. Li, and Z. Qiao.  Generalized SAV–exponential integrator schemes for Allen–Cahn type gradient flows.  SIAM J. Numer. Anal., 60(4):1905–1931, 2022.
  * [18] L. Ju, X. Li, and Z. Qiao.  Stabilized exponential–SAV schemes preserving energy dissipation law and maximum bound principle for the Allen–Cahn type equations.  J. Sci. Comput., 92(2):66, 2022.
  * [19] Q. Li, L. Mei, and B. You.  A second-order, uniquely solvable, energy stable BDF numerical scheme for the phase field crystal model.  Appl. Numer. Math., 134:46–65, 2018.
  * [20] X. Li and J. Shen.  Stability and error estimates of the SAV Fourier-spectral method for the phase field crystal equation.  Adv. Comput. Math., 46:48, 2020.
  * [21] H. Liao, B. Ji, and L. Zhang.  An adaptive BDF2 implicit time-stepping method for the phase field crystal model.  IMA J. Numer. Anal., 42(1):649–679, 2022.
  * [22] Z. Liu and X. Li.  The fast scalar auxiliary variable approach with unconditional energy stability for nonlocal Cahn–Hilliard equation.  Numer. Methods Partial Differ. Equ., 37(1):244–261, 2021.
  * [23] Z. Qiao, Z. Zhang, and T. Tang.  An adaptive time-stepping strategy for the molecular beam epitaxy models.  SIAM J. Sci. Comput., 33(3):1395–1414, 2011.
  * [24] J. Shen, T. Tang, and L. L. Wang.  Spectral methods, Algorithms, Analysis and Applications.  Springer, Berlin, 2010.
  * [25] J. Shen and J. Xu.  Convergence and error analysis for the scalar auxiliary variable(SAV) schemes to gradient flows.  SIAM J. Numer. Anal., 56(5):2895–2912, 2018.
  * [26] J. Shen, J. Xu, and J. Yang.  The scalar auxiliary variable (SAV) approach for gradient flows.  J. Comput. Phys., 353:407–416, 2018.
  * [27] J. Shen, J. Xu, and J. Yang.  A new class of efficient and robust energy stable schemes for gradient flows.  SIAM Rev., 61(3):474–506, 2019.
  * [28] C. Wang and S. M. Wise.  An energy stable and convergent finite–difference scheme for the modified phase field crystal equation.  SIAM J. Numer. Anal., 49(3):945–969, 2011.
  * [29] S. M. Wise, C. Wang, and J. S. Lowengrub.  An energy-stable and convergence finite–difference scheme for the phase field crystal equation.  SIAM J. Numer. Anal., 47(1):2269–2288, 2009.
  * [30] X. Yang.  Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends.  J. Comput. Phys., 327:294–316, 2016.
  * [31] X. Yang and D. Han.  Linearly first- and second-order, unconditionally energy stable schemes for the phase field crystal model.  J. Comput. Phys., 88(330):1116–1134, 2017.
  * [32] Z. Yang, L. Lin, and S. Dong.  A family of second-order energy-stable schemes for Cahn–Hilliard type equations.  J. Comput. Phys., 383:24–54, 2019.
  * [33] Z. Zhang, Y. Ma, and Z. Qiao.  An adaptive time-stepping strategy for solving the phase field crystal model.  J. Comput. Phys., 249:204–215, 2013.
  * [34] J. Zhao, Q. Wang, and X. Yang.  Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach.  Int. J. Numer. Meth. Eng., 110(3):279–300, 2017.


# Exploring single-song autoencoding schemes
for audio-based music structure analysis

###### Abstract

The ability of deep neural networks to learn complex data relations and
representations is established nowadays, but it generally relies on large sets
of training data. This work explores a “piece-specific” autoencoding scheme,
in which a low-dimensional autoencoder is trained to learn a latent/compressed
representation specific to a given song, which can then be used to infer the
song structure. Such a model does not rely on supervision nor annotations,
which are well-known to be tedious to collect and often ambiguous in Music
Structure Analysis. We report that the proposed unsupervised auto-encoding
scheme achieves the level of performance of supervised state-of-the-art
methods with 3 seconds tolerance when using a Log Mel spectrogram
representation on the RWC-Pop dataset.

Index Terms—  Autoencoders, Music Structure Analysis, Audio Signal Processing

##  1 Introduction

Music Structure Analysis (MSA) consists in segmenting a music piece in several
distinct parts, which represent a mid-level description. Segmentation is
usually based on criteria such as homogeneity, novelty, repetition and
regularity [1]. In practice, MSA often relies on similarity between passages
of a song summarized in an autosimilarity matrix [2, 3, 4, 5], in which each
coefficient represents an estimation of the similarity between two musical
fragments.

Recently, Deep Neural Networks (DNN) have lead to some high level of
excitement in Music Information Retrieval (MIR) research, and notably in MSA
[6, 7]. In general, DNN approaches rely on large databases which make it
possible to learn a large number of parameters, which in turn yields better
performance than previously established machine learning approaches. This is
the consequence of the ability of DNNs to learn complex nonlinear mappings
through which musical objects can be well separated [8].

However, DNNs generally learn and exploit “deep” features stemming from
multiple examples used in a training phase, and then evaluate the potential of
learned features in a test phase, where no optimization is performed.

In this work, we consider a different DNN approach: we train a song-dependent
autoencoder (AE), i.e. an AE which is specifically trained to compress a given
song. By doing so, we make the hypothesis that the ability of AE to learn
compressed representations of musical objects occurring within a given song
will provide a set of features that are useful to infer the structure of that
song.

To study the relevance of this approach, this work presents a Convolutional
Neural Network (CNN)-based technique for MSA, evaluated on the RWC-Pop and
SALAMI datasets [9, 10] in their audio form. Songs are represented with
different features.

The rest of the article is structured as follows: Section 2 presents the
formal concepts framing our work, Section 3 presents the details of our
approach for music processing, Section 4 presents the segmentation process,
and Section 5 reports on the experimental results on both datasets.

##  2 Autoencoders

###  2.1 Basics of autoencoders

Autoencoders are neural networks, which, by design, perform unsupervised
dimensionality reduction. Throughout the years, AE have received increasing
interest, notably due to their ability to extract salient latent
representations without the need of large amount of annotations. Recently, AE
also showed great results as a generation tool [11].

An autoencoder is divided in two parts: an encoder, which compresses the input
$x\in\mathbb{R}^{n}$ into a latent representation $z\in\mathbb{R}^{m}$ of
smaller dimension (generally, $m\ll n$), and a decoder, which reconstructs
$\hat{x}$ from $z$.

In this work, layers are of two types: fully-connected, or convolutional.
Convolutional layers lead to impressive results in image processing due to
their ability to discover local correlations (such as lines or edges), which
turn to higher-order features with the depth of the network [12]. While local
correlations are less obvious in spectrogram processing [13], Convolutional
Neural Networks (CNN) still perform well in MIR tasks, such as MSA [6].

###  2.2 AE architecture

The tested AE is a CNN. We use the nowadays quite standard ReLu function as
the activation function, except in the latent layer because it could lead to
null latent variables. The encoder is composed of five hidden layers: 2
convolutional/max-pooling blocks, followed by a fully-connected layer,
controlling the size of the latent space. Convolutional kernels are of size
3x3, and the pooling is of size 2x2. Convolutional layers define respectively
4 and 16 feature maps.

The decoder is composed of 3 hidden layers: a fully-connected layer (inverse
of the previous one) and 2 “transposed convolutional” layers of size 3x3 and
stride 2x2. A transposed convolution is similar to the convolution operation
taken in the backward pass: an operation which takes one scalar as input and
returns several scalars as output [14, Ch. 4]. Hence, it is well suited to
inverse the convolution process. This network is represented in Figure 1.

Fig. 1: Architecture of the autoencoder

##  3 AE for Music Processing

###  3.1 Motivations

The underlying idea of this work is that music structure can be related to
compression of information. Indeed, a common view of music structure is to
consider structural segments as internally coherent passages, and automatic
retrieval techniques generally focus on finding them by maximizing homogeneity
and repetition and/or by setting boundaries between dissimilar segments, as
points of high novelty. In the context of compressed representations, such as
AEs, each passage is transformed in a vector of small dimension, compelling
this representation to summarize the original content. From this angle,
similar passages are expected to be represented by similar representations, as
they share underlying properties (such as coherence and redundancy), while
dissimilar passages are bound to create strong discrepancies at their
boundaries. Thus, we expect that compressed representations will enhance the
original structure while reducing incidental signal-wise properties which do
not participate to the structure.

###  3.2 Barwise music processing

Following our former work in [5], we process music as barwise spectrograms,
with a fixed number of frames per bar.

Conversely to fixed-size frame analysis, barwise computation guarantees that
the information contained in each frame does not depend on the tempo, but on
the metrical positions, which is a more abstract musical notion to describe
time. As a consequence, comparisons of bars are more reliable than comparisons
of frames of arbitrarily fixed size as it allows to cope with small variations
of tempo.

In addition, pop music (i.e. our case study) is generally quite regular at the
bar level: repetitions occur rather clearly at the bar scale (or its
multiples) and motivic patterns tend to develop within a limited number of
bars, suggesting that frontiers mainly occur between bars.

Accordingly, this method relies on a consistent bar division of music. It also
requires a powerful tool to detect bars, as otherwise errors could propagate
and affect the performance. Results reported in [5] tend to show that the
madmom toolbox [15] is efficient in that respect on the RWC-Pop dataset.
Nonetheless, this processing may hinder the retrieval of boundaries based on a
change of tempo.

Following this approach, each bar is represented with 96-equally spaced
frames, selected from a high-resolution spectrogram (computed with a hop-
length of 32 samples) [5].

###  3.3 Features

Music is represented with different features throughout our experiments,
focusing on different aspects of music such as harmony or timbre.

Chromagram A chromagram represents the time-frequency aspect of music as
sequences of 12-row vectors, corresponding to the 12 semi-tones of the
classical western music chromatic scale (C, C#, …, B), which is largely used
in Pop music. Each row represents the weight of a semi-tone (and its octave
counterparts) at a particular instant.

Mel spectrogram A Mel spectrogram corresponds to the STFT representation of a
song, whose frequency bins are recast in the Mel scale. Mel spectrograms are
dimensioned following the work of [6] (80 filters, from 80Hz to 16kHz). STFT
are computed as power spectrograms.

Log Mel spectrogram A Log Mel spectrogram corresponds to the logarithmic
values of the precedent spectrogram. This representation accounts for the
exponential decay of frequency power.

MFCC spectrogram Mel-Frequency Cepstral Coefficients (MFCCs) are timbre-
related coefficients, obtained by a discrete cosine transform of the Log Mel
spectrogram. This spectrogram contains 32 coefficients, following [16].

##  4 Structural segmentation

The ability of a single-song AE to separate and group bars is evaluated on the
MSA task, as presented in [17].

Thereby, for a given song in a given representation, we compute latent
representations ${z_{b}}$ for all bars $1\leq b\leq B$ of the song. These
latent representations are at the heart of the segmentation strategy, via
their autosimilarity matrix.

Denoting $Z\in\mathbb{R}^{d_{ls}\times B}$ the matrix resulting from the
concatenation of all $d_{ls}$-dimensional $z_{b}$ vectors, its autosimilarity
is defined as $Z^{T}Z$, i.e. the $B\times B$ matrix of the dot products
between every $z_{b}$. Examples of autosimilarities are shown in Figure 2.

Fig. 2: Barwise autosimilarity of the raw Log Mel spectrogram (left) and of
the latent representations (right), computed on the song “Pop01” from RWC-Pop.

Following the dynamic programming algorithm inspired from [18] and presented
in [5], the segmentation is computed as the global maximum of the sum of all
individual segments costs. The cost of each segment is obtained by applying a
sliding convolving kernel on the diagonal of the autosimilarity matrix. This
convolving kernel is a square matrix, the size of which is that of the
potential segment.

Contrary to the novelty kernel of Foote [2], largely described in the
literature [1], the proposed kernel aims at framing square blocks of high
similarity, occurring when consecutive bars are similar. In that sense, while
Foote’s kernel focuses on novelty, ours favors homogeneity. Still, the dynamic
programming algorithm allows to find a trade-off between maximizing
homogeneity of individual blocks and novelty across them as a whole. When the
diagonal of the autosimilarity is structured in several self-similar blocks,
the algorithm frames and partitions these blocks.

210 Fig. 3: Kernel of size 10

The diagonal is equal to zero, so as to disregard perfect self-similarity of
each bar with itself (normalized to 1). The other elements of this kernel are
designed so as to emphasize the short-term similarity in the 4 contiguous
bars, and still catch longer-term similarity for long segments. This kernel
performs best in comparative
experiments111https://gitlab.inria.fr/amarmore/musicntd/-/blob/v0.2.0/Notebooks/5%20-%20Different%20kernels.ipynb
of [5].

This raw convolutional cost is combined with a(n ir-)regularity cost, favoring
frequent segments’ sizes [18], adapted as in [5].

##  5 Experiments

This segmentation pipeline is studied on the RWC-Pop dataset, which consists
in 100 Pop songs of high recording quality [9], along with the MIREX10
annotations [19]. The best autoencoding method is also evaluated on the SALAMI
dataset [10]. We focus on boundary retrieval and ignore segment labelling.

Boundaries are evaluated using the hit-rate metric, which considers a boundary
valid if it is approximately equal to an annotation, within a fixed tolerance
window. Consistently with MIREX standards [20], tolerances are equal to 0.5s
and 3s. The hit-rate is then expressed in terms of Precision, Recall, and
F-measure.

###  5.1 Related work

Table 1: Results of best-performing AEs and state-of-the-art on the RWC-Pop dataset. (*Results: 2015 MIREX contest [20].) |  | $P_{0.5}$ | $R_{0.5}$ | $F_{0.5}$ | $P_{3}$ | $R_{3}$ | $F_{3}$
---|---|---|---|---|---|---|---
Latent representations autosimilarity | Chromagram | 51.8% | 54% | 52.3% | 70.5% | 73.3% | 71.2%
Mel Spectrogram | 52.1% | 55.3% | 53.3% | 73.1% | 77.6% | 74.7%
Log Mel Spectrogram | 59.5% | 61.3% | 59.9% | 79.1% | 81.7% | 79.9%
MFCC | 54.3% | 54.4% | 54% | 76.6% | 76.4% | 76%
Barwise feature representation autosimilarity | Chromagram | 39.6% | 38.4% | 38.6% | 63.8% | 61.2% | 62%
Mel Spectrogram | 46.4% | 47% | 46.3% | 70.9% | 71.3% | 70.6%
Log Mel Spectrogram | 46.5% | 46.1% | 46% | 72.2% | 71% | 71.1%
MFCC | 41.9% | 40.6% | 40.9% | 68.1% | 65.6% | 66.4%
State-of-the-art | Foote [2] | 42,0% | 30,0% | 34,5% | 67,1% | 47,7% | 55,0%
Spectral clustering [3] | 49.2% | 45% | 45% | 65.5% | 60.6% | 60.3%
NTD [5] | 58.4% | 60.7% | 59.0% | 72.5% | 75.3% | 73.2%
Supervised CNN [6]* | 80.4% | 62.7% | 69.7% | 91.9% | 71.1% | 79.3%

Related work mainly splits in two categories: blind and learning-based
techniques.

Blind segmentation techniques, such as this work, do not use training datasets
and generally focus on autosimilarity matrices. The seminal paper in this area
is that by Foote [2], based on a novelty kernel, which detects boundaries as
points of strong dissimilarity between the near past and future of a given
instant. The strong point of this algorithm is its good performance given its
relative simplicity.

Later on, McFee and Ellis made use of spectral clustering as a segmentation
method by interpreting the autosimilarity as a graph and clustering
principally connected vertices as segments [3]. This work performed best among
blind segmentation techniques in the last structural segmentation MIREX
campaign in 2016 [20]. Recently, we used tensor decomposition (Nonnegative
Tucker Decomposition, NTD) as a way to describe music as barwise patterns,
which then served as features for the computation of an autosimilarity matrix
[5].

Still, to the best of the authors’ knowledge, the current state-of-the-art
approach for the task of RWC-Pop dataset structural segmentation is a
supervised CNN developed by Grill and Schlüter [6]. This technique also uses
signal self-similarity as input, but in the form of lag matrices, and in
conjunction with Log Mel spectrograms. Notably, this technique is trained on
the SALAMI dataset only.

We compared the song-dependent autoencoders described in the present work with
the aforementioned techniques. Results for [2] and [3] were computed with the
MSAF toolbox [21], and boundaries were aligned to the closest downbeat (as in
[5]). Results for [6] were taken from the 2015 MIREX campaign [20].

###  5.2 Practical considerations

Networks are developed with Pytorch 1.8.0 [22]. They are optimized with the
Adam optimizer [23], with a learning rate of 0.001, divided by 10 when the
loss function reaches a plateau (20 iterations without improvement) until
1e-5. The optimization stops if no progress is made during 100 consecutive
epochs, or after a total of 1000 epochs. Each song is processed in 8-size
mini-batches. All networks are initialized with the uniform distribution
defined in [24], also known as “kaiming” initialization. Bars were estimated
with madmom [15], and spectrograms computed with librosa [25]. All
segmentation scores were computed with mir_eval [26].

The entire code for this work is open-source, and contains experimental
notebooks for
reproducibility222https://gitlab.inria.fr/amarmore/musicae/-/tree/v0.1.0.
Performing 1000 epochs for a song takes between 1.5 minute (for chromagrams)
and 5 minutes (for Mel/Log Mel spectrograms).

Fig. 4: F-measure with 0.5 seconds tolerance for AE with different features
and latent space dimensions $d_{ls}$.

###  5.3 Results

Figure 4 presents results on the RWC-Pop dataset, when varying the dimension
of the latent space $d_{ls}$. The Log Mel spectrogram clearly outperfoms the
other features in this task, and achieves the state-of-the-art level of
performance with 3-seconds tolerance, as presented in Table 1.

Table 2 presents the results of autoencoding this Log Mel representation on
the SALAMI dataset [10]. More precisely, we used the test subset of the 2015
MIREX competition, using the upper-level annotations, and keeping for each
song the annotation resulting in the best $F_{3}$ among the two available. The
latent space dimension was trained on the rest of the dataset, resulting in
$d_{ls}=32$. Even though they do not reach the state of the art level, results
on SALAMI still appear competitive with a tolerance of 3 seconds.

Results on both datasets are still improvable with 0.5s tolerance window. A
wrong estimation at 0.5s can be due to an incorrect frontier estimation, but
also to an incorrect bar estimation. This could be particularly the case for
the SALAMI dataset, where music is available in very different genres and
recording conditions, such as live recordings.

Table 2: Results on the SALAMI dataset (test subset of MIREX 2015), for the Log Mel spectrogram. Note: this subset contains song on which the CNN has been trained [6]. | $F_{0.5}$ | $F_{3}$
---|---|---
Latent representation autosimilarity | 36.6% | 59.3%
Barwise feature autosimilarity | 32% | 53.9%
CNN [6] (taken from MIREX [20]) | 54.1% | 62.3%

##  6 Conclusion

Unsupervised AEs appear as competitive schemes for Music Structure Analysis.
In our experiments, the Log Mel feature outperforms the other ones, and
reaches the state-of-the-art at 3 seconds tolerance on the RWC-Pop music
dataset.

Future work will focus on improving the proposed paradigm using strategies
such as transfer learning from a song-independent AE or exploring various
network architectures. Finally, while the convolutional dynamic programming
algorithm is competitive, it is unstable to small variations in
autosimilarities and should be made more robust.

Altogether, we believe that these first results pave the way to an interesting
paradigm using single-song AEs for the description of structural elements in
music.

## References

  * [1] O. Nieto and al.,  “Audio-based music structure analysis: Current trends, open challenges, and applications,”  Trans. Int. Soc. for Music Information Retrieval, vol. 3, no. 1, 2020.
  * [2] J. Foote,  “Automatic audio segmentation using a measure of audio novelty,”  in 2000 IEEE Int. Conf. Multimedia and Expo. ICME2000. Proc. Latest Advances in the Fast Changing World of Multimedia. IEEE, 2000, vol. 1, pp. 452–455.
  * [3] B. McFee and D. Ellis,  “Analyzing song structure with spectral clustering,”  in Int. Soc. Music Information Retrieval (ISMIR), 2014, pp. 405–410.
  * [4] O. Nieto and T. Jehan,  “Convex non-negative matrix factorization for automatic music structure identification,”  in 2013 IEEE Int. Conf. Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2013, pp. 236–240.
  * [5] A. Marmoret, J.E. Cohen, N. Bertin, and F. Bimbot,  “Uncovering audio patterns in music with nonnegative tucker decomposition for structural segmentation,”  in ISMIR, 2020, pp. 788–794.
  * [6] T. Grill and J. Schlüter,  “Music boundary detection using neural networks on combined features and two-level annotations,”  in ISMIR, 2015, pp. 531–537.
  * [7] M.C. McCallum,  “Unsupervised learning of deep features for music segmentation,”  in ICASSP. IEEE, 2019, pp. 346–350.
  * [8] E.J. Humphrey, J.P. Bello, and Y. LeCun,  “Feature learning and deep architectures: New directions for music informatics,”  J. Intelligent Information Systems, vol. 41, no. 3, pp. 461–481, 2013.
  * [9] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka,  “RWC Music Database: Popular, Classical and Jazz Music Databases,”  in ISMIR, 2002, vol. 2, pp. 287–288.
  * [10] J.B.L. Smith, J.A. Burgoyne, I. Fujinaga, D. De Roure, and J.S. Downie,  “Design and creation of a large-scale database of structural annotations.,”  in ISMIR. Miami, FL, 2011, vol. 11, pp. 555–560.
  * [11] J. Engel and al.,  “Neural audio synthesis of musical notes with wavenet autoencoders,”  in Int. Conf. Machine Learning. PMLR, 2017, pp. 1068–1077.
  * [12] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner,  “Gradient-based learning applied to document recognition,”  Proc. of the IEEE, vol. 86, no. 11, pp. 2278–2324, 1998.
  * [13] G. Peeters and G. Richard,  “Deep Learning for Audio and Music,”  in Multi-faceted Deep Learning: Models and Data. Springer, 2021\.
  * [14] V. Dumoulin and F. Visin,  “A guide to convolution arithmetic for deep learning,”  arXiv preprint arXiv:1603.07285, 2016.
  * [15] S. Böck, F. Korzeniowski, J. Schlüter, F. Krebs, and G. Widmer,  “madmom: a new Python Audio and Music Signal Processing Library,”  in Proc. 24th ACM Int. Conf. Multimedia, Amsterdam, The Netherlands, 10 2016, pp. 1174–1178.
  * [16] B. McFee and D. Ellis,  “Learning to segment songs with ordinal linear discriminant analysis,”  in ICASSP. IEEE, 2014, pp. 5197–5201.
  * [17] J. Paulus, M. Müller, and A. Klapuri,  “Audio-based music structure analysis.,”  in ISMIR. Utrecht, 2010, pp. 625–636.
  * [18] G. Sargent, F. Bimbot, and E. Vincent,  “Estimating the structural segmentation of popular music pieces under regularity constraints,”  IEEE/ACM Trans. Audio, Speech, and Language Processing, vol. 25, no. 2, pp. 344–358, 2016.
  * [19] F. Bimbot, G. Sargent, E. Deruty, C. Guichaoua, and E. Vincent,  “Semiotic description of music structure: An introduction to the quaero/metiss structural annotations,”  in 53rd Int. Conf. Audio Engineering Soc., 2014.
  * [20] “MIREX wiki,” https://www.music-ir.org/mirex/wiki/MIREX_HOME,  (Accessed: Sep. 29, 2021).
  * [21] O. Nieto and J.P. Bello,  “Systematic exploration of computational music structure research.,”  in ISMIR, 2016, pp. 547–553.
  * [22] A. Paszke and al.,  “Pytorch: An imperative style, high-performance deep learning library,”  in Advances in Neural Information Processing Systems 32, pp. 8024–8035. Curran Associates, Inc., 2019.
  * [23] D.P. Kingma and J. Ba,  “Adam: A method for stochastic optimization,”  arXiv preprint arXiv:1412.6980, 2014.
  * [24] K. He, X. Zhang, S. Ren, and J. Sun,  “Delving deep into rectifiers: Surpassing human-level performance on imagenet classification,”  in Proc. IEEE Int. Conf. Computer Vision, 2015, pp. 1026–1034.
  * [25] B. McFee and al.,  “librosa/librosa: 0.8.1rc2,” May 2021.
  * [26] C. Raffel and al.,  “mir_eval: A transparent implementation of common MIR metrics,”  in ISMIR, 2014.