∎

11institutetext: Qifeng Zhang 22institutetext: Department of Mathematics,
Zhejiang Sci-Tech University, Hangzhou, 310018, China
22email<EMAIL_ADDRESS>33institutetext: Jiyuan Zhang
44institutetext: Department of Mathematics, Zhejiang Sci-Tech University,
Hangzhou, 310018, China
44email<EMAIL_ADDRESS>55institutetext: Zhimin Zhang 66institutetext:
Department of Mathematics, Wayne State University, Detroit, Michigan 48202,
USA
66email<EMAIL_ADDRESS>
# Error estimates of invariant-preserving difference schemes for the rotation-
two-component Camassa–Holm system with small energy ††thanks: Qifeng Zhang was
supported by Zhejiang Provincial Natural Science Foundation of China (Grant
No. LZ23A010007). Zhimin Zhang was supported by the National Natural Science
Foundation of China (Grant Nos. 12131005 and 11871092).

Qifeng Zhang  Jiyuan Zhang  Zhimin Zhang

(Received: date / Accepted: date)

###### Abstract

A rotation-two-component Camassa-Holm (R2CH) system was proposed recently to
describe the motion of shallow water waves under the influence of gravity.
This is a highly nonlinear and strongly coupled system of partial differential
equations. A crucial issue in designing numerical schemes is to preserve
invariants as many as possible at the discrete level. In this paper, we
present a provable implicit nonlinear difference scheme which preserves at
least three discrete conservation invariants: energy, mass, and momentum, and
prove the existence of the difference solution via the Browder theorem. The
error analysis is based on novel and refined estimates of the bilinear
operator in the difference scheme. By skillfully using the energy method, we
prove that the difference scheme not only converges unconditionally when the
rotational parameter diminishes, but also converges without any step-ratio
restriction for the small energy case when the rotational parameter is
nonzero. The convergence orders in both settings (zero or nonzero rotation
parameter) are $O(\tau^{2}+h^{2})$ for the velocity in the $L^{\infty}$-norm
and the surface elevation in the $L^{2}$-norm, where $\tau$ denotes the
temporal stepsize and $h$ the spatial stepsize, respectively. The theoretical
predictions are confirmed by a properly designed two-level iteration scheme.
Comparing with existing numerical methods in the literature, the proposed
method demonstrates its effectiveness for long-time simulation over larger
domains and superior resolution for both smooth and non-smooth initial values.

###### Keywords:

R2CH system; Invariants; Error estimate; Small energy; Long time simulation

###### MSC:

65M06 65M12 26A33 35R11

††journal:

##  1 Introduction

In this paper we propose, analyze and test a two-level invariant-preserving
nonlinear difference scheme for solving a rotation-two-component Camassa-Holm
(R2CH) system of the form

$\displaystyle u_{t}-u_{xxt}-\kappa
u_{x}+3uu_{x}=\sigma(2u_{x}u_{xx}+uu_{xxx})-\mu
u_{xxx}-(1-2\Omega\kappa)\rho\rho_{x}+2\Omega\rho(\rho u)_{x},$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad
x\in\mathbb{R},\;t\in[0,T],$ (1.1a) $\displaystyle\rho_{t}+(\rho
u)_{x}=0,\quad x\in\mathbb{R},\;t\in[0,T],$ (1.1b)

subject to the initial value conditions

$\displaystyle u(x,0)=u^{0}(x),\quad\rho(x,0)=\rho^{0}(x),\quad
x\in\mathbb{R},$ (1.2)

and periodic boundary value conditions

$\displaystyle u(x,t)=u(x+L,t),\quad\rho(x,t)=\rho(x+L,t),\quad
x\in\mathbb{R},\;t\in[0,T].$ (1.3)

The system (1.1)–(1.3) was introduced by Fan, Gao, and Liu FGL2016 to depict
the motion of shallow water waves at a free surface involving the Coriolis
force caused by the Earth’s rotation in the equatorial ocean regions. The
variable $u(x,t)$ in (1.1) denotes the horizontal fluid velocity along with
the $x$-direction, and the variable $\rho(x,t)$ is the altitude from the free
surface elevation to equilibrium. The coefficient $\kappa$ denotes a
underlying linear shear flow, and the parameter $\sigma>0$ provides the
competition/balance index. $\mu$ is a real dimensionless constant, and
$\Omega$ defines the average rotational angular velocity of the Earth.
Throughout the paper, we always assume $\Omega\in[0,1/4)$, and
$1-2\Omega\kappa>0$, see e.g., CFGL2017 .

It can be easily shown that at least three invariants for the R2CH system
(1.1) can be expressed: (see e.g., FGL2016 )

  * •

Energy:

$\displaystyle
E(u,\rho)=\frac{1}{2}\int_{\mathbb{R}}\big{[}u^{2}+u_{x}^{2}+(1-2\Omega\kappa)(\rho-1)^{2}\big{]}{\rm
d}x.$ (1.4)

  * •

Momentum:

$\displaystyle
H(u,\rho)=\int_{\mathbb{R}}\big{[}u+\Omega(\rho-1)^{2}\big{]}{\rm d}x.$ (1.5)

  * •

Mass:

$\displaystyle I(u,\rho)=\int_{\mathbb{R}}(\rho-1){\rm d}x.$ (1.6)

As of now, there have been intense theoretical studies on the R2CH system
(1.1) such as the existence of global solutions Zhang2017 ; Moon2017 , the
wave-breaking of solitary waves CFGL2017 ; LPZ2019 , the local well-posedness
ZL2018 , the peakon-delta weak solutions FY2019 , the persistence properties
YLQ2020 , the traveling waves ZW2021 , and so forth. Once the Earth’s rotation
parameter equals zero, the system (1.1) degenerates to a two-component
Camassa–Holm (2CH) system, which was introduced in OR1996 for the first time,
and was later re-derived in CI2008 . The extensive analyses with respect to
the theoretical solutions are referred to CL2011 ; GL2011 ; GY2011 ; Henry2009
. Furthermore, when the variable $\rho(x,t)$ disappears, (1.1) reduces to the
classical CH equation, which was first proposed as a bi-Hamiltonian system by
Fuchssteiner et al. in FF1981 . Subsequently, it was reformulated as a
physical model to simulate the unidirectional propagation of shallow water
waves by Camassa and Holm in CH1993 . Its theoretical studies have been very
extensive, see e.g., BC2007 ; Dan2003 ; HX2008 . In addition, some researchers
considered more general three-component, four-component, and multi-component
CH systems, see e.g., GX2011 ; KLQ2021 ; LLP2014 ; HI2010 . Several
significant research timelines are summarized in Figure 1.

Figure 1: Milestones of the research timeline

In contrast with extensive and in-depth theoretical studies, most of the
present numerical studies are limited to the special cases of the R2CH system.
For example, with regard to the classical CH equation, there are plentiful
papers devoted to investigating the numerical behavior of weak solutions and
peakons, see e.g., HR2006a ; HR2006b ; CKR2008 ; COR2008 ; CR2012 ; XS2008 ;
LX2016 . Other important advances include the particle method CL2008 , multi-
symplectic/multi-symplectic wavelet collocation schemes DBX2008 ; ZST2011 ,
the spectral/spectral element/spectral projection methods WX2015 ; YYW2021 ;
KR2006 , the Galerkin finite element method ADM2019 , and the multiple scalar
auxiliary variables/the invariant energy quadratization methods JGCW2020 ;
JWG2021 . In the meantime, for the 2CH system, Cohen et al. investigated a
multi-symplectic structure-preserving scheme CMR2014 , Liu et al. developed an
invariant-preserving difference scheme LP2016 , Li et al. proposed a multi-
symplectic compact method LQZS2017 , Yu et al. presented a three-step
iterative algorithm YFS2018 , Chertock et al. applied a finite volume particle
method CKL2020 , and Galtung et al. studied a variational discretization in
Lagrangian variables GG2021 . Almost for sure, there is in no way an
exhaustive list, but these results indicate the trend and efforts in the
development of invariant-preserving numerical methods for the CH equation and
2CH system.

So far as is known to us, the numerical methods for the R2CH system are still
lacking. The first numerical attempt for the R2CH system is attributed to
Zhang et al. ZLZ2022 , who designed a conditionally convergent, and
energy/mass-preserving finite difference scheme. However, there are several
drawbacks for the scheme including: (i) It is difficult to determine whether
the difference scheme preserves the momentum; (ii) When $\Omega=0$, a
convergent condition $\tau\lesssim h$ in the error estimate is indispensable;
(iii) When $\Omega\neq 0$, the error estimate is missing; (iv) Numerical
simulation will distort over a long-time simulation on large domain unless
utilizing a more refined grid.

The original motivation in the present paper is to construct numerical
schemes, which are able to fully preserve invariants (1.4)–(1.6) in the R2CH
system (1.1). Other motivations include overcoming the difficulties listed
above and remedying deficiencies during the numerical analysis. Nevertheless,
the R2CH system is very complicated, and developing robust and accurate
numerical schemes is a highly nontrivial task.

To achieve these goals, we have already learned that the invariant-preserving
property is an important index to judge the merits of an algorithm, see e.g.,
FQ1991 ; HPR2022 ; LV1995 . Therefore, we first ensure that the constructed
difference scheme preserves as many invariants as possible. In addition, the
fully implicit nonlinear numerical schemes usually have a better invariant-
preserving ability and long-time numerical simulation capability. As a
consequence, we start to derive the implicit time discretization, and then
establish a class of fully implicit nonlinear numerical schemes for the R2CH
system. Unfortunately, nonlinear difference schemes are usually
computationally expensive in solving nonlinear system. As a possible
compromise to balance computational effort and to preserve all invariants, we
propose a well-designed two-level linearized iterative method to solve the
nonlinear difference scheme of the R2CH system, see (5.42) in Section 5.

Additionally, numerical analysis for most of the shallow water waves is often
intractable due to complex high-order nonlinear derivative terms and coupled
nonlinear terms. As a result, much previous research work focuses on the
numerical schemes and numerical simulation, but ignore rigorous convergence
analysis, see e.g., COR2008 ; CKL2020 ; GG2021 ; LP2016 ; LX2016 . Moreover,
convergence analysis for nonlinear difference schemes is usually tougher
compared with linear ones since their error systems would encounter more
unknown information. To overcome these difficulties appeared in the nonlinear
scheme for the R2CH system, by means of a more refined observation on the
bilinear operators, and technical energy analysis, we find that some
troublesome terms in the error system can be coincidentally cancelled out by
taking a new inner product. At the same time, we find a new recursion
relation, which could skillfully deal with the discrete time difference
quotient. These, together with the nice properties of the bilinear operators
and several key techniques, make the desired convergence possible.

The main contribution of the present work is to analyze a coupled nonlinear
difference scheme for solving the R2CH system, and dedicated to establishing a
complete convergence theory. In particular:

  * •

The nonlinear difference scheme not only preserves the energy and the mass of
the original problem naturally, but also preserves the momentum explicitly;

  * •

When $\Omega=0$, the nonlinear difference scheme converges unconditionally
(without time step restriction), which allows large time steps in calculation;

  * •

When $\Omega\neq 0$, it is worth noting that the nonlinear difference scheme
converges unconditionally provided the initial energy is suitably small (or
the initial values are small). However, it does not need a time step
restriction and any boundedness hypothesis of numerical solutions;

  * •

The difference scheme has better resolution in long-time numerical simulation
on a large area even with a coarse grid comparing with that in CKL2020 ;
CMR2014 ; LP2016 ; YFS2018 ; ZLZ2022 under the smooth/nonsmooth initial
values.

  * •

The analyzing techniques based on the framework of the bilinear operator open
the door of the difference methods to solve other types of the shallow water
problems.

We begin in Section 2 by introducing several useful notations and lemmas. In
Section 3, we derive a fully implicit nonlinear difference scheme for the R2CH
system followed by detailed analyses for the invariant-preserving properties
and a rigorous proof for the existence of a numerical solution. Section 4 is
the main body of the paper, which focuses on the convergence analysis
involving two scenarios: $\Omega=0$ and $\Omega\neq 0$. In Section 5, we carry
out two specific numerical examples including several different cases to test
the theoretical results by designing an iterated scheme. We end the paper by
comparing the pros and cons of the present difference scheme with that in the
literature in Section 6. Some conclusions and outlooks are summarized in the
last Section.

##  2 Notations and lemmas

To define our finite difference method, we introduce some notations. Once the
positive integers $M$ and $N$ are selected, we let $h=L/M$ and $\tau=T/N$.
Denote $x_{i}=ih,\;i\in\mathbb{Z}$, $t_{n}=n\tau,\;0\leqslant n\leqslant N$,
$\Omega_{h}=\\{x_{i}\,|\,x_{i}=ih,\;i\in\mathbb{Z}\\}$,
$\Omega_{\tau}=\\{t_{n}\,|\,t_{n}=n\tau,\;0\leqslant n\leqslant N\\}$,
$\Omega_{h\tau}=\Omega_{h}\times\Omega_{\tau}$. For arbitrary grid functions
$v=\\{v_{i}^{n}\,|\,i\in\mathbb{Z},\;0\leqslant n\leqslant N\\}$, and
$u=\\{u_{i}^{n}\,|\,i\in\mathbb{Z},\;0\leqslant n\leqslant N\\}$ defined on
$\Omega_{h\tau}$, we bring in the notations below:

$\displaystyle
v_{i}^{n+\frac{1}{2}}=\frac{1}{2}(v_{i}^{n}+v_{i}^{n+1}),\quad\delta_{t}v_{i}^{n+\frac{1}{2}}=\frac{1}{\tau}(v_{i}^{n+1}-v_{i}^{n}),\quad\delta_{x}v_{i-\frac{1}{2}}^{n}=\frac{1}{h}(v_{i}^{n}-v_{i-1}^{n}),$
$\displaystyle\delta_{x}^{2}v_{i}^{n}=\frac{1}{h}(\delta_{x}v_{i+\frac{1}{2}}-\delta_{x}v_{i-\frac{1}{2}}),\quad\Delta_{x}v_{i}^{n}=\frac{1}{2h}(v_{i+1}^{n}-v_{i-1}^{n}),\quad(uv)_{i}^{n+\frac{1}{2}}=\frac{u_{i}^{n}v_{i}^{n}+u_{i}^{n+1}v_{i}^{n+1}}{2}.$

The grid function spaces in space and time are, respectively, denoted by

$\mathcal{V}_{h}=\big{\\{}v\,|\,v=\\{v_{i}\\},v_{i+M}=v_{i},\;i\in\mathbb{Z}\big{\\}},$

and

$\mathcal{S}_{\tau}=\big{\\{}w\,|\,w=(w^{0},w^{1},\cdots,w^{N})~{}{\rm
is~{}the~{}grid~{}function~{}defined~{}on}~{}\Omega_{\tau}\big{\\}}.$

For the integers $i$ and $n$, define the index sets
$\mathbb{I}_{M}=\\{i\,|\,1\leqslant i\leqslant M\\}$,
$\mathbb{K}_{N}=\\{n\,|\,0\leqslant n\leqslant N\\}$ and
$\mathbb{K}_{N}^{0}=\\{n\,|\,0\leqslant n\leqslant N-1\\}$. For arbitrary grid
functions $u$, $v\in\mathcal{V}_{h}$, we define the discrete $L^{2}$-inner
products by

$\displaystyle(u,v)=h\sum_{i\in\mathbb{I}_{M}}u_{i}v_{i},\quad\langle\delta_{x}u,\delta_{x}v\rangle=h\sum_{i\in\mathbb{I}_{M}}(\delta_{x}u_{{i-\frac{1}{2}}})(\delta_{x}v_{{i-\frac{1}{2}}}),$

and the corresponding norms (seminorm) by

$\displaystyle\|v\|=\sqrt{(v,v)},\quad|v|_{1}=\sqrt{\langle\delta_{x}v,\delta_{x}v\rangle},\quad\|v\|_{\infty}=\max_{i\in\mathbb{I}_{M}}\limits|v_{i}|.$

Furthermore, denote $(uv)_{i}=u_{i}v_{i}$ and define a bilinear operator (see
e.g., Guo1974 ; Guo1981 ) to facilitate numerical analysis later as follows.

$\displaystyle\psi(u,v)_{i}=\frac{1}{3}[u_{i}\Delta_{x}v_{i}+\Delta_{x}(uv)_{i}],\quad
i\in\mathbb{I}_{M}.$ (2.1)

