Arbitrarily high-order conservative schemes for the generalized Korteweg-de Vries equation

Kai Yang Florida International University

Abstract.

This paper proposes a new class of arbitrarily high-order conservative numerical schemes for the generalized Korteweg-de Vries (KdV) equation. This approach is based on the scalar auxiliary variable (SAV) method. The equation is reformulated into an equivalent system by introducing a scalar auxiliary variable, and the energy is reformulated into a sum of two quadratic terms. Therefore, the quadratic preserving Runge-Kutta method will preserve all the three invariants (momentum, mass and the reformulated energy) in the discrete time flow (assuming the spatial variable is continuous). With the Fourier pseudo-spectral spatial discretization, the scheme conserves the first and third invariant quantities (momentum and energy) exactly in the space-time full discrete sense. The discrete mass possesses the precision of the spectral accuracy. Our numerical experiment shows the great efficiency of this scheme in simulating the breathers for the mKdV equation.

Key words and phrases:

generalized KdV, SAV, energy conservation, high-order conservative scheme, breathers

1. Introduction

This paper considers the numerical approximation of the generalized Korteweg–de Vries (gKdV) equation

(1.1)

\displaystyle\begin{cases}u_{t}=-(u_{xx}+\frac{1}{p}u^{p})_{x},\quad x\in\mathbb{R},\,\,t>0\\ u(x,0)=u_{0},\end{cases}

where $p$ is a positive integer.

When $p=2$ , it is known as the KdV equation.When $p=3$ , it is the modified KdV (mKdV) equation. One appealing feature of these two cases is that they form the integrable systems and possess infinitely many integrals, which are invariant in time, see [40]. This is accounted by reformulating the KdV equation into the Lax pairs [33]. When $p\geq 4$ , it is known as the generalized KdV equation, which is not integrable. In general, the equation (1.1) conserves the following three invariant quantities:

(1.2)		$\displaystyle I[u(t)]\stackrel{{\scriptstyle\rm{def}}}{{=}}\int u(x,t)dx=I[u_{0}];$
(1.3)		$\displaystyle M[u(t)]\stackrel{{\scriptstyle\rm{def}}}{{=}}\int[u(x,t)]^{2}dx=M[u_{0}];$
(1.4)		$\displaystyle E[u(t)]\stackrel{{\scriptstyle\rm{def}}}{{=}}\int\left[\frac{1}{2}\left(u_{x}(x,t)\right)^{2}-\frac{1}{p(p+1)}\left(u(x,t)\right)^{p+1}\right]dx=E[u_{0}],$

known as the momentum, mass and energy (Hamiltonian), respectively.

These evolution equations are among the simplest of a general class of models featuring nonlinear advection (the term $(\frac{1}{p}u^{p})_{x}$ in this case). This family of equations arise as mathematical models for the propagation of physical waves in a wide variety of situations, such as shallow-water waves with weakly non-linear restoring forces, long internal waves in a density-stratified ocean, ion-acoustic waves in a plasma, acoustic waves on a crystal lattice, etc., e.g., see [3], [27], [6], [5], [14], [7]. The equations also attracted attention from the mathematical theory side. Analytical and numerical investigations include well-posedness, [30], [31]; soliton stability, [8], [39]; dispersion limit, [21]; and singularity formations, [32]. Nevertheless, there are still many open questions, and therefore, an efficient and accurate numerical algorithm would be desirable for future investigations.

Numerical studies of the gKdV equation trace back to 1970’s, ranging from finite difference methods [18], [34], [16]; finite element methods [2], [4], [28]; spectral methods [17], [23], [26], [38], [43]; operator splitting methods [25], [24] and discontinuous Galerkin methods [11], [35], [49], [52], [36]. Besides the order of accuracy and the stiffness from the term $u_{xxx}$ , preserving the conservation laws (1.2)-(1.4) is also a main concern in designing a numerical method for (1.1). Indeed, for conservative PDEs, numerical methods, which can preserve the corresponding invariants, are often advantageous: besides the accuracy, a conservative scheme can preserve good stability properties, especially in long-time simulations. On the other hand, to our best knowledge, most of the papers above are only capable of preserving one or, possibly, two quantities from (1.2)-(1.4), with the second order accuracy in time.

The purpose of this paper is to present a numerical scheme for the gKdV equation that preserves all these three invariant quantities (1.2)-(1.4) exactly in the discrete time, together with an arbitrarily high order of temporal accuracy. This is achieved by applying the scalar auxiliary variable (SAV) approach. The SAV approach, which was developed from the IEQ approach ([50] and [51]), was proposed for minimizing the free energy by gradient flows, aiming to keep the energy stable in the discrete time flow, see [46], [45] and convergence analysis in [47]. The main difference between these two methods is that while the IEQ method is introducing an additional variable function from the nonlinear terms in the equation, the SAV method is introducing a scalar function from the potential part in the energy. Hence, when considering the spatial discretized system (assume $u$ is discretized into $N$ grid points or nodes), the IEQ method results in a $2N\times 1$ system, and the SAV approach results in an $(N+1)\times 1$ system, which reduces the computational cost significantly. The SAV or IEQ methods are then applied to the Hamiltonian PDEs in [10] and [15]. Inspired by the same idea, we reformulate the equation (1.1) into an equivalent system by introducing an auxiliary variable. The reformulated system conserves the original momentum and mass, and also the modified energy, which is rewritten into the sum of two quadratic terms. Therefore, by the standard numerical ODE theory, the quadratic preserving (symplectic) Runge-Kutta methods will preserve all these three invariants exactly in the discrete time flow (assume the space variable is continuous). The standard Fourier pseudo-spectral method is chosen for the spatial discretization. This spatial discretization preserves the momentum and energy exactly in the spatial discrete sense (assuming time is continuous), and consequently, the error of the discrete mass only comes from the spatial discretization, which is of spectral accuracy and usually unnoticeable in our numerical experiments.

This paper is organized as follows. In Section 2, we give the equivalent form of the reformulated gKdV equation (1.1) and the modified energy (1.4) based on the SAV approach. In Section 3, we show that the family of the quadratic preserving (symplectic) Runge-Kutta methods will preserve the momentum, mass and modified energy in the discrete time flow. In Section 4, we describe the spatial discretization. We prove that the conservation laws (1.2) and (1.4) hold in the spatial discrete sense. We also show the error of (1.3) from this discretization. Combing with results in Section 3, we prove that the proposed scheme preserves the momentum and energy exactly in the space-time fully discretized scheme in the computation. In the meanwhile, we also give an error estimate of the discrete mass, which is $\sim C_{d}\cdot t$ for some constant $C_{d}$ coming from the error of the spatial discretization. Due to the high accuracy of the Fourier spectral method, the constant $C_{d}$ is generally on the order of $10^{-14}$ , several orders smaller than the tolerance for solving the resulting nonlinear system, and thus, the discrete mass error is almost unnoticeable. In Section 5, we first describe the numerical algorithm for solving the resulting nonlinear system from the symplectic Runge-Kutta method. Next, we briefly describe the several existing algorithms for the gKdV equations, such as the modified Crank-Nicholson method, conventional SAV with Leap-Frog (SAV-LF) approach, the Strang-Splitting (SS) method from [25], and the 4th order modified exponential time differencing (mETDRK4) method ([13] and [29]). They will be used in this paper for comparison purpose. Finally, we list our numerical examples. Numerical results show the fulfillment of conservation laws and a high accuracy of the proposed scheme, especially in simulating the breathers, which is a type of highly oscillatory non-dispersing pulse solutions for the mKdV equations.

Acknowledgment:

The author is partially supported by the NSF grant DMS-1927258 (PI: Svetlana Roudenko). The author is thankful for Dr. Roudenko’s helpful discussion, reading and remarks on the paper.

2. Model reformulation by the SAV approach

In this section, we reformulate the equation (1.1) into an equivalent system, which possesses a modified energy function of a new variable.

We first define the inner product for the real valued functions functions $f,g\in L^{2}(\mathbb{R}):$

(2.1)

\displaystyle(f,g)\stackrel{{\scriptstyle\rm{def}}}{{=}}\int_{\mathbb{R}}f(x)g(x)dx.

Assume $u^{p+1}\in L^{1}_{x}([0,T])$ . Consider a new scalar variable

v\stackrel{{\scriptstyle\rm{def}}}{{=}}v(u,t)=\sqrt{(u^{p},u)+C_{0}},

where $C_{0}$ is a large enough positive constant to prevent the expression $(u^{p},u)+C_{0}$ under the square root to become negative during our simulation time $t\in[0,T]$ . In Section 3, we will provide an adjustment process for $C_{0}$ when $v(t_{m})$ is close to $0$ at some time $t=t_{m}$ . Therefore, for now we only need to choose a constant $C_{0}$ at $t=t_{m}$ to make sure $v(t)>0$ for $t\in[t_{m},t_{m+1}]$ . This can be easily achieved, since we assume that the solution is well-posed (smooth in time) in $t\in[0,T]$ .

The equation (1.1) is then reformulated into the following system

(2.2)

\displaystyle\begin{cases}u_{t}=-\left(u_{xx}+\frac{1}{p}\frac{u^{p}v}{\sqrt{(u^{p},u)+C_{0}}}\right)_{x}\stackrel{{\scriptstyle\rm{def}}}{{=}}f(u,v),\\ v_{t}=\frac{p+1}{2\sqrt{(u^{p},u)+C_{0}}}(u^{p},u_{t})\stackrel{{\scriptstyle\rm{def}}}{{=}}g(u,v),\end{cases}

with the initial condition

\displaystyle u(x,0)=u_{0},\qquad v_{0}=\sqrt{(u_{0}^{p},u_{0})+C_{0}}.

Since the modified system (2.2) is identical to (1.1) when the spatial and time variables $(x,t)$ are continuous, the modified system (2.2) conserves the momentum, mass and modified energy in the form

(2.3)		$\displaystyle I[u(t)]\stackrel{{\scriptstyle\rm{def}}}{{=}}(u,1)\equiv I[u_{0}];$
(2.4)		$\displaystyle M[u(t)]\stackrel{{\scriptstyle\rm{def}}}{{=}}(u,u)\equiv M[u_{0}];$
(2.5)		$\displaystyle E[u(t),v(t)]\stackrel{{\scriptstyle\rm{def}}}{{=}}-\frac{1}{2}(u_{xx},u)-\frac{1}{p(p+1)}(v^{2}-C_{0})\equiv E[u_{0},v_{0}].$

3. Temporal discretization

In this section, we describe the temporal discretization of the scheme. We first show that the family of the symplectic Runge-Kutta (SRK) methods will conserve all three quantities (2.3) (2.4) and (2.5) in the discrete time flow. Then, we propose a strategy to adjust the constant $C_{0}$ in the computation, which makes the reformulated system (2.2) valid all the time and keeps the energy invariant.

3.1. Symplectic Runge-Kutta method

We first recall the $s$ -stage collocation Runge-Kutta method. Consider the IVP of the ODE in general form

{u}_{t}=f({u},t),\quad{u}(t_{0})={u}_{0}.

Define $\tau$ to be the time step. Let $c_{1},\cdots,c_{s}$ be distinct real numbers (usually $0\leq c_{i}\leq 1$ ). The collocation polynomial ${u}(t)$ is a polynomial of degree $s$ satisfying

{u}_{t}(t_{0}+c_{i}\tau)=f({u}(t_{0}+c_{i}\tau),t_{0}+c_{i}\tau),

and the numerical solution is defined by ${u}_{1}={u}(t_{0}+\tau)$ .

