Arbitrarily high-order energy-preserving schemes for the Camassa-Holm equation

TL;DR

This paper introduces a new class of high-order energy-preserving numerical schemes for the Camassa-Holm equation, combining invariant energy quadratization, Fourier pseudo-spectral methods, and symplectic Runge-Kutta integrators to ensure energy conservation.

Contribution

The paper presents a novel framework that achieves arbitrarily high-order energy-preserving schemes for the Camassa-Holm equation using invariant energy quadratization and spectral methods.

Findings

Schemes exactly preserve semi-discrete energy conservation law.

Numerical results confirm high accuracy and efficiency.

The approach is applicable to high-order discretizations.

Abstract

In this paper, we develop a novel class of arbitrarily high-order energy-preserving schemes for the Camassa-Holm equation. With the aid of the invariant energy quadratization approach, the Camassa-Holm equation is first reformulated into an equivalent system, which inherits a quadratic energy. {The new system is then discretized} by the standard Fourier pseudo-spectral method, which can exactly preserve the semi-discrete energy conservation law. Subsequently, { a symplectic Runge-Kutta method such as the Gauss collocation method is applied} for the resulting semi-discrete system to arrive at an arbitrarily high-order fully discrete scheme. {We prove that the obtained schemes can conserve the discrete energy conservation law}. Numerical results are addressed to confirm accuracy and efficiency of the proposed schemes.

Tables

Table 1. Table. 1: The numerical error and convergence rate for different second-order schemes with h = 2 ​ π 128 ℎ 2 𝜋 128 h=\frac{2\pi}{128} and different time steps at T = 1 𝑇 1 T=1 .

Scheme	$\tau$	$L^{\infty }$	Rate
IEQ-LCNS	$\frac{1}{100}$	2.083e-04	-
	$\frac{1}{200}$	5.182e-05	2.01
	$\frac{1}{400}$	1.293e-05	2.00
	$\frac{1}{800}$	3.230e-06	2.00
EPFPS	$\frac{1}{100}$	4.355e-05	-
	$\frac{1}{200}$	1.089e-05	2.00
	$\frac{1}{400}$	2.722e-06	2.00
	$\frac{1}{800}$	6.806e-07	2.00
MSFPS	$\frac{1}{100}$	4.344e-05	-
	$\frac{1}{200}$	1.086e-05	2.00
	$\frac{1}{400}$	2.716e-06	2.00
	$\frac{1}{800}$	6.790e-07	2.00

Table 2. Table. 2: The numerical error and convergence rate for different high-order schemes with h = 2 ​ π 128 ℎ 2 𝜋 128 h=\frac{2\pi}{128} and different time steps at T = 1 𝑇 1 T=1 .

Scheme	$\tau$	$L^{\infty }$	Rate
4th-order HIEQ-GM	$\frac{1}{30}$	2.817e-07	-
	$\frac{1}{60}$	1.765e-08	4.00
	$\frac{1}{120}$	1.104e-09	4.00
6th-order HIEQ-GM	$\frac{1}{30}$	2.231e-010	-
	$\frac{1}{60}$	3.523e-012	5.98
	$\frac{1}{120}$	5.534e-014	5.99

Equations

$$u_t-u_{xxt}+3u{u}_x-2u_x{u}_{xx}-u{u}_{xxx}=0,(x,t)\in \Omega \times (0,T],$$

$$u(x+L,t)=u(x,t),(x,t)\in \Omega \times [0,T],$$

$$u(x,0)=u_0(x),x\in \Omega ,$$

$$\mathcal{M}={\int }_{\Omega }udx,\mathcal{I}={\int }_{\Omega }(u^2+u_x^2)dx,\mathcal{H}=-\frac{1}{2}{\int }_{\Omega }(u^3+u{u}_x^2)dx,$$

$$\frac{\partial u}{\partial t}=\mathcal{D}\frac{\delta \mathcal{H}}{\delta u},$$

$$\mathcal{H}={\int }_{\Omega }H(u,u_x)dx,H(u,u_x)=-\frac{1}{2}(u^3+u{u}_x^2),$$