###### Lemma 2.1

ZL2020 For two arbitrary spatial grid functions $u,\,v\in\mathcal{V}_{h}$, we
have

$\displaystyle(\psi(u,v),v)=0,\quad(\Delta_{x}u,u)=0,\quad(\Delta_{x}u,v)=-(u,\Delta_{x}v),\quad(\delta_{x}^{2}u,v)=-\langle\delta_{x}u,\delta_{x}v\rangle.$

The discrete form of the classical embedded inequality is introduced as
follows, see also ZL2020 . The continuous edition is similar to the discrete
one, which is omitted here for brevity.

###### Lemma 2.2

For an arbitrary spatial grid function $v\in\mathcal{{V}}_{h}$ and
$\varepsilon>0$, we have

$\|v\|_{\infty}^{2}\leqslant\varepsilon|v|_{1}^{2}+\Big{(}\frac{1}{\varepsilon}+\frac{1}{L}\Big{)}\|v\|^{2}.$

In order to simplify the numerical analysis of the temporal direction, we give
the following lemma.

###### Lemma 2.3

For two arbitrary temporal grid functions $u,\;v\in\mathcal{{S}}_{\tau}$, we
have

$\displaystyle\;(\delta_{t}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}v^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;\frac{1}{2\tau}[(u^{n+1},u^{n+1}v^{n+1})-(u^{n},u^{n}v^{n})]$
$\displaystyle\;-\frac{1}{4}(u^{n+1}-u^{n},(u^{n+1}-u^{n})\delta_{t}v^{n+\frac{1}{2}})-\frac{1}{2}(u^{n+1}u^{n},\delta_{t}v^{n+\frac{1}{2}}).$

###### Proof

We easily have

$\displaystyle\;\quad\;(\delta_{t}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}v^{n+\frac{1}{2}})$
$\displaystyle\;=(\delta_{t}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}v^{n+\frac{1}{2}}-(uv)^{n+\frac{1}{2}}+(uv)^{n+\frac{1}{2}})$
$\displaystyle\;=(\delta_{t}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}v^{n+\frac{1}{2}}-(uv)^{n+\frac{1}{2}})+(\delta_{t}u^{n+\frac{1}{2}},(uv)^{n+\frac{1}{2}})$
$\displaystyle\;=:A^{n}+B^{n}.$ (2.2)

We calculate each term in (2.2).

$\displaystyle A^{n}=$
$\displaystyle\;(\delta_{t}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}v^{n+\frac{1}{2}}-(uv)^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;\frac{1}{4}(\delta_{t}u^{n+\frac{1}{2}},(u^{n}+u^{n+1})(v^{n+1}+v^{n})-2(u^{n+1}v^{n+1}+u^{n}v^{n}))$
$\displaystyle=$
$\displaystyle\;\frac{1}{4}(\delta_{t}u^{n+\frac{1}{2}},u^{n}v^{n+1}-u^{n}v^{n}+u^{n+1}v^{n}-u^{n+1}v^{n+1})$
$\displaystyle=$
$\displaystyle\;\frac{1}{4\tau}(u^{n+1}-u^{n},u^{n}v^{n+1}-u^{n}v^{n}+u^{n+1}v^{n}-u^{n+1}v^{n+1})$
$\displaystyle=$
$\displaystyle\;\frac{1}{4}(u^{n+1}-u^{n},u^{n}\delta_{t}v^{n+\frac{1}{2}}-u^{n+1}\delta_{t}v^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;-\frac{1}{4}(u^{n+1}-u^{n},(u^{n+1}-u^{n})\delta_{t}v^{n+\frac{1}{2}}),$
(2.3)

and

$\displaystyle B^{n}=$
$\displaystyle\;(\delta_{t}u^{n+\frac{1}{2}},(uv)^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;\frac{1}{2\tau}(u^{n+1}-u^{n},u^{n+1}v^{n+1}+u^{n}v^{n})$
$\displaystyle=$
$\displaystyle\;\frac{1}{2\tau}[(u^{n+1},u^{n+1}v^{n+1})-(u^{n},u^{n}v^{n})]+\frac{1}{2\tau}[(u^{n+1},u^{n}v^{n})-(u^{n},u^{n+1}v^{n+1})]$
$\displaystyle=$
$\displaystyle\;\frac{1}{2\tau}[(u^{n+1},u^{n+1}v^{n+1})-(u^{n},u^{n}v^{n})]-\frac{1}{2}(u^{n}u^{n+1},\delta_{t}v^{n+\frac{1}{2}}).$
(2.4)

Plugging (2.3) and (2.4) into (2.2), it completes the proof.

##  3 Numerical analysis

Assume the exact solutions to the problem (1.1)–(1.3) satisfy

$u(x,t)\in\mathbb{C}^{5,3}(\mathbb{R}\times[0,T]),\quad\rho(x,t)\in\mathbb{C}^{3,3}(\mathbb{R}\times[0,T])$

and denote

$\displaystyle c_{0}=\max_{0\leqslant x\leqslant L,0\leqslant t\leqslant
T}\\{|u_{x}(x,t)|,\,|u_{xx}(x,t)|,\,|u_{t}(x,t)|,\,|\rho(x,t)|,\,|\rho_{x}(x,t)|\\},$
(3.1) $\displaystyle c_{\rm max}=\max_{0\leqslant x\leqslant L,0\leqslant
t\leqslant T}|u(x,t)|.$ (3.2)

###  3.1 The derivation of the difference scheme

Denote the exact solutions at the grid functions in short by

$\displaystyle U_{i}^{n}=u(x_{i},t_{n}),\;\Pi_{i}^{n}=\rho(x_{i},t_{n}),\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}.$

Considering (1.1) directly at the virtual point $(x_{i},t_{n+\frac{1}{2}})$,
and combining Taylor expansion, we have

$\displaystyle\delta_{t}U_{i}^{n+\frac{1}{2}}-\delta_{t}\delta_{x}^{2}U_{i}^{n+\frac{1}{2}}-\kappa\Delta_{x}U_{i}^{n+\frac{1}{2}}+3\psi(U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})_{i}$
$\displaystyle=3\sigma\psi(\delta_{x}^{2}U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})_{i}-\mu\Delta_{x}\delta_{x}^{2}U_{i}^{n+\frac{1}{2}}-(1-2\Omega\kappa)\Pi_{i}^{n+\frac{1}{2}}\Delta_{x}\Pi_{i}^{n+\frac{1}{2}}$
$\displaystyle\quad+2\Omega\Pi_{i}^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})_{i}+Q_{i}^{n+\frac{1}{2}},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0},$ (3.3a)
$\displaystyle\delta_{t}\Pi_{i}^{n+\frac{1}{2}}+\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})_{i}=R_{i}^{n+\frac{1}{2}},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0},$ (3.3b)

there is a constant $c_{1}$ independent of $\tau$ and $h$ such that

$\displaystyle|Q_{i}^{n+\frac{1}{2}}|\leqslant c_{1}(\tau^{2}+h^{2}),\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0},$ (3.4a)
$\displaystyle|R_{i}^{n+\frac{1}{2}}|\leqslant c_{1}(\tau^{2}+h^{2}),\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0},$ (3.4b)

Taking the initial and boundary value conditions (1.2)–(1.3) into account, we
have

$\displaystyle U_{i}^{0}=u^{0}(x_{i}),\;\Pi_{i}^{0}=\rho^{0}(x_{i}),\quad
i\in\mathbb{I}_{M},$ (3.5a) $\displaystyle
U_{i}^{n}=U_{i+M}^{n},\;\Pi_{i}^{n}=\Pi_{i+M}^{n},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}.$ (3.5b)

Negelcting the local truncation errors (3.4), replacing the functions
$U_{i}^{n}$ and $\Pi_{i}^{n}$ with their numerical approximations $u_{i}^{n}$
and $\rho_{i}^{n}$ in (3.3) and (3.5), respectively, a two-level coupled
nonlinear difference scheme reads

$\displaystyle\delta_{t}u_{i}^{n+\frac{1}{2}}-\delta_{t}\delta_{x}^{2}u_{i}^{n+\frac{1}{2}}-\kappa\Delta_{x}u_{i}^{n+\frac{1}{2}}+3\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})_{i}$
$\displaystyle=3\sigma\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})_{i}-\mu\Delta_{x}\delta_{x}^{2}u_{i}^{n+\frac{1}{2}}-(1-2\Omega\kappa)\rho_{i}^{n+\frac{1}{2}}\Delta_{x}\rho_{i}^{n+\frac{1}{2}}$
$\displaystyle\quad+2\Omega\rho_{i}^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0},$ (3.6a)
$\displaystyle\delta_{t}\rho_{i}^{n+\frac{1}{2}}+\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i}=0,\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0}.$ (3.6b) $\displaystyle
u_{i}^{0}=u^{0}(x_{i}),\;\rho_{i}^{0}=\rho^{0}(x_{i}),\quad
i\in\mathbb{I}_{M},$ (3.6c) $\displaystyle
u_{i}^{n}=u_{i+M}^{n},\;\rho_{i}^{n}=\rho_{i+M}^{n},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}.$ (3.6d)

###### Remark 3.1

Note that the scheme (3.6) is exactly the standard Crank-Nicolson scheme once
we replace $3\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})_{i}$ by
$u_{i}^{n+\frac{1}{2}}\Delta_{x}u_{i}^{n+\frac{1}{2}}$ and
$3\sigma\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})_{i}$ by
$2\Delta_{x}u_{i}^{n+\frac{1}{2}}\delta_{x}^{2}u_{i}^{n+\frac{1}{2}}+u_{i}^{n+\frac{1}{2}}\Delta_{x}\delta_{x}^{2}u_{i}^{n+\frac{1}{2}}$,
respectively. In this case, the numerical analysis is challengeable and
deserves further study.

###  3.2 Conservative invariants and boundedness

###### Theorem 3.1

Suppose
$\\{u_{i}^{n},\,\rho_{i}^{n}\,|\,i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}\\}$
is the numerical solution of the scheme (3.6). Then for any
$n\in\mathbb{K}_{N}$, we have the discrete invariants:

  * •

Energy:

$\displaystyle E^{n}=E^{0},\quad{\rm where}\quad
E^{n}=\frac{1}{2}\big{[}\|u^{n}\|^{2}+|u^{n}|_{1}^{2}+(1-2\Omega\kappa)\|\rho^{n}\|^{2}\big{]};$
(3.7)

  * •

Momentum:

$\displaystyle H^{n}=H^{0},\quad{\rm where}\quad
H^{n}=(u^{n},1)+\Omega\|\rho^{n}\|^{2};$ (3.8)

  * •

Mass:

$\displaystyle I^{n}=I^{0},\quad{\rm where}\quad I^{n}=(\rho^{n},1).$ (3.9)

###### Proof

(Energy): Taking the $L^{2}$-inner product of (3.6a) with $u^{n+\frac{1}{2}}$,
we readily have

$\displaystyle(\delta_{t}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})-(\delta_{t}\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})-\kappa(\Delta_{x}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})+3(\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),u^{n+\frac{1}{2}})$
$\displaystyle=\;$ $\displaystyle
3\sigma(\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),u^{n+\frac{1}{2}})-\mu(\Delta_{x}\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})-(1-2\Omega\kappa)\cdot(\rho^{n+\frac{1}{2}}\Delta_{x}\rho^{n+\frac{1}{2}},u^{n+\frac{1}{2}})$
$\displaystyle+2\Omega(\rho^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),u^{n+\frac{1}{2}}),\quad
n\in\mathbb{K}_{N}^{0}.$

Using Lemma 2.1, for $n\in\mathbb{K}_{N}^{0}$, the above equality deduces to

$\displaystyle\frac{1}{2\tau}(\|u^{n+1}\|^{2}-\|u^{n}\|^{2})+\frac{1}{2\tau}(|u^{n+1}|_{1}^{2}-|u^{n}|_{1}^{2})=(1-2\Omega\kappa)\cdot(\rho^{n+\frac{1}{2}},\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})).$
(3.10)

Again taking the $L^{2}$-inner product of (3.6b) with $\rho^{n+\frac{1}{2}}$,
we have

$\displaystyle(\delta_{t}\rho^{n+\frac{1}{2}},\rho^{n+\frac{1}{2}})+(\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),\rho^{n+\frac{1}{2}})=0,\quad
n\in\mathbb{K}_{N}^{0}.$ (3.11)

Multiplying (3.11) by $1-2\Omega\kappa$ and summing up with (3.10), we have

$\displaystyle\frac{1}{2\tau}\big{[}(\|u^{n+1}\|^{2}-\|u^{n}\|^{2})+(|u^{n+1}|_{1}^{2}-|u^{n}|_{1}^{2})+(1-2\Omega\kappa)\cdot(\|\rho^{n+1}\|^{2}-\|\rho^{n}\|^{2})\big{]}=0,\quad
n\in\mathbb{K}_{N}^{0}.$

Rearranging the equality above, we have

$\displaystyle\frac{1}{2}\big{[}\|u^{n+1}\|^{2}+|u^{n+1}|_{1}^{2}+(1-2\Omega\kappa)\cdot\|\rho^{n+1}\|^{2}\big{]}=\frac{1}{2}\big{[}\|u^{n}\|^{2}+|u^{n}|_{1}^{2}+(1-2\Omega\kappa)\cdot\|\rho^{n}\|^{2}\big{]},\quad
n\in\mathbb{K}_{N}^{0},$

which implies (3.7).
(Momentum): Taking the $L^{2}$-inner product of (3.6a) with $1$, we have

$\displaystyle(\delta_{t}u^{n+\frac{1}{2}},1)-(\delta_{t}\delta_{x}^{2}u^{n+\frac{1}{2}},1)-\kappa(\Delta_{x}u^{n+\frac{1}{2}},1)+3(\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),1)$
$\displaystyle=\;$ $\displaystyle
3\sigma(\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),1)-\mu(\Delta_{x}\delta_{x}^{2}u^{n+\frac{1}{2}},1)-(1-2\Omega\kappa)\cdot(\rho^{n+\frac{1}{2}}\Delta_{x}\rho^{n+\frac{1}{2}},1)$
$\displaystyle+2\Omega(\rho^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),1),\quad
n\in\mathbb{K}_{N}^{0}.$

Noticing Lemma 2.1, the equality above becomes

$\displaystyle\frac{1}{\tau}\left[(u^{n+1},1)-(u^{n},1)\right]-2\Omega(\rho^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),1)=0,\quad
n\in\mathbb{K}_{N}^{0}.$ (3.12)

Taking an inner product of (3.6b) with $2\Omega\rho^{n+\frac{1}{2}}$, we have

$\displaystyle
2\Omega(\delta_{t}\rho^{n+\frac{1}{2}},\rho^{n+\frac{1}{2}})+2\Omega(\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),\rho^{n+\frac{1}{2}})=0,\quad
n\in\mathbb{K}_{N}^{0}.$ (3.13)

Adding (3.12) and (3.13) together, we have

$\displaystyle\frac{1}{\tau}\left[(u^{n+1},1)-(u^{n},1)+\Omega(\|\rho^{n+1}\|^{2}-\|\rho^{n}\|^{2})\right]=0,\quad
n\in\mathbb{K}_{N}^{0}.$

Rearranging the above equality, we have

$\displaystyle(u^{n+1},1)+\Omega\|\rho^{n+1}\|^{2}=(u^{n},1)+\Omega\|\rho^{n}\|^{2},\quad
n\in\mathbb{K}_{N}^{0}.$

By reduction, it consequently yields (3.8).
(Mass): Taking an $L^{2}$-inner product of (3.6b) directly with $1$, we have

$\displaystyle(\delta_{t}\rho^{n+\frac{1}{2}},1)+(\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),1)=0,\quad
n\in\mathbb{K}_{N}^{0}.$ (3.14)

By the summation by parts, the second term vanishes in (3.14). Thus we have
(3.9), which ends the proof.

###### Remark 3.2

By Theorem 3.1, we have

$\displaystyle\|u^{n}\|^{2}+|u^{n}|_{1}^{2}\leqslant\|u^{0}\|^{2}+|u^{0}|_{1}^{2}+(1-2\Omega\kappa)\|\rho^{0}\|^{2},\quad
n\in\mathbb{K}_{N},$ (3.15a)
$\displaystyle\|\rho^{n}\|^{2}\leqslant\frac{1}{1-2\Omega\kappa}\left(\|u^{0}\|^{2}+|u^{0}|_{1}^{2}+(1-2\Omega\kappa)\|\rho^{0}\|^{2}\right),\quad
n\in\mathbb{K}_{N}.$ (3.15b)