From [22], we know that the $s$ -stage collocation method is equivalent to the $s$ -stage Runge-Kutta method. Indeed, consider

k_{i}:={u}_{t}(t_{0}+c_{i}\tau)=f({u}(t_{0}+c_{i}\tau),t_{0}+c_{i}\tau).

Letting $l_{i}(z)=\Pi_{l\neq i}(z-c_{l})/(c_{i}-c_{l})$ be the Lagrange polynomial, the Lagrange interpolation formula gives us

{u}_{t}(t_{0}+z\tau)=\sum_{j=1}^{s}k_{j}l_{j}(z).

Integrating from $0$ to $c_{i}$ for each $i$ yields

u(t_{0}+c_{i}\tau)=u_{0}+\tau\sum_{j=1}^{s}k_{j}\int_{0}^{c_{i}}l_{j}(z)dz.

Denote $a_{ij}=\int_{0}^{c_{i}}l_{j}(z)dz$ , $b_{i}=\int_{0}^{1}l_{i}(z)dz$ and the intermediate values $u_{i}=u(t_{0}+c_{i}\tau)$ . We have the Runge-Kutta formulation

\displaystyle u_{i}=u_{0}+\sum_{j=1}^{s}k_{j}a_{ij},

and, the solution is updated by

\displaystyle u(t_{0}+\tau)=u_{0}+\sum_{i=1}^{s}k_{i}b_{i}.

We rewrite the coefficients $\mathbf{A}=(a_{ij})$ , $\mathbf{b}=(b_{1},b_{2},\cdots,b_{s})$ and $\mathbf{c}=(c_{1},c_{2},\cdots,c_{s})^{T}$ in the Butcher’s Tableaus ([9]):

\begin{array}[]{c|c}\mathbf{c}&\mathbf{A}\\ \hline\cr&\mathbf{b}\end{array}.

For example, we list three commonly used Runge-Kutta methods in the Butcher’s Tableaus in Table 1. $c_{i}$ ’s are chosen from the Gaussian-Legendre collocation points. Thus, they are known as the $s$ -stage Gaussian-Legendre Runge-Kutta method, namely IRK2, IRK4 and IRK6 from their order of accuracy when $s=1,2,3$ , respectively. There are many other types of collocation Runge-Kutta methods, we refer the interested reader to [12], [41] and [42].

\begin{array}[]{c|c}\frac{1}{2}&\frac{1}{2}\\ \hline\cr&1\\ \end{array}

(a) IRK2

\begin{array}[]{c|cc}\frac{1}{2}-\frac{1}{6}\sqrt{3}&\frac{1}{4}&\frac{1}{4}-\frac{1}{6}\sqrt{3}\\ \frac{1}{2}+\frac{1}{6}\sqrt{3}&\frac{1}{4}+\frac{1}{6}\sqrt{3}&\frac{1}{4}\\ \hline\cr&\frac{1}{2}&\frac{1}{2}\\ \end{array}

(b) IRK4

\begin{array}[]{c|ccc}\frac{1}{2}-\frac{\sqrt{15}}{10}&\frac{5}{36}&\frac{2}{9}-\frac{\sqrt{15}}{15}&\frac{5}{36}-\frac{\sqrt{15}}{30}\\ \frac{1}{2}&\frac{5}{36}+\frac{\sqrt{15}}{24}&\frac{2}{9}&\frac{5}{36}-\frac{\sqrt{15}}{24}\\ \frac{1}{2}+\frac{\sqrt{15}}{10}&\frac{5}{36}+\frac{\sqrt{15}}{30}&\frac{2}{9}+\frac{\sqrt{15}}{15}&\frac{5}{36}\\ \hline\cr&\frac{5}{18}&\frac{4}{9}&\frac{5}{18}\\ \end{array}

Table 1. Butcher’s Tableaus for the

s

-stage Gaussian-Legendre collocation Runge-Kutta methods with

s=1,2,3

Now, we adapt the collocation Runge-Kutta methods into the equation (2.2). Denote the semi-discretization in time $u^{m}\approx u(x,t_{m})$ , $v^{m}\approx v(t_{m})$ . For simplicity, we rewrite the equation (2.2) as the system

(3.1)

\displaystyle\begin{cases}u_{t}=f(u,v),\\ v_{t}=g(u,v).\end{cases}

Denote the intermediate values $U_{i}=u(t_{m}+c_{i}\tau)$ and $V_{i}=v(t_{m}+c_{i}\tau)$ to be the solution satisfying (2.2) at the time $t=t_{m}+c_{i}\tau$ . Then, $U_{i}$ and $V_{i}$ can be obtained from the relation

(3.2)

\displaystyle U_{i}={u}^{m}+\tau\sum_{j=1}^{s}a_{ij}{f}_{j},\quad V_{i}={v}^{m}+\tau\sum_{j=1}^{s}a_{ij}{g}_{j},

where ${f}_{i}=f(U_{i},V_{i})$ and ${g}_{i}=g(U_{i},V_{i})$ . Then, ${u}^{m+1}$ and $v^{m+1}$ are updated by

(3.3)

\displaystyle{u}^{m+1}={u}^{m}+\tau\sum_{i=1}^{s}b_{i}{f}_{i},\quad v^{m+1}=v^{m}+\tau\sum_{i=1}^{s}b_{i}{g}_{i}.

Denote $I^{m}=I[u^{m}],\,M^{m}=M[u^{m}]$ and $E^{m}=E[u^{m},v^{m}]$ to be the invariant quantities from (2.3)-(2.5) at the time $t=t_{m}$ . We have the following theorem.

Theorem 3.1.

The $s$ -stage collocation Runge-Kutta method preserves the momentum (2.3), i.e.,

I^{m+1}=I^{m}.

Furthermore, the $s$ -stage symplectic Runge-Kutta method (which is also called the quadratic preserving Runge-Kutta method) satisfies

(3.4)

\displaystyle b_{i}a_{ij}+b_{j}a_{ji}-b_{i}b_{j}=0,\qquad\mbox{for}\quad i,j=1,\cdots,s,

conserves the mass (2.4) and energy (2.5) exactly in the discrete time flow, i.e.,

M^{m+1}=M^{m},\quad\mbox{and}\quad E^{m+1}=E^{m}.

Proof.

The proof is standard, similar to [12], [41] and [37]. Putting (3.3) in (2.3) yields

\displaystyle I^{m+1}-I^{m}=(u^{m+1}-u^{m},1)=\tau(\sum_{i=1}^{s}b_{i}{f}_{i},1)=\tau\sum_{i=1}^{s}b_{i}({f}_{i},1)=0.

The last equality comes from the conservation law $(f_{i},1)=0$ from (2.2) at each collocation time $t=t_{m}+c_{i}\tau$ .

Similarly, putting (3.3) in (2.4), and using (3.2), yields

	$\displaystyle M^{m+1}$	$\displaystyle=({u}^{m+1},{u}^{m+1})=({u}^{m}+\tau\sum_{i=1}^{s}b_{i}{f}_{i},{u}^{m}+\tau\sum_{i=1}^{s}b_{i}{f}_{i})$
		$\displaystyle=({u}^{m},{u}^{m})-2\tau\sum_{i=1}^{s}b_{i}({f_{i}},\tau\sum_{j=1}^{s}a_{ij}{f}_{j}-U_{i})+\tau^{2}\sum_{i,j=1}^{s}b_{i}b_{j}(f_{i},f_{j})$
		$\displaystyle=M^{m}+2\tau\sum_{i=1}^{s}b_{i}({f}_{i},U_{i})-\tau^{2}\sum_{i,j=1}^{s}(b_{i}a_{ij}+b_{j}a_{ji}-b_{i}b_{j})({f}_{i},{f}_{j})=M^{m},$

since $(f_{i},U_{i})=(u_{t},u)=\frac{1}{2}\frac{d}{dt}M=0$ at each collocation time $t=t_{m}+c_{i}\tau$ from the conservation law in (2.4).

Similarly, to show $E^{m+1}=E^{m}$ , we have

	$\displaystyle E^{m+1}$	$\displaystyle=-\frac{1}{2}(u_{xx}^{m+1},{u}^{m+1})-\frac{1}{p(p+1)}((v^{m+1})^{2}-C_{0})$
		$\displaystyle=-\frac{1}{2}(({u}^{m}+\tau\sum_{i=1}^{s}b_{i}{f}_{i})_{xx},{u}^{m}+\tau\sum_{i=1}^{s}b_{i}{f}_{i})-\frac{1}{p(p+1)}((v^{m}+\tau\sum_{i=1}^{s}b_{i}{g}_{i})^{2}-C_{0})$
		$\displaystyle=-\frac{1}{2}(u_{xx}^{m},{u}^{m})-\frac{1}{p(p+1)}\left((v^{m})^{2}-C_{0}\right)+\tau\sum_{i=1}^{s}b_{i}(\partial_{x}^{2}f_{i},U_{i})+\tau\sum_{i=1}^{s}b_{i}\left(\frac{2}{p(p+1)}V_{i}g_{i}\right)$
		$\displaystyle+\tau^{2}\sum_{i,j=1}^{s}(b_{i}a_{ij}+b_{j}a_{ji}-b_{i}b_{j})\left((\partial_{x}^{2}f_{i},f_{j})+{\frac{1}{p(p+1)}}g_{i}g_{j}\right)$
		$\displaystyle=E^{m}$

by $(\partial_{x}^{2}f_{i},U_{i})+\frac{2}{p(p+1)}V_{i}g_{i}=0$ for each $i$ , which follows from differentiating the conservation law (2.5) with respect to $t$ . This completes the proof. ∎

From the theorem, we know that any arbitrary order of collocation Runge-Kutta method will conserve the momentum (2.3), the symplectic Runge-Kutta method will conserve the mass (2.4) and modified energy (2.5) in the discrete time flow. Therefore, an arbitrarily high order time integrator can be constructed. In fact, from the proof of Theorem 3.1, reformulating the original equation (1.1) into the system (2.2) is only for the purpose of the energy conservation.

3.2. A $C_{0}$ adjustment process

The key part for this SAV approach is to make sure that the term $\int u^{p+1}dx+C_{0}$ is positive during the computational time $t\in[0,T]$ . When considering the mKdV ( $p=3$ ) case, $\int u^{p+1}dx+C_{0}>0$ will automatically hold for any positive constant $C_{0}$ , since $p+1=4$ is even. This is also true for considering the nonlinear Schrödinger (NLS) equations or the Gross-Pitaevskii (GP) equations in [15] and [37], as the authors there only consider the cubic nonlinearity $|u|^{2}u$ . However, when considering the KdV ( $p=2$ ) case, a prior bound for $\int u^{p+1}(t)dx$ is needed for all $t\in[0,T]$ . Otherwise, $\int u^{p+1}(t_{m})dx+C_{0}<0$ may happen at some time $t=t_{m}$ (though $\int u^{p+1}(t_{0})dx+C_{0}>0$ is satisfied in the beginning), and result in the failure of the algorithm. Instead of choosing $C_{0}$ carefully in the beginning of the simulation, we introduce a $C_{0}$ adjustment process when the value $\int u^{p+1}(t_{m})dx+C_{0}$ is approaching $0$ , and thus, we only need to choose the constant $C_{0}$ such that $\int u^{p+1}(t_{0})dx+C_{0}>0$ in the beginning of the computation.

