⇦Chinese Physics Letters, 2017, Vol. 34, No. 9, Article code 090202 A Multi-Symplectic Compact Method for the Two-Component Camassa–Holm Equation with Singular Solutions * Xiang Li(李翔)1, Xu Qian(钱旭)1,2**, Bo-Ya Zhang(张博亚)1, Song-He Song(宋松和)1 Affiliations 1College of Science and State Key Laboratory of High Performance Computing, National University of Defense Technology, Changsha 410073 2Academy of Ocean Science and Engineering, National University of Defense Technology, Changsha 410073 Received 14 June 2017 *Supported by the National Natural Science Foundation of China under Grant Nos 11571366 and 11501570, the Open Foundation of State Key Laboratory of High Performance Computing of China, the Research Fund of the National University of Defense Technology under Grant No JC15-02-02, and the Fund from HPCL.
**Corresponding author. Email: qianxu@nudt.edu.cn Citation Text: Li X, Qian X, Zhang B Y and Song S H 2017 Chin. Phys. Lett. 34 090202 Abstract The two-component Camassa–Holm equation includes many intriguing phenomena. We propose a multi-symplectic compact method to solve the two-component Camassa–Holm equation. Based on its multi-symplectic formulation, the proposed method is derived by the sixth-order compact finite difference method in spatial discretization and the symplectic implicit midpoint scheme in temporal discretization. Numerical experiments finely describe the velocity and density variables in the two-component integrable system and distinctly display the evolvement of the singular solutions. Moreover, the proposed method shows good conservative properties during long-time numerical simulation. DOI:10.1088/0256-307X/34/9/090202 PACS:02.60.Cb, 02.60.Lj, 02.70.Bf © 2017 Chinese Physics Society Article Text In this study, we aim to focus on the following two-component Camassa–Holm (CH2) equation[1,2] $$\begin{align} v_t-v_{txx}+\mu\rho\rho_x=\,&2v_xv_{xx}-3vv_x+vv_{xxx},\\ \rho_t+(v\rho)_x=\,&0,~~ \tag {1} \end{align} $$ which is extended from the classical Camassa–Holm (CH) equation[3,4] by admitting an integrable multi-component generalization of density. The CH2 equation approximates the governing equation for the velocity and fluid density dynamics of unidirectional waves on shallow water. Naturally, we use $v(x,t)$ to represent scalar velocity at time $t\geq0$ in the spatial direction $x$ and $\rho(x,t)$ to represent density (or the free surface elevation from equilibrium), where $v(x,t),\rho(x,t)\in [0,L]\times \mathbb{R}^+$. The parameter $\mu$ is a real coefficient related to the critical wave speed, and the case $\mu>0$ ($\mu < 0$) corresponds to the situation in which the gravity acceleration points downwards (upwards). We take the periodic boundary conditions $v(0,t)=v(L,t)$ and $\rho(0,t)=\rho(L,t)$ on the given domain. The CH2 equation is a generalization of the well-known Korteweg–de Vries (KdV) equation (appearing at linear order in an asymptotic expansion for unidirectional shallow water waves) and the CH equation (if taking $\rho_0=0$ as initial condition for Eq. (1)), which is more complicated to deal with. Two kinds of noteworthy qualities of CH2,[4,5] which are spontaneous emergence of peakon (singular) solutions in velocity and corner-like solutions in fluid density,[6] have always grasped the attention of numerous scholars. To probe into the steepening mechanism, Holm et al. made a slight modification of the CH2 equation to ensure that both the velocity and the fluid density are the peakon solutions.[7] The same as the KdV and the CH equations, the CH2 system possesses the complete integrability of bi-Hamiltonian equations.[8,9] Recently, Cohen et al. firstly presented the multi-symplectic formulation for the CH2 equation.[10] To avoid non-physical oscillation,[10] we present a new combined multi-symplectic compact method to explore the singularity of the CH2 system. This compact finite difference method (CFDM) approximates to the spatial derivatives with fewer space nodes for higher theoretical accuracy.[11,12] Introducing some canonical momenta, Eq. (1) can be expressed as a system of equations containing only first-order derivatives, which can be further written as the following multi-symplectic formulation $$ {\boldsymbol M}{\boldsymbol z}_t+{\boldsymbol K}{\boldsymbol z}_x=\nabla_{\boldsymbol z}S({\boldsymbol z}),~~ \tag {2} $$ where ${\boldsymbol z}= v{\boldsymbol e}_1+\phi{\boldsymbol e}_2+q{\boldsymbol e}_3+u{\boldsymbol e}_4+\theta{\boldsymbol e}_5+\rho{\boldsymbol e}_6+\gamma{\boldsymbol e}_7+\beta{\boldsymbol e}_8$ (${\boldsymbol e}_i$ are unitary basis vectors), and the skew-symmetric matrices are $$\begin{align} {\boldsymbol M}=\,&\frac12{\boldsymbol e}_1\otimes{\boldsymbol e}_2-\frac12{\boldsymbol e}_2\otimes{\boldsymbol e}_1 -\frac12{\boldsymbol e}_5\otimes{\boldsymbol e}_1+\frac12{\boldsymbol e}_1\otimes{\boldsymbol e}_5\\ &+\frac \mu2{\boldsymbol e}_6\otimes{\boldsymbol e}_8-\frac \mu2{\boldsymbol e}_8\otimes{\boldsymbol e}_6,\\ {\boldsymbol K}=\,&-{\boldsymbol e}_1\otimes{\boldsymbol e}_4+{\boldsymbol e}_4\otimes{\boldsymbol e}_1 +{\boldsymbol e}_2\otimes{\boldsymbol e}_3-{\boldsymbol e}_3\otimes{\boldsymbol e}_2\\ &+\mu{\boldsymbol e}_7\otimes{\boldsymbol e}_8-\mu2{\boldsymbol e}_8\otimes{\boldsymbol e}_7. \end{align} $$ The Hamiltonian function is given by $S(z)=-q v-v^{3}/2-v\theta ^{2}/2-\mu v\rho ^{2}/2+\theta u+\mu \gamma \rho$ and $\otimes$ denotes the Kronecker product. Furthermore, Eq. (2) satisfies the multi-symplectic conservation law $$ \partial_t\omega+\partial_x \kappa=0,~~ \tag {3} $$ where $\omega=d{\boldsymbol z} \wedge {\boldsymbol M}d{\boldsymbol z}$ and $\kappa=d{\boldsymbol z} \wedge {\boldsymbol K} d{\boldsymbol z}$. The symbol $\wedge$ indicates the wedge product, and $d{\boldsymbol z}$ is the solution of the variational equation for Eq. (3). The first-order form of Eq. (3) is $$\begin{align} \begin{cases} \frac 12\phi_t-\frac 12\theta_t-u_x=-q-\frac12\theta^2-\frac \mu2\rho^2-\frac 32v^2,\\ \frac 12v_t=q_x,\\ \phi_x=v,\\ v_x=\theta,\\ \frac12 v_t=-v\theta+u,\\ \frac \mu2 \beta_t=\mu\gamma-\mu v\rho,\\ \mu\beta_x=\mu \rho,\\ \frac \mu 2 \rho_t=-\mu \gamma_x. \end{cases} \end{align} $$ Under the periodic boundary conditions, the first two important global conserved quantities[13,14] (Hamiltonian functions) associated with the CH2 equation are $$\begin{align} {\boldsymbol H}_1=\,&\frac 12\int (v^2+v_x^2+\rho^2)dx, \\ {\boldsymbol H}_2=\,&\frac 12\int (v^3+vv_x^2+\mu v\rho^2)dx.~~ \tag {4} \end{align} $$ As shown in Eq. (4), ${\boldsymbol H}_1$ represents the conservation law for momentum, and ${\boldsymbol H}_2$ is obtained by an action principle expressed in terms of a special potential.[1] Considering the sixth-order compact finite difference method (CFDM) in spatial discretization,[12] we discretize the double periodic spatial region $[-L,L]\times[0,T]$ with coordinates $x_j=(j\times h-L), j=0,1,\ldots,N$ and $t_n=n\times \tau, n=0,1,\ldots,\hat{N}$, where $h=2L/N$ and $\tau=T/\hat{N}$ are spatial step length and temporal step length, respectively. To ensure that the scheme would be steady-going, we set the stencil-state and the weighted factors ($\sigma$, $\lambda$ and $\epsilon$) to be symmetric. Therefore, the stated formula for the first-order derivative $(V_x)_j$ is $$\begin{align} &\sigma(V_x)_{j-1}+(V_x)_{j}+\sigma(V_x)_{j+1}\\ =\,&\lambda\frac {V_{j+2}-V_{j-2}}{4h}+\epsilon\frac {V_{j+1}-V_{j-1}}{2h}. \end{align} $$ As the theoretical precision depending upon the accuracy condition constrains on these unknown coefficients, we designate the sixth-order restrictive condition $$ \epsilon+2^4\lambda=10\sigma. $$ If we select $\sigma$ as the free parameter, we can obtain a batch of schemes, $\lambda=1/3(4\sigma-1),\epsilon =2/3(\sigma+2)$. Thus it further derives a sixth-order approximation with $\sigma=\frac 13$, $\lambda=\frac 19$, $\epsilon=\frac {14}9$, and the major truncation error is $\frac 4{7!}V^{(7)}h^6$. With appropriate periodic or homogeneous boundary conditions, the above discrete approximation can be written as the matrix form, namely, $$ {\boldsymbol A}\hat{V}_x= {\boldsymbol B} \hat{V} \Longrightarrow \hat{V}_x = {\boldsymbol D}_1 \hat{V},~~ \tag {5} $$ where $\hat{V}=(V_1,V_2,\ldots,V_j,\ldots,V_N)$, ${\boldsymbol D}_1={\boldsymbol A}^{-1}{\boldsymbol B}$, and $$\begin{align} {\boldsymbol A}=\frac 13 \begin{pmatrix} 3&1&0&\dots&0&1\\ 1&3&1&\dots&0&0\\ &&\ddots&\ddots&\ddots&\\ 0&0&\dots&1&3&1\\ 1&0&\dots&0&1&3\\ \end{pmatrix}, \end{align} $$ $$\begin{align} {\boldsymbol B}=\frac 1{36h} \left(\tiny\begin{matrix} 0&28&1&0&0\dots&0&-1&-28\\ -28&0&28&1&0\dots&0&0&-1\\ -1&-28&0&28&1\dots&0&0&0\\ &&\ddots&\ddots&\ddots&\ddots&\ddots&\\ 0&0&0\dots&-1&-28&0&28&1\\ 1&0&0\dots&0&-1&-28&0&28\\ 28&1&0\dots&0&0&-1&-28&0\\ \end{matrix}\right). \end{align} $$ Obviously, ${\boldsymbol A}$ is a real-symmetric matrix, and ${\boldsymbol B}$ is a real skew-symmetric matrix. For the theoretical analysis of non-dissipative and discrete conservation laws, one can see Ref. [11] and references therein. Then, applying CFDM for spatial discretization and the implicit midpoint scheme for temporal discretization in CH2, we obtain the scheme $$\begin{align} \begin{cases} \frac 12 D_t{\it \Phi}^n -\frac 12D_t{\it \Theta}^n-A_t {\boldsymbol D}_1U^n\\ =-A_tQ^n-\frac12 (A_t{\it \Theta}^n)^2-\frac \mu2 (A_tP^n)^2-\frac 32 (A_tV^n)^2,\\ \frac 12D_tV^n=A_t{\boldsymbol D}_1Q^n,\\ A_t{\boldsymbol D}_1{\it \Phi}^n=A_tV^n,\\ A_t{\boldsymbol D}_1V^n=A_t{\it \Theta}^n,\\ \frac12 D_t V^n=-(A_tV^n)\cdot(A_t{\it \Theta}^n)+A_tU^n,\\ \frac \mu2 D_t B^n=\mu A_t{\it \Upsilon}^n-\mu (A_tV^n)\cdot(A_t P^n),\\ \mu A_t{\boldsymbol D}_1 B^n=\mu A_t P^n,\\ \frac \mu 2 D_t P^n=-\mu A_t{\boldsymbol D}_1{\it \Upsilon}^n. \end{cases} \end{align} $$ They can be rewritten as the following multi-symplectic compact formulation $$ {\boldsymbol M}D_t {\boldsymbol z}^n + {\boldsymbol K}{({\boldsymbol D}_1A_t{\boldsymbol z}^n)}=\nabla_{\boldsymbol z}S(A_t{\boldsymbol z}^n),~~ \tag {6} $$ where $V^n=[v_1^n,v_2^n,\ldots,v_N^n]^{\rm T}$ and $v_j^i=v(x_j,t_i)$. For a fixed step $\tau$, we have the first-order difference operator $D_tV^n=\frac {(V^{n+1}-V^n)}\tau$ and the average operator $A_tV^n=\frac {(V^{n+1}+V^n)}2=V^{1/2}$. The vector operation is defined by $(V^n)\cdot({\it \Theta}^n)=[v_1^n\theta_1^n,v_2^n\theta_2^n,\ldots,v_N^n\theta_N^n]^{\rm T}$. Precisely speaking, the implicit midpoint scheme has been proved as one of the typical discrete methods for Hamiltonian PDEs, which leads to multi-symplectic algorithm.[15-18] Rewriting the variational equation of Eq. (6) reads $$\begin{align} &{\boldsymbol M}D_td{\boldsymbol z}_j^n + {\boldsymbol K} \sum^N_{l=1} {{({\boldsymbol D}_1)}_{j,l}A_td{\boldsymbol z}^n_l}\\ =\,&\nabla_{{\boldsymbol z}{\boldsymbol z}}S(A_t{\boldsymbol z}_j^n)A_td{\boldsymbol z}_j^n, \end{align} $$ where ${\boldsymbol M}$, ${\boldsymbol K}$ and ${\boldsymbol D}_1$ all are skew-symmetric. Taking the wedge product with $A_td{\boldsymbol z}_j^n$, we can obtain $$ \nabla_{{\boldsymbol z}{\boldsymbol z}}S(A_t{\boldsymbol z}_j^n)A_td{\boldsymbol z}_j^n \wedge A_td{\boldsymbol z}_j^n =0.~~ \tag {7} $$ It engenders the discrete multi-symplectic conservation laws $$ D_t\omega_j^n+\sum^N_{l=1}{({\boldsymbol D}_1)}_{j,l}A_t\kappa_{j,l}^n=0,~~ \tag {8} $$ and $\omega_j^n=d{\boldsymbol z}_j^n \wedge {\boldsymbol M}d{\boldsymbol z}_j^n$, $\kappa_{j,l}^n=d{\boldsymbol z}_j^n \wedge {\boldsymbol K} d{\boldsymbol z}_l^n + d{\boldsymbol z}_l^n \wedge {\boldsymbol K} d{\boldsymbol z}_j^n$, which are equivalent to Eq. (3). Eliminating the intermediate variables, we can obtain the full-discrete scheme $$\begin{align} &D_tA_t(V^n)-\frac 12 D_1A_t(D_1V^{1/2})^2\\ &-D_1^2A_t[V^{1/2}\cdot(D_1V^{1/2})]-D_1^2D_tA_tV^n\\ &+\frac \kappa 2 D_1A_t({\boldsymbol P}^{1/2})^2+\frac 32D_1A_t(V^{1/2})^2=0,\\ &D_tA_t{\boldsymbol P}^n+D_1A_t(V^{1/2}\cdot {\boldsymbol P}^{1/2})=0.~~ \tag {9} \end{align} $$ Two discrete Hamiltonian invariants can be derived as follows: $$\begin{align} {\boldsymbol H}_1=\,&\frac {h}2 \sum\limits_{n} ((A_tV^n)^2+(A_tD_1V^n)^2+\kappa(A_t {\boldsymbol P}^n)^2), \\ {\boldsymbol H}_2=\,&\frac {h}2 \sum\limits_{n} ((A_tV^n)^3+(A_tV^n)\cdot(A_tD_1V^n)^2 \\ &+\kappa(A_tV^n)\cdot(A_t \mathbb{{\boldsymbol P}}^n)^2).~~ \tag {10} \end{align} $$
cpl-34-9-090202-fig1.png
Fig. 1. CH2: exact and numerical profiles of traveling wave's velocity and density ($L=5$, $N=500$, $\tau=0.001$, $h=\frac LN$) at $t=5$ s.
cpl-34-9-090202-fig2.png
Fig. 2. CH2: propagation of the numerical peakon anti-peakon solution of velocity ($L=1$, $N=250$, $\tau=0.004$, $h=\frac LN$) at $t=2$ s.
cpl-34-9-090202-fig3.png
Fig. 3. CH2: propagation of the numerical peakon–anti-peakon solution of density ($L=1$, $N=250$, $\tau=0.004$, $h=\frac LN$) at $t=2$ s.
cpl-34-9-090202-fig4.png
Fig. 4. CH2: propagation of the numerical dam-break solution of velocity ($L=30$, $N=600$, $\tau=0.01$, $h=\frac {2L}N$) at $t=20$ s.
Next, we conduct three typical numerical examples to show the singularity of the CH2 system. (i) Periodic traveling wave: we consider an exact solution of the form $v(x-ct)$ and $\rho(x-ct)$ at finite speed $c=2$ and $\mu=1$ by solving the second-order differential equation $v^{\prime\prime}=\mu \frac 4 {(v-c)^3}+\frac 2 {(v-c)^2}+v$.[3,10] From Fig. 1, we notice that our numerical results agree very well with the exact ones. (ii) Peakon–anti-peakon solution: the first remarkable feature of the CH2 equation is the presence of the peaked solitary waves. We choose the following piecewise smooth initial conditions with $\mu=1$ $$\begin{align} v_0=\,&\begin{cases} \frac 1{2\sinh(\frac 14)}\sinh(x),& x\geq0 \cap x\leq \frac 14,\\\!\! \frac 1{\sinh(-\frac 12)}\sinh(x-\frac 12),& x>\frac 14 \cap x\leq \frac 34,\\\!\! \frac 1{2\sinh(\frac 14)}\sinh(x-1),& x>\frac 34 \cap x < 1, \end{cases}\\ \rho_0=\,&1.5. \end{align} $$ Figure 2 displays the evolution of the peakon–anti-peakon solution in fluid velocity, and there are the jump discontinuities at $x=\frac 14$ and $x=\frac 34$, respectively. In addition, it is clear that the peakon profile is unstable but keeps its singularity with time. Since $\rho_0>0$, we find that the corner-like solution in fluid density retains strictly to be positive and in regularity, and the fluid density focuses around the collision points ($x=\frac 14$ and $x=\frac 34$) of the velocity, as well as the water waves moving over an underlying shear flow at the corners in Fig. 3. Slight oscillations near the crests may be caused by the discretization error of numerical method. (iii) Dam-break solution: a further remarkable property is that the equation develops singularities of the wave-breaking solutions in finite time. We consider Eq. (1) with $\mu=1$ and the dam-break initial conditions are $$\begin{align} v_0=\,&0,\\ \rho_0=\,&1+\tanh(x+a)-\tanh(x-a),~~a=2. \end{align} $$
cpl-34-9-090202-fig5.png
Fig. 5. CH2: propagation of the numerical dam-break solution of density ($L=30$, $N=600$, $\tau=0.01$, $h=\frac {2L}N$) at $t=20$ s.
The dam-break solution arises when a body of water of uniform depth is retained behind a barrier at $x=\pm a$ (initial phases) over a flat bed. Figures 4 and 5 show the propagation of the wave breaking solutions at particular times. It is clear that this solution remains bounded while its slope becomes unbounded in finite time. Each of the single waveforms is a singular solution and replicates a characteristic of the travelling peakon waves with a peak at their crest. In this case, both variables admit the singular property and keep symmetry with more and more singular elements appearing as time progresses. An overhead view of this dam-break solution, shown in Figs. 6 and 7, represents the distinct wave breaking phenomena which are the most intriguing long-standing problems of water wave theory. According to the blow-up criterion,[9,13] we simulate out that, unlike the fluid velocity, the wave breaking of density emerges after a short time.
cpl-34-9-090202-fig6.png
Fig. 6. CH2: propagation of the numerical dam-break solution of velocity ($L=30$, $N=600$, $\tau=0.01$, $h=\frac {2L}N$) at $t=20$ s.
cpl-34-9-090202-fig7.png
Fig. 7. CH2: propagation of the numerical dam-break solution of density ($L=30$, $N=600$, $\tau=0.01$, $h=\frac {2L}N$) at $t=20$ s.
cpl-34-9-090202-fig8.png
Fig. 8. CH: Propagation of the numerical dam-break solution of velocity ($L=33$, $N=660$, $\tau=0.01$, $h=\frac {2L}N$) at $t=20$ s.
For comparison, we choose a spatially confined velocity profile for the single-component CH equation $$\begin{align} v_0=\,&1+\tanh(x+a)-\tanh(x-a),~~a=1,\\ \rho_0=\,&0. \end{align} $$ Figure 8 definitely indicates that the leftward peakons are broken and only rightward peakons emerge in variable velocity. For the single-component CH equation, taking $\rho_0=0$ as the initial condition smashes the inherent left–right symmetry of the CH2 system. As shown in Fig. 9, four groups of data (including three examples as above and a change of initial phase for example 2) are simulated out to measure the global errors of two Hamiltonian invariants. Apparently, all of the data are controlled in a stable range with high accuracy and we further verify the authenticity and reliability of the singular solutions and our scheme.
cpl-34-9-090202-fig9.png
Fig. 9. The global errors of two discrete Hamiltonian invariants.
In conclusion, we have shown an efficient and feasible numerical method to depict the singular solutions of CH2 equation by virtue of multi-symplectic algorithm. By means of numerical experiment, we analyze two types of notable properties of CH2 systems, and the Hamiltonian invariants are also kept within slight error for long-time simulation. In future work, it is worth attempting more effective numerical schemes for such a model of water wave theory. References An integrable shallow water equation with peaked solitonsAdvances in Applied MechanicsOn an integrable two-component Camassa–Holm shallow water systemBreakdown of a shallow water equationInverse scattering transform for the Camassa–Holm equationA Two-component Generalization of the Camassa-Holm Equation and its SolutionsSingular solutions of a modified two-component Camassa-Holm equationOn smooth traveling waves of an integrable two-component Camassa–Holm shallow water systemWell-posedness and blow-up for a modified two-component Camassa–Holm equationA multi-symplectic numerical integrator for the two-component Camassa–Holm equationA compact split-step finite difference method for solving the nonlinear Schrödinger equations with constant and variable coefficientsHigh-order compact splitting multisymplectic method for the coupled nonlinear Schrödinger equationsMulti-symplectic integration of the Camassa–Holm equationMulti-symplectic wavelet collocation method for the nonlinear Schrödinger equation and the Camassa–Holm equationTwo New Simple Multi-Symplectic Schemes for the Nonlinear Schrödinger EquationMulti-symplectic methods for the coupled 1D nonlinear Schrödinger systemMulti-symplectic method for generalized fifth-order KdV equationMulti-symplectic Runge–Kutta methods for nonlinear Dirac equations
[1] Camassa R and Holm D D 1993 Phys. Rev. Lett. 71 1661
[2] Camassa R, Holm D D and Hyman J M 1994 Adv. Appl. Mech. 31 1
[3] Constantin A and Ivanov R I 2008 Phys. Lett. A 372 7129
[4] McKean H P 1998 Asian J. Math. 2 867
[5] Constantin A, Gerdjikov V S and Ivanov R I 2006 Inverse Probl. 22 2197
[6] Chen M, Liu S and Zhang Y 2006 Lett. Math. Phys. 75 1
[7] Holm D D, Ó L and Tronci C 2009 Phys. Rev. E 79 016601
[8] Mustafa O G 2009 Wave Motion 46 397
[9] Yu S 2012 Appl. Anal. 91 1321
[10] Cohen D, Matsuo T and Raynaud X 2014 J. Nonlinear. Math. Phys. 21 442
[11] Dehghan M and Taleei A 2010 Comput. Phys. Commun. 181 43
[12] Ma Y P, Kong L H and Hong J L 2011 Comput. Math. Appl. 61 319
[13] Cohen D, Owren B and Raynaud X 2008 J. Comput. Phys. 227 5492
[14] Zhu H J, Song S H and Tang Y F 2011 Comput. Phys. Commun. 182 616
[15] Wang Y S, Li Q H and Song Y Z 2008 Chin. Phys. Lett. 25 1538
[16] Sun J Q and Qin M Z 2003 Comput. Phys. Commun. 155 221
[17] Hu W P and Deng Z C 2008 Chin. Phys. B 17 3923
[18] Hong J L and Li C 2006 J. Comput. Phys. 211 448