In combination of Lemma 2.2 with (3.15a), we further deduce that
$\|u^{n}\|_{\infty}$ is bounded.

###### Remark 3.3

According to the discrete energy (3.7), we easily have
$(\rho^{n}-1,1)=(\rho^{0}-1,1),\;n\in\mathbb{K}_{N}$, which can be viewed as a
discrete edition of the continuous conservation of mass in (1.6). Furthermore,
combining with (3.7) and (3.9), we know that the discrete energy satisfies
$\frac{1}{2}\big{[}\|u^{n}\|^{2}+|u^{n}|_{1}^{2}+(1-2\Omega\kappa)\|\rho^{n}-1\|^{2}\big{]}=\frac{1}{2}\big{[}\|u^{0}\|^{2}+|u^{0}|_{1}^{2}+(1-2\Omega\kappa)\|\rho^{n}-1\|^{2}\big{]},\;n\in\mathbb{K}_{N}$,
which recovers the energy conservation of the analytical form (1.4).
Analogously, with the help of (3.8) and (3.9), the discrete momentum
conservation becomes
$(u^{n},1)+\Omega\|\rho^{n}-1\|^{2}=(u^{n},1)+\Omega\|\rho^{n}-1\|^{2},\;n\in\mathbb{K}_{N}$.

###  3.3 Existence

The following Browder theorem is very useful for the proof of the existence.

######  Theorem 3.2 (Browder theorem AK1993 )

Suppose $(H,(\cdot,\cdot))$ is an inner product space of finite dimension.
$\|\cdot\|$ is the corresponding norm, $\Theta:H\rightarrow H$ is a continuous
map. Further assume that

$\exists\,\alpha>0,\quad\forall\,z\in
H,\quad\|z\|=\alpha,\quad(\Theta(z),z)\geqslant 0.$

Then there is an element $z^{*}\in H$ such that $\Theta(z^{*})=0$ when
$\|z^{*}\|\leqslant\alpha$.

######  Theorem 3.3 (Existence)

The difference scheme (3.6) has a solution.

###### Proof

Suppose $\\{u^{n},\,\rho^{n}\\}$ has been determined. Denote
$u_{i}=u_{i}^{n+\frac{1}{2}}$ and $\rho_{i}=\rho_{i}^{n+\frac{1}{2}}$. The
difference scheme (3.6a)–(3.6b) could be rewritten as

$\displaystyle\frac{2}{\tau}(u_{i}-u_{i}^{n})-\frac{2}{\tau}(\delta_{x}^{2}u_{i}-\delta_{x}^{2}u_{i}^{n})-\kappa\Delta_{x}u_{i}+3\psi(u,u)_{i}-3\sigma\psi(\delta_{x}^{2}u,u)_{i}$
$\displaystyle\quad+\mu\Delta_{x}\delta_{x}^{2}u_{i}+(1-2\Omega\kappa)\rho_{i}\Delta_{x}\rho_{i}-2\Omega\rho_{i}\Delta_{x}(\rho
u)_{i}=0,\quad i\in\mathbb{I}_{M},$ (3.16a)
$\displaystyle\frac{2}{\tau}(\rho_{i}-\rho_{i}^{n})+\Delta_{x}(\rho
u)_{i}=0,\quad i\in\mathbb{I}_{M}.$ (3.16b)

Suppose $u\in\mathcal{U}_{h}$, we define an operator
$\Theta:\mathcal{U}_{h}\rightarrow\mathcal{U}_{h}$ by

$\displaystyle\Theta(u)_{i}=\frac{2}{\tau}(u_{i}-u_{i}^{n})-\frac{2}{\tau}(\delta_{x}^{2}u_{i}-\delta_{x}^{2}u_{i}^{n})-\kappa\Delta_{x}u_{i}+3\psi(u,u)_{i}-3\sigma\psi(\delta_{x}^{2}u,u)_{i}$
$\displaystyle\qquad\qquad+\mu\Delta_{x}\delta_{x}^{2}u_{i}+(1-2\Omega\kappa)\rho_{i}\Delta_{x}\rho_{i}-2\Omega\rho_{i}\Delta_{x}(\rho
u)_{i},\quad i\in\mathbb{I}_{M},$ (3.17)

where $\rho$ is determined by (3.16b) when $u$ is known.

In what follows, we will show the existence of $u_{i}$. Taking a standard
$L^{2}$-inner product of $\Theta(u)$ in (3.17) with $u$, it follows

$\displaystyle(\Theta(u),u)=\frac{2}{\tau}(\|u\|^{2}-(u^{n},u))+\frac{2}{\tau}(\|\delta_{x}u\|^{2}-(\delta_{x}u^{n},\delta_{x}u))-\kappa(\delta_{x}u,u)+3(\psi(u,u),u)$
$\displaystyle\qquad\qquad\quad-3\sigma(\psi(\delta_{x}^{2}u,u),u)+\mu(\Delta_{x}\delta_{x}^{2}u,u)+(1-2\Omega\kappa)(\rho\Delta_{x}\rho,u)-2\Omega(\rho(\rho
u),u).$ (3.18)

With the aid of Lemma 2.1, it yields

$\displaystyle(\delta_{x}u,u)=0,\quad(\psi(u,u),u)=0,\quad(\psi(\delta_{x}^{2}u,u),u)=0,\quad(\Delta_{x}\delta_{x}^{2}u,u)=0.$

Using the summation by parts, we have

$\displaystyle(\rho\Delta_{x}\rho,u)$ $\displaystyle=(\Delta_{x}\rho,\rho u)$
$\displaystyle=-(\rho,\Delta_{x}(\rho
u))\overset{\eqref{Exeq1.2}}{=}\Big{(}\rho,\frac{2}{\tau}(\rho-\rho^{n})\Big{)}$
$\displaystyle=\frac{2}{\tau}(\|\rho\|^{2}-(\rho,\rho^{n}))$
$\displaystyle\geqslant\frac{2}{\tau}(-\frac{1}{4}\|\rho^{n}\|^{2})$

and

$\displaystyle(\rho\Delta_{x}(\rho u),u)=(\Delta_{x}(\rho u),\rho u)=0.$

In addition,

$\displaystyle\|u\|^{2}-(u,u^{n})+\|\delta_{x}u\|^{2}-(\delta_{x}u,\delta_{x}u^{n})$
$\displaystyle\geqslant$
$\displaystyle\|u\|^{2}-\|u\|\cdot\|u^{n}\|+\|\delta_{x}u\|^{2}-\|\delta_{x}u\|\cdot\|\delta_{x}u^{n}\|$
$\displaystyle\geqslant$
$\displaystyle\|u\|^{2}-(\frac{1}{2}\|u\|^{2}+\frac{1}{2}\|u^{n}\|^{2})+\|\delta_{x}u\|^{2}-(\|\delta_{x}u\|^{2}+\frac{1}{4}\|\delta_{x}u^{n}\|^{2})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\|u\|^{2}-\frac{1}{2}\|u^{n}\|^{2}-\frac{1}{4}\|\delta_{x}u^{n}\|^{2},$

in which the Cauchy Schwarz inequality and $ab\leqslant\varepsilon
a^{2}+\frac{1}{4\varepsilon}b^{2}$ are utilized.

Substituting the above formulas into (3.18), we have

$\displaystyle(\Theta(u),u)$
$\displaystyle\geqslant\frac{2}{\tau}\Big{(}\frac{1}{2}\|u\|^{2}-\frac{1}{2}\|u^{n}\|^{2}-\frac{1}{4}\|\delta_{x}u^{n}\|^{2}\Big{)}+(1-2\Omega\kappa)\cdot\frac{2}{\tau}(-\frac{1}{4}\|\rho^{n}\|^{2})$
$\displaystyle=\frac{1}{\tau}\Big{(}\|u\|^{2}-\|u^{n}\|^{2}-\frac{1}{2}\|\delta_{x}u^{n}\|^{2}-\frac{1-2\Omega\kappa}{2}\|\rho^{n}\|^{2}\Big{)}.$

When
$\|u\|=\Big{(}\|u^{n}\|^{2}+\frac{1}{2}\|\delta_{x}u^{n}\|^{2}+\frac{1-2\Omega\kappa}{2}\|\rho^{n}\|^{2}\Big{)}^{1/2}$,
we have $(\Theta(u),u)\geqslant 0$. By the Browder theorem 3.2, there is a
$u^{\ast}\in\mathcal{U}_{h}$ satisfying

$\|u^{\ast}\|\leqslant\Big{(}\|u^{n}\|^{2}+\frac{1}{2}\|\delta_{x}u^{n}\|^{2}+\frac{1-2\Omega\kappa}{2}\|\rho^{n}\|^{2}\Big{)}^{1/2}$

such that $\Theta(u^{\ast})=0$.

This finishes the proof.

###### Remark 3.4

In fact, once $u$ has been determined, we can calculate $\rho$ by the
homogeneous system

$\displaystyle\frac{2}{\tau}\rho_{i}+\Delta_{x}(\rho u)_{i}=0,\quad
i\in\mathbb{I}_{M}.$

Taking an inner product of the above equality with $\rho$, we have

$\displaystyle\|\rho\|^{2}+\frac{\tau}{2}(\rho,\Delta_{x}(\rho u))=0.$

Rearranging the above inequality and by means of the Cauchy-Schwarz
inequality, it follows

$\displaystyle\|\rho\|^{2}=-\frac{\tau}{2}(\rho,\Delta_{x}(\rho
u))=\frac{\tau}{2}(\rho
u,\Delta_{x}\rho)\leqslant\frac{\tau}{2}\|u\|_{\infty}\cdot\|\rho\|\cdot\|\Delta_{x}\rho\|\leqslant\frac{\tau}{2h}\|u\|_{\infty}\cdot\|\rho\|^{2}.$

Therefore, there is a unique solution $\rho$ to the difference scheme (3.16b)
when $\frac{\tau}{2h}\|u\|_{\infty}<1$.

##  4 Convergence

###  4.1 Convergence for the difference scheme with zero rotation parameter
($\Omega=0$)

In this case, the problem (1.1) deduces to the familiar generalized two-
component Dullin-Gottwald-Holm system, see e.g., HGG2013 . Denote

$\displaystyle e_{i}^{n}=U_{i}^{n}-u_{i}^{n},\quad
f_{i}^{n}=\Pi_{i}^{n}-\rho_{i}^{n},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N},$ $\displaystyle
c_{2}=\frac{1}{2}\max\\{1+c_{\rm max}+3c_{0}+\sigma c_{0},\;c_{\rm
max}+4\sigma c_{0}\\},\quad
c_{3}=c_{1}\exp\Big{(}\frac{3}{2}c_{2}T\Big{)}\sqrt{\frac{L}{c_{2}}},$

then the convergence with $\Omega=0$ is followed.

###### Theorem 4.1

Suppose
$\\{U_{i}^{n},\,\Pi_{i}^{n}\,|\,i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}\\}$ is
the solution of the problem (1.1)–(1.3),
$\\{u_{i}^{n},\,\rho_{i}^{n}\,|\,i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}\\}$
is the approximation solution of the difference scheme (3.6). When $\Omega=0$
and $3c_{2}\tau\leq 1$, we have

$\displaystyle\|e^{n}\|\leqslant
c_{3}(\tau^{2}+h^{2}),\quad|e^{n}|_{1}\leqslant
c_{3}(\tau^{2}+h^{2}),\quad\|f^{n}\|\leqslant c_{3}(\tau^{2}+h^{2}),\quad
n\in\mathbb{K}_{N}.$

###### Proof

When $\Omega=0$, subtracting (3.6) from (3.3) and (3.5), we get the following
error system

$\displaystyle\delta_{t}e_{i}^{n+\frac{1}{2}}-\delta_{t}\delta_{x}^{2}e_{i}^{n+\frac{1}{2}}-\kappa\Delta_{x}e_{i}^{n+\frac{1}{2}}+3[\psi(U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})_{i}-\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})_{i}]$
$\displaystyle=3\sigma[\psi(\delta_{x}^{2}U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})_{i}-\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})_{i}]-\mu\Delta_{x}\delta_{x}^{2}e_{i}^{n+\frac{1}{2}}$
$\displaystyle\quad-[\Pi_{i}^{n+\frac{1}{2}}\Delta_{x}\Pi_{i}^{n+\frac{1}{2}}-\rho_{i}^{n+\frac{1}{2}}\Delta_{x}\rho_{i}^{n+\frac{1}{2}}]+Q_{i}^{n+\frac{1}{2}},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0},$ (4.19a)
$\displaystyle\delta_{t}f_{i}^{n+\frac{1}{2}}+[\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})_{i}-\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i}]=R_{i}^{n+\frac{1}{2}},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0},$ (4.19b) $\displaystyle
e_{i}^{0}=0,\;f_{i}^{0}=0,\quad i\in\mathbb{I}_{M},$ (4.19c) $\displaystyle
e_{i}^{n}=e_{i+M}^{n},\;f_{i}^{n}=f_{i+M}^{n},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}.$ (4.19d)

An inner product of (4.19a) with $e^{n+\frac{1}{2}}$ is carried out, we have

$\displaystyle\;(\delta_{t}e^{n+\frac{1}{2}},e^{n+\frac{1}{2}})-(\delta_{t}\delta_{x}^{2}e^{n+\frac{1}{2}},e^{n+\frac{1}{2}})-\kappa(\Delta_{x}e^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle\;+3(\psi(U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;3\sigma(\psi(\delta_{x}^{2}U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})-\mu(\Delta_{x}\delta_{x}^{2}e^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle-(\Pi^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}}-\rho^{n+\frac{1}{2}}\Delta_{x}\rho^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(Q^{n+\frac{1}{2}},e^{n+\frac{1}{2}}),\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0}.$

Employing Lemma 2.1, the equality above becomes

$\displaystyle\frac{1}{2\tau}\left[(\|e^{n+1}\|^{2}-\|e^{n}\|^{2})+(|e^{n+\frac{1}{2}}|_{1}^{2}-|e^{n}|_{1}^{2})\right]+3(\psi(U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;3\sigma(\psi(\delta_{x}^{2}U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})-(\Pi^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}}-\rho^{n+\frac{1}{2}}\Delta_{x}\rho^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle+(Q^{n+\frac{1}{2}},e^{n+\frac{1}{2}}),\quad
n\in\mathbb{K}_{N}^{0}.$ (4.20)

With the application of Lemma 2.1, we have

$\displaystyle-(\psi(U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle-(\psi(U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(U^{n+\frac{1}{2}}-e^{n+\frac{1}{2}},U^{n+\frac{1}{2}}-e^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle-(\psi(U^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+\psi(e^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(e^{n+\frac{1}{2}},e^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle-(\psi(e^{n+\frac{1}{2}},U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle-\frac{1}{3}(e^{n+\frac{1}{2}}\Delta_{x}U^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+\frac{1}{3}(e^{n+\frac{1}{2}}U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}}).$
(4.21)

Analogously, the first term and second term on the right-hand side of (4.20)
becomes

$\displaystyle(\psi(\delta_{x}^{2}U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle(\psi(\delta_{x}^{2}U^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+\psi(\delta_{x}^{2}e^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(\delta_{x}^{2}e^{n+\frac{1}{2}},e^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle(\psi(\delta_{x}^{2}e^{n+\frac{1}{2}},U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\frac{1}{3}(\delta_{x}^{2}e^{n+\frac{1}{2}},\Delta_{x}U^{n+\frac{1}{2}}\cdot
e^{n+\frac{1}{2}})-\frac{1}{3}(\delta_{x}^{2}e^{n+\frac{1}{2}}\cdot
U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}}),$ (4.22)

and

$\displaystyle-(\Pi^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}}-\rho^{n+\frac{1}{2}}\Delta_{x}\rho^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle-(\Pi^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}}-(\Pi^{n+\frac{1}{2}}-f^{n+\frac{1}{2}})\Delta_{x}(\Pi^{n+\frac{1}{2}}-f^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle-(\Pi^{n+\frac{1}{2}}\Delta_{x}f^{n+\frac{1}{2}}+f^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}}-f^{n+\frac{1}{2}}\Delta_{x}f^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle-(\Pi^{n+\frac{1}{2}}\Delta_{x}f^{n+\frac{1}{2}},e^{n+\frac{1}{2}})-(f^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(f^{n+\frac{1}{2}}\Delta_{x}f^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle(f^{n+\frac{1}{2}},\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}}))-(f^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}},e^{n+\frac{1}{2}})-(f^{n+\frac{1}{2}},\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}})).$
(4.23)

Substituting (4.21)–(4.23) into (4.20), we have

$\displaystyle\frac{1}{2\tau}\left[(\|e^{n+1}\|^{2}-\|e^{n}\|^{2})+(|e^{n+\frac{1}{2}}|_{1}^{2}-|e^{n}|_{1}^{2})\right]$
$\displaystyle=$
$\displaystyle\;-(e^{n+\frac{1}{2}}\Delta_{x}U^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(e^{n+\frac{1}{2}}U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})+\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}},\Delta_{x}U^{n+\frac{1}{2}}\cdot
e^{n+\frac{1}{2}})$ $\displaystyle-\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}}\cdot
U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})+(f^{n+\frac{1}{2}},\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}}))-(f^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle-(f^{n+\frac{1}{2}},\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}}))+(Q^{n+\frac{1}{2}},e^{n+\frac{1}{2}}),\quad
n\in\mathbb{K}_{N}^{0}.$ (4.24)

