<EMAIL_ADDRESS>(T. Yan<EMAIL_ADDRESS>(J. Zhang),
<EMAIL_ADDRESS>(Q. Zhang)

#  Fully Conservative Difference Schemes for the Rotation-Two-Component
Camassa-Holm System with Smooth/Nonsmooth Initial Data111The research was
supported by Zhejiang Provincial Natural Science Foundation of China (Grant
No. LZ23A010007).

Tong Yan  1  Jiwei Zhang  2  and Qifeng Zhang  1  11affiliationmark:
Department of Mathematics, Zhejiang Sci-Tech University, Hangzhou, 310018,
P.R. China.
22affiliationmark:  School of Mathematics and Statistics, and Hubei Key
Laboratory of Computational Science, Wuhan University, Wuhan, 430072, P.R.
China

###### Abstract

The rotation-two-component Camassa–Holm system, which possesses strongly
nonlinear coupled terms and high-order differential terms, tends to have
continuous nonsmooth solitary wave solutions, such as peakons, stumpons,
composite waves and even chaotic waves. In this paper an accurate semi-
discrete conservative difference scheme for the system is derived by taking
advantage of its Hamiltonian invariants. We show that the semi-discrete
numerical scheme preserves at least three discrete conservative laws: mass,
momentum and energy. Furthermore, a fully discrete finite difference scheme is
proposed without destroying anyone of the conservative laws. Combining a
nonlinear iteration process and an efficient threshold strategy, the accuracy
of the numerical scheme can be guaranteed. Meanwhile, the difference scheme
can capture the formation and propagation of solitary wave solutions with
satisfying long time behavior under the smooth/nonsmooth initial data. The
numerical results reveal a new type of asymmetric wave breaking phenomenon
under the nonzero rotational parameter.

###### keywords:

R2CH system, semi-discrete scheme, discrete conservation law, peakon solution

65M06, 65M22, 35Q35, 35G25

##  1 Introduction

As one of the typical representatives of the nonlinearly dispersive partial
differential equations, the Camassa–Holm (CH) equation

$\displaystyle m_{t}+um_{x}+2mu_{x}=0,\quad{\rm with}\quad m=u-u_{xx}$ (1)

was proposed to simulate the evolution of shallow water waves [1]. Here $m$
represents the momentum related to the fluid velocity $u$. The CH equation
instead of the famous KdV equation allows traveling wave solutions in the
explicit form of

$u(x,t)=c^{\prime}{\rm e}^{-|x-ct|}$

with a sharp peak (peakons), where $c^{\prime}$ is the amplitude and $c$ is
the wave velocity. The CH equation has at least the following remarkable
properties: complete integrability, infinitely conservation laws and bi-
Hamiltonian structure [5, 10]. In this paper, we focus on the generalized two-
component case of the CH equation (1) given by

$\displaystyle m_{t}+\sigma(um_{x}+2mu_{x})=-3(1-\sigma)uu_{x}+Au_{x}-\mu
u_{xxx}-(1-2\Omega A)\rho\rho_{x}+2\Omega\rho(\rho u)_{x},\quad t\in(0,T],$
(2) $\displaystyle\rho_{t}+(\rho u)_{x}=0,\quad t\in(0,T]$ (3)

with the periodic boundary condition for $x\in{R}$. The system (2)–(3) is
called the rotation-two-component Camassa–Holm (R2CH) system.

The R2CH system was firstly proposed by Fan et al. [7] in 2016, which not only
inherits most of the properties of the solutions of the Camassa–Holm equation,
but also adds a new variable $\rho$ to describe the height of the water waves.
Moreover, it also introduces a constant rotational speed $\Omega$
($\Omega\in[0,1/4)$) to characterize the effect of the Earth’s rotation on the
shallow water waves. Therefore, the R2CH system can depict the propagation
phenomenon of shallow water waves more accurately.

Recalling $m=u-u_{xx}$, the R2CH system (2)–(3) can be rewritten as

$\displaystyle u_{t}-u_{xxt}-Au_{x}+3uu_{x}=\sigma(2u_{x}u_{xx}+uu_{xxx})-\mu
u_{xxx}-(1-2\Omega A)\rho\rho_{x}+2\Omega\rho(\rho u)_{x},\quad t\in(0,T],$
(4) $\displaystyle\rho_{t}+(\rho u)_{x}=0,\quad t\in(0,T].$ (5)

Several significant works have been developed in theory after its first
appearance. For example, blow-up scenarios of strong solutions are discussed
in [20, 14], blow-up criteria and wave-breaking phenomena are established in
[13], solitary waves with singularities, like peakons and cuspons are
introduced in [2], peakon weak solutions in distribution sense are considered
in [6], the persistence properties of the system in weighted $L^{p}$-spaces
are investigated in [17] and the effect of the Coriolis force on traveling
waves is studied in [21].

As we know, the system (4)–(5) has at least three conservative laws

$\displaystyle I_{1}(u,\rho)=\int_{R}\rho{\rm d}x,$ (6) $\displaystyle
I_{2}(u,\rho)=\int_{R}(u+\Omega\rho^{2}){\rm d}x,$ (7) $\displaystyle
E(u,\rho)=\frac{1}{2}\int_{R}(u^{2}+u_{x}^{2}+(1-2\Omega A)\rho^{2}){\rm d}x.$
(8)

In addition, for a special case of $\sigma=1$ and $\Omega=0$, the system
(4)–(5) has another three-order conservative law, which is described by

$\displaystyle H=\int_{R}(u^{3}+uu_{x}^{2}-Au^{2}-\mu u_{x}^{2}+u\rho^{2}){\rm
d}x,$

see also [3, 23]. The detailed derivation is postponed in the appendix.

While theoretical studies have been well carried out in many aspects,
numerical works for the R2CH system (4)–(5) are still few in the literature
and in urgent need for development. At present, the only numerical scheme is
finite difference discretization based on a bilinear operator, see e.g., [19].
However, the numerical simulation in [19] only covers the smooth initial
values with conditional stability. In the case of nonsmooth initial data, it
is not only theoretically difficult to analyze, but also difficult to capture
the nonsmooth solutions. In order to develop more accurate numerical
algorithms, this paper extends the excellent ideas in [16] to solve the system
(4)–(5).

Li and Vu-Quoc pointed out in [11]: “ _in some areas, the ability to preserve
some invariant properties of the original differential equation is a criterion
to judge the success of a numerical simulation_.” Hence, one of our principal
targets in deriving numerical scheme is to preserve the specially intrinsic
conservative structure of the system. Other targets include the study of novel
wave-breaking phenomena during the collision and evolution under the nonsmooth
initial data in a long time simulation.

The study of the high-order discrete conservation law is a very challenging
subject. As Liu and Xing pointed out in [12]: “ _it appears a rather difficult
task to preserve all conservation laws_.” In this paper, we make a tentative
attempt to simulate a third-order discrete quantity $H^{n}$ in the last
example in Section 5.3. It seems that $H^{n}$ is very close to the conserved
quantity $H$. However, a rigorous theoretical interpretation is still lacking
to determine whether the quantity $H^{n}$ is a discrete conservation law.

The main contributions of the paper are concretely summarized as follows.

  * •

The present difference scheme possesses at least three discrete conservative
laws, enabling it to accurately capture the behaviors of solutions to the R2CH
system under smooth/nonsmooth initial conditions.

  * •

The difference scheme demonstrates several new phenomena with nice resolution
for the first time, including the phenomena of the short-wave breaking and
interaction of the long-wave propagation.

  * •

The numerical accuracy of the difference scheme can be guaranteed and
oscillations of nonsmooth solutions can also be effectively eliminated by a
threshold technique.

  * •

The difference scheme improves the numerical results in the literature and has
potential impacts for predicting the propagation of solutions of other shallow
water wave equations.

The rest of the paper is arranged as follows. In Section 2, a semi-discrete
finite difference scheme is derived, which preserves the semi-discrete mass,
momentum and the energy. After that, a fully discrete finite difference scheme
is established in Section 3, which preserves all the conservative laws in the
discrete counterpart. Algorithm implementation in combination of a fixed point
iteration method is carried out in Section 4. Numerical examples including the
dam-break problem, singe peakon, multipeakon and peakon anti-peakon
interaction problems demonstrate the discrete conservative laws and nice
evolution of the solutions in Section 5 followed by a brief conclusion in
Section 6.

##  2 Semi-discrete scheme

To establish a conservative finite difference scheme for solving (2)–(3) with
the periodic boundary condition on the computational domain $[0,\,L]$, we
first discretize the system in space and derive a spatial semi-discrete
scheme, which is shown to preserve conservation laws. Let $h=L/M$ be the
spatial stepsize for a given integer $M$ and define the spatial variable by
$x_{i}=ih$, $i=0,1,\cdots,M$. Denote $u_{i}(t)=u(x_{i},t)$ and
$m_{i}(t)=m(x_{i},t)$. Put forward on this basis, we propose the semi-discrete
conservative finite difference scheme by

$\displaystyle\frac{{\rm d}}{{\rm
d}t}m_{i}(t)+\frac{\sigma}{2h}\big{[}(m_{i+1}(t)u_{i+1}(t)-m_{i-1}(t)u_{i-1}(t))+m_{i}(t)(u_{i+1}(t)-u_{i-1}(t))\big{]}$
$\displaystyle\quad=-\frac{3(1-\sigma)}{4h}\big{[}u_{i+1}^{2}(t)-u_{i-1}^{2}(t)\big{]}+\frac{A}{2h}\big{[}u_{i+1}(t)-u_{i-1}(t)\big{]}$
$\displaystyle\quad\quad-\frac{\mu}{{2h^{3}}}\big{[}u_{i+2}(t)-2u_{i+1}(t)+2u_{i-1}(t)-u_{i-2}(t)\big{]}-\frac{(1-2\Omega
A)}{{4h}}\big{[}\rho_{i+1}^{2}(t)-\rho_{i-1}^{2}(t)\big{]}$
$\displaystyle\quad\quad+\frac{\Omega}{2h}\rho_{i}(t)\big{[}(u_{i+1}(t)+u_{i}(t))(\rho_{i+1}(t)+\rho_{i}(t))-(u_{i-1}(t)+u_{i}(t))(\rho_{i-1}(t)+\rho_{i}(t))\big{]},$
(1) $\displaystyle\frac{{\rm d}}{{\rm
d}t}\rho_{i}(t)+\frac{1}{4h}\big{[}(u_{i+1}(t)+u_{i}(t))(\rho_{i+1}(t)+\rho_{i}(t))-(u_{i}(t)+u_{i-1}(t))(\rho_{i}(t)+\rho_{i-1}(t))\big{]}=0,$
(2) $\displaystyle
m_{i}(t)=u_{i}(t)-\frac{1}{2h}\big{[}v_{i+1}(t)-v_{i-1}(t)\big{]},$ (3)

where $v_{i}(t)=\frac{1}{2h}\big{[}u_{i+1}(t)-u_{i-1}(t)\big{]}$ and
$i=1,2,\cdots,M$.

With the periodicity of the boundary conditions for discrete mesh functions
$u_{i}$ and $\rho_{i}$, we have

$\displaystyle u_{i}=u_{i+M},\quad\rho_{i}=\rho_{i+M},\quad i=0,1,\cdots,M-1.$

Moreover, the periodicity for $u$ ensures the periodicity of the intermediate
variable $m$. Throughout the whole paper, we omit the spatial stepsize $h$ in
the semi-discrete and discrete conservative laws. On this basis, we define
semi-discrete mass, momentum, and energy as follows

$\displaystyle I_{1}(t)=\sum_{i=1}^{M}\rho_{i}(t),$ (4) $\displaystyle
I_{2}(t)=\sum_{i=1}^{M}\left[u_{i}(t)+\Omega\rho_{i}^{2}(t)\right],$ (5)
$\displaystyle
E(t)=\frac{1}{2}\sum_{i=1}^{M}\bigg{[}u_{i}^{2}(t)+\Big{(}\frac{u_{i+1}(t)-u_{i-1}(t)}{2h}\Big{)}^{2}+(1-2\Omega
A)\rho_{i}^{2}(t)\bigg{]}.$ (6)

Therefore, we show that the semi-discrete scheme (1)–(3) is conservative in
the following sense.

######  Theorem 2.1.

Consider the semi-discrete scheme (1)–(3) of the R2CH system (2)–(3). Then the
semi-discrete mass, momentum, and the energy with $\sigma=1$ are conservative
for the R2CH system (2)–(3) in the following sense

$\displaystyle\frac{{\rm d}}{{\rm d}t}I_{1}(t)=0,\quad\frac{{\rm d}}{{\rm
d}t}I_{2}(t)=0,\quad\frac{{\rm d}}{{\rm d}t}E(t)=0.$

######  Proof 2.2.

(I). In combination of (2) and (4), we have

$\displaystyle\frac{{\rm d}}{{\rm d}t}I_{1}(t)=\sum_{i=1}^{M}\frac{{\rm
d}}{{\rm
d}t}\rho_{i}(t)=\sum_{i=1}^{M}\frac{(u_{i-1}+u_{i})(\rho_{i-1}+\rho_{i})-(u_{i+1}+u_{i})(\rho_{i+1}+\rho_{i})}{4h}=0,$

where the periodicity of the boundary conditions for $u$ and $\rho$ is used.

(II). Similarly, we have

$\displaystyle\;\sum_{i=1}^{M}\frac{m_{i+1}(t)u_{i+1}(t)-m_{i-1}(t)u_{i-1}(t)}{2h}=0,$
$\displaystyle\;\sum_{i=1}^{M}m_{i}(t)\frac{u_{i+1}(t)-u_{i-1}(t)}{2h}=0,$
$\displaystyle\;\sum_{i=1}^{M}\frac{u_{i+1}^{2}(t)-u_{i-1}^{2}(t)}{4h}=0,$
$\displaystyle\;\sum_{i=1}^{M}\frac{u_{i+1}(t)-u_{i-1}(t)}{2h}=0,$
$\displaystyle\;\sum_{i=1}^{M}\frac{u_{i+2}(t)-2u_{i+1}(t)+2u_{i-1}(t)-u_{i-2}(t)}{2h^{3}}=0,$
$\displaystyle\;\sum_{i=1}^{M}\frac{\rho_{i+1}^{2}(t)-\rho_{i-1}^{2}(t)}{4h}=0.$

Thus, summing over $i$ from $1$ to $M$ on both sides of (1), we have

$\displaystyle\sum_{i=1}^{M}\frac{{\rm d}}{{\rm
d}t}m_{i}(t)=2\Omega\sum_{i=1}^{M}\rho_{i}(t)\frac{(u_{i+1}(t)+u_{i}(t))(\rho_{i+1}(t)+\rho_{i}(t))-(u_{i-1}(t)+u_{i}(t))(\rho_{i-1}(t)+\rho_{i}(t))}{4h}.$
(7)

Taking the derivative on both sides of (3) with respect to $t$ and summing
over $i$ from 1 to $M$, we have

$\displaystyle\sum_{i=1}^{M}\frac{{\rm d}}{{\rm
d}t}m_{i}(t)=\sum_{i=1}^{M}\left(\frac{{\rm d}}{{\rm
d}t}u_{i}(t)-\frac{v^{\prime}_{i+1}(t)-v^{\prime}_{i-1}(t)}{2h}\right)=\sum_{i=1}^{M}\frac{{\rm
d}}{{\rm d}t}u_{i}(t).$ (8)

Using (2), (7) and (8), we have

$\displaystyle\frac{{\rm d}}{{\rm d}t}I_{2}(t)=$
$\displaystyle\;\sum_{i=1}^{M}\frac{{\rm d}}{{\rm
d}t}u_{i}(t)+2\Omega\sum_{i=1}^{M}\rho_{i}(t)\frac{{\rm d}}{{\rm
d}t}\rho_{i}(t)$ $\displaystyle=$
$\displaystyle\;2\Omega\sum_{i=1}^{M}\rho_{i}(t)\frac{(u_{i+1}(t)+u_{i}(t))(\rho_{i+1}(t)+\rho_{i}(t))-(u_{i-1}(t)+u_{i}(t))(\rho_{i-1}(t)+\rho_{i}(t))}{4h}$
$\displaystyle\;+2\Omega\sum_{i=1}^{M}\rho_{i}(t)\frac{(u_{i-1}(t)+u_{i}(t))(\rho_{i-1}(t)+\rho_{i}(t))-(u_{i+1}(t)+u_{i}(t))(\rho_{i+1}(t)+\rho_{i}(t))}{4h}$
$\displaystyle=$ $\displaystyle\;0.$

(III). Finally, we show that $\frac{{\rm d}}{{\rm d}t}E(t)=0.$ By observing
(3), we have

$\displaystyle\frac{{\rm d}}{{\rm d}t}E(t)=$
$\displaystyle\;\sum_{i=1}^{M}\big{[}u_{i}(t)u^{\prime}_{i}(t)+v_{i}(t)v^{\prime}_{i}(t)+(1-2\Omega
A)\rho_{i}(t)\rho^{\prime}_{i}(t)\big{]}$ $\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\Big{[}u_{i}(t)\big{(}u^{\prime}_{i}(t)-\frac{v^{\prime}_{i+1}(t)-v^{\prime}_{i-1}(t)}{2h}\big{)}+(1-2\Omega
A)\rho_{i}(t)\rho^{\prime}_{i}(t)\Big{]}$ $\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}u_{i}(t)m^{\prime}_{i}(t)+(1-2\Omega
A)\sum_{i=1}^{M}\rho_{i}(t)\rho^{\prime}_{i}(t).$ (9)