Suppose at $t=t_{m}$ , $\int(u^{m})^{p+1}dx+C_{0}<Tol$ , where $Tol$ is a given positive number (e.g., $Tol=5$ ). Then, we choose another constant $\tilde{C_{0}}$ such that $\int(u^{m})^{p+1}dx+\tilde{C_{0}}>Tol$ . For example, we can take $\tilde{C_{0}}=10-\int(u^{m})^{p+1}dx$ , which leads to our new $\tilde{v}^{m}\approx\sqrt{10}$ . Then, by using $E[u^{m},v^{m}]=E[u^{m},\tilde{v}^{m}]$ from (2.5), we have our new $\tilde{v}^{m}$

(3.5)

\displaystyle\tilde{v}^{m}=\sqrt{(v^{m})^{2}+\tilde{C_{0}}-C_{0}}.

Finally, we substitute the $v^{m}$ and $C_{0}$ in (4.5) with $\tilde{v}^{m}$ and $\tilde{C_{0}}$ , and then, continue with the time evolution for $t=t_{m+1},t_{m+2},\cdots$ .

Remark 3.1.

Note that $v^{2}=\int u^{p+1}dx+C_{0}$ holds only at the collocation points $t=t_{m}+\tau c_{i}$ for each $i=1,2,\cdots,s$ in $t\in[t_{m},t_{m+1}]$ . However, the constant $c_{s}$ may not necessarily equal to $0$ or $1$ , e.g., see Table 1, which means $v^{2}=\int u^{p+1}dx+C_{0}$ does not hold at $t_{m}$ in the discrete time flow. Therefore, the new $\tilde{v}^{m}$ can only be evaluated by (3.5) to keep the discrete energy (2.5) invariant. In the actual computation, due to the dispersive nature of the negative data (e.g., Example 3), the $C_{0}$ adjustment process has not been executed yet.

4. Spatial discretization

In this section, we describe the spatial discretization. We show that the Fourier pseudo-spectral method will keep the momentum and energy invariant under such spatial discretization. Together with the temporal discretization in the previous section, the proposed scheme will preserve the discrete momentum and energy exactly in the discrete time flow. Moreover, we give an upper bound for the error of the discrete mass, which is caused by the spatial discretization.

Without loss of generality, we truncate the whole space into a bounded domain $x\in[-L,L]$ for sufficiently large $L$ with periodic boundary conditions at the boundary. We use the Fourier pseudo-spectral method to discretize the space because of the high-order accuracy and the application of the fast Fourier transform (FFT) (see. e.g., [44] and [48, Chapter 3]).

Now, we briefly introduce the spatial discretization strategy. Let $N$ be the number of nodes and $h=2L/N$ to be the spatial step size. Denote $x_{j}=hj$ and $u_{j}\approx u(x_{j})$ to be the spatial discretization for $j=-N/2,\cdots,N/2-1$ . By applying the standard discrete Fourier expansion, one obtains

(4.1)

\displaystyle u_{j}=\frac{1}{\sqrt{N}}\sum_{k=-N/2}^{N/2-1}\hat{u}_{k}e^{i2k\pi x_{j}/L},\quad\mbox{where}\quad\hat{u}_{k}=\frac{1}{\sqrt{N}}\sum_{j=-N/2}^{N/2-1}u_{j}e^{-i2k\pi x_{j}/L}.

By introducing vectors $\mathbf{u}=(u_{-N/2},\cdots,u_{N/2-1})^{T}$ , $\mathbf{\hat{u}}=(\hat{u}_{-N/2},\cdots,\hat{u}_{N/2-1})^{T}$ , and the discrete Fourier transform matrices

\mathbf{F}_{k,j}=\frac{1}{\sqrt{N}}e^{-i2k\pi x_{j}/L},\quad\mathbf{F}^{-1}_{j,k}=\frac{1}{\sqrt{N}}e^{i2k\pi x_{j}/L},\quad-N/2\leq j,k\leq N/2-1,

one can write $\mathbf{u}$ and $\mathbf{\hat{u}}$ in a matrix form as $\mathbf{u}=\mathbf{F}^{-1}\mathbf{\hat{u}}$ and $\mathbf{\hat{u}}=\mathbf{Fu}$ . Also note that $\mathbf{F}^{-1}=(\mathbf{\bar{F}})^{T}.$ Using the properties of the Fourier transform, one can easily obtain

\displaystyle u^{\prime}_{j}=\frac{i2\pi}{L}\sum_{k=-N/2}^{N/2-1}k\hat{u}_{k}e^{i2k\pi x_{j}/L},\quad u^{\prime\prime}_{j}=-\left(\frac{2\pi}{L}\right)^{2}\sum_{k=-N/2}^{N/2-1}k^{2}\hat{u}_{k}e^{i2k\pi x_{j}/L}.

Thus, we have the differential matrices $\mathbf{D_{1}}$ and $\mathbf{D_{2}}$ :

(4.2)

\displaystyle\mathbf{u}^{\prime}=\mathbf{D_{1}u}=\mathbf{\bar{F}^{T}\Lambda_{1}Fu},\quad\mbox{and}\quad\mathbf{u}^{\prime\prime}=\mathbf{D_{2}u}=\mathbf{\bar{F}^{T}\Lambda_{2}Fu},

where $\mathbf{\Lambda_{1}}=i\mathbf{\Lambda}$ , $\mathbf{\Lambda_{2}}=\mathbf{\Lambda^{2}_{1}}$ , and $\mathbf{\Lambda}=\frac{2\pi}{L}\mbox{diag}(\frac{-N}{2},\cdots,\frac{N-1}{2})$ . We note that $\mathbf{D_{1,2}}$ are real matrices. Moreover, the matrix $\mathbf{D_{2}}$ is symmetric. The matrix $\mathbf{D_{1}}$ is antisymmetric ( $\mathbf{D_{1}}=-\mathbf{D_{1}^{T}}$ ) and circulant, i.e., $\mathbf{D_{1}}_{(j,k)}=d_{j-k}$ . Also note that both the row sum and the column sum of the matrix $\mathbf{D_{1}}$ equal to $0$ , i.e., for each fixed $k$ or $j$ ,

(4.3)

\displaystyle\sum_{j=-N/2}^{N/2-1}d_{j,k}=\sum_{k=-N/2}^{N/2-1}d_{j,k}=0.

This has been studied in many literatures, see e.g., [19] and [48, Chapter 3]. We remark here that in the actual computation, $\mathbf{F}$ and $\mathbf{F}^{-1}$ are not explicitly assembled, instead, the fast Fourier transform (FFT) is typically used in the computation.

For simplicity, if $\mathbf{u}$ is a vector, we denote $\mathbf{u^{p}}=(u_{-N/2}^{p},u_{-N/2+1}^{p},\cdots u_{N/2-1}^{p})^{T}$ to be the pointwise power. Then, we define the discrete norms. Given two vectors $\mathbf{u},\mathbf{v}$ , define

(4.4)

\displaystyle(\mathbf{u},\mathbf{v})_{h}=h\mathbf{\bar{v}}^{T}\mathbf{u}=h\sum_{j=-N/2}^{N/2-1}u_{j}\bar{v}_{j},\quad\|\mathbf{u}\|_{h}=\sqrt{(\mathbf{u},\mathbf{u})_{h}}.

Note that when $\mathbf{u},\mathbf{v}\in\mathbb{R}$ , $(\mathbf{u},\mathbf{v})_{h}=(\mathbf{v},\mathbf{u})_{h}$ . It is easy to see that $(\mathbf{D_{2}u},\mathbf{u})_{h}=(\mathbf{u},\mathbf{D_{2}u})_{h}\in\mathbb{R}$ , since $\mathbf{D_{2}}$ is real and symmetric, and $(\mathbf{D_{1}u},\mathbf{u})_{h}=(\mathbf{u},\mathbf{D_{1}u})_{h}=0$ , since $\mathbf{D_{1}}$ is antisymmetric.

Denote $v_{h}$ to be the spatial semi-discretization for $v(t)$ , i.e., $v_{h}\approx v(t)$ . Then, the system of equations (2.2) is discretized into the following system

(4.5)

\displaystyle\begin{cases}\mathbf{u}_{t}=-\mathbf{D_{1}}\left(\mathbf{D_{2}u}+\frac{1}{p}\frac{\mathbf{u^{p}}v_{h}}{\sqrt{(\mathbf{u^{p}},\mathbf{u})_{h}+C_{0}}}\right)\equiv f(\mathbf{u},v),\\ (v_{h})_{t}=\frac{p+1}{2\sqrt{(\mathbf{u^{p}},\mathbf{u})_{h}+C_{0}}}(\mathbf{u^{p}},\mathbf{u}_{t})_{h}\equiv g(\mathbf{u},v).\end{cases}

Note that $\mathbf{u}$ and $f(\mathbf{u},v_{h})$ are an $N\times 1$ vectors, and $v_{h}$ and $g(\mathbf{u},v_{h})$ are scalars.

Now, the spatial discrete momentum, mass and energy are defined as follows

(4.6)		$\displaystyle I_{h}(\mathbf{u})=(\mathbf{u},\mathbf{1})_{h};$
(4.7)		$\displaystyle M_{h}(\mathbf{u})=(\mathbf{u},\mathbf{u})_{h};$
(4.8)		$\displaystyle E_{h}(\mathbf{u},v_{h})=-\frac{1}{2}(\mathbf{D_{2}u},\mathbf{u})_{h}-\frac{1}{p(p+1)}(v_{h}^{2}-C_{0}),$

where $\mathbf{1}=(1,1,\cdots,1)^{T}$ is an $N\times 1$ vector. We next obtain the conservation of the discrete momentum and energy, together with the error on the discrete mass.

Theorem 4.1.

The Fourier pseudo-spectral discretization to the equation (4.5) leads to the following properties:

(4.9)

\displaystyle\frac{d}{dt}I_{h}(\mathbf{u})=0;\quad\frac{d}{dt}E_{h}(\mathbf{u},v_{h})=0.

Moreover,

(4.10)

\displaystyle\frac{d}{dt}M_{h}(\mathbf{u})=-\frac{2h}{p}\mathbf{u^{T}D_{1}u^{p}}.

Proof.

Note that $I_{h}\in\mathbb{R}$ , thus,

\displaystyle\frac{d}{dt}I_{h}(\mathbf{u})=(\mathbf{u}_{t},\mathbf{1})_{h}=(\mathbf{1},\mathbf{u}_{t}).

We note that from (4.3), we have

\mathbf{1^{T}D_{1}}=\left(\sum_{k=-N/2}^{N/2-1}d_{-N/2,k},\sum_{j=-N/2}^{N/2-1}d_{-N/2+1,k},\cdots,\sum_{j=-N/2}^{N/2-1}d_{N/2-1,k}\right)=\mathbf{0^{T}}.

Therefore, substitute $\mathbf{u}_{t}$ by (4.5) yields

(4.11)

\displaystyle(\mathbf{u}_{t},\mathbf{1})_{h}=h\mathbf{1^{T}}\mathbf{u}_{t}=-h\mathbf{1^{T}}\mathbf{D_{1}}\left(\mathbf{D_{2}u}+\frac{1}{p}\frac{\mathbf{u^{p}}v_{h}}{\sqrt{(\mathbf{u^{p}},\mathbf{u})_{h}+C_{0}}}\right)=0.

Next, we prove $\frac{d}{dt}E_{h}=0$ , provided by $v_{h}(v_{h})_{t}=\frac{p+1}{2}(\mathbf{u^{p}},\mathbf{u}_{t})_{h}$ from (4.5), i.e.,