A similar inner product of (4.19b) with $f^{n+\frac{1}{2}}$ is made, we have

$\displaystyle\frac{1}{2\tau}(\|f^{n+1}\|^{2}-\|f^{n}\|^{2})+(\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})-\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})=(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}}).$
(4.25)

Similar to the derivation in (4.21), we have

$\displaystyle-(\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})-\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;-(\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})-\Delta_{x}(\Pi^{n+\frac{1}{2}}-f^{n+\frac{1}{2}})(U^{n+\frac{1}{2}}-e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;-(\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}})+\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}})-\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;-(\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})-(\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})+(\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}}).$
(4.26)

Plugging (4.26) into (4.25), and simple calculation yields

$\displaystyle\;\frac{1}{2\tau}(\|f^{n+1}\|^{2}-\|f^{n}\|^{2})$
$\displaystyle=$
$\displaystyle\;-(\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})-(\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})$
$\displaystyle\;+(\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})+(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}}).$
(4.27)

Summing up (4.24) with (4.27), we have

$\displaystyle\;\frac{1}{2\tau}\left[\|e^{n+1}\|^{2}-\|e^{n}\|^{2})+(|e^{n+\frac{1}{2}}|_{1}^{2}-|e^{n}|_{1}^{2})+(\|f^{n+1}\|^{2}-\|f^{n}\|^{2})\right]$
$\displaystyle=$
$\displaystyle\;-(e^{n+\frac{1}{2}}\Delta_{x}U^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(e^{n+\frac{1}{2}}U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})+\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}},\Delta_{x}U^{n+\frac{1}{2}}\cdot
e^{n+\frac{1}{2}})$
$\displaystyle\;-\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}}\cdot
U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})-(f^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}},e^{n+\frac{1}{2}})-(\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})$
$\displaystyle\;+(Q^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}}),$
$\displaystyle=$ $\displaystyle\;\sum_{i=1}^{8}J_{i},\quad
n\in\mathbb{K}_{N}^{0}.$ (4.28)

With the help of the Cauchy-Schwarz inequality and combining assumption
conditions (3.1) with (3.2), we easily have the estimates as

$\displaystyle
J_{1}=-(e^{n+\frac{1}{2}}\Delta_{x}U^{n+\frac{1}{2}},e^{n+\frac{1}{2}})\leqslant
c_{0}\|e^{n+\frac{1}{2}}\|^{2},$ $\displaystyle
J_{2}=(e^{n+\frac{1}{2}}U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})\leqslant
c_{\rm max}\|e^{n+\frac{1}{2}}\|\cdot|e^{n+\frac{1}{2}}|_{1},$ $\displaystyle
J_{5}=-(f^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}},e^{n+\frac{1}{2}})\leqslant
c_{0}\|f^{n+\frac{1}{2}}\|\cdot\|e^{n+\frac{1}{2}}\|,$ $\displaystyle
J_{7}=(Q^{n+\frac{1}{2}},e^{n+\frac{1}{2}})\leqslant\|Q^{n+\frac{1}{2}}\|\cdot\|e^{n+\frac{1}{2}}\|,$
$\displaystyle
J_{8}=(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}})\leqslant\|R^{n+\frac{1}{2}}\|\cdot\|f^{n+\frac{1}{2}}\|.$

The remaining terms ($J_{3}$, $J_{4}$ and $J_{6}$) will be estimated
respectively. Firstly, we estimate $J_{3}$ as

$\displaystyle J_{3}=$
$\displaystyle\;\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}},\Delta_{x}U^{n+\frac{1}{2}}\cdot
e^{n+\frac{1}{2}})$ $\displaystyle=$
$\displaystyle\;-\sigma\langle\delta_{x}e^{n+\frac{1}{2}},\delta_{x}(\Delta_{x}U^{n+\frac{1}{2}}\cdot
e^{n+\frac{1}{2}})\rangle$ $\displaystyle=$
$\displaystyle\;-\sigma\sum_{i\in\mathbb{I}_{M}}\delta_{x}e_{i}^{n+\frac{1}{2}}(\Delta_{x}U_{i+1}^{n+\frac{1}{2}}\cdot
e_{i+1}^{n+\frac{1}{2}}-\Delta_{x}U_{i}^{n+\frac{1}{2}}\cdot
e_{i}^{n+\frac{1}{2}})$ $\displaystyle=$
$\displaystyle\;-\sigma\sum_{i\in\mathbb{I}_{M}}\delta_{x}e_{i}^{n+\frac{1}{2}}\Big{[}\Delta_{x}U_{i+1}^{n+\frac{1}{2}}\cdot(e_{i+1}^{n+\frac{1}{2}}-e_{i}^{n+\frac{1}{2}})+(\Delta_{x}U_{i+1}^{n+\frac{1}{2}}-\Delta_{x}U_{i}^{n+\frac{1}{2}})\cdot
e_{i}^{n+\frac{1}{2}}\Big{]}$ $\displaystyle=$ $\displaystyle\;-\sigma
h\sum_{i\in\mathbb{I}_{M}}\delta_{x}e_{i}^{n+\frac{1}{2}}\Big{[}\Delta_{x}U_{i+1}^{n+\frac{1}{2}}\cdot\delta_{x}e_{i+\frac{1}{2}}^{n+\frac{1}{2}}+\delta_{x}(\Delta_{x}U_{i+\frac{1}{2}}^{n+\frac{1}{2}})\cdot
e_{i}^{n+\frac{1}{2}}\Big{]}$ $\displaystyle\leqslant$
$\displaystyle\;\sigma|e^{n+1}|_{1}^{2}\cdot\|\Delta_{x}U^{n+\frac{1}{2}}\|_{\infty}+\sigma|e^{n+\frac{1}{2}}|_{1}\cdot\|\delta_{x}(\Delta_{x}U^{n+\frac{1}{2}})\|_{\infty}\cdot\|e^{n+\frac{1}{2}}\|$
$\displaystyle\leqslant$ $\displaystyle\;\sigma c_{0}|e^{n+1}|_{1}^{2}+\sigma
c_{0}|e^{n+\frac{1}{2}}|_{1}\cdot\|e^{n+\frac{1}{2}}\|.$

Then we directly estimate $J_{4}$ as

$\displaystyle J_{4}=$
$\displaystyle\;-\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}}\cdot
U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})$ $\displaystyle=$
$\displaystyle\;-\sigma
h\sum_{i\in\mathbb{I}_{M}}(\delta_{x}^{2}e_{i}^{n+\frac{1}{2}})U_{i}^{n+\frac{1}{2}}\cdot\Delta_{x}e_{i}^{n+\frac{1}{2}}$
$\displaystyle=$
$\displaystyle\;\frac{\sigma}{2}\sum_{i\in\mathbb{I}_{M}}\Big{[}(\delta_{x}e_{i-\frac{1}{2}}^{n+\frac{1}{2}})^{2}-(\delta_{x}e_{i+\frac{1}{2}}^{n+\frac{1}{2}})^{2}\Big{]}U_{i}^{n+\frac{1}{2}}$
$\displaystyle=$ $\displaystyle\;\frac{\sigma
h}{2}\sum_{i\in\mathbb{I}_{M}}(\delta_{x}e_{i+\frac{1}{2}}^{n+\frac{1}{2}})^{2}\delta_{x}U_{i+\frac{1}{2}}^{n+\frac{1}{2}}$
$\displaystyle\leqslant$ $\displaystyle\;\frac{\sigma
c_{0}}{2}|e^{n+\frac{1}{2}}|_{1}^{2}.$

In light of the summation by parts, $J_{6}$ is calculated as

$\displaystyle J_{6}=$
$\displaystyle\;-(\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}},\Delta_{x}f^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;h\sum_{i\in\mathbb{I}_{M}}f_{i}^{n+\frac{1}{2}}U_{i}^{n+\frac{1}{2}}\Delta_{x}f_{i}^{n+\frac{1}{2}}$
$\displaystyle=$
$\displaystyle\;\frac{1}{2}\sum_{i\in\mathbb{I}_{M}}f_{i}^{n+\frac{1}{2}}U_{i}^{n+\frac{1}{2}}(f_{i+1}^{n+\frac{1}{2}}-f_{i-1}^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;\frac{1}{2}\sum_{i\in\mathbb{I}_{M}}(f_{i}^{n+\frac{1}{2}}f_{i+1}^{n+\frac{1}{2}}U_{i}^{n+\frac{1}{2}}-f_{i}^{n+\frac{1}{2}}f_{i-1}^{n+\frac{1}{2}}U_{i}^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;\frac{h}{2}\sum_{i\in\mathbb{I}_{M}}f_{i}^{n+\frac{1}{2}}f_{i+1}^{n+\frac{1}{2}}\cdot\frac{U_{i}^{n+\frac{1}{2}}-U_{i+1}^{n+\frac{1}{2}}}{h}$
$\displaystyle\leqslant$
$\displaystyle\;\frac{c_{0}}{2}\|f^{n+\frac{1}{2}}\|^{2}.$

Substituting $J_{i}\;(1\leqslant i\leqslant 8)$ into (4.28), we have

$\displaystyle\;\frac{1}{2\tau}\left[\|e^{n+1}\|^{2}-\|e^{n}\|^{2})+(|e^{n+\frac{1}{2}}|_{1}^{2}-|e^{n}|_{1}^{2})+(\|f^{n+1}\|^{2}-\|f^{n}\|^{2})\right]$
$\displaystyle\leqslant$ $\displaystyle\;c_{0}\|e^{n+\frac{1}{2}}\|^{2}+c_{\rm
max}\|e^{n+\frac{1}{2}}\|\cdot|e^{n+\frac{1}{2}}|_{1}+\sigma
c_{0}|e^{n+1}|_{1}^{2}+\sigma
c_{0}|e^{n+\frac{1}{2}}|_{1}\cdot\|e^{n+\frac{1}{2}}\|+\frac{1}{2}\sigma
c_{0}|e^{n+\frac{1}{2}}|_{1}^{2}$
$\displaystyle\;+c_{0}\|f^{n+\frac{1}{2}}\|\cdot\|e^{n+\frac{1}{2}}\|+\frac{c_{0}}{2}\|f^{n+\frac{1}{2}}\|^{2}+\|Q^{n+\frac{1}{2}}\|\cdot\|e^{n+\frac{1}{2}}\|+\|R^{n+\frac{1}{2}}\|\cdot\|f^{n+\frac{1}{2}}\|$
$\displaystyle\leqslant$
$\displaystyle\;c_{0}\|e^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}c_{\rm
max}\|e^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}c_{\rm
max}|e^{n+\frac{1}{2}}|_{1}^{2}+\sigma
c_{0}|e^{n+\frac{1}{2}}|_{1}^{2}+\frac{1}{2}\sigma
c_{0}|e^{n+\frac{1}{2}}|_{1}^{2}+\frac{1}{2}\sigma
c_{0}\|e^{n+\frac{1}{2}}\|^{2}$ $\displaystyle\;+\frac{1}{2}\sigma
c_{0}|e^{n+\frac{1}{2}}|_{1}^{2}+\frac{c_{0}}{2}\|f^{n+\frac{1}{2}}\|^{2}+\frac{c_{0}}{2}\|e^{n+\frac{1}{2}}\|^{2}+\frac{c_{0}}{2}\|f^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}\|Q^{n+\frac{1}{2}}\|^{2}$
$\displaystyle\;+\frac{1}{2}\|e^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}\|R^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}\|f^{n+\frac{1}{2}}\|^{2}$
$\displaystyle=$ $\displaystyle\;(\frac{1}{2}+\frac{c_{\rm
max}}{2}+\frac{3}{2}c_{0}+\frac{\sigma}{2}c_{0})\|e^{n+\frac{1}{2}}\|^{2}+(\frac{c_{\rm
max}}{2}+2\sigma
c_{0})|e^{n+\frac{1}{2}}|_{1}^{2}+(c_{0}+\frac{1}{2})\|f^{n+\frac{1}{2}}\|^{2}$
$\displaystyle\;+\frac{1}{2}(\|Q^{n+\frac{1}{2}}\|^{2}+\|R^{n+\frac{1}{2}}\|^{2})$
$\displaystyle\leqslant$
$\displaystyle\;c_{2}(\|e^{n+\frac{1}{2}}\|^{2}+|e^{n+\frac{1}{2}}|_{1}^{2}+\|f^{n+\frac{1}{2}}\|^{2})+\frac{1}{2}(\|Q^{n+\frac{1}{2}}\|^{2}+\|R^{n+\frac{1}{2}}\|^{2})$
$\displaystyle\leqslant$
$\displaystyle\;\frac{c_{2}}{2}(\|e^{n+1}\|^{2}+\|e^{n}\|^{2}+|e^{n+1}|_{1}^{2}+|e^{n}|_{1}^{2}+\|f^{n+1}\|^{2}+\|f^{n}\|^{2})$
$\displaystyle\;+\frac{1}{2}(\|Q^{n+\frac{1}{2}}\|^{2}+\|R^{n+\frac{1}{2}}\|^{2}),\quad
n\in\mathbb{K}_{N}^{0}.$ (4.29)

Denote

$\displaystyle F^{n}=\|e^{n}\|^{2}+|e^{n}|_{1}^{2}+\|f^{n}\|^{2},\quad
n\in\mathbb{K}_{N}.$

Therefore, (4.29) becomes

$\displaystyle\frac{1}{2\tau}(F^{n+1}-F^{n})\leqslant\frac{c_{2}}{2}(F^{n+1}+F^{n})+\frac{1}{2}(\|Q^{n+\frac{1}{2}}\|^{2}+\|R^{n+\frac{1}{2}}\|^{2}),\quad
n\in\mathbb{K}_{N}^{0}.$

Multiplying the inequality above by $2\tau$ and then noticing (3.4), we have

$\displaystyle(1-c_{2}\tau)F^{n+1}\leqslant(1+c_{2}\tau)F^{n}+2Lc_{1}^{2}\tau(\tau^{2}+h^{2})^{2},\quad
n\in\mathbb{K}_{N}^{0}.$

When $c_{2}\tau\leqslant\frac{1}{3}$, we have

$\displaystyle
F^{n+1}\leqslant(1+3c_{2}\tau)F^{n}+3Lc_{1}^{2}\tau(\tau^{2}+h^{2})^{2},\quad
n\in\mathbb{K}_{N}^{0}.$

Employing the discrete Gronwall inequality, we have

$\displaystyle
F^{n}\leqslant\exp(3c_{2}T)\cdot\frac{Lc_{1}^{2}}{c_{2}}(\tau^{2}+h^{2})^{2}=c_{3}^{2}(\tau^{2}+h^{2})^{2},\quad
n\in\mathbb{K}_{N}.$

According to the definition of $F^{n}$, Theorem 4.1 holds, which completes the
proof.

###  4.2 Convergence for the difference scheme with nonzero rotation parameter
($\Omega\neq 0$)

Next we further analyze the convergence of the difference scheme (3.6) when
$\Omega\neq 0$ by taking an inner product of (3.6b) twice but with different
quantities, and based on the technical energy analysis. For this purpose,
denote