Noticing that

$\displaystyle\;\sum_{i=1}^{M}u_{i}(t)\Big{(}\frac{m_{i+1}(t)u_{i+1}(t)-m_{i-1}(t)u_{i-1}(t)}{2h}+m_{i}(t)\frac{u_{i+1}(t)-u_{i-1}(t)}{2h}\Big{)}$
$\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\Big{(}\frac{m_{i+1}(t)u_{i+1}(t)u_{i}(t)-m_{i-1}(t)u_{i-1}(t)u_{i}(t)}{2h}+m_{i}(t)u_{i}(t)v_{i}(t)\Big{)}$
$\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\Big{(}\frac{m_{i}(t)u_{i}(t)u_{i-1}(t)-m_{i}(t)u_{i}(t)u_{i+1}(t)}{2h}+m_{i}(t)u_{i}(t)v_{i}(t)\Big{)}$
$\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\Big{(}-m_{i}(t)u_{i}(t)v_{i}(t)+m_{i}(t)u_{i}(t)v_{i}(t)\Big{)}=0.$
(10)

Eq. (9) follows from (1), (2) and (10) that

$\displaystyle\;\sum_{i=1}^{M}u_{i}(t)m^{\prime}_{i}(t)+(1-2\Omega
A)\sum_{i=1}^{M}\rho_{i}(t)\rho^{\prime}_{i}(t)$ $\displaystyle=$
$\displaystyle\;-\frac{\sigma}{2h}\sum_{i=1}^{M}u_{i}(t)\Big{[}\big{(}m_{i+1}(t)u_{i+1}(t)-m_{i-1}(t)u_{i-1}(t)\big{)}+m_{i}(t)\big{(}u_{i+1}(t)-u_{i-1}(t)\big{)}\Big{]}$
$\displaystyle\;-\frac{3(1-\sigma)}{4h}\sum_{i=1}^{M}\Big{[}u_{i}(t)(u_{i+1}^{2}(t)-u_{i-1}^{2}(t))\Big{]}+\frac{A}{2h}\sum_{i=1}^{M}\Big{[}u_{i}(t)(u_{i+1}(t)-u_{i-1}(t))\Big{]}$
$\displaystyle\;-\frac{\mu}{2h^{3}}\sum_{i=1}^{M}u_{i}(t)\Big{[}u_{i+2}(t)-2u_{i+1}(t)+2u_{i-1}(t)-u_{i-2}(t)\Big{]}$
$\displaystyle\;-\frac{(1-2\Omega
A)}{4h}\sum_{i=1}^{M}u_{i}(t)\Big{[}\rho_{i+1}^{2}(t)-\rho_{i-1}^{2}(t)\Big{]}$
$\displaystyle\;+2\Omega\sum_{i=1}^{M}u_{i}(t)\rho_{i}(t)\frac{(u_{i+1}(t)+u_{i}(t))(\rho_{i+1}(t)+\rho_{i}(t))-(u_{i-1}(t)+u_{i}(t))(\rho_{i-1}(t)+\rho_{i}(t))}{4h}$
$\displaystyle\;+(1-2\Omega
A)\sum_{i=1}^{M}\rho_{i}(t)\frac{(u_{i-1}(t)+u_{i}(t))(\rho_{i-1}(t)+\rho_{i}(t))-(u_{i+1}(t)+u_{i}(t))(\rho_{i+1}(t)+\rho_{i}(t))}{4h}$
$\displaystyle=$ $\displaystyle\;0,$