(4.12)		$\displaystyle\frac{d}{dt}E_{h}$	$\displaystyle=-(\mathbf{D_{2}u},\mathbf{u}_{t})_{h}-\frac{2}{p(p+1)}v_{h}(v_{h})_{t}=-\left((\mathbf{D_{2}u},\mathbf{u}_{t})_{h}+\frac{1}{p}(\mathbf{u^{p}},\mathbf{u}_{t})_{h}\right)$
		$\displaystyle=\left(\mathbf{D_{2}u}+\frac{1}{p}\mathbf{u^{p}},\mathbf{D_{1}}(\mathbf{D_{2}u}+\frac{1}{p}\mathbf{u^{p}})\right)_{h}=0,$

since $\mathbf{D_{1}}$ is antisymmetric.

For proving (4.10), we first note that $\mathbf{D_{1}D_{2}}=\mathbf{D_{2}D_{1}}$ from the decomposition in (4.2). Then, $(\mathbf{D_{1}D_{2}u},\mathbf{u})_{h}=-(\mathbf{D_{2}u},\mathbf{D_{1}u})_{h}$ from the antisymmetry of $\mathbf{D_{1}}$ . On the other hand, $(\mathbf{D_{1}D_{2}u},\mathbf{u})_{h}=(\mathbf{D_{2}D_{1}u},\mathbf{u})_{h}=(\mathbf{D_{2}u},\mathbf{D_{1}u})_{h}$ . Thus, $(\mathbf{D_{1}D_{2}u},\mathbf{u})_{h}=0$ . Also, note that $\frac{v_{h}}{\sqrt{(\mathbf{u^{p}},\mathbf{u})_{h}+C_{0}}}=1$ in the continuous time. Then,

(4.13)		$\displaystyle\frac{d}{dt}M_{h}$	$\displaystyle=2(\mathbf{u_{t}},\mathbf{u})_{h}=-2(\mathbf{D_{1}D_{2}u},\mathbf{u})_{h}-\frac{2}{p}(\mathbf{D_{1}u^{p}},\mathbf{u})_{h}$
		$\displaystyle=-\frac{2}{p}(\mathbf{D_{1}u^{p}},\mathbf{u})_{h}=-\frac{2h}{p}\mathbf{u^{T}D_{1}u^{p}},$

which completes the proof.

∎

Theorem 4.1 shows that in the continuous time flow, the spatial discretization will keep the momentum (2.3) and energy (2.5) invariant in their corresponding discrete sense. However, the time derivative on the discrete mass (4.10) is not necessarily $0$ . We next show that the proposed scheme conserves the discrete momentum and energy, and also give an error estimate of the discrete mass. From Theorem 4.1, we have the following little corollary.

Corollary 4.1.

Suppose $(\mathbf{u},v_{h})$ is the solution to (4.5), then we have

(4.14)		$\displaystyle({f}(\mathbf{u},v_{h}),\mathbf{1})_{h}=0,$
(4.15)	$\displaystyle-$	$\displaystyle(\mathbf{D_{2}u},{f}(\mathbf{u},v_{h}))_{h}-\frac{2}{p(p+1)}v_{h}g(\mathbf{u},v_{h})=0,$
(4.16)		$\displaystyle({f}(\mathbf{u},v_{h}),\mathbf{u})_{h}=-\frac{2h}{p}\mathbf{u^{T}D_{1}u^{p}}.$

Proof.

Substituting $\mathbf{u}_{t}={f}(\mathbf{u},v_{h})$ in the equation (4.11) and (4.13) yields (4.14) and (4.16), respectively. From (4.12), we have

\displaystyle 0=\frac{d}{dt}E_{h}=-(\mathbf{D_{2}u},\mathbf{u}_{t})_{h}-\frac{2}{p(p+1)}v_{h}(v_{h})_{t}=-(\mathbf{D_{2}u},\mathbf{f}(\mathbf{u},v_{h}))_{h}-\frac{2}{p(p+1)}v_{h}g(\mathbf{u},v_{h}),

which completes the proof. ∎

Define the full discretization in space and time $u_{j}^{m}\approx u(x_{j},t_{m})$ and $v_{h}^{m}\approx v(t_{m})$ . Then, denote $\mathbf{u}^{m}$ and $v_{h}^{m}$ to be the numerical solution to (4.5). Define the discrete momentum, mass and energy as follows:

(4.17)		$\displaystyle I_{h}^{m}=(\mathbf{u}^{m},\mathbf{1})_{h};$
(4.18)		$\displaystyle M_{h}^{m}=(\mathbf{u}^{m},\mathbf{u}^{m})_{h};$
(4.19)		$\displaystyle E_{h}^{m}=-\frac{1}{2}(\mathbf{D_{2}u}^{m},\mathbf{u}^{m})_{h}-\frac{1}{p(p+1)}\left((v_{h}^{m})^{2}-C_{0}\right).$

The following theorem shows the conservation of the discrete momentum (4.17) and the discrete energy (4.19) for the proposed scheme. It also shows the upper bound of the error for the discrete mass (4.18).

Theorem 4.2.

Suppose the discretized system (4.5) is solved by the $s$ -stage symplectic Runge-Kutta method which satisfies (3.4), then the solution $\mathbf{u}^{m}$ and $v_{h}^{m}$ to the equation (4.5) satisfies

(4.20)

\displaystyle I_{h}^{m+1}=I_{h}^{m};\quad E_{h}^{m+1}=E_{h}^{m}.

In the meanwhile, denote $U_{i}^{m}$ and $V_{i}^{m}$ to be the intermediate step for the SRK method at $t=t_{m}+c_{i}\tau$ . Also note that $U_{i}^{m}\approx\mathbf{u}^{m}$ , then

(4.21)

\displaystyle|M_{h}^{m}-M_{h}^{0}|\leq t_{m}\frac{4h}{p}\max_{1\leq\nu\leq m,1\leq i\leq s}|(U^{\nu}_{i})^{\mathbf{T}}\mathbf{D_{1}}(U^{\nu}_{i})^{\mathbf{p}}|\approx t_{m}\frac{4h}{p}\max_{1\leq\nu\leq m}|(\mathbf{u^{T}})^{\nu}\mathbf{D_{1}}(\mathbf{u^{p}})^{\nu}|.

Proof.

The conservation (4.20) are obtained from substituting the continuous variables and inner products in Theorem 3.1 with the spatial discrete sense (4.4), and substitute $\mathbf{u}$ in the equations (4.14) and (4.15) in Corollary 4.1 by $U_{i}^{m}$ and $V_{i}^{m}$ , since they are intermediate steps which satisfy the equation (4.5).

For (4.21), by substituting the continuous variables and inner products into the spatial discrete sense and use (4.16) yields

	$\displaystyle M_{h}^{m}-M_{h}^{m-1}$	$\displaystyle=2\tau\sum_{i=1}^{s}b_{i}(f_{i},U_{i}^{m-1})_{h}=-2\tau\sum_{i=1}^{s}b_{i}\frac{2h}{p}(U^{m-1}_{i})^{\mathbf{T}}\mathbf{D_{1}}(U^{m-1}_{i})^{\mathbf{p}}$
		$\displaystyle\leq\frac{4\tau h}{p}\max_{1\leq i\leq s}\|(U^{m-1}_{i})^{\mathbf{T}}\mathbf{D_{1}}(U^{m-1}_{i})^{\mathbf{p}}\|,$

since $\sum_{i=1}^{s}b_{i}=1$ . Therefore,

	$\displaystyle\|M_{h}^{m}-M_{h}^{0}\|$	$\displaystyle=\|\sum_{\nu=1}^{m}(M_{h}^{\nu}-M_{h}^{\nu-1})\|\leq\sum_{\nu=1}^{m}\|M_{h}^{\nu}-M_{h}^{\nu-1}\|\leq\frac{4\tau h}{p}\sum_{\nu=1}^{m}\max_{1\leq i\leq s}\|(U^{\nu}_{i})^{\mathbf{T}}\mathbf{D_{1}}(U^{\nu}_{i})^{\mathbf{p}}\|$
		$\displaystyle\leq\frac{4\tau h}{p}\sum_{\nu=1}^{m}\max_{1\leq\nu\leq m,1\leq i\leq s}\|(U^{\nu}_{i})^{\mathbf{T}}\mathbf{D_{1}}(U^{\nu}_{i})^{\mathbf{p}}\|=\frac{4t_{m}h}{p}\max_{1\leq\nu\leq m,1\leq i\leq s}\|(U^{\nu}_{i})^{\mathbf{T}}\mathbf{D_{1}}(U^{\nu}_{i})^{\mathbf{p}}\|,$

the last equality comes from $t_{m}=m\tau$ and completes the prove. ∎

Theorem 4.2 shows that the discrete momentum and energy will be preserved exactly in the proposed scheme. In the meanwhile, the inequality (4.21) shows the upper bound of the error for the discrete mass, which grows linearly with respect to time. However, note that the term $h\mathbf{u^{T}}\mathbf{D_{1}}\mathbf{u^{p}}$ is the Fourier spectral approximation of the integral $\frac{1}{p+1}\int(u^{p+1})_{x}dx$ , i.e.,

h\mathbf{u^{T}}\mathbf{D_{1}}\mathbf{u^{p}}\approx((u^{p})_{x},u)=\frac{1}{p+1}((u^{p+1})_{x},1)=0.

Hence, instead of decreasing in the 1st order with respect to $h$ , the term $|h\mathbf{u^{T}}\mathbf{D_{1}}\mathbf{u^{p}}|$ is of spectral accuracy and very small (generally on the order of $10^{-14}$ ), which is usually several order lower than the tolerance of the fixed point solver for the resulting nonlinear system. Moreover, similar to previous results observed in [36] for the discontinuous Galerkin scheme, despite of the failure of the mass conservation, the error of discrete mass remains at the same level of the conserved quantities (discrete momentum and energy) in the simulation. These conserved discrete momentum and energy are likely to imply the cancellations of the mass error. As a consequence, the error of the discrete mass becomes unnoticeable during our simulation even up to the time $t_{m}=200$ , see details in the next section.

5. Numerical results

5.1. A fast solver

In this section, we show numerical examples for the proposed scheme. Before illustrating the examples, we describe a fast solver for solving the resulting nonlinear system from the IRK methods. This fast solver is similar to the one in [15], which is based on the fixed point iteration. It is on the same order of the computational cost in solving the original equation, despite of the new auxiliary variable being introduced. To be specific, we consider the IRK4 ( $s=2$ ) case. The other cases can be easily generalized. We introduce the auxiliary variables and rewrite the system as

(5.1)

\displaystyle\begin{cases}U_{i}=\mathbf{u}^{m}+\tau\sum_{j=1}^{2}a_{ij}\mathbf{f}_{j},\\ \Phi_{i}=\frac{(U_{i})^{p}}{\sqrt{((U_{i})^{p},U_{i})_{h}+C_{0}}},\quad g_{i}=\frac{p+1}{2}(\Phi_{i},\mathbf{f}_{i})_{h},\\ V_{i}={v_{h}}^{m}+\tau\sum_{j=1}^{2}a_{ij}{g}_{j},\quad i=1,2,\end{cases}

and

(5.2)

\displaystyle\begin{cases}\mathbf{f}_{1}=-\mathbf{D_{1}}\left(\mathbf{D_{2}}(\mathbf{u}^{m}+\tau\sum_{j=1}^{2}a_{1j}\mathbf{f}_{j})+\frac{1}{p}\Phi_{1}V_{1})\right),\\ \mathbf{f}_{2}=-\mathbf{D_{1}}\left(\mathbf{D_{2}}(\mathbf{u}^{m}+\tau\sum_{j=1}^{2}a_{2j}\mathbf{f}_{j})+\frac{1}{p}\Phi_{2}V_{2})\right).\end{cases}