$\displaystyle c_{4}=\frac{1}{2}(3c_{0}+\sigma c_{0}+1),\qquad
c_{5}=\frac{c_{0}}{2}(4\sigma+1-2\Omega\kappa+c_{0}),$ $\displaystyle
c_{6}=\frac{1}{2}\big{[}(2c_{0}+1)(1-2\Omega\kappa)+2(c_{0}+c_{\rm
max})\Omega\big{]},\qquad c_{7}=\frac{1}{2}\big{(}2-\Omega\kappa+2\Omega
c_{0}\big{)},$ $\displaystyle
c_{8}=\frac{1}{2}\max\left\\{c_{4},\,c_{5},\,\frac{c_{6}+3c_{0}\Omega}{1-2\Omega(\kappa+c_{0})}\right\\},\qquad
c_{9}=c_{1}\exp(3c_{8}T)\sqrt{\frac{L}{2}\cdot\frac{c_{7}}{c_{8}}},$
$\displaystyle E(0)=E(u^{0}(x),\rho^{0}(x)),\qquad
c_{10}=\sqrt{\frac{1+\sqrt{1+4L^{2}}}{L}\cdot E(0)}.$

Then the convergence under the restriction of small initial energy is listed
as follows.

###### Theorem 4.2

Suppose
$\\{U_{i}^{n},\,\Pi_{i}^{n}\,|\,i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}\\}$ is
the solution of the problem (1.1)–(1.3),
$\\{u_{i}^{n},\,\rho_{i}^{n}\,|\,i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}\\}$
is the solution of the difference scheme (3.6). When
$c_{10}\leqslant\frac{1}{2\Omega}-\kappa$ and $6c_{8}\tau<1$, we have the
error estimates

$\displaystyle\|e^{n}\|\leqslant
c_{9}(\tau^{2}+h^{2}),\quad|e^{n}|_{1}\leqslant
c_{9}(\tau^{2}+h^{2}),\quad\|f^{n}\|\leqslant\frac{c_{9}}{\sqrt{1-2\Omega(\kappa+c_{10})}}(\tau^{2}+h^{2}),\quad
n\in\mathbb{K}_{N}.$

###### Proof

When $\Omega\neq 0$, subtracting (3.3) and (3.5) from (3.6) directly, we get
the following error system of equations

$\displaystyle\delta_{t}e_{i}^{n+\frac{1}{2}}-\delta_{t}\delta_{x}^{2}e_{i}^{n+\frac{1}{2}}-\kappa\Delta_{x}e_{i}^{n+\frac{1}{2}}+3[\psi(U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})_{i}-\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})_{i}]$
$\displaystyle=3\sigma[\psi(\delta_{x}^{2}U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})_{i}-\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}})_{i}]-\mu\Delta_{x}\delta_{x}^{2}e_{i}^{n+\frac{1}{2}}$
$\displaystyle\quad-(1-2\Omega\kappa)[\Pi_{i}^{n+\frac{1}{2}}\Delta_{x}\Pi_{i}^{n+\frac{1}{2}}-\rho_{i}^{n+\frac{1}{2}}\Delta_{x}\rho_{i}^{n+\frac{1}{2}}]$
$\displaystyle\quad+2\Omega[\Pi_{i}^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})_{i}-\rho_{i}^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i}]+Q_{i}^{n+\frac{1}{2}},$
$\displaystyle\qquad\qquad\qquad\qquad\qquad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0},$ (4.30a)
$\displaystyle\delta_{t}f_{i}^{n+\frac{1}{2}}+[\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})_{i}-\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i}]=R_{i}^{n+\frac{1}{2}},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0},$ (4.30b) $\displaystyle
e_{i}^{0}=0,\;f_{i}^{0}=0,\quad i\in\mathbb{I}_{M},$ (4.30c) $\displaystyle
e_{i}^{n}=e_{i+M}^{n},\;f_{i}^{n}=f_{i+M}^{n},\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}.$ (4.30d)

Denote

$\displaystyle
G^{n}=\|e^{n}\|^{2}+|e^{n}|_{1}^{2}+(1-2\Omega\kappa)\|f^{n}\|^{2}-2\Omega(f^{n},f^{n}U^{n}).$
(4.31)

By employing Lemma 2.2 in continuous counterpart and the energy conservation
(1.4), and then taking $\varepsilon=\frac{1+\sqrt{1+4L^{2}}}{2L}$, we have

$\displaystyle c^{2}_{\rm max}=\|u(\cdot,t)\|_{\infty}^{2}$
$\displaystyle\leqslant\varepsilon|u(\cdot,t)|_{1}^{2}+\Big{(}\frac{1}{L}+\frac{1}{\varepsilon}\Big{)}\|u(\cdot,t)\|^{2}$
$\displaystyle=\varepsilon(|u(\cdot,t)|_{1}^{2}+\|u(\cdot,t)\|^{2})$
$\displaystyle\leqslant 2\varepsilon E(u,\rho)=2\varepsilon E(0)=c_{10}^{2}.$

When $c_{10}=\sqrt{2\varepsilon E(0)}\leqslant\frac{1}{2\Omega}-\kappa$, we
have $1-2\Omega\kappa-2\Omega c_{10}\geqslant 0$. The nonnegativity of $G^{n}$
at the moment is guaranteed because of

$\displaystyle G^{n}\geqslant$
$\displaystyle\|e^{n}\|^{2}+|e^{n}|_{1}^{2}+(1-2\Omega\kappa-2\Omega c_{\rm
max})\|f^{n}\|^{2}$ $\displaystyle\geqslant$
$\displaystyle\|e^{n}\|^{2}+|e^{n}|_{1}^{2}+(1-2\Omega\kappa-2\Omega
c_{10})\|f^{n}\|^{2}\geqslant 0.$ (4.32)

Taking a $L^{2}$-inner product of (4.30a) with $e^{n+\frac{1}{2}}$ and
noticing Lemma 2.1, we have

$\displaystyle\frac{1}{2\tau}\left[(\|e^{n+1}\|^{2}-\|e^{n}\|^{2})+(|e^{n+\frac{1}{2}}|_{1}^{2}-|e^{n}|_{1}^{2})\right]+3(\psi(U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;3\sigma(\psi(\delta_{x}^{2}U^{n+\frac{1}{2}},U^{n+\frac{1}{2}})-\psi(\delta_{x}^{2}u^{n+\frac{1}{2}},u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle-(1-2\Omega\kappa)(\Pi^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}}-\rho^{n+\frac{1}{2}}\Delta_{x}\rho^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(Q^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle+2\Omega(\Pi^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})-\rho^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}}),\quad
i\in\mathbb{I}_{M},\;n\in\mathbb{K}_{N}^{0}.$ (4.33)

Noticing that expression in the second term on the right-hand side of the
equation above can be written as

$\displaystyle\;\Pi_{i}^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})_{i}-\rho^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i}$
$\displaystyle=$
$\displaystyle\;\Pi_{i}^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}}-\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i}+f_{i}^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i}$
$\displaystyle=$
$\displaystyle\;\Pi_{i}^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}}+f^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i}+f_{i}^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}})_{i},$

and using Lemma 2.1 again, the second term in the right-hand side of (4.33)
becomes

$\displaystyle\;(\Pi^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}})-\rho^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;(\Pi^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})+(f^{n+\frac{1}{2}}\Delta_{x}(\rho^{n+\frac{1}{2}}u^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;(\Pi^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}(U^{n+\frac{1}{2}}-e^{n+\frac{1}{2}})),e^{n+\frac{1}{2}})$
$\displaystyle\;+(f^{n+\frac{1}{2}}\Delta_{x}[(\Pi^{n+\frac{1}{2}}-f^{n+\frac{1}{2}})(U^{n+\frac{1}{2}}-e^{n+\frac{1}{2}})],e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;(\Pi^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})-(\Pi^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle\;+(f^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})-(f^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle\;-(f^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})+(f^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;(\Pi^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})+(f^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle\;-(f^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}}).$
(4.34)

Substituting (4.21)–(4.23) and (4.34) into (4.33), for any
$n\in\mathbb{K}_{N}^{0}$, we have

$\displaystyle\;\frac{1}{2\tau}\left[(\|e^{n+1}\|^{2}-\|e^{n}\|^{2})+(|e^{n+\frac{1}{2}}|_{1}^{2}-|e^{n}|_{1}^{2})\right]$
$\displaystyle=$
$\displaystyle\;-(e^{n+\frac{1}{2}}\Delta_{x}U^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(e^{n+\frac{1}{2}}U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})+\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}},\Delta_{x}U^{n+\frac{1}{2}}\cdot
e^{n+\frac{1}{2}})$
$\displaystyle\;-\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}}\cdot
U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})+(1-2\Omega\kappa)[(f^{n+\frac{1}{2}},\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}}))$
$\displaystyle\;-(f^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}},e^{n+\frac{1}{2}})-(f^{n+\frac{1}{2}},\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}}))]+2\Omega[(\Pi^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle\;+(f^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})-(f^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})]+(Q^{n+\frac{1}{2}},e^{n+\frac{1}{2}}).$
(4.35)

Taking an inner product of (4.19b) with $(1-2\Omega\kappa)f^{n+\frac{1}{2}}$
and noticing (4.27), we can derive that

$\displaystyle\;\frac{(1-2\Omega\kappa)}{2\tau}(\|f^{n+1}\|^{2}-\|f^{n}\|^{2})$
$\displaystyle=$
$\displaystyle\;(1-2\Omega\kappa)[-(\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})-(\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})$
$\displaystyle\;+(\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})]+(1-2\Omega\kappa)(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}}).$
(4.36)

Summing up (4.35) and (4.36), we have

$\displaystyle\;\frac{1}{2\tau}\left[(\|e^{n+1}\|^{2}-\|e^{n}\|^{2})+(|e^{n+\frac{1}{2}}|_{1}^{2}-|e^{n}|_{1}^{2})+(1-2\Omega\kappa)(\|f^{n+1}\|^{2}-\|f^{n}\|^{2})\right]$
$\displaystyle=$
$\displaystyle\;-(e^{n+\frac{1}{2}}\Delta_{x}U^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(e^{n+\frac{1}{2}}U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})+\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}},\Delta_{x}U^{n+\frac{1}{2}}\cdot
e^{n+\frac{1}{2}})$
$\displaystyle\;-\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}}\cdot
U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})+(1-2\Omega\kappa)[-(f^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle\;-(\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})]+2\Omega[(\Pi^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
$\displaystyle\;+(f^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})-(f^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})]$
$\displaystyle\;+(Q^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(1-2\Omega\kappa)(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}}).$
(4.37)

In order to eliminate the difficulty of estimating the term
$-(f^{n+\frac{1}{2}}\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})$
in (4.37), we further taking an inner product of (4.19b) with $-2\Omega
f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}$. Similar to the derivation in (4.27), we
have

$\displaystyle\;-2\Omega(\delta_{t}f^{n+\frac{1}{2}},f^{n+\frac{1}{2}}U^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;2\Omega[(\Delta_{x}(\Pi^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}}U^{n+\frac{1}{2}})+(\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),f^{n+\frac{1}{2}}U^{n+\frac{1}{2}})$
$\displaystyle\;-(\Delta_{x}(f^{n+\frac{1}{2}}e^{n+\frac{1}{2}}),f^{n+\frac{1}{2}}U^{n+\frac{1}{2}})]-2\Omega(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}).$
(4.38)

Directly applying Lemma 2.3 to the first term in (4.38), it becomes

$\displaystyle\;-2\Omega(\delta_{t}f^{n+\frac{1}{2}},f^{n+\frac{1}{2}}U^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;-2\Omega(\delta_{t}f^{n+\frac{1}{2}},f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}-(fU)^{n+\frac{1}{2}}+(fU)^{n+\frac{1}{2}})$
$\displaystyle=$
$\displaystyle\;-\frac{\Omega}{\tau}[(f^{n+1},f^{n+1}U^{n+1})-(f^{n},f^{n}U^{n})]+\Omega(f^{n}f^{n+1},\delta_{t}U^{n+\frac{1}{2}})$
$\displaystyle\;+\frac{\Omega}{2}(f^{n+1}-f^{n},(f^{n+1}-f^{n})\delta_{t}U^{n+\frac{1}{2}}).$
(4.39)

Substituting (4.39) into (4.38) and adding with (4.37) together, we have

$\displaystyle\;\frac{1}{2\tau}(G^{n+1}-G^{n})$ $\displaystyle=$
$\displaystyle\;-(e^{n+\frac{1}{2}}\Delta_{x}U^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(e^{n+\frac{1}{2}}U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})+\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}},\Delta_{x}U^{n+\frac{1}{2}}\cdot
e^{n+\frac{1}{2}})$
$\displaystyle\;-\sigma(\delta_{x}^{2}e^{n+\frac{1}{2}}\cdot
U^{n+\frac{1}{2}},\Delta_{x}e^{n+\frac{1}{2}})+(1-2\Omega\kappa)[-(f^{n+\frac{1}{2}}\Delta_{x}\Pi^{n+\frac{1}{2}},e^{n+\frac{1}{2}})$
$\displaystyle\;-(\Delta_{x}(f^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),f^{n+\frac{1}{2}})]+(Q^{n+\frac{1}{2}},e^{n+\frac{1}{2}})+(1-2\Omega\kappa)(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}})$
$\displaystyle\;+2\Omega(f^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})-2\Omega(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}}U^{n+\frac{1}{2}})$
$\displaystyle\;-\Omega(f^{n}f^{n+1},\delta_{t}U^{n+\frac{1}{2}})-\frac{\Omega}{2}(f^{n+1}-f^{n},(f^{n+1}-f^{n})\delta_{t}U^{n+\frac{1}{2}})$
$\displaystyle=$ $\displaystyle\;\sum_{i=1}^{12}P_{i},\quad
n\in\mathbb{K}_{N}^{0},$

in which four terms diminish. The estimates of $P_{i}\;(1\leqslant i\leqslant
8)$ are similar to $J_{i}\;(1\leqslant i\leqslant 8)$. The remaining terms
$P_{i}\;(9\leqslant i\leqslant 12)$ will be estimated respectively as

$\displaystyle
P_{9}=2\Omega(f^{n+\frac{1}{2}}\Delta_{x}(\Pi^{n+\frac{1}{2}}U^{n+\frac{1}{2}}),e^{n+\frac{1}{2}})\leqslant
2c_{0}^{2}\Omega\|f^{n+\frac{1}{2}}\|\cdot\|e^{n+\frac{1}{2}}\|,$
$\displaystyle
P_{10}=-2\Omega(R^{n+\frac{1}{2}},f^{n+\frac{1}{2}}U^{n+\frac{1}{2}})\leqslant
2c_{\rm max}\Omega\|R^{n+\frac{1}{2}}\|\cdot\|f^{n+\frac{1}{2}}\|,$
$\displaystyle
P_{11}=-\Omega(f^{n}f^{n+1},\delta_{t}U^{n+\frac{1}{2}})\leqslant
c_{0}\Omega\|f^{n}\|\cdot\|f^{n+1}\|,$ $\displaystyle
P_{12}=-\frac{\Omega}{2}(f^{n+1}-f^{n},(f^{n+1}-f^{n})\delta_{t}U^{n+\frac{1}{2}})\leqslant
c_{0}\Omega(\|f^{n}\|^{2}+\|f^{n+1}\|^{2}).$

Therefore,