which completes the proof.

##  3 Fully discrete difference scheme

In the previous section, a semi-discrete finite difference scheme is
established for the R2CH system (2)–(3). Below, a new temporal discretization
is presented to guarantee the long time accurate calculation without
destroying the conservation properties of the original system. To this end,
let $\tau=T/N$ the time stepsize and $t^{n}=n\tau$ $(n=0,1,\cdots,N)$ the
partition of time variable with a given integer $N$. Denoting
$m_{i}^{n}=m(x_{i},t^{n})$, $u_{i}^{n}=u(x_{i},t^{n})$,
$\rho_{i}^{n}=\rho(x_{i},t^{n})$, then (1)–(3) can be discretized implicitly
by

$\displaystyle\frac{m_{i}^{*}-m_{i}^{n}}{\tau/2}+\frac{\sigma}{2h}\big{[}(m_{i+1}^{*}u_{i+1}^{*}-m_{i-1}^{*}u_{i-1}^{*})+m_{i}^{*}(u_{i+1}^{*}-u_{i-1}^{*})\big{]}$
$\displaystyle\quad=-\frac{3(1-\sigma)}{4h}\big{[}(u_{i+1}^{*})^{2}-(u_{i-1}^{*})^{2}\big{]}+\frac{A}{2h}\big{(}u_{i+1}^{*}-u_{i-1}^{*}\big{)}$
$\displaystyle\quad\quad-\frac{\mu}{2h^{3}}\big{(}u_{i+2}^{*}-2u_{i+1}^{*}+2u_{i-1}^{*}-u_{i-2}^{*}\big{)}-\frac{(1-2\Omega
A)}{4h}\big{[}(\rho_{i+1}^{*})^{2}-(\rho_{i-1}^{*})^{2}\big{]}$
$\displaystyle\quad\quad+\frac{\Omega}{2h}\rho_{i}^{*}\big{[}(u_{i+1}^{*}+u_{i}^{*})(\rho_{i+1}^{*}+\rho_{i}^{*})-(u_{i-1}^{*}+u_{i}^{*})(\rho_{i-1}^{*}+\rho_{i}^{*})\big{]},$
(1)
$\displaystyle\frac{\rho_{i}^{*}-\rho_{i}^{n}}{\tau/2}+\frac{1}{4h}\big{[}(u_{i+1}^{*}+u_{i}^{*})(\rho_{i+1}^{*}+\rho_{i}^{*})-(u_{i-1}^{*}+u_{i}^{*})(\rho_{i-1}^{*}+\rho_{i}^{*})\big{]}=0,$
(2) $\displaystyle
m_{i}^{n}=u_{i}^{n}-\frac{1}{2h}\big{(}v_{i+1}^{n}-v_{i-1}^{n}\big{)},\quad
v_{i}^{n}=\frac{1}{2h}(u_{i+1}^{n}-u_{i-1}^{n}),$ (3) $\displaystyle
u_{i}^{n+1}=2u_{i}^{*}-u_{i}^{n},\quad
m_{i}^{n+1}=2m_{i}^{*}-m_{i}^{n},\quad\rho_{i}^{n+1}=2\rho_{i}^{*}-\rho_{i}^{n},$
(4)

where $i=1,\cdots,M$ and $n=0,1,\cdots,N-1$. Thus, the discrete mass, momentum
and energy are defined as

$\displaystyle I_{1}^{n}=\sum_{i=1}^{M}\rho_{i}^{n},$ (5) $\displaystyle
I_{2}^{n}=\sum_{i=1}^{M}\left[u_{i}^{n}+\Omega(\rho_{i}^{n})^{2}\right],$ (6)
$\displaystyle
E^{n}=\frac{1}{2}\sum_{i=1}^{M}\bigg{[}(u_{i}^{n})^{2}+\Big{(}\frac{u_{i+1}^{n}-u_{i-1}^{n}}{2h}\Big{)}^{2}+(1-2\Omega
A)(\rho_{i}^{n})^{2}\bigg{]}.$ (7)

The newly established discrete scheme will be shown to preserve the above
three conserved quantities for the original system (4)–(5). Then we have the
following theorem.

######  Theorem 3.1.

Consider the difference scheme (1)–(4) of the R2CH system (2)–(3). Then the
discrete mass, momentum, and energy with $\sigma=1$ are conservative for the
R2CH system (2)–(3) in the following sense

$\displaystyle I_{1}^{n+1}=I_{1}^{n},\quad I_{2}^{n+1}=I_{2}^{n},\quad
E^{n+1}=E^{n}.$

######  Proof 3.2.

(I). Firstly, we show that $I_{1}^{n+1}=I_{1}^{n}$. Noticing (2) and (4), it
readily knows that

$\displaystyle\frac{\rho_{i}^{n+1}-\rho_{i}^{n}}{\tau}+\frac{(u_{i+1}^{*}+u_{i}^{*})(\rho_{i+1}^{*}+\rho_{i}^{*})-(u_{i-1}^{*}+u_{i}^{*})(\rho_{i-1}^{*}+\rho_{i}^{*})}{4h}=0.$
(8)

Summing (8) with respect to $i$ from 1 to $M$ and utilizing the periodicity,
we have

$\displaystyle\;\sum_{i=1}^{M}\Big{[}\frac{\rho_{i}^{n+1}-\rho_{i}^{n}}{\tau}+\frac{(u_{i+1}^{*}+u_{i}^{*})(\rho_{i+1}^{*}+\rho_{i}^{*})-(u_{i-1}^{*}+u_{i}^{*})(\rho_{i-1}^{*}+\rho_{i}^{*})}{4h}\Big{]}$
$\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\frac{\rho_{i}^{n+1}-\rho_{i}^{n}}{\tau}=\;\frac{I_{1}^{n+1}-I_{1}^{n}}{\tau}=\;0,$

which implies $I_{1}^{n+1}=I_{1}^{n}$.

(II). Next we show $I_{2}^{n+1}=I_{2}^{n}$. In combination of (1) and (3), and
noticing the periodicity, we have

$\displaystyle\;\sum_{i=1}^{M}\frac{u_{i}^{n+1}-u_{i}^{n}}{\tau}$
$\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\frac{m_{i}^{n+1}-m_{i}^{n}}{\tau}+\sum_{i=1}^{M}\frac{(v_{i+1}^{n+1}-v_{i-1}^{n+1})-(v_{i+1}^{n}-v_{i-1}^{n})}{2h\tau}$
$\displaystyle=$
$\displaystyle\;\Omega\sum_{i=1}^{M}\rho_{i}^{*}\frac{(u_{i+1}^{*}+u_{i}^{*})(\rho_{i+1}^{*}+\rho_{i}^{*})-(u_{i-1}^{*}+u_{i}^{*})(\rho_{i-1}^{*}+\rho_{i}^{*})}{2h}.$
(9)

Again using (8), we find that

$\displaystyle\;\sum_{i=1}^{M}\frac{(\rho_{i}^{n+1})^{2}-(\rho_{i}^{n})^{2}}{\tau}=2\sum_{i=1}^{M}\frac{\rho_{i}^{n+1}+\rho_{i}^{n}}{2}\cdot\frac{\rho_{i}^{n+1}-\rho_{i}^{n}}{\tau}$
$\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\rho_{i}^{*}\frac{(u_{i-1}^{*}+u_{i}^{*})(\rho_{i-1}^{*}+\rho_{i}^{*})-(u_{i+1}^{*}+u_{i}^{*})(\rho_{i+1}^{*}+\rho_{i}^{*})}{2h}.$
(10)

Combining (9) with (10), we have

$\displaystyle\frac{I_{2}^{n+1}-I_{2}^{n}}{\tau}=\sum_{i=1}^{M}\frac{u_{i}^{n+1}-u_{i}^{n}}{\tau}+\Omega\sum_{i=1}^{M}\frac{(\rho_{i}^{n+1})^{2}-(\rho_{i}^{n})^{2}}{\tau}=0,$

which indicates that $I_{2}^{n+1}=I_{2}^{n}$.

(III). Finally, we prove $E^{n+1}=E^{n}$. Using (10), we have

$\displaystyle\sum_{i=1}^{M}u_{i}^{*}\Big{(}\frac{m_{i+1}^{*}u_{i+1}^{*}-m_{i-1}^{*}u_{i-1}^{*}}{2h}+m_{i}^{*}\frac{u_{i+1}^{*}-u_{i-1}^{*}}{2h}\Big{)}=0.$

Using the periodicity, we have

$\displaystyle\sum_{i=1}^{M}\rho_{i}^{*}\frac{(u_{i-1}^{*}+u_{i}^{*})(\rho_{i-1}^{*}+\rho_{i}^{*})-(u_{i+1}^{*}+u_{i}^{*})(\rho_{i+1}^{*}+\rho_{i}^{*})}{4h}=\sum_{i=1}^{M}(\rho_{i}^{*})^{2}\frac{u_{i-1}^{*}-u_{i+1}^{*}}{4h}.$
(11)

Multiplying (1) with $u_{i}^{*}$ and summing over $i$ from 1 to $M$, and
rearranging the corresponding result, we have