Then, $\mathbf{u}^{m+1}$ and $v_{h}^{m+1}$ can be updated by (3.3).

The system (5.1) and (5.2) can be solved by the fixed point iteration. At the $l$ th iteration, we have

(5.3)

\displaystyle\begin{cases}\left(\mathbf{I}+\tau a_{11}\mathbf{D_{3}}\right)\mathbf{f}_{1}^{l+1}+\tau a_{12}\mathbf{D_{3}f}_{2}^{l+1}=-(\mathbf{D_{3}u}^{m}-\frac{1}{p}\mathbf{D_{1}}(\Phi_{1}^{l}V_{1}^{l})),\\ \tau a_{21}\mathbf{D_{3}f}_{1}^{l+1}+\left(\mathbf{I}+\tau a_{22}\mathbf{D_{3}}\right)\mathbf{f}_{2}^{l+1}=-(\mathbf{D_{3}u}^{m}-\frac{1}{p}\mathbf{D_{1}}(\Phi_{2}^{l}V_{2}^{l})),\end{cases}

where $\mathbf{I}$ is the identity matrix and $\mathbf{D_{3}}=\mathbf{D_{1}D_{2}}$ .

Note that from the system (5.3), each block in the matrix

\mathbf{M}=\begin{pmatrix}\mathbf{I}+\tau a_{11}\mathbf{D_{3}}&\tau a_{12}\mathbf{D_{3}}\\ \tau a_{21}\mathbf{D_{3}}&\mathbf{I}+\tau a_{22}\mathbf{D_{3}}\end{pmatrix}

is commute to each other (i.e., $\mathbf{M}_{i,j}\mathbf{M}_{k,l}=\mathbf{M}_{k,l}\mathbf{M}_{i,j}$ ). Therefore, $\mathbf{M}$ can be easily inverted:

\mathbf{M}^{-1}=\begin{pmatrix}\mathbf{J}^{-1}(\mathbf{I}+\tau a_{22}\mathbf{D_{3}})&-\mathbf{J}^{-1}(\tau a_{12}\mathbf{D_{3}})\\ -\mathbf{J}^{-1}(\tau a_{21}\mathbf{D_{3}})&\mathbf{J}^{-1}(\mathbf{I}+\tau a_{11}\mathbf{D_{3}})\end{pmatrix},

where $\mathbf{J}=(\mathbf{I}+\tau a_{11}\mathbf{D_{3}})(\mathbf{I}+\tau a_{22}\mathbf{D_{3}})-\tau^{2}a_{12}a_{21}\mathbf{D_{3}D_{3}}$ .

Next, note that the matrix $\mathbf{M}^{-1}$ can be further decomposed by FFT. Therefore, the system (5.3) will be solved efficiently with the leading order of computational cost $\mathcal{O}(2N\log(N))$ . After obtain $\mathbf{f}_{i}^{l+1}$ , we can use (5.1) to obtain the values of other variables $U_{i}^{l+1},\Phi_{i}^{l+1},g_{i}^{l+1},$ and $V_{i}^{l+1}$ . The total computational cost remains on the same order.

We set $\mathbf{f}_{1}^{0}=\mathbf{f}_{2}^{0}=f(\mathbf{u}^{m},v^{m}_{h})$ to start the iteration, and the iteration terminates when

\max_{i}\|\mathbf{f}_{i}^{l+1}-\mathbf{f}_{i}^{l}\|_{\infty}<\epsilon,

where $i=1,2$ and $\epsilon$ is set to be $10^{-11}$ for Example 1, or $10^{-12}$ for Example 2 and 3 in our simulations.

In general, for the $s$ -stage Runge-Kutta method, we can always find the explicit $\mathbf{M}^{-1}$ for the system similar to (5.3). The computational cost is consequently on the order of $\mathcal{O}(s\times N\log(N))$ . The $s$ -stage Gaussian-Legendre collocation methods will obtain the maximum order of $2s$ on temporal accuracy.

5.2. Brief description of existing time integrators

Since we are going to compare our proposed methods to other time integrators, we briefly describe the other time integrators we used in this paper. For simplicity, we only discuss the temporal discretization in this subsection. The space-time full discretization can be easily adapted by applying the Fourier pseudo-spectral spatial discretization as described in the previous section.

5.2.1. Modified Crank-Nicholson method

We first review the widely used 2nd order conservative modified Crank-Nicolson (MCN) method as follows:

(5.4)

\displaystyle\frac{u^{m+1}-u^{m}}{\tau}+\left(\frac{1}{2}(u^{m+1}_{xx}+u^{m}_{xx})+\frac{1}{p(p+1)}\frac{(u^{m+1})^{p+1}-(u^{m})^{p+1}}{(u^{m+1})^{2}-(u^{m})^{2}}(u^{m+1}+u^{m})\right)_{x}=0.

The nonlinear system (5.4) can be solved by the fixed point iteration in the same spirit from (5.3). We omit the details for the purpose of conciseness.

Proposition 5.1.

The modified Crank-Nicolson scheme (5.4) preserves the following two invariants:

I_{h}^{m+1}=I^{m}\quad E_{h}^{m+1}=E^{m}.

Consequently, if the the scheme is discretized by the Fourier pseudo-spectral method (section 4), then the discrete momentum and energy will be preserved ( $I_{h}^{m+1}=I_{h}^{m}$ and $E_{h}^{m+1}=E_{h}^{m}$ ).

Proof.

By taking the inner product with $1$ and the term

\frac{1}{2}(u^{m+1}_{xx}+u^{m}_{xx})+\frac{1}{p(p+1)}\frac{(u^{m+1})^{p+1}-(u^{m})^{p+1}}{(u^{m+1})^{2}-(u^{m})^{2}}(u^{m+1}+u^{m}),

we obtain that three invariants $I^{m}$ and $E^{m}$ are preserved in the temporal semi-discretized sense (consider the spatial variable is continuous), respectively. Furthermore, by using the argument in the previous section, we can obtain $I_{h}^{m+1}=I_{h}^{m}$ and $E_{h}^{m+1}=E_{h}^{m}$ . ∎

5.2.2. Semi-implicit Scalar auxiliary variable method

We next consider the original semi-implicit SAV approach from [46]. The main advantage is that solving the resulting nonlinear algebraic system can be avoided. Other than the conventional Crank-Nicholson-Adam-Bashforth (CNAB) approach, we use the Leap-Frog scheme here, since it is time reversible. To be specific, we discretize the equation (2.2) as follows

(5.5)

\displaystyle\begin{cases}\frac{u^{m+1}-u^{m-1}}{2\tau}+\tilde{u}^{m}_{xxx}+\left(\frac{1}{p}q^{m}\tilde{v}^{m}\right)_{x}=0,\\ v^{m+1}-v^{m-1}=\frac{p+1}{2}(q^{m},u^{m+1}-u^{m-1}),\end{cases}

where $\tilde{u}^{m}=\frac{1}{2}(u^{m+1}+u^{m-1})$ , $q^{m}=q(u^{m})=\frac{(u^{m})^{p}}{\sqrt{((u^{m})^{p},u^{m})+C_{0}}}$ , and $\tilde{v}^{m}=\frac{1}{2}(v^{m+1}+v^{m-1})$ .

In order to solve this equation system, we rewrite the first equation in (5.5) as

(5.6)

\displaystyle\frac{\tilde{u}^{m}-u^{m-1}}{\tau}+\tilde{u}^{m}_{xxx}+\left(\frac{1}{p}q^{m}\tilde{v}^{m}\right)_{x}=0,

Assuming $\tilde{u}^{m}=u_{1}^{m}+\tilde{v}^{m}u_{2}^{m}$ , putting it in (5.6) yields the linear decoupled system

(5.7)

\displaystyle\begin{cases}(1+\tau\partial_{x}^{3})u_{1}=u^{m-1},\\ (1+\tau\partial_{x}^{3})u_{2}=\frac{1}{p}(q^{m})_{x}.\end{cases}

Hence, $u_{1}$ and $u_{2}$ can be obtained by inverting the operator $(1+\tau\partial_{x}^{3})$ , which will be easy if the Fourier spectral spatial discreitzation is considered.

Then, substituting $v^{m+1}=2\tilde{v}^{n}-v^{m-1}$ and putting it into the second equation of (5.5), using the relation $u^{m+1}=2\tilde{u}^{m}-u^{m-1}$ and $\tilde{u}^{m}=u_{1}^{m}+\tilde{v}^{m}u_{2}^{m}$ , we obtain the formula for $\tilde{v}^{m}$ :

(5.8)

\displaystyle\tilde{v}^{m}=\left({1-\frac{p+1}{2}(q^{m},u_{2}^{m})}\right)^{-1}\left(\frac{p+1}{2}(q^{m},u_{1}^{m}-u^{m-1})+v^{m-1}\right).

The scheme (5.5) is a semi-implicit multistep method. The first step $u^{1}$ can be obtained from the MCN method (5.4) above. Moreover, we can have the modified energy conserved and the error of the mass is approximated by the third order in time ( $\mathcal{O}(\tau^{3})$ ) in the following proposition.

Proposition 5.2.

The scheme (5.5) preserves the momentum $I^{m+1}=I^{m}$ and the modified energy $E^{m+1}=E^{m}$ ; and the local error of the mass is at the third order in time, i.e.,

(5.9)

\displaystyle M^{m+1}=M^{m}+\mathcal{O}(\tau^{3}).

Proof.

Taking the inner product of the first equation (5.5) with $1$ will yield the conservation of the momentum $I^{m+1}=I^{m-1}$ , which implies $I^{m+1}=I^{m}$ .

For the conservation of the modified energy, we take the inner product of the first equation (5.5) with the term $\tilde{u}^{m}_{xx}+\frac{1}{p}q^{m}\tilde{v}^{m}$ . Then, substituting $\tilde{u}^{m}=\frac{1}{2}(u^{m+1}+u^{m-1})$ yields

\frac{1}{2}(u^{m+1}-u^{m-1},\tilde{u}^{m}_{xx})+\frac{1}{p}(u^{m+1}-u^{m-1},q^{m})\cdot\tilde{v}^{m}=0.

Note that $(u^{m+1}-u^{m-1},q^{m})=\frac{2}{p+1}(v^{m+1}-v^{m-1})$ from the second equation of (5.5), and thus, we obtain the energy conservation relation $E^{m+1}=E^{m-1}$ .

For the local error of mass, we first notice that $\tilde{u}^{m}=u^{m}+\mathcal{O}(\tau^{2})$ , and consequently $\tilde{q}^{m}:=q(\tilde{u}^{m})=q^{m}+\mathcal{O}(\tau^{2})$ . Taking the inner product of the first equation (5.5) with $\tilde{u}^{m}$ and substitute $q^{m}=\tilde{q}^{m}+\mathcal{O}(\tau^{2})$ , this yields

\frac{1}{2\tau}(u^{m+1}-u^{m-1},\tilde{u}^{m})+((\tilde{q}^{m})_{x},\tilde{u}^{m})\cdot\tilde{v}^{m}+\mathcal{O}(\tau^{2})=0.

Note that the term $((\tilde{q}^{m})_{x},\tilde{u}^{m})=0$ . Thus, we obtain (5.9) and complete the proof. ∎