$\displaystyle\;\frac{1}{2\tau}(G^{n+1}-G^{n})$ $\displaystyle\leqslant$
$\displaystyle\;c_{0}\|e^{n+\frac{1}{2}}\|^{2}+c_{0}\|e^{n+\frac{1}{2}}\|\cdot|e^{n+\frac{1}{2}}|_{1}+c_{0}\sigma|e^{n+1}|_{1}^{2}+c_{0}\sigma|e^{n+\frac{1}{2}}|_{1}\cdot\|e^{n+\frac{1}{2}}\|$
$\displaystyle\;+\frac{1}{2}c_{0}\sigma|e^{n+\frac{1}{2}}|_{1}^{2}+c_{0}(1-2\Omega\kappa)\|f^{n+\frac{1}{2}}\|\cdot\|e^{n+\frac{1}{2}}\|+\frac{1}{2}c_{0}(1-2\Omega\kappa)\|f^{n+\frac{1}{2}}\|^{2}$
$\displaystyle\;+\|Q^{n+\frac{1}{2}}\|\cdot\|e^{n+\frac{1}{2}}\|+(1-2\Omega\kappa)\|R^{n+\frac{1}{2}}\|\cdot\|f^{n+\frac{1}{2}}\|+2c_{0}^{2}\Omega\|e^{n+\frac{1}{2}}\|\cdot\|f^{n+\frac{1}{2}}\|$
$\displaystyle\;+2c_{\rm
max}\Omega\|R^{n+\frac{1}{2}}\|\cdot\|f^{n+\frac{1}{2}}\|+c_{0}\Omega\|f^{n}\|\cdot\|f^{n+1}\|+c_{0}\Omega(\|f^{n}\|^{2}+\|f^{n+1}\|^{2})$
$\displaystyle\leqslant$
$\displaystyle\;c_{0}\|e^{n+\frac{1}{2}}\|^{2}+\frac{c_{0}^{2}}{2}\|e^{n+\frac{1}{2}}\|^{2}+\frac{c_{0}^{2}}{2}|e^{n+\frac{1}{2}}|_{1}^{2}+c_{0}\sigma|e^{n+1}|_{1}^{2}+\frac{1}{2}c_{0}\sigma|e^{n+1}|_{1}^{2}+\frac{1}{2}c_{0}\sigma\|e^{n+\frac{1}{2}}\|^{2}$
$\displaystyle\;+\frac{c_{0}\sigma}{2}|e^{n+\frac{1}{2}}|_{1}^{2}+\frac{1}{2}c_{0}(1-2\Omega\kappa)\|f^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}c_{0}(1-2\Omega\kappa)|e^{n+\frac{1}{2}}|_{1}^{2}$
$\displaystyle\;+\frac{1}{2}c_{0}(1-2\Omega\kappa)\|f^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}\|Q^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}\|e^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}(1-2\Omega\kappa)\|R^{n+\frac{1}{2}}\|^{2}$
$\displaystyle\;+\frac{1}{2}(1-2\Omega\kappa)\|f^{n+\frac{1}{2}}\|^{2}+c_{0}^{2}\Omega\|e^{n+\frac{1}{2}}\|^{2}+c_{\rm
max}\Omega\|f^{n+\frac{1}{2}}\|^{2}+c_{\rm
max}\Omega\|R^{n+\frac{1}{2}}\|^{2}$
$\displaystyle\;+c_{0}\Omega\|f^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}c_{0}\Omega\|f^{n}\|^{2}+\frac{1}{2}c_{0}\Omega\|f^{n+1}\|^{2}+c_{0}\Omega(\|f^{n}\|^{2}+\|f^{n+1}\|^{2})$
$\displaystyle=$
$\displaystyle\;c_{4}\|e^{n+\frac{1}{2}}\|^{2}+c_{5}|e^{n+\frac{1}{2}}|_{1}^{2}+c_{6}\|f^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}3c_{0}\Omega\|f^{n}\|^{2}+\frac{1}{2}3c_{0}\Omega\|f^{n+1}\|^{2}$
$\displaystyle\;+\frac{1}{2}\|Q^{n+\frac{1}{2}}\|^{2}+\frac{1}{2}(1-2\Omega\kappa+2c_{\rm
max}\Omega)\|R^{n+\frac{1}{2}}\|^{2}$ $\displaystyle\leqslant$
$\displaystyle\;\frac{c_{4}}{2}\|e^{n+1}\|^{2}+\frac{c_{4}}{2}\|e^{n}\|^{2}+\frac{c_{5}}{2}|e^{n+1}|_{1}^{2}+\frac{c_{5}}{2}|e^{n}|_{1}^{2}+\frac{1}{2}(c_{6}+3c_{0}\Omega)\|f^{n+1}\|^{2}$
$\displaystyle\;+\frac{1}{2}(c_{6}+3c_{0}\Omega)\|f^{n}\|^{2}+c_{7}(\|Q^{n+\frac{1}{2}}\|^{2}+\|R^{n+\frac{1}{2}}\|^{2})$
$\displaystyle\leqslant$
$\displaystyle\;c_{8}(G^{n}+G^{n+1})+c_{7}(\|Q^{n+\frac{1}{2}}\|^{2}+\|R^{n+\frac{1}{2}}\|^{2}).$

Multiplying the inequality above by $2\tau$ on both sides and noticing (3.4),
we have

$\displaystyle(1-2c_{8}\tau)G^{n+1}\leqslant(1+2c_{8}\tau)G^{n}+4c_{7}c_{1}^{2}L\tau(\tau^{2}+h^{2})^{2},\quad
n\in\mathbb{K}_{N}^{0}.$

When $2c_{8}\tau\leqslant\frac{1}{3}$, we have

$\displaystyle
G^{n+1}\leqslant(1+6c_{8}\tau)G^{n}+6c_{7}c_{1}^{2}L\tau(\tau^{2}+h^{2})^{2},\quad
n\in\mathbb{K}_{N}^{0}.$

which implies

$\displaystyle
G^{n}\leqslant\exp(6c_{8}T)\cdot\frac{c_{7}c_{1}^{2}L}{c_{8}}(\tau^{2}+h^{2})^{2}\equiv
c_{9}^{2}(\tau^{2}+h^{2})^{2},\quad n\in\mathbb{K}_{N}$

by using the discrete Gronwall inequality. According to the definition of
$G^{n}$ in (4.31), Theorem 4.2 holds. We finish the proof.

###### Remark 4.5

The restriction condition $c_{10}<\frac{1}{2\Omega}-\kappa$ implies that the
initial energy satisfies

$E(0)<\Big{(}\frac{1-2\Omega\kappa}{2\Omega}\Big{)}^{2}\cdot\frac{\sqrt{1+4L^{2}}-1}{4L}.$
(4.40)

The inequality (4.40) further indicates that Theorem 4.2 holds only for small
initial value or small initial energy. In addition, we know from the proof
that

$c_{\rm max}<\frac{1+\sqrt{1+4L^{2}}}{L}\cdot E(0)=\frac{1}{2\Omega}-\kappa.$
(4.41)

In physics, the inequality (4.41) implies that Theorem 4.2 holds when the
maximum fluid velocity is no more than $\frac{1}{2\Omega}-\kappa$.

###### Remark 4.6

Theorem 4.2 states that there is no grid ratio restriction for the spatial and
temporal stepsizes. It implies that the numerical scheme (3.6) is
unconditionally convergent. Advantages of the difference scheme (3.6) compared
with those in literature ZLZ2022 are shown in Table 1. We see that the current
scheme performs much better.

Table 1: The pros and cons of the scheme (3.6) at present with that in
literature ZLZ2022 .

$\Omega=0$ $\Omega\neq 0$ Scheme in ZLZ2022 Scheme (3.6) Scheme in ZLZ2022
Scheme (3.6) Convergence condition $\tau\lesssim h$ $\tau\lesssim 1$
$\tau\lesssim h$? $\tau\lesssim 1$, small initial energy Convergence order
$O(h^{2}+\tau^{2})$ $O(h^{2}+\tau^{2})$ $O(h^{2}+\tau^{2})?$
$O(h^{2}+\tau^{2})\;$ Invariants Energy, Mass Energy, Momentum, Mass Energy,
Mass Energy, Momentum, Mass Solvability $\checkmark$ $\checkmark$ $\checkmark$
$\checkmark\;$

“$\checkmark$” denotes that the result has been proved; “$?$” denotes that the
results are numerically correct, but lack of a theoretical proof.

##  5 Numerical results

We will present two benchmark problems (includes several cases in each one) to
detect numerical theories including the unconditional convergence, invariant-
preserving properties (energy/momentum/mass), and long-time simulation on
large domain.

To facilitate calculation, we propose a subsequent two-level iteration method
to solve the nonlinear numerical scheme (3.6). Suppose we have known
$\\{u_{i}^{n},\,\rho_{i}^{n}\,|\,i\in\mathbb{I}_{M}\\}$, the following linear
system of equations will be used to approximate the solution of the difference
scheme (3.6),

$\displaystyle\frac{2}{\tau}(u_{i}^{(l+1)}-u_{i}^{n})-\frac{2}{\tau}(\delta_{x}^{2}u_{i}^{(l+1)}-\delta_{x}^{2}u_{i}^{n})-\kappa\Delta_{x}u_{i}^{(l+1)}+3\psi(u^{(l)},u^{(l+1)})_{i}-3\sigma\psi(\delta_{x}^{2}u^{(l)},u^{(l+1)})_{i}$
$\displaystyle\quad+\mu\Delta_{x}\delta_{x}^{2}u_{i}^{(l+1)}+(1-2\Omega\kappa)\rho_{i}^{(l)}\Delta_{x}\rho_{i}^{(l+1)}-2\Omega\rho_{i}^{(l)}\Delta_{x}(\rho^{(l)}u^{(l+1)})_{i}=0,\quad
i\in\mathbb{I}_{M},$ (5.42a)
$\displaystyle\frac{2}{\tau}(\rho_{i}^{(l+1)}-\rho_{i}^{n})+\Delta_{x}(\rho^{(l)}u^{(l+1)})_{i}=0,\quad
i\in\mathbb{I}_{M},$ (5.42b) $\displaystyle
u_{i}^{(0)}=u_{i}^{n},\;\rho_{i}^{(0)}=\rho_{i}^{n},\quad i\in\mathbb{I}_{M},$
(5.42c) $\displaystyle
u_{i}^{(l)}=u_{i+M}^{(l)},\;\rho_{i}^{(l)}=\rho_{i+M}^{(l)},\quad
i\in\mathbb{I}_{M},$ (5.42d)

until

$\max\big{\\{}\max_{i\in\mathbb{I}_{M}}|u_{i}^{(l+1)}-u_{i}^{(l)}|,\max_{i\in\mathbb{I}_{M}}|\rho_{i}^{(l+1)}-\rho_{i}^{(l)}|\big{\\}}\leqslant\varepsilon,\quad
l=0,\;1,\;\ldots.$

Then we have

$u_{i}^{n+1}=2u_{i}^{(l+1)}-u_{i}^{n},\quad\rho_{i}^{n+1}=2\rho_{i}^{(l+1)}-\rho_{i}^{n},\quad
i\in\mathbb{I}_{M}.$

In numerical implementation, we fix the tolerance error $\varepsilon=1e-12$
for each iteration unless otherwise specified. The formulas

$\displaystyle\mathrm{Ord_{\infty}^{h}}=\mathbf{\log_{2}}\frac{\|\mathrm{F_{u}}(2h,\tau)\|_{\infty}}{\|\mathrm{F_{u}}(h,\tau)\|_{\infty}},\qquad\quad\mathrm{Ord_{2}^{h}}=\mathbf{\log_{2}}\frac{\|\mathrm{F_{\rho}}(2h,\tau)\|}{\|\mathrm{F_{\rho}}(h,\tau)\|},$
$\displaystyle\mathrm{Ord_{\infty}^{\tau}}=\mathbf{\log_{2}}\frac{\|\mathrm{G_{u}}(h,2\tau)\|_{\infty}}{\|\mathrm{G_{u}}(h,\tau)\|_{\infty}},\qquad\quad\mathrm{Ord_{2}^{\tau}}=\mathbf{\log_{2}}\frac{\|\mathrm{G_{\rho}}(h,2\tau)\|}{\|\mathrm{G_{\rho}}(h,\tau)\|}$

are borrowed to test the convergence orders in space and time, where the
$L^{\infty}$-norm and the $L^{2}$-norm of numerical errors are defined by

$\displaystyle\|\mathrm{F_{u}}(h,\tau)\|_{\infty}=\max_{i\in\mathbb{I}_{M},k\in\mathbb{K}_{N}}\limits\Big{|}u_{i}^{k}(h,\tau)-u_{2i}^{k}(h/2,\tau)\Big{|},$
$\displaystyle\|\mathrm{G_{u}}(h,\tau)\|_{\infty}=\max_{i\in\mathbb{I}_{M},k\in\mathbb{K}_{N}}\limits\Big{|}u_{i}^{k}(h,\tau)-u_{i}^{2k}(h,\tau/2)\Big{|},$
$\displaystyle\|\mathrm{F_{\rho}}(h,\tau)\|=\max_{k\in\mathbb{K}_{N}}\limits\Big{(}{h\sum_{i\in\mathbb{I}_{M}}|\rho_{i}^{k}(h,\tau)-\rho_{2i}^{k}(h/2,\tau)|^{2}}\Big{)}^{1/2},$
$\displaystyle\|\mathrm{G_{\rho}}(h,\tau)\|=\max_{k\in\mathbb{K}_{N}}\limits\Big{(}{h\sum_{i\in\mathbb{I}_{M}}|\rho_{i}^{k}(h,\tau)-\rho_{i}^{2k}(h,\tau/2)|^{2}}\Big{)}^{1/2},$

respectively.

###### Example 1

(Dam break problem) The initial fluid velocity and free surface elevation from
equilibrium are taken, respectively, as

$\displaystyle
u^{0}(x)=0,\quad\rho^{0}(x)=1+\tanh(x+\emph{a})-\tanh(x-\emph{a}),$

with a being a dam-breaking parameter. Six groups of parameters are selected
in Table 6.

Table 2: Selected parameters and the calculated domain in the numerical tests.

${\rm~{}~{}~{}~{}~{}\qquad\qquad Parameters}$ ${\rm Domain}~{}{\rm
of}~{}(x,t)$ Case A $\emph{a}=0.1$, $\kappa=\mu=\Omega=0$, $\sigma=1$, CMR2014
; HI2010 ; YFS2018 $[-6,6]\times[0,20];\;$ Case B $\emph{a}=4$,
$\kappa=\mu=\Omega=0$, $\sigma=1$, ZLZ2022 $[-12\pi,12\pi]\times[0,2];\;$ Case
C $\emph{a}=0.2$, $\kappa=0,\,\mu=\sigma=1,\,\Omega=73\times 10^{-6}$, ZLZ2022
$[-8,8]\times[0,1];\;$ Case D $\emph{a}=\kappa=\mu=\sigma=1,\,\Omega=73\times
10^{-6}$, ZLZ2022 $[-8,8]\times[0,1];\;$ Case E
$\emph{a}=\mu=\sigma=1,\kappa=0,\,\Omega=73\times 10^{-6}$,ZLZ2022
$[-12\pi,12\pi]\times[0,50];\;$ Case F
$\emph{a}=\kappa=\mu=\sigma=1,\,\Omega=73\times 10^{-6}$,
$[-100,100]\times[0,1000].\;$

We report the numerical data in Tables 3–5 and display the numerical behavior
in Figures 2–4.

Table 3: Numerical errors against $h$-grid size reduction with the fixed
temporal stepsizes.