$$\frac{\delta \mathcal{H}}{\delta u}=\frac{\partial H}{\partial u}-\frac{\partial }{\partial x}\frac{\partial H}{\partial u_x}=-\frac{3}{2}{u}^2-\frac{1}{2}{u}_x^2+(u{u}_x{)}_x.$$

$$\frac{d}{dt}\mathcal{H}=(\frac{\delta \mathcal{H}}{\delta u},\frac{\partial u}{\partial t})=(\frac{\delta \mathcal{H}}{\delta u},\mathcal{D}\frac{\delta \mathcal{H}}{\delta u})=0,$$

$$q=g(u,u_x)=-\frac{1}{2}(u^2+u_x^2),$$

$$\mathcal{H}={\int }_{\Omega }uqdx.$$

$$\displaystyle\left\{\begin{aligned} &\partial_{t}u=\mathcal{D}\Bigg{(}q+u\frac{\partial g}{\partial u}-\partial_{x}\big{(}u\frac{\partial g}{\partial u_{x}}\big{)}\Bigg{)},\\ &\partial_{t}q=\frac{\partial g}{\partial u}u_{t}+\frac{\partial g}{\partial u_{x}}\partial_{x}u_{t},\end{aligned}\right.\ \$$

$$\displaystyle u(x,0)=u_{0}(x),\ q(x,0)=-\frac{1}{2}\big{(}u(x,0)^{2}+u_{x}(x,0)^{2}\big{)},$$

$$\frac{\partial g}{\partial u}=-u,\frac{\partial g}{\partial u_x}=-u_x.$$

$$\frac{d}{dt}\mathcal{H}=0,\mathcal{H}={\int }_{\Omega }uqdx.$$

$$\frac{d}{dt}\mathcal{H}$$

$$\displaystyle=(u_{t},q)+\Big{(}u,\frac{\partial g}{\partial u}u_{t}+\frac{\partial g}{\partial u_{x}}\partial_{x}u_{t}\Big{)}$$

$$\displaystyle=(q,u_{t})+(u\frac{\partial g}{\partial u},u_{t})-\Big{(}\partial_{x}(u\frac{\partial g}{\partial u_{x}}),u_{t}\Big{)}$$

$$\displaystyle=\Big{(}q+u\frac{\partial g}{\partial u}-\partial_{x}\big{(}u\frac{\partial g}{\partial u_{x}}\big{)},u_{t}\Big{)}$$

$$\displaystyle=\Bigg{(}q+u\frac{\partial g}{\partial u}-\partial_{x}\big{(}u\frac{\partial g}{\partial u_{x}}\big{)},\mathcal{D}\Big{(}q+u\frac{\partial g}{\partial u}-\partial_{x}\big{(}u\frac{\partial g}{\partial u_{x}}\big{)}\Big{)}\Bigg{)}$$

$$=0,$$

$$U_j=U_{j+N}.$$

$$⟨\mathbf{U},\mathbf{V}{⟩}_h=h\sum _{j=0}^{N-1}{U}_j{V}_j,\forall \mathbf{U},\mathbf{V}\in {\mathbb{V}}_h.$$

$$\Vert \mathbf{U}{\Vert }_{h,\infty }=\underset{0\le j\le N-1}{\max }|U_j|.$$

$$\mathbf{U}\cdot \mathbf{V}=$$

$$g_j(x)=\frac{1}{N}\sum _{l=-N/2}^{N/2}\frac{1}{a_l}{e}^{\text{i}l\mu (x-x_j)},$$

$$I_{N}U(x)=\sum _{j=0}^{N-1}{U}_j{g}_j(x),$$

$$\frac{{\partial }^s{I}_{N}U(x_j)}{\partial x^s}$$

$$({\mathbf{D}}_s{)}_{j,k}=\frac{d^s{g}_k(x_j)}{d{x}^s},j,k=0,1,\cdots ,N-1.$$

$$\frac{\partial I_{N}U(x_j)}{\partial x}=[{\mathbf{D}}_1\mathbf{U}{]}_j,\frac{{\partial }^2{I}_{N}U(x_j)}{\partial x^2}=[{\mathbf{D}}_2\mathbf{U}{]}_j,j=0,1,\cdots ,N-1,$$

$$\displaystyle{\bm{D}}_{1}={\bm{F}}_{N}^{H}\Lambda_{1}{\bm{F}}_{N},\ \Lambda_{1}=\text{i}\mu\text{diag}\big{(}0,1,\cdots,\frac{N}{2}-1,0,-\frac{N}{2}+1,\cdots,-1\big{)},$$

$$\displaystyle{\bm{D}}_{2}={\bm{F}}_{N}^{H}\Lambda_{2}{\bm{F}}_{N},\ \Lambda_{2}=\big{[}\text{i}\mu\text{diag}(0,1,\cdots,\frac{N}{2}-1,\frac{N}{2},-\frac{N}{2}+1,\cdots,-1)\big{]}^{2},$$

$$\displaystyle\left\{\begin{aligned} &\frac{d}{dt}{\bm{U}}={\bm{D}}\Bigg{(}{\bm{Q}}-{\bm{U}}^{2}+{\bm{D}}_{1}\Big{(}({\bm{D}}_{1}{\bm{U}}\big{)}\cdot{\bm{U}}\Big{)}\Bigg{)},\\ &\frac{d}{dt}{\bm{Q}}=-{\bm{U}}\cdot\frac{d}{dt}{\bm{U}}-\big{(}{\bm{D}}_{1}{\bm{U}}\big{)}\cdot\Big{(}{\bm{D}}_{1}\frac{d}{dt}{\bm{U}}\Big{)},\end{aligned}\right.\$$

$$\displaystyle\frac{d}{dt}E=0,\ E=\Big{\langle}{{\bm{U}},\bm{Q}}\Big{\rangle}_{h}.$$

Full text

AMS subject classification: 65M20, 65M70, 65P10

Keywords: Invariant energy quadratization approach, Camassa-Holm equation, high-order energy-preserving scheme.

1 Introduction

In this paper, we consider the following Camassa-Holm (CH) equation [6, 7]

[TABLE]

with periodic boundary condition

[TABLE]

and initial condition

[TABLE]

where $\Omega=[a,a+L]$ is a bounded domain. The CH equation arises as a model for the propagation of unidirectional shallow water waves, with $u(x,t)$ representing the height of the fluid free surface above a flat bottom. Under the periodic boundary condition, the CH equation possesses the following conserved quantities

[TABLE]

where $\mathcal{M}$, $\mathcal{I}$ and $\mathcal{H}$ correspond to mass, momentum and Hamiltonian energy of the original problem, respectively.

The analysis and numerical solution of the CH equation has been widely investigated. Constantin and Escher [14] showed that the solution of the CH equation could develop singularities at a finite time, even for smooth initial data value with compact support. Li and Olver [29] established the local well-posedness of the CH equation in the nonhomogeneous Sobolev space $H^{s}$ with $s>3/2$. Numerical strategies for solving the CH equation include finite difference methods [10, 24], pseudo-spectral methods [28, 25], local discontinuous Galerkin method [43], operator splitting methods [5, 17], multi-symplectic methods [12, 13, 47] and other effective methods (e.g., see Refs. [8, 18]). However, the mentioned methods cannot exactly preserve the energy of the original system.

It is well known that the energy conservation is an important property of Hamiltonian partial differential equations (PDEs). In Ref. [33], Matsuo et al. presented an energy-conserving Galerkin scheme for the CH equation. Gong and Wang constructed an energy-preserving wavelet collocation scheme for the CH equation [19]. Other energy-preserving schemes can be found in Refs. [16, 32]. However, the existing energy-preserving schemes are only second order in time, which can’t provide satisfactory solutions in long time simulations with a given large time step. Therefore, in order to compute long time accurate solutions, the time step has to be refined, leading to expensive costs. To remove such obstacle, the most commonly used approach is to construct higher-order energy-preserving schemes, which make large marching steps practical while preserving the accuracy. To our best knowledge, there has been no reference considering a high-order energy-preserving scheme for the CH equation.

In [35], Quispel and McLaren derived third- and fourth-order averaged vector flied (AVF) methods. Subsequently, Li, Wang and Qin [30] have extended the fourth-order AVF method to sixth order. At present, high-order AVF methods have been successfully applied to develop high-order energy-preserving schemes for Hamiltonian PDEs (e.g., see [4, 27]). For linear problems, the proposed high-order schemes are concise and a customized fast solver has been presented to solve the resulting discrete linear equations efficiently [27]. However, since high-order AVF methods require to calculate high-order derivatives of the vector field, the resulting schemes are tedious for a nonlinear problem such as the CH equation. Based on the method of discrete line integral, Brugnano et al. [3] proposed Hamiltonian Boundary Value Methods (HBVMs), which can be recast as a multistage Runge-Kutta (RK) method. HBVMs are of arbitrarily high order and can exactly preserve energy of polynomial Hamiltonian systems. For non-polynomial cases, a practical energy-preserving scheme is usually gained in the sense that the energy error remains bounded within machine precision, but one has to increase RK stages [2]. Recently, energy-preserving continuous stage Runge-Kutta (CSRK) methods have been attracting a lot of interest [22, 31, 34, 41]. This kind of methods can eliminate the limit of HBVMs to cover non-polynomial Hamiltonian systems. However, the implementation of the CSRK method requires the computation of integrals. If the appearing integrals are replaced by a fixed high-order quadrature, the resulting method is closely related to HBVMs [22]. In addition, the above mentioned high-order energy-preserving methods are only valid for Hamiltonian systems with constant skew-symmetric structural matrix. For Hamiltonian systems with non-canonical structure matrix, these methods should be further discussed (e.g., see [1, 11, 42]).

In this paper, we present a novel strategy for efficiently developing arbitrarily high-order energy-preserving schemes for the CH equation. We first utilize the invariant energy quadratization (IEQ) approach [21, 44, 45, 46] to transform the CH equation into an equivalent system, which admits a quadratic energy. Then, the reformulated system is discretized in space by the standard Fourier pseudo-spectral method. Subsequently, arbitrarily high-order fully discrete schemes are derived by applying in time a symplectic RK method such as the Gauss collocation method. The newly proposed schemes have several desired advantages: (i) they are energy-preserving for the reformulated Camassa-Holm equation and can reach arbitrarily high-order; (ii) the required RK stages are optimal; (iii) they can be directly applied for efficiently solving non-canonical Hamiltonian PDEs where the energies are quadratic. It is worth noting that the proposed methods preserve the quadratic energy of the modified system, but not the Hamiltonian energy of the original CH equation.

The rest of the paper is organized as follows. In Section 2, an equivalent reformulation of the CH equation is presented based on the IEQ approach. In Section 3, a semi-discrete system, which inherits a quadratic energy, is obtained by using the standard Fourier pseudo-spectral method for the spatial discretization. In Section 4, we apply the Gauss collocation method for the semi-discrete system to arrive at a fully discrete scheme, which is proved to be energy-preserving. Several numerical experiments are presented in Section 5. We draw some conclusions in Section 6.

2 Model reformulation using the IEQ approach

In this section, we reformulate the CH equation into an equivalent form with a quadratic energy functional, which is called the IEQ reformulation. The IEQ reformulation for the CH equation provides an elegant platform for efficiently developing arbitrarily high-order energy-preserving schemes.

Firstly, Eq. (1.1) can be rewritten equivalently into the following Hamiltonian system

[TABLE]

where $\mathcal{D}=(1-\partial_{xx})^{-1}\partial_{x}$ is a skew-adjoint operator, $\mathcal{H}$ is the Hamiltonian energy

[TABLE]

and $\frac{\delta\mathcal{H}}{\delta u}$ denotes the variational derivative of $\mathcal{H}$ with respect to $u$

[TABLE]

One intrinsic property of (2.1) is energy-preserving, i.e.,

[TABLE]

where $(\cdot,\cdot)$ is the $L^{2}$-inner product defined by $(f,g)=\int_{\Omega}fgdx$.

Then, we formulate the idea of the IEQ approach for the Hamiltonian system (2.1). Let

[TABLE]

the Hamiltonian energy functional is then rewritten as

[TABLE]

According to energy variational, the Hamiltonian system (2.1) can be reformulated into the following equivalent form

[TABLE]

with consistent initial conditions

[TABLE]

where

[TABLE]

Theorem 2.1**.**

The system (2.5) satisfies the following quadratic energy

[TABLE]

Proof**.**

By some calculations, we obtain from the system (2.5)

[TABLE]

where the last equality follows from the skew-adjoint of $\mathcal{D}$. This completes the proof.∎

3 Energy-preserving spatial discretization

Many energy-preserving schemes have been designed and investigated for solving the CH equation in the literature, but little attention is paid to the energy-preserving properties brought by spatial discretization. In this section, the Fourier pseudo-spectral method is applied for the reformulated Camassa-Holm equation (2.5) to derive a spatial semi-discrete scheme, which is shown to preserve the semi-discrete quadratic energy (2.4).

Let $\Omega_{h}=\{x_{j}|x_{j}=a+jh,\ 0\leq j\leq N\}$ be a partition of $\Omega=[a,a+L]$ with mesh size $h=L/N$, where $N$ is an even number. A discrete mesh function ${U}_{j}=U(x_{j}),\ j\in\mathbb{Z}$ satisfies the periodic boundary condition if and only if

[TABLE]

Let \mathbb{V}_{h}=\big{\{}{\bm{U}}|{\bm{U}}=(U_{0},U_{1},\cdots,U_{N-1})^{T}\big{\}} be the space of mesh functions on $\Omega_{h}$ that satisfy the periodic boundary condition (3.1). We define the discrete inner product as follows

[TABLE]

The discrete $L^{\infty}$-norm of ${\bm{U}}\in\mathbb{V}_{h}$ is defined as

[TABLE]

In addition, we denote $`\cdot$’ as the componentwise product of vectors ${\bm{U}},{\bm{V}}\in\mathbb{V}_{h}$, that is,

[TABLE]

For brevity, we denote $\underbrace{{\bm{U}}\cdot...\cdot{\bm{U}}}_{p}$ as ${\bm{U}}^{p}$.

Let $S_{N}=\text{span}\{g_{j}(x),\ 0\leq j\leq N-1\}$ be the interpolation space, where $g_{j}(x)$ is the trigonometric polynomial of degree $N/2$ given by

[TABLE]

with $a_{l}=\left\{\begin{aligned} &1,\ |l|<\frac{N}{2},\\ &2,\ |l|=\frac{N}{2},\end{aligned}\right.$ and $\mu=\frac{2\pi}{L}$. We define the interpolation operator $I_{N}:C(\Omega)\to S_{N}$ as follows [9]

[TABLE]

where $U_{j}=U(x_{j})$. Taking the derivative with respect to $x$, and then evaluating the resulting expression at the collocation point $x_{j}$, we have

[TABLE]

where $j=0,\cdots,N-1$ and ${\bm{D}}_{s}$ is an $N\times N$ matrix with elements given by

[TABLE]

In particular, for first and second derivatives, we have, respectively

[TABLE]

where ${\bm{D}}_{1}$ is a real skew-symmetric matrix, and ${\bm{D}}_{2}$ is a real symmetric matrix. We note that [38]

[TABLE]

where ${\bm{F}}_{N}$ is the discrete Fourier transform matrix with elements \big{(}{\bm{F}}_{N}\big{)}_{j,k}=\frac{1}{\sqrt{N}}e^{-\text{\rm i}jk\frac{2\pi}{N}}, ${\bm{F}}_{N}^{H}$ is the conjugate transpose matrix of ${\bm{F}}_{N}$.

Applying the standard Fourier pseudo-spectral method to the system (2.5) in space, we have

[TABLE]

where ${\bm{D}}=({\bm{I}}-{\bm{D}}_{2})^{-1}{\bm{D}}_{1}$, and ${\bm{U}},{\bm{Q}}\in\mathbb{V}_{h}.$

Theorem 3.1**.**

The system (3.4) preserves the following semi-discrete quadratic energy

[TABLE]

Proof**.**

The proof strictly follows that done for Theorem 2.1, thus, for brevity, we omit it.

Remark 3.1**.**

If the Fourier pseudo-spectral method is employed for the original system (2.1), we can obtain a new semi-discrete scheme

[TABLE]

which can also be proved to preserve a semi-discrete Hamiltonian energy

[TABLE]

We note that the quadratic energy (2.4) is only equivalent to the Hamiltonian energy (2.2) in continuous sense, but not for the semi-discrete sense.

4 Energy-preserving fully discretized schemes

In this section, we first derive a class of high-order energy-preserving schemes by using the Gauss collocation method in time for the IEQ reformulation (3.4). Then, we show that, together with the other time integrators [20, 26], the IEQ reformulation (3.4) also provides an elegant platform for efficiently developing linear-implicitly energy-preserving schemes.

4.1 High-order energy-preserving schemes

Applying an $s$-stage collocation method to the system (3.4) in time, we obtain the following scheme.

Scheme 4.1**.**

Let $c_{1},\cdots,c_{s}$ be distinct real numbers (usually $0\leq c_{i}\leq 1$). For given ${\bm{U}}^{n},{\bm{Q}}^{n}\in\mathbb{V}_{h}$, ${\bm{u}}(t)$ and ${\bm{v}}(t)$ are two $N$ dimensional vector polynomials of degree $s$ satisfying, respectively,

[TABLE]

where $t_{n}^{i}=t_{n}+c_{i}\tau,i=1,\cdots,s$. And the numerical solution is defined by ${\bm{U}}^{n+1}={\bm{u}}(t_{n}+\tau)$ and ${\bm{Q}}^{n+1}={\bm{v}}(t_{n}+\tau)$, respectively.

Theorem 1.4 on page 31 of Ref. [23] indicates that the collocation method is equivalent to a RK method. If we take $c_{1},\cdots,c_{s}$ as the zeros of the $s$th shifted Legendre polynomial

[TABLE]

Scheme 4.1 is called Gauss collocation method and has order $2s$, and the zeros are called Gauss collocation points. Collocation points for Gauss collocation methods of order 4 and 6 are given explicitly in Ref. [23].

Remark 4.1**.**

As pointed out above, the proposed scheme can not preserve the Hamiltonian energy (2.2) of the original Camassa-Holm equation, but only the quadratic energy of the modified system (2.5).

It is well known that symplectic Runge-Kutta schemes preserve all the quadratic invariants of the ODE problem (see Refs. [15, 36, 37]). Since Gauss collocation schemes of any order are symplectic (see Ref. [36] and references therein), the following theorem is straightforward.

Theorem 4.1**.**

The $s$-stage Gauss collocation Scheme 4.1 is energy-preserving, i.e., it satisfies the following quadratic energy

[TABLE]

Remark 4.2**.**

Since any other symplectic Runge-Kutta method would preserve the quadratic energy, other arbitrarily high-order schemes which preserve the discrete quadratic energy (4.4) can be easily obtained.

4.2 Linearly-implicit energy-preserving schemes

In this subsection, a novel, linearly-implicit and energy-preserving scheme for the CH equation is obtained by utilizing in time the linearized Crank-Nicolson method for the semi-discrete system (3.4). The resulting scheme is denoted by IEQ-LCNS.

Scheme 4.2**.**

Applying the linearized Crank-Nicolson method to discretize the semi-discrete system (3.4) in time, we obtain a fully discretized scheme, as follows:

[TABLE]

where ${\delta_{t}^{+}{\bm{U}}^{n}}=\frac{{\bm{U}}^{n+1}-{\bm{U}}^{n}}{\tau},\hat{\bm{U}}^{n+\frac{1}{2}}=\frac{3{\bm{U}}^{n}-{\bm{U}}^{n-1}}{2}$ and ${\bm{Q}}^{n+\frac{1}{2}}=\frac{{\bm{Q}}^{n+1}+{\bm{Q}}^{n}}{2}$ and ${\bm{U}}^{1},{\bm{Q}}^{1}$ is the solution of the following equation

[TABLE]

Theorem 4.2**.**

The IEQ-LCNS (4.5)-(4.6) can exactly preserve the discrete quadratic energy (4.4).

Proof**.**

We can deduce from (4.5) that

[TABLE]

where

[TABLE]

and the last equality follows from the skew-symmetry of ${\bm{D}}$. Thus, according to (Proof), we have

[TABLE]

An argument similar to (4.6) used in (Proof) shows that

[TABLE]

This completes the proof.∎

Subsequently, we show that the above scheme can be solved efficiently. Indeed, we first rewrite (4.5) as

[TABLE]

where

[TABLE]

Then, by eliminating ${\bm{Q}}^{n+\frac{1}{2}}$ from (4.8), we have

[TABLE]

where

[TABLE]

Finally, we obtain ${\bm{U}}^{n+\frac{1}{2}}$ from (4.9) by using the following iteration method for linear equations (4.9)

[TABLE]

where we take the initial iteration vector ${\bm{U}}^{n+\frac{1}{2},0}={\bm{U}}^{n}$ and each iteration will terminate if the infinity norm of the error between two adjacent iterative steps is less than $10^{-14}$. Then, ${\bm{Q}}^{n+\frac{1}{2}}$ is obtained from the second equality of (4.8). Subsequently, we have ${\bm{U}}^{n+1}=2{\bm{U}}^{n+\frac{1}{2}}-{\bm{U}}^{n}$ and ${\bm{Q}}^{n+1}=2{\bm{Q}}^{n+\frac{1}{2}}-{\bm{Q}}^{n}$.

Remark 4.3**.**

We should note from (4.9) that, the IEQ approach need introduce an auxiliary variable, but the auxiliary variable can be eliminated in practical computations.

5 Numerical examples

In this section, we will investigate the accuracy, CPU time and invariants-preservation of the proposed schemes. Theoretically, the newly proposed high-order schemes 4.1 could reach arbitrarily high-order accuracy in time (with proper choice of the Gauss collocation points), and they all can exactly preserve the discrete quadratic energy (4.4). For simplicity, in the rest of this paper, the Gauss methods of order 4 (denoted by 4th-order HIEQ-GM) and 6 (denoted by 6th-order HIEQ-GM) are only used for demonstration purposes. Also, the results are compared with the energy-preserving Fourier pseudo-spectral scheme (denoted by EPFPS) and the multi-symplectic Fourier pseudo-spectral scheme (denoted by MSFPS), where we substitute the Fourier pseudo-spectral method into the wavelet collocation method for space directions in Refs. [19] and [47], respectively. For the convergence rate, we use the formula

[TABLE]

where $\tau_{l},error_{l},(l=1,2)$ are step sizes and errors with the step size $\tau_{l}$, respectively.

5.1 Accuracy test

We consider the periodic smooth solution with initial condition

[TABLE]

and the periodic boundary condition. The exact solution is obtained numerically by 6th-order HIEQ-GM under a very small time step $\tau=0.001$ and spatial step $h=\frac{2\pi}{128}$ at $T=1$. In Tables 1 and 2, we display the numerical error in discrete $L^{\infty}$-norm and the convergence rate for different schemes at $T=1$, respectively. As illustrated in Table 1, all of the schemes have second-order convergence rate in time, and from Table 2, it is clear to see that 4th-order HIEQ-GM and 6th-order HIEQ-GM can arrive at fourth-order and sixth-order convergence rates in time, respectively. In Fig. 1, we plot the global numerical error in discrete $L^{\infty}$-norm versus the CPU time for different schemes. The plot shows that, for a given global error, the sixth-order scheme is computationally cheapest. The IEQ-LCNS admits larger numerical errors than the ones provided by EPFPS and MSFPS, however, it is computationally cheaper.

5.2 Peakon solution

We consider the periodic peaked traveling wave with initial condition [43]

[TABLE]

where $c$ is the wave speed, $L$ is the period, and $x_{0}$ is the position of the trough. In the numerical experiment, the parameters are chosen as $c=1$, $L=1$, and $x_{0}=0$ and the periodic boundary condition is considered. The errors of invariants are plotted in Fig. 2. In Fig. 2 (a), we can see that IEQ-LCNS and MSFPS can only preserve the Hamiltonian energy approximately and the error provided by IEQ-LCNS is largest. In theory, the proposed high-order schemes cannot exactly preserve the discrete Hamiltonian energy, however, from Fig. 2 (a), we can observe that the resulting errors provided by 4th-order HIEQ-GM and 6th-order HIEQ-GM, respectively, can be preserved up to the machine accuracy and are much smaller than the one provided by EPFPS. In Figs. 2 (b)-(c) , it is clear to see that the errors of the momentum are bounded and all of the schemes can exactly preserve the mass conservation law. Fig. 2 (d) show that the proposed schemes can exactly preserve the discrete quadratic energy (4.4), which conforms the theoretical analysis.

5.3 Three-peakon interaction

In this example, we consider the three-peakon interaction of the CH equation with initial condition [43]

[TABLE]

where

[TABLE]

The parameters are $c_{1}=2,c_{2}=1,c_{3}=0.8,x_{1}=-5,x_{2}=-3,x_{3}=-1$ and $L=30$, and the computational domain is $\Omega=[0,L]$ with the periodic boundary condition. In Fig. 3, we display the interaction of three peakons by 4th-order HIEQ-GM at $t=0,1,2,3,4,6,8$ and $10$, respectively. We can see clearly that the moving peak interaction is resolved very well. The three-peakon interaction obtained by other schemes are not presented since they are close to Fig. 3. The errors of invariants are plotted in Fig. 4, which shows that all of the schemes can exactly preserve the mass conservation law and the momentum errors provided by the schemes are bounded. The Hamiltonian energy errors provided by the high-order schemes are smallest and the newly proposed schemes can exactly preserve the discrete quadratic energy (4.4).

6 Concluding remarks

In this paper, we combine the idea of the IEQ approach with symplectic Runge-Kutta methods to propose a new class of energy-preserving methods for the CH equation. The proposed schemes could reach arbitrarily high-order accuracy while exactly preserving the discrete quadratic energy of the modified system. Numerical examples are addressed to illustrate the accuracy and energy-preserving property of the proposed schemes. Compared with some existing low-order structure-preserving schemes, the proposed high-order schemes show remarkable efficiency and the advantage in preserving the discrete Hamiltonian energy.

We conclude this paper with some remarks. First, the presented strategy can be directly extended to propose high-order energy-preserving methods for the Hamiltonian PDEs where the Hamiltonian functionals are quadratic. For a general case (including the non-canonical Hamiltonian PDE), we dealt with it, as follows: we first utilize the idea of the IEQ approach to transform Hamiltonian energy as a quadratic form and then, following the energy variational, the original system is reformulated into an equivalent system, which inherits such quadratic energy. Finally, the resulting system is solved by using in time a symplectic Runge-Kutta method. Second, compared with a existing energy-preserving method such as HBVMs, the proposed method can not preserve the Hamiltonian energy of the original system. Thus, such trade-offs among methods should be more carefully investigated. Finally, we can also apply the idea of the scalar auxiliary variable (SAV) approach [39, 40] to reformulate the CH equation into a new equivalent system which inherits a quadratic energy. Thus, a possible future work is to develop high-order energy-preserving methods for Hamiltonian PDEs by combine the idea of the SAV approach with the symplectic Runge-Kutta method.

Acknowledgments

The authors would like to express sincere gratitude to the referees for their insightful comments and suggestions. Chaolong Jiang’s work is partially supported by the National Natural Science Foundation of China (Grant No. 11901513), the Yunnan Provincial Department of Education Science Research Fund Project (Grant No. 2019J0956) and the Science and Technology Innovation Team on Applied Mathematics in Universities of Yunnan. Yushun Wang’s work is partially supported by the National Natural Science Foundation of China (Grant No. 11771213) and the National Key Research and Development Project of China (Grant Nos. 2016YFC0600310, 2018YFC0603500, 2018YFC1504205). Yuezheng Gong’s work is partially supported by the Natural Science Foundation of Jiangsu Province (Grant No. BK20180413) and the National Natural Science Foundation of China (Grant No. 11801269).