One can see that unlike the SAV-IRK methods in Theorem 3.1, the semi-implicit SAV-Leap Frog (SAV-LF) method will not preserve the mass. By following the same argument, similar results can be obtained by for the conventional SAV-CNAB method from [46]. Later in our numerical examples, we will show that the failure of mass conservation will lead to the failure of simulating the breathers.

5.2.3. Strang splitting method

Another widely used method for the gKdV equations is the operator splitting method, e.g., see [25] and [24]. We briefly describe the 2nd order Strang splitting (SS) method here. We consider the ODE (or PDE) of the form

(5.10)

\displaystyle u_{t}=\mathbf{L}u+\mathbf{N}(u),

where $\mathbf{L}$ is the linear part ( $\partial_{x}^{3}$ for gKdV), and $\mathbf{N}(u)$ is the nonlinear term ( $\frac{1}{p}(u^{p})_{x}$ for gKdV). The spirit of operator splitting method is to split the gKdV equation (1.1) into two equations. One is the linear dispersive equation

(5.11)

\displaystyle u_{t}+u_{xxx}=0.

We denote the time evolution from $t=t_{m}$ to $t_{m+1}$ for this linear dispersive equation as $u^{m+1}_{\tau}=A_{\tau}u^{m}$ .

The other is the conservation law equation

(5.12)

\displaystyle u_{t}+(\frac{1}{p}u^{p})_{x}=0.

We denote the time evolution as $u^{m+1}_{\tau}=S_{\tau}u^{m}$ .

We use the strategy in [25]. The time evolution $u^{m+1}_{\tau}=A_{\tau}u^{m}$ is approximated by the standard Crank-Nicolson scheme:

\displaystyle(u^{m+1}-u^{m})+\frac{\tau}{2}\left(u^{m+1}_{xxx}+u^{m}_{xxx}\right)=0,

and the time evolution for $u^{m+1}_{\tau}=S_{\tau}u^{m}$ is approximated by the IRK2 (mid-point) method:

\displaystyle(u^{m+1}-u^{m})+\frac{\tau}{p}\left(\left(\frac{u^{m+1}+u^{m}}{2}\right)^{p}\right)_{x}=0.

Again, in the same spirit of (5.4), the fixed point iteration can be used to solve the resulting nonlinear system after the spatial discretization. Then, the solution can be updated by $u^{m+1}=S_{\frac{\tau}{2}}A_{\tau}S_{\frac{\tau}{2}}u^{m}.$

The main advantage of the operator splitting method is that we can treat the linear and nonlinear part via different kinds of approaches, and thus, leads to great flexibility. We refer the interest readers to [25] and [24] for the detailed study of operator splitting methods for the gKdV equations.

5.2.4. Exponential Time Differencing method

The other widely used efficient numerical method for the stiff ODE’s are the Exponential Time Differencing methods. Recall the Duhammel’s form of the equation (5.10):

(5.13)

\displaystyle u^{m+1}=e^{\mathbf{L}\tau}u^{m}+\int_{0}^{\tau}e^{\mathbf{L}(\tau-s)}\mathbf{N}(u(t_{m}+s))ds.

The linear stiff part $\mathbf{L}u$ could be solved exactly as $e^{\mathbf{L}\tau}u^{m}$ from proper spatial discretizations (e.g., Fourier pseudo-spectral method), and the numerical quadrature can be constructed to approximate the nonlinear part in (5.13). One possible fourth order scheme is proposed by Cox and Matthews (known as modified ETDRK4 or mETDRK4), which gives the following formulas:

	$\displaystyle a_{m}=e^{\mathbf{L}\tau/2}u^{m}+\mathbf{L}^{-1}(e^{\mathbf{L}\tau/2}-\mathbf{I})\mathbf{N}(u^{m});$
	$\displaystyle b_{m}=e^{\mathbf{L}\tau/2}u^{m}+\mathbf{L}^{-1}(e^{\mathbf{L}\tau/2}-\mathbf{I})\mathbf{N}(a_{m});$
	$\displaystyle c_{m}=e^{\mathbf{L}\tau/2}a_{m}+\mathbf{L}^{-1}(e^{\mathbf{L}\tau/2}-\mathbf{I})(2\mathbf{N}(b_{m})-\mathbf{N}(u^{m}));$
	$\displaystyle\mathbf{g}_{1}=\tau^{-2}\mathbf{L}^{-3}[-4-\tau\mathbf{L}+e^{\tau\mathbf{L}}(4-3\tau\mathbf{L}+(\tau\mathbf{L})^{2})];$
	$\displaystyle\mathbf{g}_{2}=2\tau^{-2}\mathbf{L}^{-3}[2+\tau\mathbf{L}+e^{\tau\mathbf{L}}(-2+\tau\mathbf{L})];$
	$\displaystyle\mathbf{g}_{3}=\tau^{-2}\mathbf{L}^{-3}[-4-3\tau\mathbf{L}-(\tau\mathbf{L})^{2}+e^{\tau\mathbf{L}}(4-\tau\mathbf{L})];$
(5.14)		$\displaystyle u^{m+1}=e^{\tau\mathbf{L}}u^{m}+\mathbf{g}_{1}\mathbf{N}(u^{m})+\mathbf{g}_{2}[\mathbf{N}(a_{m})+\mathbf{N}(b_{m})]+\mathbf{g}_{3}\mathbf{N}(c_{m}),$

where $\mathbf{I}$ is the identity operator. The potential singularity of the term $\mathbf{L}^{-1}(e^{\mathbf{L}\Delta t/2}-\mathbf{I})$ at the $0$ th Fourier node can be resolved by using the contour integrals from [29].

The mETDRK scheme (5.2.4) is an explicit scheme with 4th order accuracy, and thus, very efficient. For the gKdV case in this paper, we notice that this scheme (5.2.4) still has the CFL condition restriction $\tau\leq Const\cdot L/N$ , since the first order derivative is shown in the nonlinear term.

5.3. Numerical examples

In this subsection, we list our numerical examples by using the proposed schemes. We denote the $s$ -stage implicit Gaussian-Legendre Runge-Kutta methods in Table 3.1 as SAV-IRK2, SAV-IRK4 and SAV-IRK6, for $s=1,2,3$ , respectively.

We first compare the proposed SAV-IRK4 method with the existing (MCN, SAV-LF, SS and mETDRK4) methods by simulating the breathers, which is non-disperse oscillatory pulse, for the modified KdV (mKdV) equations. We find that only the conservative methods (SAV-IRK4 and MCN) are capable of simulating these type of solutions for long time, and the SAV-IRK4 method is the most efficient one.

Then, we switch to the widely considered soliton interaction solutions and scattering solutions in the next two examples. In these cases, while the invariant quantities $I_{h}^{m}$ , $M_{h}^{m}$ and ${E}_{h}^{m}$ could be preserved well as shown in Theorem 3.1, the 4th order mETDRK4 method is still the most efficient one, since it is an explicit method and no iteration is needed from solving the nonlinear system.

We track the the $L^{\infty}$ error $\mathcal{E}^{m}$ , error of discrete momentum $\mathcal{E}_{I}^{m}$ , discrete mass $\mathcal{E}^{m}_{M}$ and discrete (modified) energy $\mathcal{E}^{m}_{E}$ at $t=t_{m}$ as follows:

	$\displaystyle\mathcal{E}^{m}=\\|u_{exact}^{m}-\mathbf{u}^{m}\\|_{\infty};\quad\mathcal{E}_{I}^{m}=\|I_{h}^{m}-I_{h}^{0}\|;$
	$\displaystyle\mathcal{E}^{m}_{M}=\max_{j<m}\|M_{h}^{j}-M_{h}^{0}\|;\quad\mathcal{E}^{m}_{E}=\max_{j<m}\|E_{h}^{j}-E_{h}^{0}\|.$

Example $1$

We first consider the breathers for the modified KdV equation ( $p=3$ ), which describe the non-dispersing oscillatory pulses. The two parameter family of exact solutions can be written as

	$\displaystyle B_{\alpha,\beta}(x,t)=$	$\displaystyle 2\sqrt{6}\beta\operatorname{sech}(\beta(x+\gamma t))\times$