Case A, $\tau=1/50$ Case B, $\tau=1/50$ $h$
$\|\mathrm{F_{u}}(h,\tau)\|_{\infty}$ $\mathrm{Ord_{\infty}^{h}}$
$\|\mathrm{F_{\rho}}(h,\tau)\|$ $\mathrm{Ord_{2}^{h}}$
$\|\mathrm{F_{u}}(h,\tau)\|_{\infty}$ $\mathrm{Ord_{\infty}^{h}}$
$\|\mathrm{F_{\rho}}(h,\tau)\|$ $\mathrm{Ord_{2}^{h}}$ $\;0.6\;$ $3.1656{\rm
e}-02$ $*$ $8.2588{\rm e}-02$ $*$ $1.1779{\rm e}-01$ $*$ $3.1404{\rm e}-01$
$*\;$ $\;0.3\;$ $8.0761{\rm e}-03$ $1.9707$ $3.1629{\rm e}-02$ $1.3847$
$4.6146{\rm e}-02$ $1.3519$ $1.2761{\rm e}-01$ $1.2992\;$ $\;0.15\;$
$2.2533{\rm e}-03$ $1.8416$ $7.2516{\rm e}-03$ $2.1249$ $1.3660{\rm e}-02$
$1.7562$ $4.1642{\rm e}-02$ $1.6156\;$ $\;0.075\;$ $5.7025{\rm e}-04$ $1.9824$
$1.8181{\rm e}-03$ $1.9959$ $3.5951{\rm e}-03$ $1.9259$ $1.1627{\rm e}-02$
$1.8406\;$ $\;0.0375\;$ $1.4320{\rm e}-04$ $1.9936$ $4.5442{\rm e}-04$
$2.0003$ $9.0983{\rm e}-04$ $1.9824$ $3.0100{\rm e}-03$ $1.9496\;$ Case C,
$\tau=1/1000$ Case D, $\tau=1/1000$ $h$ $\|\mathrm{F_{u}}(h,\tau)\|_{\infty}$
$\mathrm{Ord_{\infty}^{h}}$ $\|\mathrm{F_{\rho}}(h,\tau)\|$
$\mathrm{Ord_{2}^{h}}$ $\|\mathrm{F_{u}}(h,\tau)\|_{\infty}$
$\mathrm{Ord_{\infty}^{h}}$ $\|\mathrm{F_{\rho}}(h,\tau)\|$
$\mathrm{Ord_{2}^{h}}$ $\;0.4\;$ $3.7012{\rm e}-03$ $*$ $8.9008{\rm e}-03$ $*$
$1.6694{\rm e}-02$ $*$ $3.8179{\rm e}-02$ $*\;$ $\;0.2\;$ $1.0564{\rm e}-03$
$1.8088$ $2.1868{\rm e}-03$ $2.0251$ $4.9151{\rm e}-03$ $1.7640$ $9.8260{\rm
e}-03$ $1.9581\;$ $\;0.1\;$ $2.7672{\rm e}-04$ $1.9327$ $5.4452{\rm e}-04$
$2.0058$ $1.3122{\rm e}-03$ $1.9053$ $2.4850{\rm e}-03$ $1.9834\;$ $\;0.05\;$
$7.0471{\rm e}-05$ $1.9733$ $1.3601{\rm e}-04$ $2.0012$ $3.3245{\rm e}-04$
$1.9808$ $6.2326{\rm e}-04$ $1.9953\;$ $\;0.025\;$ $1.7670{\rm e}-05$ $1.9957$
$3.3996{\rm e}-05$ $2.0003$ $8.3398{\rm e}-05$ $1.9951$ $1.5595{\rm e}-04$
$1.9988\;$

Table 4: Numerical errors against $\tau$-grid size reduction with the fixed
spatial stepsizes.

Case A, $h=6/25$ Case B, $h=6/25$ $\tau$ $\|\mathrm{F_{u}}(h,\tau)\|_{\infty}$
$\mathrm{Ord_{\infty}^{\tau}}$ $\|\mathrm{F_{\rho}}(h,\tau)\|$
$\mathrm{Ord_{2}^{\tau}}$ $\|\mathrm{F_{u}}(h,\tau)\|_{\infty}$
$\mathrm{Ord_{\infty}^{\tau}}$ $\|\mathrm{F_{\rho}}(h,\tau)\|$
$\mathrm{Ord_{2}^{\tau}}$ $\;0.25\;$ $1.2391{\rm e}-03$ $*$ $4.2968{\rm e}-03$
$*$ $5.3429{\rm e}-02$ $*$ $1.2532{\rm e}-01$ $*\;$ $\;0.125\;$ $3.1403{\rm
e}-04$ $1.9803$ $1.0730{\rm e}-03$ $2.0016$ $1.6057{\rm e}-02$ $1.7345$
$3.8919{\rm e}-02$ $1.6871\;$ $\;0.0625\;$ $7.8767{\rm e}-05$ $1.9952$
$2.6815{\rm e}-04$ $2.0006$ $4.1766{\rm e}-03$ $1.9428$ $1.0397{\rm e}-02$
$1.9043\;$ $\;0.03125\;$ $1.9708{\rm e}-05$ $1.9988$ $6.7032{\rm e}-05$
$2.0001$ $1.0538{\rm e}-03$ $1.9867$ $2.6455{\rm e}-03$ $1.9746\;$
$\;0.015625\;$ $4.9280{\rm e}-06$ $1.9997$ $1.6757{\rm e}-05$ $2.0001$
$2.6404{\rm e}-04$ $1.9968$ $6.6432{\rm e}-04$ $1.9936\;$ Case C, $h=4/25$
Case D, $h=4/25$ $\tau$ $\|\mathrm{F_{u}}(h,\tau)\|_{\infty}$
$\mathrm{Ord_{\infty}^{\tau}}$ $\|\mathrm{F_{\rho}}(h,\tau)\|$
$\mathrm{Ord_{2}^{\tau}}$ $\|\mathrm{F_{u}}(h,\tau)\|_{\infty}$
$\mathrm{Ord_{\infty}^{\tau}}$ $\|\mathrm{F_{\rho}}(h,\tau)\|$
$\mathrm{Ord_{2}^{\tau}}$ $\;0.0125\;$ $2.7005{\rm e}-06$ $*$ $3.4715{\rm
e}-06$ $*$ $3.4991{\rm e}-05$ $*$ $4.1005{\rm e}-05$ $*\;$ $\;0.00625\;$
$6.7526{\rm e}-07$ $1.9997$ $8.6798{\rm e}-07$ $1.9998$ $8.7511{\rm e}-06$
$1.9994$ $1.0254{\rm e}-05$ $1.9996\;$ $\;0.003125\;$ $1.6883{\rm e}-07$
$1.9999$ $2.1700{\rm e}-07$ $2.0000$ $2.1878{\rm e}-06$ $2.0000$ $2.5638{\rm
e}-06$ $1.9999\;$ $\;0.0015625\;$ $4.2309{\rm e}-08$ $1.9965$ $5.4154{\rm
e}-08$ $2.0025$ $5.4701{\rm e}-07$ $1.9998$ $6.4092{\rm e}-07$ $2.0000\;$
$\;0.00078125\;$ $1.0528{\rm e}-08$ $2.0066$ $1.3626{\rm e}-08$ $1.9907$
$1.3677{\rm e}-07$ $1.9998$ $1.6025{\rm e}-07$ $1.9998\;$

Table 5: Numerical conservative invariants of $E^{n}$, $H^{n}$, and $I^{n}$ at
time $t_{n}$.

Case A, $(h,\tau)=(1/5,1/256)$ $t_{n}\quad$ $\qquad E^{n}\qquad$ $\qquad
H^{n}\qquad$ $\qquad I^{n}\qquad$ $0\quad$ $6.426590811396586\quad$ $0\quad$
$12.39999498602724\quad$ $2\quad$ $6.426590811396584\quad$
$0.000000000000006\quad$ $12.39999498602724\quad$ $4\quad$
$6.426590811396586\quad$ $0.000000000000024\quad$ $12.39999498602725\quad$
$6\quad$ $6.426590811396584\quad$ $0.000000000000055\quad$
$12.39999498602725\quad$ $8\quad$ $6.426590811396582\quad$
$0.000000000000073\quad$ $12.39999498602724\quad$ $10\quad$
$6.426590811396588\quad$ $0.000000000000023\quad$ $12.39999498602725\quad$
Case B, $(h,\tau)=(1/5,1/256)$ $t_{n}\quad$ $\qquad E^{n}\qquad$ $\qquad
H^{n}\qquad$ $\qquad I^{n}\qquad$ $0\quad$ $67.70071037403311\quad$ $0\quad$
$91.40037694550523\quad$ $2\quad$ $67.70067103780045\quad$
$-0.000000000000133\quad$ $91.40037694550517\quad$ $4\quad$
$67.70064579547363\quad$ $-0.000000000000383\quad$ $91.40037694550520\quad$
$6\quad$ $67.70062829372516\quad$ $-0.000000000000820\quad$
$91.40037694550522\quad$ $8\quad$ $67.70061909000481\quad$
$-0.000000000000772\quad$ $91.40037694550519\quad$ $10\quad$
$67.70061445300891\quad$ $-0.000000000000737\quad$ $91.40037694550520\quad$
Case C, $(h,\tau)=(1/10,1/256)$ $t_{n}\quad$ $\qquad E^{n}\qquad$ $\qquad
H^{n}\qquad$ $\qquad I^{n}\qquad$ $0\quad$ $8.905545767953516\quad$
$0.001300209682121\quad$ $16.79999981448777\quad$ $2\quad$
$8.905545767953521\quad$ $0.001300209690680\quad$ $16.79999981448777\quad$
$4\quad$ $8.905545767953528\quad$ $0.001300209671917\quad$
$16.79999981448778\quad$ $6\quad$ $8.905545767953543\quad$
$0.001300209652087\quad$ $16.79999981448778\quad$ $8\quad$
$8.905545767953550\quad$ $0.001300209648438\quad$ $16.79999981448778\quad$
$10\quad$ $8.905545767953560\quad$ $0.001300209621099\quad$
$16.79999981448778\quad$ Case D, $(h,\tau)=(1/5,1/256)$ $t_{n}\quad$ $\qquad
E^{n}\qquad$ $\qquad H^{n}\qquad$ $\qquad I^{n}\qquad$ $0\quad$
$14.14719145331662\quad$ $0.002065791557752\quad$ $19.99999836196486\quad$
$2\quad$ $14.14719145331671\quad$ $0.002065792012461\quad$
$19.99999836196486\quad$ $4\quad$ $14.14719145331683\quad$
$0.002065792010409\quad$ $19.99999836196486\quad$ $6\quad$
$14.14719145331692\quad$ $0.002065792072498\quad$ $19.99999836196486\quad$
$8\quad$ $14.14719145331738\quad$ $0.002065794300180\quad$
$19.99999836196487\quad$ $10\quad$ $14.14719145331774\quad$
$0.002065796011144\quad$ $19.99999836196487\quad$

(a) Case A-$u$ (b) Case A-$\rho$ (c) Case D-$u$ (d) Case D-$\rho$

Figure 2: Unconditional convergence tests against spatial grid with the fixed
temporal stepsize $\tau$. Case A: (a) and (b); Case D: (c) and (d).

(a) Case E-$u$, view(45,70)  (b) Case E-$\rho$, view(45,70)  (c) Case E-$u$,
view(45,70)  (d) Case E-$\rho$, view(45,70)

(e) Case E-$u$, contour  (f) Case E-$\rho$, contour  (g) Case E-$u$, contour
(h) Case E-$\rho$, contour

Figure 3: The portraits of the velocity and altitude in long-time simulation
with same stepsizes $h=1/16$ and $\tau=1/20$ in Case E. Upper: (a) and (b)
calculated by the difference scheme (3.6), (c) and (d) calculated by the
difference scheme in ZLZ2022 ; Lower: corresponding contours.

(a) $u$, view(45,70) (b) $\rho$, view(45,70) (c) $u$, view(45,70) (d) $\rho$,
view(45,70)

(e) $u$-contour (f) $\rho$-contour (g) $u$-contour (h) $\rho$-contour

(i) Numerical invariants (j) Errors of the numerical invariants

Figure 4: The portraits of the velocity and altitude in long-time simulation
calculated by the difference scheme (3.6) (a)–(b), and three-level linearized
difference scheme (c)–(d) with same stepsizes $h=1/10$, $\tau=1/50$ in Case F;
(e)–(h) denote the corresponding contours; (i)–(j): Numerical invariants, and
relevant error curves with the same grid stepsizes for the difference scheme
(3.6).

Convergent accuracy test in Cases A–F. Table 3 displays the numerical errors
and spatial convergence orders when the temporal stepsize $\tau=1/50$ is fixed
for Case A and Case B. It is observed that the finite difference scheme is
quadratically convergent in the temporal direction. These data are consistent
with the theoretical result in Theorem 4.1 for case of $\Omega=0$. The second
part of Table 3 displays the numerical errors in space and convergence orders
in Cases C and D with a fixed and tiny temporal stepsize $\tau=1/1000$ such
that the temporal error may be negligible. Similar results are observed, which
are consistent with Theorem 4.2 for case of $\Omega\neq 0$.

Tables 5 displays the temporal error behavior and convergence orders with the
fixed spatial stepsize $h=6/25$ for Cases A and B, and $h=4/25$ for Cases C
and D, respectively. Clearly, the temporal convergence orders approach to two
whether $\Omega=0$ or $\Omega\neq 0$. These data are also consistent with
Theorems 4.1 and 4.2.

Invariants test in Cases A–D. We next study the long-time behaviors of
solutions for the difference scheme (3.6). Table 5 shows the discrete energy,
momentum, and mass for case of $\Omega=0$, in which the computations are
performed with $h=1/5$ and $\tau=1/256$. We observe that all the invariants
are conserved very well, even for long-time dynamics. Compared with the
results in ZLZ2022 , the new method preserves one more invariant: momentum. It
should be noticed that the computed energy is approximately half of that in
ZLZ2022 due to the different definition of the discrete energy.

When $\Omega=73\times 10^{-6}$, Table 5 also shows similar discrete invariants
for Cases C and D, which are consistent with Theorems 3.1. We also observe
from Table 5 that the discrete momentum is close to zero when $\Omega=0$,
while it is non-zero when $\Omega=73\times 10^{-6}$.

Unconditional convergence test in Cases A–D. To explore the unconditional
convergence of (3.6), we take Case A and Case D as examples for illustration.
Figures 2(a) and 2(b) demonstrate the convergence behavior of the velocity and
height in Case A with the spatial grid size reduction when the fixed temporal
stepsizes $\tau=5/16$, $5/32$, $5/64$, $5/128$, $5/256$ are used,
respectively. We see that no matter how the spatial stepsize varies, the
numerical errors always approach fixed values, which confirms that the
difference scheme (3.6) converges unconditionally (namely, no grid ratio
restriction. Otherwise, the numerical error would increase steeply with the
refined spatial grid). The similar phenomenon except minor differences is
observed in Case D with the fixed temporal stepsizes $\tau=1/16$, $1/32$,
$1/64$, $1/128$, $1/256$. Cases B and C also show a similar phenomenon, in
which unconditionally convergent figures are omitted for brevity.

Long time behavior test on large domain in Cases E–F. We first compare the
simulation effect of the difference scheme (3.6) with that in ZLZ2022 under
the same coarser grids $h=1/16,\;\tau=1/20$ (at this point,
${\tau}/{h}=0.8<1$, see Figure 3), and we see that the solution portraits in
the present paper are better than that generated by the difference scheme in
ZLZ2022 . To further study the performance of (3.6) and that in ZLZ2022 , we
conduct numerical simulations for Case F on a larger area for a long-time (at
the moment, $L=200$ and $T=1000$). As illustrated in Figure 4, the numerical
results show amazing periodicity for both schemes. However, the difference
scheme (3.6) owns better resolution even using refined stepsizes
$h=1/10,\;\tau=1/50$ (${\tau}/{h}=0.2$). Though the smooth initial data are
utilized, we observe that the occurrence of the bounded peakon solutions.
Figure 4(i) and Figure 4(j) denote numerical invariants of the difference
scheme (3.6) and their numerical errors, which demonstrates superior
simulation performance.

###### Example 2

(Peakon anti-peakon interaction) We consider a R2CH system with nonsmooth
initial values defined by

$u^{0}(x)={\rm e}^{-|x-5|}-{\rm e}^{-|x+5|},\quad\rho^{0}(x)=0.5.$

The problem depicts the interaction of peakons and anti-peakons. Special cases
of the problem have been studied in the early literature on the computational
domain $[-20,20]$ with $\Omega=0$, see e.g., CKL2020 ; CMR2014 ; LP2016 ;
YFS2018 . In the circumstances, we have $L=40$. We will solve the underlying
problem based on the difference scheme (3.6) for the following five cases,
especially the case of $\Omega\neq 0$, see Table 6 for the details. The
numerical results are shown in Figure 5 and Figure 6.

Table 6: The selected parameters in the numerical tests.

${\rm\qquad\qquad\qquad Parameters}$ Case A $\sigma=1$,
$\kappa=\mu=\Omega=0;\qquad$CKL2020 ; CMR2014 ; LP2016 ; YFS2018 Case B
$\sigma=1$, $\kappa=\mu=0$, $\Omega=0.2;\qquad$ Case C $\sigma=1$,
$\kappa=\mu=0$, $\Omega=73\times 10^{-6};\qquad$ Case D $\sigma=1$,
$\kappa=1$, $\mu=0$, $\Omega=73\times 10^{-6};\qquad$ Case E $\sigma=1$,
$\kappa=1$, $\mu=1$, $\Omega=73\times 10^{-6}.\qquad$

(a) Case A, $t=1$ (b) Case A, $t=3$ (c) Case A, $t=6$ (d) Case A, $t=8$

(e) Case B, $t=1$ (f) Case B, $t=3$ (g) Case B, $t=6$ (h) Case B, $t=8$