$\displaystyle 0=$
$\displaystyle\;\sum_{i=1}^{M}\frac{m_{i}^{n+1}-m_{i}^{n}}{\tau}u_{i}^{*}+(1-2\Omega
A)\sum_{i=1}^{M}\frac{(\rho_{i+1}^{*})^{2}-(\rho_{i-1}^{*})^{2}}{4h}u_{i}^{*}$
$\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\bigg{(}\frac{u_{i}^{n+1}-u_{i}^{n}}{\tau}u_{i}^{*}-\frac{v_{i+1}^{n+1}-v_{i-1}^{n+1}-v_{i+1}^{n}+v_{i-1}^{n}}{2h\tau}u_{i}^{*}\bigg{)}+(1-2\Omega
A)\sum_{i=1}^{M}(\rho_{i}^{*})^{2}\frac{u_{i-1}^{*}-u_{i+1}^{*}}{4h}$
$\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\bigg{[}\frac{u_{i}^{n+1}-u_{i}^{n}}{\tau}u_{i}^{*}+\Big{(}\frac{v_{i}^{n+1}-v_{i}^{n}}{\tau}\Big{)}\Big{(}\frac{u_{i+1}^{*}-u_{i-1}^{*}}{2h}\Big{)}\bigg{]}$
$\displaystyle\;+(1-2\Omega
A)\sum_{i=1}^{M}\rho_{i}^{*}\frac{(u_{i-1}^{*}+u_{i}^{*})(\rho_{i-1}^{*}+\rho_{i}^{*})-(u_{i+1}^{*}+u_{i}^{*})(\rho_{i+1}^{*}+\rho_{i}^{*})}{4h}$
$\displaystyle=$
$\displaystyle\;\sum_{i=1}^{M}\bigg{[}\Big{(}\frac{u_{i}^{n+1}-u_{i}^{n}}{\tau}\Big{)}\Big{(}\frac{u_{i}^{n+1}+u_{i}^{n}}{2}\Big{)}+\Big{(}\frac{v_{i}^{n+1}-v_{i}^{n}}{\tau}\Big{)}\Big{(}\frac{v_{i}^{n+1}+v_{i}^{n}}{2}\Big{)}\bigg{]}$
$\displaystyle\;+(1-2\Omega
A)\sum_{i=1}^{M}\frac{\rho_{i}^{n+1}+\rho_{i}^{n}}{2}\cdot\frac{\rho_{i}^{n+1}-\rho_{i}^{n}}{\tau}$
$\displaystyle=$ $\displaystyle\;\frac{E^{n+1}-E^{n}}{\tau},$

where the periodicity, (4) and $\sigma=1$ are used in the first equality, the
periodicity and (3) are used in the second equality, the periodicity and (11)
are used in the third equality, and (8) is used in the penultimate equality.
This completes the proof.

##  4 Algorithm implementation

Denote $\boldsymbol{u}=(u_{1},u_{2},\cdots,u_{M})^{T}$,
$\boldsymbol{m}=(m_{1},m_{2},\cdots,m_{M})^{T}$ and
$\boldsymbol{\rho}=(\rho_{1},\rho_{2},\cdots,\rho_{M})^{T}$. Recalling (3), we
have a linear system of equations $\boldsymbol{m}=\boldsymbol{B}\textbf{u}$,
where $\boldsymbol{B}$ is a symmetric circulant matrix defined by

$\textbf{B}={\rm\textbf{circ}}(c[0],c[1],\cdots,c[M-1])$

with $c[0]=1+\frac{1}{2h^{2}}$, $c[1]=0$, $c[2]=-\frac{1}{4h^{2}}$,
$c[3]=\cdots=c[M-3]=0$, $c[M-2]=-\frac{1}{4h^{2}}$ and $c[M-1]=0$. Obviously,
$\boldsymbol{B}$ is a dominant diagonal matrix, which allows us to express the
vector $\boldsymbol{u}$ in terms of $\boldsymbol{m}$ as
$\boldsymbol{u}=\boldsymbol{B}^{-1}\boldsymbol{m}$. To avoid ambiguity we let
$(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i}$ denote the $i$th element of the
vector. Consequently, (1)–(4) can be rewritten equivalently in terms of
$m^{*}$ as follows

$\displaystyle\frac{m_{i}^{*}-m_{i}^{n}}{\tau/2}+\frac{\sigma}{2h}\big{[}(m_{i+1}^{*}(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i+1}-m_{i-1}^{*}(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i-1})+m_{i}^{*}((\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i+1}-(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i-1})\big{]}$
$\displaystyle\quad=-\frac{3(1-\sigma)}{4h}\big{[}(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i+1}^{2}-(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i-1}^{2}\big{]}+\frac{A}{2h}\big{[}(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i+1}-(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i-1}\big{]}$
$\displaystyle\quad\quad-\frac{\mu}{2h^{3}}\big{[}(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i+2}-2(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i+1}+2(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i-1}-(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i-2}\big{]}$
$\displaystyle\quad\quad-\frac{(1-2\Omega
A)}{4h}\big{[}(\rho_{i+1}^{*})^{2}-(\rho_{i-1}^{*})^{2}\big{]}+\frac{\Omega}{2h}\rho_{i}^{*}\big{[}((\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i+1}+(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i})(\rho_{i+1}^{*}+\rho_{i}^{*})$
$\displaystyle\quad\quad-((\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i-1}+(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i})(\rho_{i-1}^{*}+\rho_{i}^{*})\big{]},$
(1)
$\displaystyle\frac{\rho_{i}^{*}-\rho_{i}^{n}}{\tau/2}+\frac{1}{4h}\big{[}((\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i+1}+(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i})(\rho_{i+1}^{*}+\rho_{i}^{*})$
$\displaystyle\quad\quad-((\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i-1}+(\boldsymbol{B}^{-1}\boldsymbol{m}^{*})_{i})(\rho_{i-1}^{*}+\rho_{i}^{*})\big{]}=0,$
(2)
$\displaystyle\;m_{i}^{n+1}=2m_{i}^{*}-m_{i}^{n},\quad\rho_{i}^{n+1}=2\rho_{i}^{*}-\rho_{i}^{n}.$
(3)

At present, we have a nonlinear system of equations with variables
$\boldsymbol{m}^{*}$ and $\boldsymbol{\rho}^{*}$ only, which can be solved by
a fixed point iteration method.

The algorithm flow chart of the scheme (1)–(3) is given as follows.
Specifically, to solve the solutions at the $(n+1)$th time level, we suppose
that
$\boldsymbol{u}^{n,l}=\Big{(}u_{1}^{n,l},u_{2}^{n,l},\cdots,u_{M}^{n,l}\Big{)}^{T}$,
$\boldsymbol{\rho}^{n,l}=\Big{(}\rho_{1}^{n,l},\rho_{2}^{n,l},\cdots,\rho_{M}^{n,l}\Big{)}^{T}$
and
$\boldsymbol{m}^{n,l}=\Big{(}m_{1}^{n,l},m_{2}^{n,l},\cdots,m_{M}^{n,l}\Big{)}^{T}$
have been determined, then Algorithm 1 is utilized to solve the scheme
(1)–(3).

Algorithm 1: The iterative process for solving (1)–(3)

1. | Set the tolerance error _tol_.
---|---
2. | Give the initial guess $\boldsymbol{u}^{*,l}:=\boldsymbol{u}^{n,l}$.
3. | For $n=1,2,\cdots,N$ Do:
4. |  Compute $\boldsymbol{\rho}^{*,l+1}$, $\boldsymbol{m}^{*,l+1}$, then obtain $\boldsymbol{u}^{*,l+1}$
5. | EndDo
6. | If $\|\boldsymbol{u}^{*,l}-\boldsymbol{u}^{*,l+1}\|_{\infty}>\emph{tol}$, $\boldsymbol{u}^{*,l}:=\boldsymbol{u}^{*,l+1}$, GoTo 2;
7. | else, $\boldsymbol{u}^{n+1}:=2\boldsymbol{u}^{*,l+1}-\boldsymbol{u}^{n,l}$; $\boldsymbol{\rho}^{n+1}:=2\boldsymbol{\rho}^{*,l+1}-\boldsymbol{\rho}^{n,l}$, $\boldsymbol{m}^{n+1}=\boldsymbol{Bu}^{n+1}$ .

######  Remark 4.1.

Provided the solution of the system (4)–(5) is less regular, then numerical
oscillation will occur during the calculation. To eliminate this phenomenon,
an efficient strategy is utilized by adding local numerical viscosity. The
following two viscous terms

$\displaystyle\;R_{i}^{u}=\frac{\epsilon_{i}^{u}}{2h}(u_{i+1}^{*}-2u_{i}^{*}+u_{i-1}^{*}),\quad\text{and}\quad\;R_{i}^{\rho}=\frac{\epsilon_{i}^{\rho}}{2h}(\rho_{i+1}^{*}-2\rho_{i}^{*}+\rho_{i-1}^{*}),$

where