(5.15)			$\displaystyle\left[\frac{\cos(\alpha(x+\delta t))-(\beta/\alpha)\sin(\alpha(x+\delta t))\tanh(\beta(x+\gamma t)}{1+(\beta/\alpha)^{2}\sin^{2}(\alpha(x+\delta t))\operatorname{sech}^{2}(\beta(x+\gamma t))}\right],$

with $\gamma=3\alpha^{2}-\beta^{2}$ , and $\delta=\alpha^{2}-3\beta^{2}$ . The phase velocity of the pulse is given by $\delta$ , and the group velocity by $\gamma$ traveling to the left. The authors in [1] showed that breathers have the orbital stability. In [20], the numerical experiments show that the numerical instability may easily occur as time evolves and the numerical error accumulates, which eventually leads to the failure of simulation. Thus, extra care is needed in simulating this type of solutions, especially for the long time simulation.

In our numerical experiment, we take $u_{0}=B_{\alpha,\beta}(x,0)$ , on the periodic domain $[-10\pi,10\pi]$ with $N=2^{10}$ nodes. We take $\alpha=3$ and $\beta=1$ . The time step $\tau$ is taking to be as large as possible such that the breathers will still behave good at the stopping time $T=1000$ , but will collapse if we take twice of it. This can be an indication of the stability for different numerical schemes. Surprisingly, the SAV-LF scheme does not work for simulating the breathers. In our numerical experiment, when we take $\tau=1e-4$ , the solution from the SAV-LF scheme collapse around time $t=0.55$ , and the solution collapse around $t=0.76$ if we take $\tau=1e-5$ . Hence, the increasing mass error (see Proposition 5.2) from this scheme will lead to the instability to the breathers.

The top left subplot in Figure 1 illustrates the solution profiles of the breather at different time on the periodic domain. The rest of the subplots in Figure 1 show the conservation of the discrete momentum, mass and energy (modified energy for SAV-IRK4). These results justify the conservation properties of the schemes proposed in Theorem 3.1 and Proportion 5.1. Moreover, the SS and mETDRK4 schemes also preserve the discrete momentum from our simulation.

Refer to caption — Figure 1. Top left: Breather profiles. The solution travels periodically on the domain $[-L,L]$ . Top right: error of discrete momentum for from different time integrators. Bottom left: error of discrete mass. Bottom right: error of discrete energy (modified energy for SAV-IRK4).

From [1], we have the relation $M[B]=2\beta M[Q]$ , and $E[B]=2\beta\gamma|E[Q]|$ , where $Q=\sqrt{6}\operatorname{sech}(x)$ being the soliton solution from the equation

-Q_{xx}+Q-\frac{1}{3}Q^{3}=0.

The consistency can be checked by tracking the numerical quantities of $\beta$ and $\gamma$ with respect to time. That is, tracking

{\beta}_{num}^{m}=\frac{E^{m}_{h}}{2M[Q]}\quad\mbox{and}\quad{\gamma}_{num}^{m}=\frac{E^{m}_{h}}{2{\beta}_{num}^{m}|E[Q]|},

at different time $t=t_{m}$ and comparing with the inital given $\beta$ and $\gamma$ . Note that the $E_{h}^{m}$ here for SAV-IRK scheme is the energy from the form (1.4), not the modified energy (2.5).

In Figure 2, there is no visual difference for the values of $\beta_{num}$ between SS and mETDRK4 methods, and so are for the conservative MCN and SAV-IRK4 methods, see the left subplot. When tracking $\gamma_{num}$ , we observe that the SS method has the most deviation from the accurate value of $\gamma=26$ . In the meanwhile, $\gamma_{num}$ from the mETDRK4 method shows a good agreement to $\gamma$ , but the deviation increases linearly (see the purple circle line in the right subplot). However, the conservative methods, MCN and SAV-IRK4, keep close to the accurate value $\gamma$ all the time.

Table 2 gives more details on our simulation. The simulation is processed via Matlab 2019b on the workstation with 18 cores Intel processor i9-10980XE. One can see that the SAV-IRK4 scheme is the most efficient one, since we can take the largest time step size $\tau=0.02$ , and the quantities $\max_{m}(|\beta-\beta_{num}^{m}|)$ and $\max_{m}(|\gamma-\gamma_{num}^{m}|)$ still remain small ( $\sim 10^{-11}$ ) at the end of simulation. While the mETDRK4 scheme (with small $\tau=0.0005$ ) takes less CPU time than the SAV-IRK4 scheme, as shown in Figure 2, the $\max_{m}(|\gamma-\gamma_{num}^{m}|)=0.061$ , which is relatively large and potential collapse may happen if extending the simulation time. Indeed, when we extend our simulation time to $T=2000$ with the same time step size in Table 2, the breather collapse from the mETDRK4 scheme (red dash), but still keeps its shape from the SAV-IRK4 method (blue solid) (see the left subplot in Figure 3). In the right subplot of Figure 3, the oscillation of $|\beta-\beta_{num}^{m}|$ and $|\gamma-\gamma_{num}^{m}|$ are within the range of $10^{-9}$ , which suggests the possible cancellation of the numerical errors from the SAV-IRK4 scheme, other than the accumulation of the numerical errors from the mETDRK4 scheme which leads to the failure of the simulation.

Method	$\tau$	$\max_{m}(\|\beta-\beta_{num}^{m}\|)$	$\max_{m}(\|\gamma-\gamma_{num}\|)$	CPU time
MCN	$2e-3$	$8e-7$	$2e-5$	$739.09$ s
SAV-LF	NA	NA	NA	NA
SS	$2e-3$	$7e-4$	$1.016$	$882.78$ s
SAV-IRK4	$2e-2$	$3e-11$	$2e-11$	$389.96$ s
mETDRK4	$5e-4$	$7e-4$	$0.061$	$287.68$ s

Table 2. Error and computational time to

T=1000

for simulation the breathers from different schemes.

Our remaining examples consider the comparison between the IRK methods directly from (1.1) and the SAV-IRK methods, showing that the conservative SAV-IRK approach will preserve the discrete momentum and energy as desired, regardless for the soliton type of solutions or the fast oscillating disperse soltuions. The results from the mETDRK4 method is also included as a comparison.

Example $2$

. We consider the interaction of two soliton solutions to the KdV ( $p=2$ ) equation. The exact solution satisfies

(5.16)

\displaystyle u(x,t)=12\frac{\gamma_{1}^{2}e^{\theta_{1}}+\gamma_{2}^{2}e^{\theta_{2}}+2(\gamma_{2}-\gamma_{1})^{2}e^{\theta_{1}+\theta_{2}}+a^{2}(\gamma_{2}^{2}e^{\theta_{1}}+\gamma_{1}^{2}e^{\theta_{2}})e^{\theta_{1}+\theta_{2}}}{(1+e^{\theta_{1}}+e^{\theta_{2}}+a^{2}e^{\theta_{1}+\theta_{2}})^{2}},

where we take the parameters

	$\displaystyle\gamma_{1}=0.4,\,\,\gamma_{2}=0.6,\,\,a^{2}=\left(\frac{\gamma_{1}-\gamma_{2}}{\gamma_{1}+\gamma_{2}}\right)^{2}=\frac{1}{25},$
	$\displaystyle\theta_{1}=\gamma_{1}x-\gamma_{1}^{3}t+x_{1},\,\,\theta_{2}=\gamma_{2}x-\gamma_{2}^{3}t+x_{2},\,\,x_{1}=10,\,\,x_{2}=25.$

This initial condition represents two solitons, a larger one is on the left of the smaller one (see the left subplot in Figure 4). Both of the solitons travel to the right. The larger soliton travels with faster speed, and thus, will catch up with the smaller one. Then, they are supposed to merge and split again. We take $L=30\pi$ with $N=2048$ in space. The time step is taken to be $\tau=0.1$ and the stopping time $T=200$ . Figure 4 shows the solution profile obtained by the SAV-IRK4 method. The left subplot is the initial condition $u_{0}$ , the right subplot is the time evolution. One can see that the solitons travel to the right with different speeds, intersect, and then split again with their initial speeds. This is similar to the result in [52].

Figure 5 tracks the results obtained by different time integrators. The top left subplot in Figure 5 shows $\|\mathbf{u}^{m}-u_{exact}\|_{\infty}$ with respect to time, where $u_{exact}=u(\mathbf{x},t_{m})$ is the exact solution. One can see that the error by using the mETDRK4 method is between the error by IRK2 type methods (IRK2 and SAV-IRK4) and IRK4 type methods (IRK4 and SAV-IRK4). This is reasonable, since IRK2 type methods are of the second order accuracy in time, whereas the rest three are of the 4th order accuracy in time. The IRK4 type methods are more accurate in simulating the soliton interactions at the same time step level. However, IRK4 type methods are implicit methods which require iterations at each time step. For example, given the same time step size, the mETDRK4 method is usually $10$ to $40$ times faster than the SAV-IRK4 scheme. Hence, for this soliton type of solutions, choosing smaller time step for the mETDRK4 method will generate more accurate results than the IRK4 methods in shorter actual computational time. The two IRK2 methods (IRK2 and SAV-IRK2), and the two IRK4 methods (IRK4 and SAV-IRK4) are indistinguishable from each other in the subplot. In fact, their $L^{\infty}$ difference are on the order of $10^{-4}$ for IRK2 methods, and $10^{-9}$ for IRK4 methods.

The rest of the subplots, from the top right to the bottom right in Figure 5, track the errors for the discrete momentum ( $I_{h}^{m}$ ), discrete mass ( $M_{h}^{m}$ ) and discrete energy ( $E_{h}^{m}$ from (1.4)) from IRK and mETDRK4 methods or the modified discrete energy ( $E_{h}^{m}$ from (2.5)) from SAV-IRK methods. The conservation of the $I_{h}^{m}$ and modified energy ${E}_{h}^{m}$ show the good agreement to the Theorem 3.1. We also observe that in this case, error of discrete mass in Theorem 4.2 is very small, even smaller than our fixed point solver tolerance ( $10^{-12}$ ), where the theoretical conservation could not be provided in Theorem 3.1.

We also list the $L^{\infty}$ error $\mathcal{E}^{m}$ for the SAV-IRK2 and SAV-IRK4 methods from different time steps $\tau$ in Table 3. The rate stands for the ratio of the error $\mathcal{E}^{m}$ in the previous row divide the error $\mathcal{E}^{m}$ in the current row, and etc. One can see that the error for SAV-IRK2 method decreases at the second order speed, and the error for SAV-IRK4 method decreases at the 4th order speed, except for $\tau=\frac{1}{40}$ , where the temporal error is too close to the error of solving the nonlinear system, and thus, affects the final results. Therefore, the SAV-IRK2 and SAV-IRK4 are still of the 2nd order and the 4th order temporal accuracy, respectively. Moreover, by taking $\tau=\frac{1}{160}$ for the mETDRK4 scheme, we obtain the more accurate result ( $\mathcal{E}^{m}\approx 5e-12$ ) in shorter computational time ( $\sim 9.27$ s).

	SAV-IRK2			SAV-IRK4
$\tau$	error	rate	CPU time	error	rate	CPU time
$\frac{1}{5}$	$1.7e-3$	NA	$6.579$ s	$6.40e-8$	NA	$4.949$ s
$\frac{1}{10}$	$4.3e-4$	$3.998$	$13.966$ s	$3.89e-9$	$16.47$	$12.016$ s
$\frac{1}{20}$	$1.1e-4$	$3.998$	$22.787$ s	$2.47e-10$	$15.75$	$21.944$ s
$\frac{1}{40}$	$2.72e-5$	$4.001$	$34.330$ s	$1.78e-11$	$13.82$	$38.257$ s

Table 3. The convergence rates of SAV-IRK2 and SAV-IRK4 in Example 2. We note that by using the mETDRK4, we can obtain the error

\approx 5e-12

(more accurate) in

\sim 9.27

seconds (faster).

Example $3$

. We next consider the scattering solution for the KdV equation with the initial condition $u_{0}=-\operatorname{sech}^{2}(x)$ . This example is a soliton type data with negative sign. The KdV evolution will lead $u_{0}$ to the dispersion, and possible negative value for $\int u^{3}dx$ . Here, the $C_{0}$ adjustment process in (3.5) will make $v(t)$ stay positive, and thus, will keep the algorithm applicable for all times. The exact solution is not explicitly given, since it is not exactly the soliton due to the negative sign and coefficients chosen. We compute the reference solution $u_{ref}$ by both mETDRK4 and SAV-IRK4 methods separately with an ultimately small time step ( $\tau=1/25600$ ). The results from these two different approaches differ at the level of $10^{-12}$ . Therefore, we can take the reference solution as the “exact” solution. In the simulation, we take $\tau=0.01$ , with $L=30\pi$ and $N=2048$ , the same as in the previous example.

The left subplot in Figure 6 shows the solution profile evolving in time. One can see that it disperses to the left. The right subplot shows the solution profile at the stopping time $T=1$ . The fast oscillations on the left can be observed. Figure 7 tracks the data $\mathcal{E}^{m}$ , $\mathcal{E}_{I}^{m}$ , $\mathcal{E}_{M}^{m}$ and $\mathcal{E}_{E}^{m}$ that we are interested in. We conclude that for this kind of scattering data, the mETDRK4 method is still the most efficient since it produces the smallest $L^{\infty}$ error in the shortest computational time.

The top right subplot in Figure 7 shows that the momentum is conserved for all these five approaches. The bottom subplots track the errors of mass and energy. The results are as expected and similar to the Example 2.

Table 4 shows the convergence rates for mETDRK4, SAV-IRK2 and SAV-IRK4. From Table 4, we find that the convergence rates is lower than the expected rate for all these approaches. However, as the time step $\tau$ decreases, the rates are approaching their theoretical values ( $4$ for the 2nd order methods and $16$ for the 4th order methods). The sub-convergence rate is most likely caused by the fast oscillations of the solution as well as insufficiently small $\tau$ .

	mETDRK4		SAV-IRK2		SAV-IRK4
$\tau$	error	rate	error	rate	error	rate
$\frac{1}{100}$	$1.85e-6$	NA	$1.76e-2$	NA	$1.76e-3$	NA
$\frac{1}{200}$	$1.34e-7$	$13.81$	$6.38e-3$	$2.76$	$2.85e-4$	$6.18$
$\frac{1}{400}$	$1.12e-8$	$11.96$	$2.09e-3$	$3.04$	$3.74e-5$	$7.63$
$\frac{1}{800}$	$7.74e-10$	$14.45$	$6.01e-4$	$3.49$	$3.63e-6$	$10.28$

Table 4. The convergence rates of mETDRK4, SAV-IRK2 and SAV-IRK4 in Example 3.

6. Conclusion

We propose a numerical scheme for the generalized KdV equations, which can preserve all the three invariant quantities in the discrete time flow with an arbitrarily high order of temporal accuracy from the the combination of the symplectic Runge-Kutta method and the scalar auxiliary variable reformulation. In fact, instead of the power nonlinearity $\frac{1}{p}u^{p}$ , the general flux $f(u)$ can also be adapted to this algorithm, i.e., one can create the auxiliary variable $v=\sqrt{\mathcal{F}+C_{0}}$ such that $f=\frac{\delta\mathcal{F}}{\delta u}$ . Then, the modified energy ${E}=\frac{1}{2}(u_{x},u_{x})-v^{2}-C_{0}$ will be conserved in the discrete SRK time flow.

In our numerical experiments, besides the fulfillment of conservation laws for our proposed numerical methods, we also show that among the popular existing time integrators, the mETDRK4 method is likely to be the most efficient one for the soliton type and the scattering type of solutions. However, for the long time simulation on the oscillatory non-dispersing solutions, such as the breathers for the mKdV equations, the conservative methods are recommended. Our numerical simulations demonstrate the high accuracy and efficiency of the proposed schemes in these cases.

References

[1] M. A. Alejo and C. Muñoz. Nonlinear stability of mkdv breathers. Communications in Mathematical Physics, 324(1):233–262, 2013.
[2] D. N. Arnold and R. Winther. A superconvergent finite element method for the Korteweg-de-Vries equation. Mathematics of Computation, 38(157):23–36, 1982.
[3] T. B. Benjamin, J. L. Bona, and J. J. Mahony. Model equations for long waves in nonlinear dispersive systems. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 272(1220):47–78, 1972.
[4] J. Bona, V. Dougalis, O. Karakashian, W. McKinney, and F. Smith. Conservative, high-order numerical schemes for the generalized Korteweg-de-Vries equation. Philosophical Transactions of the Royal Society of London. Series A: Physical and Engineering Sciences, 351(1695):107–164, 1995.
[5] J. Bona, W. Pritchard, and L. Scott. An evaluation of a model equation for water waves. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 302(1471):457–510, 1981.
[6] J. L. Bona. On solitary waves and their role in the evolution of long waves. Applications of nonlinear analysis in the physical sciences, pages 183–205, 1981.
[7] J. L. Bona, T. Colin, and D. Lannes. Long wave approximations for water waves. Archive for rational mechanics and analysis, 178(3):373–410, 2005.
[8] J. L. Bona, P. E. Souganidis, and W. A. Strauss. Stability and instability of solitary waves of Korteweg-de Vries type. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 411(1841):395–412, 1987.
[9] J. C. Butcher. Implicit Runge-Kutta processes. Mathematics of Computation, 18(85):50–64, 1964.
[10] J. Cai and J. Shen. Two classes of linearly implicit local energy-preserving approach for general multi-symplectic Hamiltonian PDEs. Journal of Computational Physics, 401:108975, 2020.
[11] Y. Cheng and C.-W. Shu. A discontinuous Galerkin finite element method for time dependent partial differential equations with higher order derivatives. Mathematics of Computation, 77(262):699–730, 2008.
[12] G. Cooper. Stability of Runge-Kutta methods for trajectory problems. IMA journal of numerical analysis, 7(1):1–13, 1987.
[13] S. M. Cox and P. C. Matthews. Exponential time differencing for stiff systems. J. Comput. Phys., 176(2):430–455, 2002.
[14] W. Craig. An existence theory for water waves and the Boussinesq and Korteweg-de Vries scaling limits. Communications in Partial Differential Equations, 10(8):787–1003, 1985.
[15] J. Cui, Y. Wang, and C. Jiang. Arbitrarily high-order structure-preserving schemes for the Gross–Pitaevskii equation with angular momentum rotation. Computer Physics Communications, 261:107767, 2021.
[16] Y. Cui and D.-K. Mao. Numerical method satisfying the first two conservation laws for the Korteweg-de Vries equation. Journal of Computational Physics, 227(1):376–399, 2007.
[17] B. Fornberg and G. B. Whitham. A numerical and theoretical study of certain nonlinear wave phenomena. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 289(1361):373–404, 1978.
[18] K. Goda. Numerical studies on recurrence of the Korteweg-de-Vries equation. Journal of the Physical Society of Japan, 42(3):1040–1046, 1977.
[19] Y. Gong, J. Cai, and Y. Wang. Multi-symplectic fourier pseudospectral method for the Kawahara equation. Communications in Computational Physics, 16(1):35–55, 2014.
[20] C. Gorria, M. A. Alejo, and L. Vega. Discrete conservation laws and the convergence of long time simulations of the mkdv equation. Journal of Computational Physics, 235:274–285, 2013.
[21] T. Grava and C. Klein. A numerical study of the small dispersion limit of the Korteweg-de Vries equation and asymptotic solutions. Physica D: Nonlinear Phenomena, 241(23-24):2246–2264, 2012.
[22] A. Guillou and J. L. Soulé. La résolution numérique des problèmes différentiels aux conditions initiales par des méthodes de collocation. ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, 3(R3):17–44, 1969.
[23] B.-Y. Guo and J. Shen. On spectral approximations using modified Legendre rational functions: application to the Korteweg-de-Vries equation on the half line. Indiana University Mathematics Journal, pages 181–204, 2001.
[24] H. Holden, K. Karlsen, N. Risebro, and T. Tao. Operator splitting for the KdV equation. Mathematics of Computation, 80(274):821–846, 2011.
[25] H. Holden, K. H. Karlsen, and N. H. Risebro. Operator splitting methods for generalized Korteweg-de-Vries equation. Journal of Computational Physics, 153(1):203–222, 1999.
[26] W. Huang and D. M. Sloan. The pseudospectral method for third-order differential equations. SIAM Journal on Numerical Analysis, 29(6):1626–1647, 1992.
[27] A. Jeffrey and T. Kakutani. Weak nonlinear dispersive waves: A discussion centered around the Korteweg–de Vries equation. SIAM Review, 14(4):582–643, 1972.
[28] O. Karakashian and W. McKinney. On optimal high-order in time approximations for the Korteweg-de-Vries equation. Mathematics of Computation, pages 473–496, 1990.
[29] A.-K. Kassam and L. N. Trefethen. Fourth-order time-stepping for stiff PDEs. SIAM J. Sci. Comput., 26(4):1214–1233, 2005.
[30] T. Kato. On the cauchy problem for the (generalized) Korteweg-de Vries equation. Studies in Appl. Math. Adv. in Math. Suppl. Stud., 8:93–128, 1983.
[31] C. E. Kenig, G. Ponce, and L. Vega. Well-posedness and scattering results for the generalized Korteweg-de Vries equation via the contraction principle. Communications on Pure and Applied Mathematics, 46(4):527–620, 1993.
[32] C. Klein and R. Peter. Numerical study of blow-up and dispersive shocks in solutions to generalized Korteweg-de Vries equations. Phys. D, 304/305:52–78, 2015.
[33] P. D. Lax. Integrals of nonlinear equations of evolution and solitary waves. Communications on Pure and Applied Mathematics, 21(5):467–490, 1968.
[34] J. Li and M. R. Visbal. High-order compact schemes for nonlinear dispersive waves. Journal of Scientific Computing, 26(1):1–23, 2006.
[35] H. Liu and J. Yan. A local discontinuous Galerkin method for the Korteweg-de Vries equation with boundary effect. Journal of Computational Physics, 215(1):197–218, 2006.
[36] H. Liu and N. Yi. A Hamiltonian preserving discontinuous Galerkin method for the generalized Korteweg–de Vries equation. Journal of Computational Physics, 321:776–796, 2016.
[37] Z. Liu, H. Zhang, X. Qian, and S. Song. Mass and energy conservative high order diagonally implicit Runge-Kutta schemes for nonlinear Schrödinger equation in one and two dimensions. arXiv preprint arXiv:1910.13700, 2019.
[38] H. Ma and W. Sun. A Legendre–Petrov–Galerkin and Chebyshev collocation method for third-order differential equations. SIAM Journal on Numerical Analysis, 38(5):1425–1438, 2000.
[39] Y. Martel and F. Merle. Instability of solitons for the critical generalized Korteweg—de Vries equation. Geometric & Functional Analysis GAFA, 11(1):74–123, 2001.
[40] R. M. Miura, C. S. Gardner, and M. D. Kruskal. Korteweg‐de Vries equation and generalizations. II. Existence of conservation laws and constants of motion. Journal of Mathematical Physics, 9(8):1204–1209, 1968.
[41] J. Sanz-Serna. Runge-Kutta schemes for Hamiltonian systems. BIT Numerical Mathematics, 28(4):877–883, 1988.
[42] J. Sanz-Serna and L. Abia. Order conditions for canonical Runge-Kutta schemes. SIAM Journal on Numerical Analysis, 28(4):1081–1096, 1991.
[43] J. Shen. A new dual-Petrov-Galerkin method for third and higher odd-order differential equations: application to the KdV equation. SIAM Journal on Numerical Analysis, 41(5):1595–1619, 2003.
[44] J. Shen, T. Tang, and L.-L. Wang. Spectral methods, volume 41 of Springer Series in Computational Mathematics. Springer, Heidelberg, 2011. Algorithms, analysis and applications.
[45] J. Shen and J. Xu. Convergence and error analysis for the scalar auxiliary variable (SAV) schemes to gradient flows. SIAM Journal on Numerical Analysis, 56(5):2895–2912, 2018.
[46] J. Shen, J. Xu, and J. Yang. The scalar auxiliary variable (SAV) approach for gradient flows. Journal of Computational Physics, 353:407–416, 2018.
[47] J. Shen, J. Xu, and J. Yang. A new class of efficient and robust energy stable schemes for gradient flows. SIAM Review, 61(3):474–506, 2019.
[48] L. N. Trefethen. Spectral methods in MATLAB, volume 10 of Software, Environments, and Tools. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2000.
[49] Y. Xu and C.-W. Shu. Error estimates of the semi-discrete local discontinuous Galerkin method for nonlinear convection–diffusion and KdV equations. Computer methods in applied mechanics and engineering, 196(37-40):3805–3822, 2007.
[50] X. Yang. Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends. Journal of Computational Physics, 327:294–316, 2016.
[51] X. Yang, J. Zhao, and Q. Wang. Numerical approximations for the molecular beam epitaxial growth model based on the invariant energy quadratization method. Journal of Computational Physics, 333:104–127, 2017.
[52] N. Yi, Y. Huang, and H. Liu. A direct discontinuous Galerkin method for the generalized Korteweg-de Vries equation: energy conservation and boundary effect. Journal of Computational Physics, 242:351–366, 2013.

Arbitrarily high-order conservative schemes for the generalized Korteweg-de Vries equation

Abstract.

Key words and phrases:

1. Introduction

2. Model reformulation by the SAV approach

3. Temporal discretization

3.1. Symplectic Runge-Kutta method

Theorem 3.1.

Proof.

3.2. A C0C_{0} adjustment process

Remark 3.1.

4. Spatial discretization

Theorem 4.1.

Proof.

Corollary 4.1.

Proof.

Theorem 4.2.

Proof.

5. Numerical results

5.1. A fast solver

5.2. Brief description of existing time integrators

5.2.1. Modified Crank-Nicholson method

Proposition 5.1.

Proof.

5.2.2. Semi-implicit Scalar auxiliary variable method

Proposition 5.2.

Proof.

5.2.3. Strang splitting method

5.2.4. Exponential Time Differencing method

5.3. Numerical examples

6. Conclusion

References

3.2. A $C_{0}$ adjustment process