(i) Case C, $t=1$ (j) Case C, $t=3$ (k) Case C, $t=6$ (l) Case C, $t=8$

(m) Case D, $t=1$ (n) Case D, $t=3$ (o) Case D, $t=6$ (p) Case D, $t=8$

(q) Case E, $t=1$ (r) Case E, $t=3$ (s) Case E, $t=6$ (t) Case E, $t=8$

Figure 5: Velocities of the waves at different instants of times calculated by
the difference scheme (3.6) for five different Cases; the spatial grid
stepsize is fixed as $h=0.02$.

(a) Case A, $t=1$ (b) Case A, $t=3$ (c) Case A, $t=6$ (d) Case A, $t=8$

(e) Case B, $t=1$ (f) Case B, $t=3$ (g) Case B, $t=6$ (h) Case B, $t=8$

(i) Case C, $t=1$ (j) Case C, $t=3$ (k) Case C, $t=6$ (l) Case C, $t=8$

(m) Case D, $t=1$ (n) Case D, $t=3$ (o) Case D, $t=6$ (p) Case D, $t=8$

(q) Case E, $t=1$ (r) Case E, $t=3$ (s) Case E, $t=6$ (t) Case E, $t=8$

Figure 6: The altitudes of wave propagation at different instance of times
calculated by the difference scheme (3.6) in five different Cases; the spatial
grid stepsize is fixed as $h=0.02$

(a) Case A-$u$ (b) Case A-$\rho$ (c) Case A-Invariants (d) Errors

(e) Case B-$u$ (f) Case B-$\rho$ (g) Case B-Invariants (h) Errors

(i) Case C-$u$ (j) Case C-$\rho$ (k) Case C-Invariants (l) Errors

(m) Case D-$u$ (n) Case D-$\rho$ (o) Case D-Invariants (p) Errors

(q) Case E-$u$ (r) Case E-$\rho$ (s) Case E-Invariants (t) Errors

Figure 7: The evolution surfaces of the velocity $u$ (Column one) and the
altitude $\rho$ (Column two), and numerical invariants (Column three) and
corresponding numerical error curves of the invariants (Column four) for five
Cases with the stepsizes $h=1/10$ and $\tau=1/50$.

Figure 5 and Figure 6 for the benchmark problem in Case A depict the evolution
of the velocity and altitude, respectively, at $t=1,\,3,\,6,\,8$ for
$\Omega=0$, which recover the elastic collision, the case studied in the
literature, see also CKL2020 ; CMR2014 ; LP2016 ; YFS2018 . The second row in
Figure 5 and Figure 6 shows the numerical evolution when the rotation
parameter $\Omega=0.2$, respectively with the same spatial interval $[-20,20]$
for Case B. Compared the results of Case A with those of Case B, we see that
the nonzero rotation parameter has an important effect on the evolution of the
solution, i.e., the solution is no longer symmetric. Next, we decrease the
rotation parameter into the practical value $\Omega=73\times 10^{-6}$ in Case
C, which is a very tiny value. We observe that the solution is very similar to
that of Case A (the zero rotation parameter) in such a short time horizon.
Moreover, let the linear underlying shear flow $\kappa=1$ in Case D, we find
that the solution also loses the symmetry compared with the result in Case C.
Analogously, we change the dimensionless parameter $\mu$ into unit one in Case
E, the collision of two waves is also no longer elastic. Furthermore, all the
invariants are still preserved very well under the above five cases, see e.g.,
Figure 7. In addition, the evolution graphs for the velocity $u$ and the
altitude $\rho$ are clearly displayed in the first and second columns in
Figure 7. In a word, these numerical results indicate that if any of the
parameters $\Omega$, $\kappa$, and $\mu$ is nonzero, the solution will be
asymmetrical, which could be qualitatively verified from the perspective of
theoretical analysis.

##  6 Conclusions and outlooks

To summarize, we propose and analyze a nonlinear difference scheme for the
R2CH system based on a framework of the bilinear operator, and obtain several
new numerical results such as unconditional convergence and invariant-
preserving properties including energy, momentum, and mass. These ensure that
the numerical scheme provides an accurate long-time evolution of solitary
waves both in smooth and nonsmooth initial values. Regrettably, the present
paper also leaves some loose ends, which are addressed as follows.

  * •

When $\Omega\neq 0$, the convergence result is valid only for small initial
data ($c_{\rm max}<\frac{1}{2\Omega}-\kappa$) due to technical reasons.
Considering the practical case $\Omega\approx 73\cdot 10^{-6}$, this is a mild
restriction. Nevertheless, it would be desirable to remove this restriction.
Indeed, we performed some numerical tests with $c_{\rm
max}\geqslant\frac{1}{2\Omega}-\kappa$, and our numerical scheme is still
working. On the other hand, the numerical theory does break down for very
large $c_{\rm max}$.

  * •

From long-time numerical simulation (see e.g., Cases E–F in Example 1), we
observe clearly that even if the initial values of the R2CH system are smooth,
the solutions may evolve into rough or cuspidal in finite time. This
phenomenon is worthy of further study.

  * •

Our preliminarily numerical tests in Example 2 capture the evolution of the
R2CH system with a nonsmooth initial velocity based on the difference scheme
(3.6), the theoretical analysis is necessary to cover this case.

  * •

It is worth applying the framework of the bilinear operator to solve and
analyze other types of the shallow water wave problems.

## Conflict of interest

The authors of this paper have no conflict of interest to declare.

## Data availability

Data will be made available on reasonable request.

###### Acknowledgements.

The authors would like to thank Prof. Zhi-zhong Sun for most helpful
discussions and suggestions. Part of the work was finished during Qifeng’s
visit in École Polytechnique Fédérale de Lausanne, and he would like to thank
Prof. Jan S. Hesthaven for his hospitality in 2021-2022.

## References

  * (1) Akrivis, G.D.: Finite difference discretization of the cubic Schrödinger equation. IMA J. Numer. Anal. 13, 115–124 (1993)
  * (2) Antonopoulos, D.C., Dougalis, V.A., Mitsotakis, D.E.: Error estimates for Galerkin finite element methods for the Camassa–Holm equation. Numer. Math. 142, 833–862 (2019)
  * (3) Bressan, A., Constantin, A.: Global conservative solutions of the Camassa-Holm equation. Arch. Rational Mech. Anal. 183, 215–239 (2007)
  * (4) Camassa, R., Holm, D.: An integrable shallow water equation with peaked solitons. Phys. Rev. Lett. 71, 1661–1664 (1993)
  * (5) Camassa, R., Lee, L.: Complete integrable particle methods and the recurrence of initial states for a nonlinear shallow-water wave equation. J. Comput. Phys. 227, 7206–7221 (2008)
  * (6) Chen. R., Fan, L., Gao, H., Liu, Y.: Breaking waves and solitary waves to the rotation-two-component Camassa-Holm system. SIAM J. Math. Anal. 49, 3573–3602 (2017)
  * (7) Chen. R., Liu, Y.: Wave breaking and global existence for a generalized two-component Camassa-Holm system. Int. Math. Res. Notices 268, 45–66 (2011)
  * (8) Chertock, A., Kurganov, A., Liu, Y.: Finite-volume-particle methods for the two-component Camassa-Holm system. Commun. Comput. Phys. 27, 480–502 (2020)
  * (9) Coclite, G., Karlsen, K., Risebro, N.: A convergent finite difference scheme for the Camassa-Holm equation with general $H^{1}$ initial data. SIAM J. Numer. Anal. 46, 1554–1579 (2008)
  * (10) Cohen, D., Matsuo, T., Raynaud, X.: A multi-symplectic numerical integrator for the two-component Camassa-Holm equation. J. Nonlinear Math. Phys. 21, 442–453 (2014)
  * (11) Cohen. D., Raynaud, X.: Convergent numerical schemes for the compressible hyperelastic rod wave equation. Numer. Math. 122, 1–59 (2012)
  * (12) Cohen, D., Owren, B., Raynaud, X.: Multi-symplectic integration of the Camassa-Holm equation. J. Comput. Phys. 227, 5492–5512 (2008)
  * (13) Constantin, A., Ivanov, R.: On an integrable two-component Camassa-Holm shallow water system. Phys. Lett. A. 372, 7129–7132 (2008)
  * (14) Danchin, R.: A note on well-posedness for Camassa-Holm equation. J. Differential Equations 192, 429–444 (2003)
  * (15) David, C., Brynjulf, O., Xavier, R.: Multi-symplectic integration of the Camassa-Holm equation. J. Comput. Phys. 227, 5492–5512 (2008)
  * (16) Fan, E., Yuen, M.: Peakon weak solutions for the rotation-two-component Camassa-Holm system. Appl. Math. Lett. 97, 53–59 (2019)
  * (17) Fan, L., Gao, H., Liu, Y.: On the rotation-two-component Camassa-Holm system modelling the equatorial water waves. Adv. Math. 291, 59–89 (2016)
  * (18) Feng, K., Qin, M.: Hamiltonian algorithms for Hamiltonian dynamical systems. Progress in Natural Science, 1(2), 105–116 (1991)
  * (19) Fuchssteinert, B., Fokas, A.: Symplectic structures, their backlund transformations and hereditary symmetries. Phys. D Nonlinear Phenomena 4, 47–66 (1981)
  * (20) Galtung, S.T., Grunert, K.: A numerical study of variational discretizations of the Camassa–Holm equation. BIT Numer. Math. 61, 1271–1309 (2021)
  * (21) Geng, X., Xue, B.: A three-component generalization of Camassa-Holm equation with $N$-peakon solutions. Adv. Math. 226, 827–839 (2011)
  * (22) Guan, C., Yin, Z.: Global weak solutions for a two-component Camassa-Holm shallow water system. J. Funct. Anal. 260, 1132–1154 (2011)
  * (23) Gui, G., Liu, Y.: On the Cauchy problem for the two-component Camassa-Holm system. Math. Z. 268, 45–66 (2011)
  * (24) Guo, B.: A class of difference scheme for two-dimensional vorticity equations with viscous fluids. Acta Mathematics Sinica 17, 242–258 (1974)
  * (25) Han, Y., Guo, F., Gao, H.: On solitary waves and wave-breaking phenomena for a generalized two-component integrable Dullin-Gottwald-Holm system. J. Nonlinear Sci. 23, 617–656 (2013)
  * (26) Hesthaven, J.S., Pagliantini, C., Rozza, G.: Reduced basis methods for time-dependent problems. Acta. Numer. 31, 265–345 (2022)
  * (27) Henry, D.: Infinite propagation speed for a two component Camassa-Holm equation. Discrete Contin. Dyn. Syst. Ser. B 12, 597–606 (2009)
  * (28) Holden, H., Raynaud, X.: A convergent numerical scheme for the Camassa-Holm equation based on multipeakons. Discrete Cont. Dyn. A 14, 505–523 (2006)
  * (29) Holden, H., Raynaud, X.: Convergence of a finite difference scheme for the Camassa-Holm equation. SIAM J. Numer. Anal. 44, 1655–1680 (2006)
  * (30) Holden, H., Raynaud, X.: Periodic conservative solutions of the Camassa-Holm equation. Ann. Inst. Fourier (Grenoble) 58, 945–988 (2008)
  * (31) Holm, D., Ivanov, R.: Multi-component generalizations of the CH equation: geometrical aspects, peakons and numerical examples. J. Phys. A 43, 492001 (2010)
  * (32) Jiang, C., Gong, Y., Cai, W., Wang, Y.: A linearly implicit structure-preserving scheme for the Camassa-Holm equation based on multiple scalar auxiliary variables approach. J. Sci. Comput. 83, Article 20 (2020)
  * (33) Jiang, C., Wang, Y., Gong, Y.: Arbitrarily high-order energy-preserving schemes for the Camassa–Holm equation. Appl. Numer. Math. 151, 85–97 (2020)
  * (34) Kalisch, H., Raynaud, X.: Convergence of a spectral projection of the Camassa-Holm equation. Numer. Meth. Part. Differ. Equ. 22, 1197–1215 (2006)
  * (35) Kang, J., Liu, X., Qu, C.: On an integrable multi-component Camassa-Holm system arising from Möbius geometry. P. Roy. Soc. A Math. Phy. 477 (2021) https://doi.org/10.1098/rspa.2021.0164
  * (36) Kuo, P., Sanz-Serna, J.M.: Convergence of methods for the numerical solution of the Korteweg-de Vries equation. IMA J. Numer. Anal. 1, 215–221 (1981)
  * (37) Li, S., Vu-Quoc, L.: Finite difference calculus invariant structure of a class of algorithms for the nonlinear Klein-Gordon equation. SIAM J. Numer. Anal. 32, 1839–1875 (1995)
  * (38) Li, X., Qian, X., Zhang, B-Y., Song, S.: A multi-symplectic compact method for the two-component Camassa-Holm equation with singular solutions. Chin. Phys. Lett. 34, 090202 (2017)
  * (39) Li, N., Liu, Q., Popowicz, Z.: A four-component Camassa-Holm type hierarchy. J. Geom. Phys. 85, 29–39 (2014)
  * (40) Liu, H., Pendleton, T.: On invariant-preserving finite difference schemes for the Camassa-Holm equation and the two-component Camassa-Holm system. Commun. Comput. Phys. 19, 1015–1041 (2016)
  * (41) Liu, H., Xing, Y.: An invariant preserving discontinuous Galerkin method for the Camassa-Holm equation. SIAM J. Sci. Comput. 38, A1919–A1934 (2016)
  * (42) Liu, J., Pucci, P., Zhang, Q.: Wave breaking analysis for the periodic rotation-two-component Camassa-Holm system. Nonlinear Anal. 187, 214–228 (2019)
  * (43) Moon, B.: On the wave-breaking phenomena and global existence for the periodic rotation-two-component Camassa-Holm system. J. Math. Anal. Appl. 451, 84–101 (2017)
  * (44) Olver, P., Rosenau, P.: Tri-Hamiltonian duality between solitons and solitary-wave solutions having compact support. Phys. Rev. E. 53, 1900–1906 (1996)
  * (45) Wang, Z. Xiang, X.: Generalized Laguerre approximations and spectral method for the Camassa–Holm equation. IMA J. Numer. Anal. 35, 1456–1482 (2015)
  * (46) Xu, Y., Shu, C-W.: A local Discontinuous Galerkin method for the Camassa-Holm equation. SIAM J. Numer. Anal. 46, 1998–2021 (2008)
  * (47) Yang, M., Li, Y., Qiao, Z.: Persistence properties and wave-breaking criteria for a generalized two-component rotational b-family system. Discret. Contin. Dyn. Syst. A 40, 2475–2493 (2020)
  * (48) Yu, C-H., Feng, B-F., Sheu, T.W.H.: Numerical solutions to a two-component Camassa-Holm equation. J. Comput. Appl. Math. 336, 317–337 (2018)
  * (49) Yu, X., Ye, X., Wang, Z.: A fast solver of Legendre-Laguerre spectral element method for the Camassa–Holm equation. Numer. Algor. 88, 1–23 (2021)
  * (50) Zhang, L., Liu, B.: Well-posedness, blow-up criteria and gevrey regularity for a rotation-two-component Camassa-Holm system. Discret. Contin. Dyn. Syst. A 38, 2655–2685 (2018)
  * (51) Zhang, Q., Liu, L.: Convergence and stability in maximum norms of linearized fourth-order conservative compact scheme for Benjamin-Bona-Mahony-Burgers’ equation. J. Sci. Comput. 87, Article 59 (2021)
  * (52) Zhang, Q., Liu, L., Zhang, Z.: Linearly implicit invariant-preserving decoupled difference scheme for the rotation-two-component Camassa-Holm system. SIAM J. Sci. Comput. 44, A2226–A2252 (2022)
  * (53) Zhang, Y.: Wave breaking and global existence for the periodic rotation-Camassa-Holm system. Discret. Contin. Dyn. Syst. A 37, 2243–2257 (2017)
  * (54) Zhao, K., Wen, Z.: Effect of the Coriolis force on bounded traveling waves of the rotation-two-component Camassa-Holm system. Appl. Anal. (2021) https://doi.org/10.1080/00036811.2021.1965587
  * (55) Zhu, H., Song, S., Tang, Y.: Multi-symplectic wavelet collocation method for the nonlinear Schrödinger equation and the Camassa-Holm equation. Comput. Phys. Commun. 182, 616–627 (2011)