$\displaystyle\epsilon_{i}^{u}=\left\\{\begin{array}[]{ll}1,\quad|u_{i+1}^{*}-2u_{i}^{*}+u_{i-1}^{*}|\geqslant\varepsilon
h,\\\ 0,\quad
otherwise,\end{array}\right.\quad\text{and}\quad\epsilon_{i}^{\rho}=\left\\{\begin{array}[]{ll}1,\quad|\rho_{i+1}^{*}-2\rho_{i}^{*}+\rho_{i-1}^{*}|\geqslant\varepsilon
h,\\\ 0,\quad otherwise.\end{array}\right.$

are added to the right-hand side of (1) and (2), respectively. The above
factor $\varepsilon$ is a threshold (usually very small, e.g.,
$\varepsilon=10^{-5}$) which can determine where the slope of the solution
becomes unbounded. It is worth mentioning that the factor $\varepsilon$ varies
from case to case in the simulation of the discontinuous solutions in the
section 5.2 below.

##  5 Numerical tests

In this section, several examples for the benchmark problems are provided to
verify the convergence, conservation laws and the performance of the proposed
scheme (1)–(4). To test the convergence orders, the posterior error estimate
is utilized in temporal direction and spatial direction. To be more specific,
for sufficient small $\tau$, we denote

$\displaystyle\|{\rm F}_{u}(h)\|_{\infty}=\max\limits_{1\leqslant i\leqslant
M,1\leqslant n\leqslant
N}\big{|}u_{i}^{n}(h,\tau)-u_{2i}^{n}(h/2,\tau)\big{|},\quad{\rm
Ord}_{\infty}^{h}=\log_{2}\frac{\|{\rm F}_{u}(h)\|_{\infty}}{\|{\rm
F}_{u}(h/2)\|_{\infty}},$

and for sufficient small $h$, denote

$\displaystyle\|{\rm G}_{u}(\tau)\|_{\infty}=\max\limits_{1\leqslant
i\leqslant M,1\leqslant n\leqslant
N}\big{|}u_{i}^{n}(h,\tau)-u_{i}^{2n}(h,\tau/2)\big{|},\quad{\rm
Ord}_{\infty}^{\tau}=\log_{2}\frac{\|{\rm G}_{u}(\tau)\|_{\infty}}{\|{\rm
G}_{u}(\tau/2)\|_{\infty}}.$

###  5.1 Part I: smooth initial data

######  Example 5.1 (Dam-break problem [16, 19]).

The R2CH system with the smooth initial data

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

are considered, where $a$ is a dam-breaking parameter. The exact solution for
the problem is unknown.

To test the convergence orders and conservation, we consider the following
four cases:

  * •

Case (I). $a=0.1$, $A=0$, $\mu=0$, $\sigma=1$, $\Omega=0$, on the domain
$[-6,6]\times[0,20]$.

  * •

Case (II). $a=4$, $A=0$, $\mu=0$, $\sigma=1$, $\Omega=0$, on the domain
$[-12\pi,12\pi]\times[0,2]$.

  * •

Case (III). $a=0.1$, $A=0.1$, $\mu=0.1$, $\sigma=1$, $\Omega=73\times
10^{-6}$, on the domain $[-8,8]\times[0,1]$.

  * •

Case (IV). $a=4$, $A=1$, $\mu=1$, $\sigma=1$, $\Omega=73\times 10^{-6}$, on
the domain $[-12\pi,12\pi]\times[0,2]$.

(Convergence) We first verify the convergence orders for the above four cases.
The temporal convergence orders of velocity $u(x,t)$ and height $\rho(x,t)$
are respectively listed in Table 1 and Table 2, which show the second-order
convergence by fixing spatial grid $M=100$ and refining $N$. Similarly, Tables
3–4 illustrate the second-order spatial convergence order by fixing $N=4000$
and refining $M$.

(Conservation) The conservations in Theorem 3.1 are demonstrated in Table 5,
which shows our scheme indeed preserves three conserved quantities defined in
(5)–(7). It is worth noting that the total momentum in Cases (I) and (II) is
zero, while it does not vanish in Cases (III) and (IV) since $\Omega$ is
nonzero in (6).

(Portraits of solutions at different instants) In addition, an applicable
reference solution is necessary to evaluate the behavior of the numerical
solution. To generate a reference solution, we take a refined grid $M=3200$ in
space. Figures 1–2 show the behaviors of the predicted velocities $u(x,t)$ and
heights $\rho(x,t)$ at time $t=2$ for Cases (II) and (IV). We observe that the
solutions in Case (II) are symmetric, while the solutions in Case (IV) are
asymmetric. Actually the symmetry depends heavily on the selection of
parameters. It can be clearly seen that the conservative scheme (1)–(4)
performs well in depicting the dam break solution of the R2CH system (2)–(3)
even using a relatively rough grid. Figures 3–4 show the evolution of the
predicted dam break solutions of the R2CH system over a long period of time.
We find that even up to $t=50$, the solution in Figure 3 still stretches
leftward and rightward in a symmetric form, while the evolution trend of the
solutions in Figure 4 is clearly less symmetric.

Table 1: Numerical errors, temporal convergence orders of velocity $u(x,t)$ of
the scheme (1)–(4) with $M=100$ for Case (I), Case (II), Case (III) and Case
(IV).

Case (I) Case (II) Case (III) Case (IV) $\quad N$ $\|{\rm
F}_{u}(\tau)\|_{\infty}$ ${{\rm Ord}_{\infty}^{\tau}}$ $\|{\rm
F}_{u}(\tau)\|_{\infty}$ ${{\rm Ord}_{\infty}^{\tau}}$ $\|{\rm
F}_{u}(\tau)\|_{\infty}$ ${{\rm Ord}_{\infty}^{\tau}}$ $\|{\rm
F}_{u}(\tau)\|_{\infty}$ ${{\rm Ord}_{\infty}^{\tau}}$ $\quad 100$ $9.0981{\rm
e}-04$ $*$ $3.3262{\rm e}-04$ $*$ $1.7801{\rm e}-07$ $*$ $8.3334{\rm e}-04$
$*$ $\quad 200$ $2.2919{\rm e}-04$ $1.9890$ $8.3267{\rm e}-05$ $1.9982$
$4.4503{\rm e}-08$ $2.0000$ $2.0881{\rm e}-04$ $1.9967$ $\quad 400$
$5.7402{\rm e}-05$ $1.9973$ $2.0827{\rm e}-05$ $1.9999$ $1.1126{\rm e}-08$
$2.0000$ $5.2245{\rm e}-05$ $1.9988$ $\quad 800$ $1.4357{\rm e}-05$ $1.9993$
$5.1994{\rm e}-06$ $1.9997$ $2.7815{\rm e}-09$ $2.0000$ $1.3056{\rm e}-05$
$2.0006$ $\quad 1600$ $3.5897{\rm e}-06$ $1.9998$ $1.3006{\rm e}-06$ $1.9995$
$6.9536{\rm e}-10$ $2.0000$ $3.2636{\rm e}-06$ $2.0001$

Table 2: Numerical errors, temporal convergence orders in height $\rho(x,t)$
of the scheme (1)–(4) with $M=100$ for Case (I), Case (II), Case (III) and
Case (IV).

Case (I) Case (II) Case (III) Case (IV) $\quad N$ $\|{\rm
F}_{\rho}(\tau)\|_{\infty}$ ${{\rm Ord}_{\infty}^{\tau}}$ $\|{\rm
F}_{\rho}(\tau)\|_{\infty}$ ${{\rm Ord}_{\infty}^{\tau}}$ $\|{\rm
F}_{\rho}(\tau)\|_{\infty}$ ${{\rm Ord}_{\infty}^{\tau}}$ $\|{\rm
F}_{\rho}(\tau)\|_{\infty}$ ${{\rm Ord}_{\infty}^{\tau}}$ $\quad 100$
$1.4094{\rm e}-03$ $*$ $3.0830{\rm e}-04$ $*$ $5.2605{\rm e}-07$ $*$
$3.7110{\rm e}-04$ $*$ $\quad 200$ $3.5345{\rm e}-04$ $1.9955$ $7.7121{\rm
e}-05$ $1.9991$ $1.3151{\rm e}-07$ $2.0000$ $9.2969{\rm e}-05$ $1.9970$ $\quad
400$ $8.8512{\rm e}-05$ $1.9976$ $1.9283{\rm e}-05$ $1.9997$ $3.2879{\rm
e}-08$ $2.0000$ $2.3258{\rm e}-05$ $1.9990$ $\quad 800$ $2.2133{\rm e}-05$
$1.9997$ $4.8210{\rm e}-06$ $1.9999$ $8.2193{\rm e}-09$ $2.0001$ $5.8135{\rm
e}-06$ $2.0002$ $\quad 1600$ $5.5335{\rm e}-06$ $1.9999$ $1.2060{\rm e}-06$
$1.9990$ $2.0549{\rm e}-09$ $2.0000$ $1.4532{\rm e}-06$ $2.0001$

Table 3: Numerical errors, spatial convergence orders in velocity $u(x,t)$ of
the scheme (1)–(4) with $N=4000$ for Case (I), Case (II), Case (III) and Case
(IV).

Case (I) Case (II) Case (III) Case (IV) $\quad M$ $\|{\rm
F}_{u}(h)\|_{\infty}$ ${{\rm Ord}_{\infty}^{h}}$ $\|{\rm F}_{u}(h)\|_{\infty}$
${{\rm Ord}_{\infty}^{h}}$ $\|{\rm F}_{u}(h)\|_{\infty}$ ${{\rm
Ord}_{\infty}^{h}}$ $\|{\rm F}_{u}(h)\|_{\infty}$ ${{\rm Ord}_{\infty}^{h}}$
$\quad 100$ $6.2259{\rm e}-04$ $*$ $1.2876{\rm e}-01$ $*$ $1.9975{\rm e}-04$
$*$ $1.1547{\rm e}-01$ $*$ $\quad 200$ $1.5547{\rm e}-04$ $2.0016$ $3.7311{\rm
e}-02$ $1.7871$ $4.9396{\rm e}-05$ $2.0157$ $4.6363{\rm e}-02$ $1.3165$ $\quad
400$ $3.8856{\rm e}-05$ $2.0004$ $1.2058{\rm e}-02$ $1.6296$ $1.2327{\rm
e}-05$ $2.0025$ $1.4034{\rm e}-02$ $1.7240$ $\quad 800$ $9.7152{\rm e}-06$
$1.9998$ $3.3317{\rm e}-03$ $1.8557$ $3.0826{\rm e}-06$ $1.9997$ $3.7075{\rm
e}-03$ $1.9205$ $\quad 1600$ $2.4288{\rm e}-06$ $2.0000$ $8.5597{\rm e}-04$
$1.9606$ $7.7051{\rm e}-07$ $2.0003$ $9.3181{\rm e}-04$ $1.9923$

Table 4: Numerical errors, spatial convergence orders in height $\rho(x,t)$ of
the scheme (1)–(4) with $N=4000$ for Case (I), Case (II), Case (III) and Case
(IV).

Case (I) Case (II) Case (III) Case (IV) $\quad M$ $\|{\rm
F}_{\rho}(h)\|_{\infty}$ ${{\rm Ord}_{\infty}^{h}}$ $\|{\rm
F}_{\rho}(h)\|_{\infty}$ ${{\rm Ord}_{\infty}^{h}}$ $\|{\rm
F}_{\rho}(h)\|_{\infty}$ ${{\rm Ord}_{\infty}^{h}}$ $\|{\rm
F}_{\rho}(h)\|_{\infty}$ ${{\rm Ord}_{\infty}^{h}}$ $\quad 100$ $8.9486{\rm
e}-04$ $*$ $1.4496{\rm e}-01$ $*$ $2.4458{\rm e}-04$ $*$ $2.4052{\rm e}-01$
$*$ $\quad 200$ $2.2263{\rm e}-04$ $2.0070$ $8.6727{\rm e}-02$ $0.7411$
$6.0340{\rm e}-05$ $2.0191$ $1.8170{\rm e}-01$ $0.4046$ $\quad 400$
$5.5586{\rm e}-05$ $2.0018$ $3.4910{\rm e}-02$ $1.3129$ $1.5037{\rm e}-05$
$2.0046$ $7.7859{\rm e}-02$ $1.2226$ $\quad 800$ $1.3894{\rm e}-05$ $2.0003$
$1.0946{\rm e}-02$ $1.6732$ $3.7620{\rm e}-06$ $1.9990$ $2.5178{\rm e}-02$
$1.6289$ $\quad 1600$ $3.4736{\rm e}-06$ $1.9999$ $2.9327{\rm e}-03$ $1.9001$
$9.4032{\rm e}-07$ $2.0003$ $7.0786{\rm e}-03$ $1.8305$

Table 5: Numerical conserved quantities $E^{n}$, $I_{1}^{n}$ and $I_{2}^{n}$
at time $t^{n}$.

Case (I), $h=0.06$, $\tau=0.01$ $\quad t_{n}\quad$ $E^{n}$ $I_{1}^{n}$
$I_{2}^{n}$ $\quad 0\quad$ $107.1098478543294\quad$ $206.6665840981688\quad$
$0$ $\quad 2\quad$ $107.1098478543293\quad$ $206.6665840981688\quad$
$-0.000000000000005\;$ $\quad 4\quad$ $107.1098478543292\quad$
$206.6665840981687\quad$ $-0.000000000000002\;$ $\quad 6\quad$
$107.1098478543292\quad$ $206.6665840981687\quad$
$~{}~{}~{}0.000000000000022\;$ $\quad 8\quad$ $107.1098478543292\quad$
$206.6665840981686\quad$ $~{}~{}~{}0.000000000000013\;$ $\quad 10\quad$
$107.1098478543292\quad$ $206.6665840981686\quad$
$~{}~{}~{}0.000000000000017\;$ Case (II), $h=0.5$, $\tau=0.005$ $\quad
t_{n}\quad$ $E^{n}$ $I_{1}^{n}$ $I_{2}^{n}$ $\quad 0\quad$
$134.6831098503824\quad$ $181.8309885794262\quad$ $~{}~{}~{}0$ $\quad 2\quad$
$134.6831098503929\quad$ $181.8309885794291\quad$
$~{}~{}~{}0.000000000000020\;$ $\quad 4\quad$ $134.6831098503961\quad$
$181.8309885794310\quad$ $~{}~{}~{}0.000000000000039\;$ $\quad 6\quad$
$134.6831098503970\quad$ $181.8309885794320\quad$
$~{}~{}~{}0.000000000000074\;$ $\quad 8\quad$ $134.6831098503970\quad$
$181.8309885794321\quad$ $~{}~{}~{}0.000000000000009\;$ $\quad 10\quad$
$134.6831098503968\quad$ $181.8309885794320\quad$
$~{}~{}~{}0.000000000000006\;$ Case (III), $h=0.16$, $\tau=0.0025$ $\quad
t_{n}\quad$ $E^{n}$ $I_{1}^{n}$ $I_{2}^{n}$ $\quad 0\quad$
$52.66545441047811\quad$ $102.4999994287402\quad$
$~{}~{}~{}0.007689268607251\;$ $\quad 1\quad$ $52.66545441059597\quad$
$102.4999994287401\quad$ $~{}~{}~{}0.007689268554501\;$ $\quad 2\quad$
$52.66545441119767\quad$ $102.4999994287401\quad$
$~{}~{}~{}0.007689262987543\;$ $\quad 3\quad$ $52.66545441165016\quad$
$102.4999994287401\quad$ $~{}~{}~{}0.007689144120514\;$ $\quad 4\quad$
$52.66545441384896\quad$ $102.4999994287401\quad$
$~{}~{}~{}0.007687960991156\;$ $\quad 5\quad$ $52.66545444157879\quad$
$102.4999994287402\quad$ $~{}~{}~{}0.007680802034650\;$ Case (IV), $h=0.45$,
$\tau=0.0025$ $\quad t_{n}\quad$ $E^{n}$ $I_{1}^{n}$ $I_{2}^{n}$ $\quad
0\quad$ $152.6185732846562\quad$ $206.0751204356079\quad$
$~{}~{}~{}0.022285565392107\;$ $\quad 1\quad$ $152.6185714348163\quad$
$206.0751204356105\quad$ $~{}~{}~{}0.022285565391962\;$ $\quad 2\quad$
$152.6185618398274\quad$ $206.0751204356114\quad$
$~{}~{}~{}0.022285565214827\;$ $\quad 3\quad$ $152.6185527930804\quad$
$206.0751204356117\quad$ $~{}~{}~{}0.022285564293643\;$ $\quad 4\quad$
$152.6185466151580\quad$ $206.0751204356121\quad$
$~{}~{}~{}0.022285563402993\;$ $\quad 5\quad$ $152.6185421518980\quad$
$206.0751204356124\quad$ $~{}~{}~{}0.022285563022157\;$

(a) $M=200$ (b) $M=400$ (c) $M=800$ (d) $M=200$

(e) $M=400$ (f) $M=800$

Figure 1: The profiles of the predicted dam-break solutions of velocities
$u(x,t)$ and heights $\rho(x,t)$ using different spatial grid points at $t=2$
in Case (II).

(a) $M=200$ (b) $M=400$ (c) $M=800$ (d) $M=200$

(e) $M=400$ (f) $M=800$

Figure 2: The profiles of the predicted dam-break solutions of velocities
$u(x,t)$ and heights $\rho(x,t)$ using different spatial grid points at $t=2$
in Case (IV).

(a) Velocity variable, view(45,70) (b) Height variable, view(45,70)

Figure 3: The predicted dam-break solutions for the R2CH system in Case (II)
show evolution of the velocity $u(x,t)$ and the height $\rho(x,t)$ with
$t=50$.

(a) Velocity variable, view(45,70) (b) Height variable, view(45,70)

Figure 4: The predicted dam-break solutions for the R2CH system in Case (IV)
show evolution of the velocity $u(x,t)$ and the height $\rho(x,t)$ with
$t=50$.

###  5.2 Part II: nonsmooth initial data

In this part, we conduct numerical simulation for the R2CH system with
nonsmooth initial data. The parameters of the numerical examples come from the
following three cases

  * •

Case (I). $A=0$, $\mu=0$, $\sigma=1$, $\Omega=0$;

  * •

Case (II). $A=0$, $\mu=0$, $\sigma=1$, $\Omega=0.1$;

  * •

Case (III). $A=1$, $\mu=1$, $\sigma=1$, $\Omega=73\times 10^{-6}$.

######  Example 5.2 (Three-peakon interaction for the CH equation[9], [16],
[18], [22]).

We now test a degenerate R2CH system (4)–(5) by taking $\rho=A=\mu=0$ and
$\sigma=1$, which will reduce to the classical CH equation. Consider the
three-peakon interaction with the following initial condition

$\displaystyle u(x,0)=\phi_{1}(x)+\phi_{2}(x)+\phi_{3}(x),$

where

$\phi_{i}(x)=\left\\{\begin{array}[]{ll}\displaystyle\frac{c_{i}}{\cosh(L/2)}\cosh(x-x_{i}),\quad|x-x_{i}|\leqslant
L/2,\\\\[12.0pt]
\displaystyle\frac{c_{i}}{\cosh(L/2)}\cosh(L-(x-x_{i})),\quad|x-x_{i}|>L/2,\\\
\end{array}\right.\quad i=1,\,2,\,3.$

The parameters are given by $c_{1}=2$, $c_{2}=1$, $c_{3}=0.8$, $x_{1}=-5$,
$x_{2}=-3$, $x_{3}=-1$ and the computational domain is $[0,L]$ with $L=30$.

In the calculation, we take the spatial stepsize $h=L/2048$ and temporal
stepsize $\tau=1/10000$ to simulate this interaction at
$t=0,\,1,\,2,\,3,\,4,\,6,\,8,\,10$, respectively. Figure 5 shows the moving
peak interaction at different instants of time. It is clear to see that the
proposed scheme performs well in resolving the complex interaction among
multiple peakons for the CH equation. The results obtained here are comparable
to those obtained by a fourth-order HIEQ-GM method in [9] and a MSWCM-AD30
method in [22].

Figure 5: Three-peakon interaction of the CH equation at
$t=0,\,1,\,2,\,3,\,4,\,6,\,8,\,10$, respectively.

######  Example 5.3 (Single peakon[4]).

We now consider a single peakon solution of the R2CH system (2)–(3) with the
initial conditions

$\displaystyle u(x,0)=e^{-|x-x_{0}|},\quad\rho(x,0)=0.5,$

where $x_{0}=10$ is used and the computational domain is prescribed as
$[0,20]$.

The parameters in Case (I) and Case (III) are taken respectively to show the
behavior of the single peakon. We take spatial stepsize $h=0.025$ and temporal
stepsize $\tau=0.0005$ in the calculation. Figures 6 and 7 depict the
predicted solution profiles of the velocity variable $u(x,t)$ and height
variable $\rho(x,t)$ at $t=1$, $3$ and $5$, respectively, under Case (I) and
Case (III). In contrast, we observe that the selection of parameters has
remarkable impact on the profiles of solutions to the system (2)–(3). Figure 8
shows the computed results of discrete conserved quantities (5)–(7) at time
$t=5$ for Case (I) and Case (III). One can easily observe that the scheme
(1)–(4) guarantees the conservation of mass, momentum and energy for different
parameters. Moreover, One can also see that the rotational parameter $\Omega$
indeed changes the conserved quantities $E^{n}$ and $I_{2}^{n}$, which is
consistent with (6)–(7).

(a) Velocities (b) Heights

Figure 6: Velocities $u(x,t)$ and heights $\rho(x,t)$ for the R2CH system in
Case (I) computed by scheme (1)–(4) at three different times with stepsizes
$h=0.025$ and $\tau=0.0005$.

(a) Velocities (b) Heights

Figure 7: Velocities $u(x,t)$ and heights $\rho(x,t)$ for the R2CH system in
Case (III) computed by scheme (1)–(4) at three different times with stepsizes
$h=0.025$ and $\tau=0.0005$.

(a) Case (I) (b) Case (III)

Figure 8: Conserved quantities defined in (5)–(7) for the R2CH system in Case
(I) and Case (III) with stepsizes $h=0.025$ and $\tau=0.0005$.

######  Example 5.4 (Peakon anti-peakon interaction–I[15]).

We consider a piecewise smooth initial conditions

$\displaystyle
u(x,0)=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{2\sinh(\frac{1}{4})}\sinh(x),\quad
0\leqslant x\leqslant\frac{1}{4},\\\\[9.0pt]
\displaystyle\frac{1}{\sinh(-\frac{1}{2})}\sinh(x-\frac{1}{2}),\quad\frac{1}{4}<x\leqslant\frac{3}{4},\qquad\rho(x,0)=1.5.\\\\[9.0pt]
\displaystyle\frac{1}{2\sinh(\frac{1}{4})}\sinh(x-1),\quad\frac{3}{4}<x<1,\\\
\end{array}\right.$

To test the influence of $\Omega$ on the solution of the R2CH system (2)–(3),
the parameters in Cases (I) and (II) are selected for calculation in this
example. Figure 9 and Figure 10 respectively display the profiles of solutions
of the velocity variable $u(x,t)$ and the height variable $\rho(x,t)$ at
different instants of time, which enables us to observe how the solutions
evolve over time. In addition, it can be observed from Figure 9(a) and Figure
10(a) that $u(x,t)$ at initial time is jump discontinuous at $x=1/4$ and
$x=3/4$. Compared with Case (I), the symmetry is broken in Case (II). Figure
11 shows the results of discrete conserved quantities defined in (5)–(7) for
Cases (I) and (II). It is observed that the total momentum is zero when the
solution is symmetric, while it is nonzero when the solution is asymmetric.
Moreover, we portray the evolution of the peakon-antipeakon for the velocity
$u(x,t)$ and the height $\rho(x,t)$ in Figures 12–13. It can be seen from
Figure 12 that the evolution of solutions presents a certain periodicity in
the long time simulation, while the evolution of solutions is drawn in a short
time in Figure 13 since it will blow up in a slightly long time.

(a) Velocities (b) Heights

Figure 9:  Velocities $u(x,t)$ and heights $\rho(x,t)$ for the R2CH system in
Case (I) computed by scheme (1)–(4) at five different times with stepsizes
$h=0.002$ and $\tau=0.001$.

(a) Velocities (b) Heights

Figure 10:  Velocities $u(x,t)$ and heights $\rho(x,t)$ for the R2CH system in
Case (II) computed by scheme (1)–(4) at five different times with stepsizes
$h=0.002$ and $\tau=0.0005$.

(a) Case (I) (b) Case (II)

Figure 11: Conserved quantities defined in (5)–(7) for the R2CH system in Case
(I) and Case (II) with stepsizes $h=0.02$ and $\tau=0.0005$.

(a) Velocity variable, view(-20,60) (b) Height variable, view(-20,60)

Figure 12: The predicted peakon solutions for the R2CH system in Case (I) show
the evolution of the velocity $u(x,t)$ and the height $\rho(x,t)$ with
stepsizes $h=0.002$ and $\tau=0.001$.

(a) Velocity variable, view(-20,60) (b) Height variable, view(-20,60)

Figure 13: The predicted peakon solutions for the R2CH system in Case (II)
show the evolution of the velocity $u(x,t)$ and the height $\rho(x,t)$ with
stepsizes $h=0.002$ and $\tau=0.001$.

######  Example 5.5 (Peakon anti-peakon interaction–II[4], [16]).

Here we study the case of a peakon anti-peakon interaction over a long time.
For this purpose, we consider the following initial data

$\displaystyle
u(x,0)=p_{1}e^{-|x-x_{1}|}+p_{2}e^{-|x-x_{2}|},\quad\rho(x,0)=0.5,$

where $x_{1}=-5$ and $x_{2}=5$ respectively represent the position of the peak
and the trough, $p_{1}=1$, $p_{2}=-1$ and the computation domain is set to be
$[-20,20]$.

We depict the evolution of the solutions for the R2CH system (2)–(3), where
the parameters are selected as Cases (I) and (III) by taking $h=0.05$ and
$\tau=0.0005$ for calculation. Figures 14–17 portray the evolution behaviors
of velocity $u(x,t)$ and height $\rho(x,t)$ under these two cases. It easily
observes that for sufficiently long periods of time, the peak solutions will
have an elastic collision, so that we obtain a dissipative solution for the
R2CH system (2)–(3). By comparing Figures 14–15 with Figures 16–17, we find
that Case (I) still maintains symmetric, while the symmetry in Case (III) is
obviously broken. Similar to Example 5.4, we assert that the total momentum
remains zero in Case (I) because the two peak solutions have the same
magnitude but opposite signs. The momentum in Case (III) is nonzero because
the solution in Case (III) loses its symmetry. We verify this through the
calculation, which is not listed here for brevity. It is worth mentioning
that, since $\rho(x,0)>0$, it can be proved that $\rho(x,t)$ remains strictly
positive and that the solution retains the regularity of the initial data, see
e.g., [8]. Particularly, $\rho(x,t)$ remains bounded. We observe that the
numerical solution for the height $\rho(x,t)$ still keep positive. Ultimately,
the evolution of the velocity $u(x,t)$ and the height $\rho(x,t)$ for these
two cases in a longer time ($t=35$) is depicted in Figures 18–19. We clearly
see the symmetry in Case (I) and asymmetry in Case (III). These phenomena for
the asymmetrical case are displayed firstly in current paper.

(a) $t=1$ (b) $t=3$ (c) $t=5$ (d) $t=6$

(e) $t=8$ (f) $t=10$

Figure 14: The velocities $u(x,t)$ for the R2CH system in Case (I) at
different times with stepsizes $h=0.05$ and $\tau=0.0005$.

(a) $t=1$ (b) $t=3$ (c) $t=5$ (d) $t=6$

(e) $t=8$ (f) $t=10$

Figure 15: The heights $\rho(x,t)$ for the R2CH system in Case (I) at
different times with stepsizes $h=0.05$ and $\tau=0.0005$.

(a) $t=1$ (b) $t=3$ (c) $t=5$ (d) $t=6$

(e) $t=8$ (f) $t=10$

Figure 16: The velocities $u(x,t)$ for the R2CH system in Case (III) at
different times with stepsizes $h=0.05$ and $\tau=0.0005$.

(a) $t=1$ (b) $t=3$ (c) $t=5$ (d) $t=6$

(e) $t=8$ (f) $t=10$

Figure 17: The heights $\rho(x,t)$ for the R2CH system in Case (III) at
different times with stepsizes $h=0.05$ and $\tau=0.0005$.

(a) Velocity variable, view(0,90) (b) Height variable, view(0,90)

Figure 18: The predicted peakon solutions for the R2CH system in Case (I) show
the evolution of the velocity $u(x,t)$ and the height $\rho(x,t)$ with $t=35$.

(a) Velocity variable, view(0,90) (b) Height variable, view(0,90)

Figure 19: The predicted peakon solutions for the R2CH system in Case (III)
show the evolution of the velocity $u(x,t)$ and the height $\rho(x,t)$ with
$t=35$.

###  5.3 Part III: the conserved quantity $H$

Recalling the conserved quantity $H$ mentioned in the introduction, it can be
defined discretely as

$\displaystyle
H^{n}=\sum_{i=1}^{M}\bigg{[}(u_{i}^{n})^{3}+u_{i}^{n}\Big{(}\frac{u_{i+1}^{n}-u_{i-1}^{n}}{2h}\Big{)}^{2}-A(u_{i}^{n})^{2}-\mu\Big{(}\frac{u_{i+1}^{n}-u_{i-1}^{n}}{2h}\Big{)}^{2}+u_{i}^{n}(\rho_{i}^{n})^{2}\bigg{]},\quad
n=1,2,\cdots,N.$ (1)

We expect that the proposed scheme (1)–(4) also can preserve the conserved
quantity. Toward this end we conduct a conservation test for $H$ through two
different types of numerical examples including a smooth initial data problem
and two nonsmooth initial data problem. For the calculation of Example 5.1,
the parameters are selected from Case (I) and Case (II). For the nonsmooth
initial data problem, Example 5.4 and Example 5.5 are selected. All the
parameters in Case (I) are used. Due to the symmetry of initial values, it
immediately obtains that $H^{0}=0$ according to (1). Figure 20 depicts the
errors between $H^{n}$ and $H^{0}$ at different times, and we can find that
$H^{n}$ can approximate to $H$ numerically.

(a) $H^{n}-H^{0}$ in Example 5.1 (b) $H^{n}-H^{0}$ in Example 5.4 (c)
$H^{n}-H^{0}$ in Example 5.5

Figure 20: (a) shows the computed errors of the conserved quantity $H$ with
$h=0.5$ and $\tau=0.005$ in Case (I) and Case (II) in subsection 5.1 for
Example 5.1; (b) and (c) show the computed errors of the conserved quantity
$H$ with $h=0.2$ and $\tau=0.0025$ in Case (I) for Example 5.4 and Example 5.5
in subsection 5.2, respectively.

##  6 Concluding remarks

In this paper, we have carried out extensively numerical study on the R2CH
system under both the smooth and nonsmooth initial data based on a well-
designed conservative finite difference discretization. The evolution of
several asymmetric and non-smooth solitary wave solutions are depicted for the
first time. The numerical conservation laws enable the calculation schemes to
accurately capture the evolution of the smooth/nonsmooth solutions of the
benchmark problems. In particular, the numerical threshold technique plays a
key role in grasping drastic change of peakon solutions. Last but surely not
the least, there is still plenty of rooms for further study. Especially, the
strict error estimate of the difference scheme is a challenging topic, which
is not covered in this paper.

## Appendix

Noticing when $\sigma=1$ and $\Omega=0$, the system (4)–(5) can be written as
a system of hyperbolic type

$\displaystyle
u_{t}+uu_{x}+\partial_{x}G*\Big{(}u^{2}+\frac{1}{2}u_{x}^{2}-Au+\mu
u_{xx}+\frac{1}{2}\rho^{2}\Big{)}=0,$ (A.1)
$\displaystyle\rho_{t}+u\rho_{x}=-u_{x}\rho,\quad x\in R,t>0,$ (A.2)

then we show the system above has the conservation law

$\displaystyle H=\int_{R}(u^{3}+uu_{x}^{2}-Au^{2}-\mu u_{x}^{2}+u\rho^{2}){\rm
d}x.$

######  Proof 6.1.

For simplicity, we set

$\displaystyle g(x)=G*\Big{(}u^{2}+\frac{1}{2}u_{x}^{2}-Au+\mu
u_{xx}\Big{)},\quad h(x)=G*\Big{(}\frac{1}{2}\rho^{2}\Big{)}.$ (A.3)

Then (A.1) takes the equivalent form

$\displaystyle u_{t}+uu_{x}+g_{x}+h_{x}=0.$ (A.4)

Using (A.4) and noticing

$\int_{R}u\rho\rho_{t}\mathrm{d}x=-\int_{R}(u\rho)(u\rho)_{x}\mathrm{d}x=0,$

we have

$\displaystyle\frac{{\rm d}}{{\rm d}t}\int_{R}(u^{3}+uu_{x}^{2}-Au^{2}-\mu
u_{x}^{2}+u\rho^{2})\mathrm{d}x$
$\displaystyle\quad=\int_{R}(3u^{2}u_{t}+u_{t}u_{x}^{2}+2uu_{x}u_{xt}-2Auu_{t}-2\mu
u_{x}u_{xt}+u_{t}\rho^{2})\mathrm{d}x$
$\displaystyle\quad=\int_{R}(3u^{2}+u_{x}^{2}-2Au)u_{t}\mathrm{d}x+\int_{R}(2uu_{x}-2\mu
u_{x})u_{xt}\mathrm{d}x+\int_{R}u_{t}\rho^{2}\mathrm{d}x$
$\displaystyle\quad=-\int_{R}(3u^{2}+u_{x}^{2}-2Au)(uu_{x}+g_{x}+h_{x})\mathrm{d}x$
$\displaystyle\qquad\,-\int_{R}(2uu_{x}-2\mu
u_{x})(u_{x}^{2}+uu_{xx}+g_{xx}+h_{xx})\mathrm{d}x$
$\displaystyle\qquad\,-\int_{R}(uu_{x}+g_{x}+h_{x})\rho^{2}\mathrm{d}x$
$\displaystyle\quad\triangleq K_{1}+K_{2}+K_{3}.$ (A.5)

Below, we calculate each term at the right hand of (A.5). First, it easily
obtains that

$\displaystyle K_{1}$
$\displaystyle=-\int_{R}(3u^{2}+u_{x}^{2}-2Au)(uu_{x}+g_{x}+h_{x})\mathrm{d}x$
$\displaystyle=-\int_{R}uu_{x}^{3}\mathrm{d}x-\int_{R}(3u^{2}+u_{x}^{2}-2Au)g_{x}\mathrm{d}x-\int_{R}(3u^{2}+u_{x}^{2}-2Au)h_{x}\mathrm{d}x.$

Using (A.3), we have

$\displaystyle K_{2}$ $\displaystyle=-\int_{R}(2uu_{x}-2\mu
u_{x})(u_{x}^{2}+uu_{xx}+g_{xx}+h_{xx})\mathrm{d}x$
$\displaystyle=-\int_{R}(2uu_{x}^{3}+2u^{2}u_{x}u_{xx})\mathrm{d}x+\int_{R}(2\mu
u_{x}^{3}+2\mu uu_{x}u_{xx})\mathrm{d}x$
$\displaystyle\quad\,-\int_{R}(2uu_{x}-2\mu
u_{x})\Big{(}g-u^{2}-\frac{1}{2}u_{x}^{2}+Au-\mu u_{xx}\Big{)}\mathrm{d}x$
$\displaystyle\quad\,-\int_{R}(2uu_{x}-2\mu
u_{x})\Big{(}h-\frac{1}{2}\rho^{2}\Big{)}\mathrm{d}x$
$\displaystyle=-\int_{R}2uu_{x}g\mathrm{d}x+\int_{R}uu_{x}^{3}\mathrm{d}x+\int_{R}2\mu
u_{x}g\mathrm{d}x-\int_{R}\mu u_{x}^{3}\mathrm{d}x$
$\displaystyle\quad\;-\int_{R}2uu_{x}h\mathrm{d}x+\int_{R}uu_{x}\rho^{2}\mathrm{d}x+\int_{R}2\mu
u_{x}h\mathrm{d}x-\int_{R}\mu u_{x}\rho^{2}\mathrm{d}x.$

Similarly, we have

$\displaystyle
K_{3}=-\int_{R}\rho^{2}(uu_{x}+g_{x}+h_{x})\mathrm{d}x=-\int_{R}uu_{x}\rho^{2}\mathrm{d}x-\int_{R}2(h-h_{xx})g_{x}\mathrm{d}x.$

Adding $K_{1}$, $K_{2}$ and $K_{3}$ together, we obtain

$\displaystyle\frac{\mathrm{d}H}{\mathrm{d}t}$
$\displaystyle=-\int_{R}(2u^{2}+u_{x}^{2}-2Au+2\mu
u_{xx})g_{x}\mathrm{d}x+\int_{R}2\mu u_{xx}g_{x}\mathrm{d}x+\int_{R}2\mu
u_{x}g\mathrm{d}x$ $\displaystyle\quad\,-\int_{R}(2u^{2}+u_{x}^{2}-2Au+2\mu
u_{xx})h_{x}\mathrm{d}x+\int_{R}2\mu u_{xx}h_{x}\mathrm{d}x+\int_{R}2\mu
u_{x}h\mathrm{d}x$ $\displaystyle\quad\,-\int_{R}\mu
u_{x}^{3}\mathrm{d}x-\int_{R}2\mu
u_{x}(h-h_{xx})\mathrm{d}x-\int_{R}2(h-h_{xx})g_{x}\mathrm{d}x$
$\displaystyle=-\int_{R}2(g-g_{xx})g_{x}\mathrm{d}x-\int_{R}2\mu
u_{x}g_{xx}\mathrm{d}x+\int_{R}2\mu u_{x}g\mathrm{d}x$
$\displaystyle\quad\,-\int_{R}2(g-g_{xx})h_{x}\mathrm{d}x-\int_{R}2\mu
u_{x}h_{xx}\mathrm{d}x+\int_{R}2\mu u_{x}h\mathrm{d}x$
$\displaystyle\quad\,-\int_{R}\mu u_{x}^{3}\mathrm{d}x-\int_{R}2\mu
u_{x}(h-h_{xx})\mathrm{d}x-\int_{R}2(h-h_{xx})g_{x}\mathrm{d}x$
$\displaystyle=-\int_{R}2\mu u_{x}g_{xx}\mathrm{d}x+\int_{R}2\mu
u_{x}g\mathrm{d}x-\int_{R}\mu u_{x}^{3}\mathrm{d}x$
$\displaystyle=-\int_{R}2\mu u_{x}\Big{(}g-u^{2}-\frac{1}{2}u_{x}^{2}+Au-\mu
u_{xx}\Big{)}\mathrm{d}x+\int_{R}2\mu u_{x}g\mathrm{d}x-\int_{R}\mu
u_{x}^{3}\mathrm{d}x$ $\displaystyle=0,$

in which

$\displaystyle-\int_{R}u^{2}u_{x}u_{xx}\mathrm{d}x=\int_{R}uu_{x}^{3}\mathrm{d}x\quad\text{and}\quad-\int_{R}2\mu
uu_{x}u_{xx}\mathrm{d}x=\int_{R}\mu u_{x}^{3}\mathrm{d}x$

have been used several times during the calculation.

## References

  * [1] R. Camassa and D. Holm, An integrable shallow water equation with peaked solitons, Phys. Rev. Lett., 71 (1993), 1661–1664.
  * [2] R. Chen, L. Fan, H. Gao and Y. Liu. Breaking waves and solitary waves to the rotation-two-component Camassa–Holm system, SIAM J. Math. Anal., 49 (2017), 3573–3602.
  * [3] R. Chen, Y. Liu and Z. Qiao, Stability of solitary waves and global existence of a generalized two-component Camassa–Holm system, Commun. Part. Diff. Eq., 36 (2011), 2162–2188.
  * [4] A. Chertock, A. Kurganov and Y. Liu, Finite-volume-particle methods for the two-component Camassa–Holm system, Commun. Comput. Phys., 27(2) (2020), 480–502.
  * [5] A. Constantin, The Hamiltonian structure of the Camassa–Holm equation, Expo. Math., 15 (1997), 53–85.
  * [6] E. Fan and M. Yuen, Peakon weak solutions for the rotation-two-component Camassa–Holm system, Appl. Math. Lett., 97 (2019), 53–59.
  * [7] L. Fan, H. Gao and Y. Liu, On the rotation-two-component Camassa–Holm system modelling the equatorial water waves, Adv. Math., 291 (2016), 59–89.
  * [8] K. Grunert, H. Holden and X. Raynaud, Global solutions for the two-component Camassa–Holm system, Commun. Part. Diff. Eq., 37(12) (2012), 2245–2271.
  * [9] C. Jiang, Y. Wang and Y. Gong, Arbitrarily high-order energy-preserving schemes for the Camassa–Holm equation, Appl. Numer. Math., 151 (2020), 85–97.
  * [10] J. Lenells, Conservation laws of the Camassa–Holm equation, J. Phys. A: Math. Gen., 38 (2005), 869–880.
  * [11] S. Li and L. Vu-Quoc, Finite difference calculus invariant structure of a class of algorithms for the nonlinear Klein–Gordon equation, SIAM J. Numer. Anal., 32(6) (1995), 1839–1875.
  * [12] H. Liu and Y. Xing, An invariant preserving discontinuous Galerkin method for the Camassa–Holm equation, SIAM J. Sci. Comput., 38 (2016), A1919–A1934.
  * [13] B. Moon, On the wave-breaking phenomena and global existence for the periodic rotation-two-component Camassa–Holm system, J. Math. Anal. Appl., 451 (2017), 84–101.
  * [14] J. Liu, P. Pucci and Q. Zhang, Wave breaking analysis for the periodic rotation-two-component Camassa–Holm system, Nonlinear Anal., 187 (2019), 214–228.
  * [15] X. Li, X. Qian, B. Zhang and S. Song, A multi-symplectic compact method for the two-component Camassa–Holm equation with singular solutions, Chinese. Phys. Lett., 34(9) (2017), 090202.
  * [16] H. Liu and T. Pendleton, On invariant-preserving finite difference schemes for the Camassa–Holm equation and the two-component Camassa–Holm system, Commun. Comput. Phys., 19(4) (2016), 1015–1041.
  * [17] M. Yang, Y. Li and Z. Qiao, Persistence properties and wave-breaking criteria for a generalized two-component rotational b-family system, Discret. Contin. Dyn. Syst. A., 40 (2020), 2475–2493.
  * [18] Y. Xu and C.-W. Shu, A local Discontinuous Galerkin method for the Camassa-Holm equation, SIAM J. Numer. Anal., 46 (2008), 1998–2021.
  * [19] Q. Zhang, L. Liu and Z. Zhang, Linearly implicit invariant-preserving decoupled difference scheme for the rotation-two-component Camassa–Holm system, SIAM J. Sci. Comput., 44(4) (2022), A2226–A2252.
  * [20] Y. Zhang, Wave breaking and global existence for the periodic rotation-Camassa–Holm system, Discret. Contin. Dyn. Syst. A., 37 (2017), 2243–2257.
  * [21] K. Zhao and Z. Wen, 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

  * [22] H. Zhu, S. Song and Y. Tang, Multi-symplectic wavelet collocation method for the nonlinear Schröinger equation and the Camassa–Holm equation, Comput. Phys. Commun., 182(3) (2011), 616–627.
  * [23] M. Zhu and J. Xu, Wave-breaking phenomena and global solutions for periodic two-component Dullin–Gottwald–Holm systems, Electron J. Differ. Eq., 44 (2013), 1–27.