Locally conservative finite difference schemes for the modified KdV equation
Gianluca Frasca-Caccia Peter E. Hydon














margin=1cm
[table]aboveskip=0pt
[table]belowskip=10pt


lemLemma
corCorollary
theoTheorem
remRemark
ℝ
𝐳
𝐱
𝐪
𝐩
𝐛
𝐜
𝐞
𝐠
𝐲
𝐮
𝐰
ψ
κ
q̃
p̃
H̃
φ
ϕ
σ
γ
η
ℋ
ℒ
𝒫
ℐ
𝒜
𝒬
𝒜
𝒬
Δx 
γ̂
([ )
d
i





/t
/t



J̃
ζ
γ̈
η̈
ω
Proof  

γ̇
η̇
sech
atan
μ_m
μ_n
sech















Locally conservative finite difference schemes for the modified KdV equation
    
G. Frasca-Caccia,   P. E. HydonCorresponding author. Email: P.E.Hydon@kent.ac.uk 

School of Mathematics, Statistics and Actuarial Science

University of Kent, Canterbury, CT2 7FS, UK
    December 30, 2023
============================================================================================================================================================================================



 

Finite difference schemes that preserve two conservation laws of a given partial differential equation can be found directly by a recently-developed symbolic approach. Until now, this has been used only for equations with quadratic nonlinearity.

In principle, a simplified version of the direct approach also works for equations with polynomial nonlinearity of higher degree. For the modified Korteweg-de Vries equation,
whose nonlinear term is cubic, this approach yields several new families of second-order accurate schemes that preserve mass and either energy or momentum. Two of these families contain Average Vector Field schemes of the type developed by Quispel and co-workers. Numerical tests show that each family includes schemes that are highly accurate compared to other mass-preserving methods that can be found in the literature.




Keywords: Finite difference methods; discrete conservation laws; modified KdV equation; energy conservation; momentum conservation.









§ INTRODUCTION


Conservation laws are among the most fundamental properties of a given system of partial differential equations (PDEs); all solutions must satisfy every conservation law. Consequently, there is a need to understand how to construct finite difference schemes that preserve multiple conservation laws of a given system.

For simplicity, we will limit our discussion to scalar PDEs for u(x,t), where x is a spatial variable and t denotes time. (The generalization to systems with two or more independent variables is straightforward, but the notation becomes messier.) A local conservation law of a given PDE,

    𝒜(x,t,[u])=0,

is a divergence expression,

    Div 𝐅=D_x{F(x,t,[u])}+D_t{G(x,t,[u])},

that is zero on the set of all solutions of the PDE. Here and henceforth, square brackets around a differentiable expression denote the expression and a finite number of its differential consequences using the total derivatives D_x and D_t. In this notation,

    Div 𝐅=0   when   [𝒜=0].

The functions F and G are, respectively, the flux and density of the conservation law. 

For suitable boundary conditions, the existence of such a conservation law implies that ∫ G x is an invariant, which means that it is constant on any solution. Some classes of PDEs have conservation laws and/or invariants built in as part of their structure. In particular, Hamiltonian PDEs are of the form

    u_t=𝒟δ/δ uℋ,    ℋ=∫ H([u]) x,

where 𝒟 is a skew-adjoint differential operator, δ/δ u denotes the variational derivative and H is a given function. For every Hamiltonian PDE, the Hamiltonian functional ℋ is an invariant.

Any Hamiltonian PDE can be embedded in a multisymplectic system; all such systems have a conservation law that expresses the invariance of a multisymplectic form <cit.>. Multisymplectic finite difference schemes, which are extensions of symplectic methods for Hamiltonian ordinary differential equations (ODEs), preserve this conservation law <cit.>. 


Sanz-Serna, de Frutos and Durán were among the first to study the benefits of using finite difference methods to preserve invariants of a Hamiltonian PDE <cit.>, using what is effectively a semi-discretization in time. An alternative approach is to discretize in space first, then apply an invariant-conserving method in time to preserve the approximated Hamiltonian functional. With this approach, invariants can be conserved by symplectic methods <cit.>, discrete line integral methods <cit.> or discrete gradient methods <cit.>. 

Among the most useful discrete gradient methods for approximating Hamiltonian PDEs is the Average Vector Field (AVF) method introduced in <cit.>, which generalizes Quispel and McLaren's energy-preserving approach for Hamiltonian ODEs <cit.>. According to <cit.>, the AVF method “is distinguished by its features of linear covariance, automatic preservation of linear symmetries, reversibility with respect to linear reversing symmetries and often by its simplicity.” Here is an outline of this method for a given PDE of the form (<ref>).

Relative to a generic lattice point 𝐧=(m,n), the uniform grid points are 

    x_i=x(m+i)=x(m)+iΔ x,   t_j=t(n+j)=t(n)+jΔ t.

In the following, u is approximated by either the semidiscretization U_i≈  u(x_i,t) or the full discretization u_i,j≈ u(x_i,t_j). The vectors whose i-th components are U_i and u_i,j are denoted by 𝐔 and 𝐮_j respectively, and ∇ is the differential operator whose i-th component is ∂/∂ U_i. 
Let ℋ(𝐔)Δ x be a semidiscretization of ℋ and let 𝒟 be a discrete skew-adjoint operator (typically depending on 𝐮_0 and 𝐮_1) that approximates 𝒟.
The AVF method approximates (<ref>) by 

    D_n 𝐮_0=𝒟(∫_0^1 ∇ℋ(ξ𝐮_1+(1-ξ)𝐮_0)dξ).


In <cit.>, McLachlan and Quispel proved that if the spatial discretization preserves a semidiscrete conservation law of the Hamiltonian PDE, the application of any discrete gradient method in time <cit.> yields a fully discrete (local) conservation law of the Hamiltonian. This is a powerful result, because conservation laws contain far more information than integral invariants do.


The current paper focuses on numerical methods that preserve two conservation laws of the modified Korteweg-de Vries (mKdV) equation,

    ≡ u_t+u^2u_x+u_xxx=0.

This PDE possesses a bi-Hamiltonian structure: there are two different ways to write it in the Hamiltonian form (<ref>), namely with

    ℋ=-∫112u^4+12uu_xx  x,    𝒟=D_x,

and with

    ℋ=-∫12u^2 dx,    𝒟=D_x^3+12(u^2D_x+uu_x).



The mKdV equation has infinitely many independent conservation laws <cit.>. Those for mass, momentum and energy, respectively, are <cit.>:

    D_xF_1+D_tG_1 = D_x(13 u^3 + u_xx)+D_tu,
        D_xF_2+D_tG_2 =D_x(14u^4+uu_xx-12u_x^2)+D_t(12 u^2),
        D_xF_3+D_tG_3 = D_x(12(13u^3+u_xx)^2+12u_xu_t-12uu_xt)+D_t(112u^4+12uu_xx).

For suitable (e.g. zero) boundary conditions, integrating (<ref>)–(<ref>) over the spatial domain yields the invariants

    ∫ G_1  x=∫ u  x,  ∫ G_2  x=12∫ u^2  x,   ∫ G_3  x=∫112u^4+12uu_xx  x.


For discretizations of the mKdV equation, the following notation is helpful. The forward shift operators in space and time are defined, for all functions f on the grid, by

    S_m(f(x_i,t_j))=f(x_i+1,t_j),     S_n(f(x_i,t_j))=f(x_i,t_j+1).

The forward difference operators D_m, D_n and forward average operators μ_m,μ_n are

    D_m=1Δ x(S_m-I),     D_n=1Δ t(S_n-I),    μ_m=12(S_m+I),    μ_n=12(S_n+I),

where I is the identity operator.


The AVF method (<ref>) for (<ref>) uses 
    ℋ=-∑_j112U_j^4+12U_jD_m^2U_j-1,    𝒟=D_mμ_mS_m^-1.

This gives the following 10-point energy-conserving scheme (and its shifts):

    AVF_EC=D_mμ_m{13(μ_nu_-1,0) μ_n(u_-1,0^2)+D_m^2μ_nu_-2,0}+D_nu_0,0.


To apply the AVF method (<ref>) to (<ref>), let ℋ=-∑_j12U_j^2 and let 𝒟 be the diagonal matrix whose entry 𝒟_j in the row and column indexed by j acts on functions f(m,n,[u]) as

    𝒟_jf=    D_m^3μ_mS_m^-2f+13S_m^-1μ_m{((μ_nu_j,0)^2+(μ_nu_j,0)(μ_nu_j+1,0)+(μ_nu_j+1,0)^2)D_mf}
       +13(μ_nu_j-1,0+μ_nu_j,0+μ_nu_j+1,0)(D_mμ_nu_j-1,0)f.

This yields the 10-point momentum-conserving scheme

    AVF_MC=D_m{13(μ_mμ_nu_-1,0)μ_m((μ_nu_-1,0)^2)+D_m^2μ_mμ_nu_-2,0}+D_nu_0,0.


A new direct approach that enables multiple conservation laws to be preserved was introduced recently in <cit.> and greatly simplified in <cit.>. This approach, which does not exploit either Hamiltonian structures or integrability, can be applied (at least, in principle) to any system of PDEs whose conservation laws are polynomial in [u]; see 4 of <cit.> for a non-Hamiltonian, non-integrable example. In the next sections we show that the two AVF schemes AVF_EC and AVF_MC are particular members of two parametrized families of methods that can be found by using this approach.

The rest of this paper is organized as follows. Section <ref> describes the simplified direct approach to finding finite difference schemes that preserve two local conservation laws, for mass and either momentum or energy. In <ref>, we apply this strategy to the mKdV equation and obtain several families of conservative schemes. In <ref>, numerical tests confirm the theoretical results and demonstrate the effectiveness of the new schemes compared with multisymplectic and narrow box schemes, each of which preserves mass, but not momentum or energy.



§ A STRATEGY FOR PRESERVING CONSERVATION LAWS


Given a (not necessarily Hamiltonian) PDE,

    𝒜(x,t,[u])=0,     (x,t)∈ V⊂ℝ^2,

a conservation law is in characteristic form if it satisfies

    Div 𝐅=.

The function  is called the characteristic of the conservation law. For simplicity, we assume throughout that the domain V is contractible.

 The vector space of total divergences is the kernel of the Euler operator, 

    ℰ=∑_i,j(-D_x)^i(-D_t)^j∂/∂ u_,ij ,    where    u_,ij=D_x^iD_t^j(u).


From Remark <ref>, if 
    ℰ()=0,

there exists 𝐅 such that Div 𝐅= is a conservation law.

A finite difference approximation of (<ref>) on the grid is a partial difference equation,

    (m,n,[u_0,0]) =0,

where square brackets [ ] around an expression defined on the grid denote the expression and a finite number of its shifts. From here on, tildes denote approximations to the corresponding continuous quantities. We will abbreviate u_0,0 to u wherever this does not cause confusion.

The aim is to find finite difference schemes that satisfy a discrete analogue of each preserved conservation law,

    𝐅=D_m{F(m,n,[u])}+D_n{G(m,n,[u])}=0,  when   [=0].

This difference conservation law is in characteristic form if 

    𝐅=(m,n,[u])(m,n,[u]);

just as for the continuous case, the multiplier function  is called the characteristic (see <cit.> for more details).
The following result, due to Kuperschmidt <cit.> and generalized in <cit.>, is a discrete version of Remark <ref> and plays a pivotal role in our approach.

The kernel of the difference Euler operator

    𝖤=∑_i,jS_m^-iS_n^-j∂/∂ u_i,j,

is the vector space of all difference divergences (<ref>).


Given a PDE (<ref>) with p conservation laws in characteristic form  (<ref>) that one wishes to preserve, our strategy is to seek approximations _i and  such that 

    𝖤(_i)= 0,     i=1,… p.

From Remark <ref>, there exist 𝐅_i such that each _i=Div 𝐅_i is a discretization of the corresponding continuous conservation law. This strategy can be implemented efficiently as follows.


  * Choose a rectangular stencil that is large enough to contain second-order approximations of  and all _i.

  * On the given stencil, the most general finite difference approximations  and _1 will depend on a large number of coefficients.

  * Impose consistency conditions giving second-order accuracy of the approximations at the centre of the stencil (which need not be a grid point).

  * Reduce the number of free parameters remaining by making some key terms as compact as possible; typically, these include highest derivatives and highest-order nonlinear terms.    

  * Some of the remaining parameters are determined by solving

    𝖤(_1)=0

symbolically, whenever a solution exists.
The discrete flux F_1 and density G_1, which satisfy _1=D_m F_1+D_nG_1, can be reconstructed from the characteristics <cit.>.

  * Steps 2 onward are iterated to preserve more conservation laws. If 𝖤(_i)=0 has no solutions for some i, the corresponding conservation law cannot be preserved on the chosen stencil without violating at least one of the earlier conservation laws.

Solving (<ref>) is a crucial step in the procedure; it is made tractable by the compactness conditions and the restriction to second-order approximations. Without these constraints, the only practical approach is to seek a Groebner basis (see <cit.>) for the polynomial system that determines the coefficients. This can take weeks of symbolic computation, even for scalar equations whose conservation laws have nothing worse than quadratic nonlinearities, approximated on the smallest possible stencil. The complexity increases exponentially with the degree of the polynomial nonlinearity and the size of the stencil.

By contrast, the compactness and second-order conditions simplify the calculations to the point that schemes preserving multiple conservation laws can be found in just a few minutes; a Groebner basis may not even be needed. These simplifications were introduced in <cit.> and used to obtain several parametrized families of schemes, some of which are highly accurate, for the KdV equation and a nonlinear wave equation with quadratic nonlinearity.

Recently, Frasca-Caccia  <cit.> outlined a new one-parameter family of mass- and energy-conserving 10-point schemes for the mKdV equation, which has cubic nonlinearity. In the current paper, the simplified strategy is used to extend this result to schemes that preserve mass and either momentum or energy, using 8-point and 10-point stencils. This demonstrates that the simplified strategy is not limited to quadratic nonlinearity. 



§ CONSERVATIVE METHODS FOR THE MKDV EQUATION


Each of the conservation laws (<ref>)–(<ref>) is in characteristic form; their characteristics are

    _1=1,    _2=u,    _3=13u^3+u_xx.


For simplicity, we consider only one-step schemes for the mKdV equation (<ref>). Therefore, our stencils are as shown in Figure <ref>, with B-A≥ 4.

On such stencils, we construct schemes for the mKdV equation of the form

    =D_mF_1+D_nG_1=0,

so that the mass conservation law (<ref>) is preserved.

From this starting-point, the strategy described in <ref> is used to preserve either (<ref>) or (<ref>). Linear terms in F_1, G_1 and _2 or _3 are approximated by linear combinations of the values of u at the stencil points, with undetermined coefficients:

    ∂^r+s/∂ x^r∂ t^su ≈1/Δ x^r1/Δ t^s∑_i=A^B∑_j=0^1α_i,ju_i,j.

Similarly, the cubic terms in _3 and F_1 are approximated by cubics:

    u^3≈∑_i=A^B   (∑_k=i^B∑_j=k^B(β_i,j,ku_i,0u_j,0u_k,0+γ_i,j,ku_i,1u_j,1u_k,1).
           .+∑_k=A^B∑_j=k^B(δ_i,j,ku_i,1u_j,0u_k,0+ε_i,j,ku_i,0u_j,1u_k,1)).

The undetermined coefficients satisfy consistency conditions to ensure that the schemes are second-order accurate at the centre of the stencil.


An approximation of a conservation law in the form (<ref>) is second-order accurate at the centre (x, t) if the approximations of F and G are second-order accurate at the points (x-Δ x/2, t) and (x, t-Δ t/2), respectively.



We now use the simplified strategy to obtain families of schemes, all of which depend on parameters that are 𝒪(Δ x^2,Δ t^2). It is possible to find values of the parameters that partially remove the leading terms of the local truncation error. However, there is no choice of the parameters that eliminates the second-order error terms identically; the optimal parameter values depend on the particular problem being approximated. No scheme in these families preserves all three conservation laws, though some come fairly close to doing so.
 


 §.§ 8-point schemes

The most compact stencil for the mKdV equation has eight points. We choose A=-2, B=1 in Figure <ref>, and seek second-order approximations of characteristics, densities and fluxes at (-1/2,1/2), (-1/2,0) and (-1,1/2), respectively.


  §.§.§ Energy-conserving methods.

To simplify the symbolic calculation giving energy-conserving schemes on the 8-point stencil, we use the following approximations to the second derivatives in F_1 and Q_3, respectively:

    u_xx≈ D_m^2μ_nu_-2,0+λ D_mD_nμ_m u_-2,0 ,   u_xx≈ D_m^2μ_mμ_nu_-2,0+λ D_mD_nμ_m^2u_-2,0,
 where λ=𝒪(Δ x^2,Δ t^2). All other terms in  and Q_3 are of the form (<ref>) or (<ref>), subject only to second-order consistency and

    𝐄(_3)= 0.

In this way, we find a family of schemes that depends on two parameters. One of these parameters is 𝒪(Δ x^4,Δ t^4); as this is negligible on fine grids, we set this parameter to zero, obtaining a family of schemes (λ) with

    F_1=13(μ_nμ_m^2u_-2,0)μ_n((μ_m^2u_-2,0)^2)+D_m^2μ_n u_-2,0+λ D_nD_mμ_mu_-2,0, G_1=μ_mu_-1,0.

Each of these schemes preserves the following discrete version of the conservation law (<ref>):

    _3= D_mF_3+ D_nG_3,

with _3=μ_mF_1 and

    F_3=    12F_1^2+12{(D_mμ_mμ_n u_-2,0) D_nμ_m^2 u_-2,0-(μ_m^2μ_nu_-2,0)D_mD_nμ_mu_-2,0}
       +12λ(D_nμ_mu_-2,0)D_nμ_mu_-1,0+148Δ x^2{(μ_m^2u_-2,0)^2+(μ_m^2u_-2,1)^2+(μ_m^2u_-2,0)(μ_m^2u_-2,1)}
       ×{(D_nμ_m^2u_-2,0)D_mμ_mμ_nu_-2,0-(μ_m^2μ_nu_-2,0)D_mD_nμ_mu_-2,0},
    G_3=    112(μ_mu_-1,0)(μ_m^3u_-2,0){μ_m^2((μ_mu_-2,0)^2)+14Δ x^2(D_mμ_mu_-2,0)D_mμ_mu_-1,0}
       +12(μ_mu_-1,0)D_m^2μ_mu_-2,0 .

In both  G_3 and F_3 there are terms that do not have any counterpart in the continuous density and flux. These all vanish as the stepsizes tend to zero.

Assuming that Δ t<Δ x, the leading error terms are 𝒪(Δ x^2). We use the notation

    EC_8(λ_1)=(λ_1Δ x^2);

the parameter λ_1 may be chosen to minimise the local truncation error.


  §.§.§ Momentum-conserving methods.

The momentum-conserving schemes introduced in this section are obtained by specifying that

    u_xx=D_m^2μ_nu_-2,0,

and choosing the following compact approximations of G_1 and Q_2:

    G_1=μ_m u_-1,0,    Q_2=μ_mμ_nu_-1,0 .

Then F_1 is determined from the general forms (<ref>) and (<ref>) by requiring consistency and 

    𝐄(_2)≡ 0.

This gives a three-parameter family of schemes that preserve mass and momentum.[Without the constraint on G_1 in (<ref>), there is another parameter. Removing (instead) the constraint on Q_2 yields two other families of schemes.] Two of the parameters are 𝒪(Δ x^4,Δ t^4), so we set these to zero and obtain the one-parameter family 𝒜(λ) with G_1 given in (<ref>) and

    F_1=    16(μ_n(u_-2,0+u_0,0))(μ_n u_-1,0)^2+D_m^2μ_nu_-2,0+λ{2(μ_nu_-1,0)μ_m((D_mμ_nu_-2,0)^2).
       .+2(D_m^2μ_nu_-2,0)μ_m^2((μ_nu_-2,0)^2)-Δ xΔ t(D_nD_mμ_mu_-2,0)μ_n((D_mμ_mu_-2,0)^2)}

where λ=𝒪(Δ x^2,Δ t^2).
For any value of λ, the discrete momentum conservation law 

    _2= D_mF_2+ D_nG_2=0

is preserved, with 𝒬_2 given in (<ref>) and

    F_2=    112(μ_nu_-1,0)^2{2(μ_nu_-2,0)(μ_mμ_nu_-1,0)+(μ_nu_-1,0)(μ_nu_0,0)}+(μ_nu_-1,0)D_m^2μ_nu_-2,0
       -12(D_mμ_nu_-2,0)D_mμ_nu_-1,0+λ(μ_mμ_nu_-2,0)(μ_mμ_nu_-1,0){(D_mμ_nu_-2,0)D_mμ_nu_-1,0.
       .+(μ_n(u_-2,0+u_0,0))D_m^2μ_nu_-2,0}+14λΔ xΔ t{(D_mμ_mμ_nu_-2,0)D_nμ_m^2u_-2,0.
       .-(μ_m^2μ_nu_-2,0)D_mD_nμ_mu_-2,0}{2μ_n((D_mμ_mu_-2,0)^2)-(D_mμ_mu_-2,0)D_mμ_mu_-2,1},
    G_2=    12(μ_mu_-1,0)^2+λΔ tΔ x(μ_mu_-1,0)(D_m^2μ_mu_-2,0){14(D_mμ_mu_-1,0)D_mμ_mu_-2,0-(D_mμ_m^2u_-2,0)^2}.

All terms in F_1, F_2 and G_2 that do not approximate any quantity in the continuous expressions vanish as the stepsizes tend to zero.
They are identically zero if λ=0. 

If Δ t <Δ x, the leading term in the local truncation error is 𝒪(Δ x^2); by choosing λ = λ_2 Δ x^2 optimally, one may be able to remove at least part of this error. In the numerical tests section, we use the notation

    _8(λ_2)=(λ_2Δx^2).




 §.§ 10-point schemes

We now develop schemes that preserve two conservation laws of the mKdV equation on the 10-point stencil with A=-2, B=2 in Figure <ref>. Approximations of characteristics, densities and fluxes of the conservation laws are second-order accurate at (0,1/2), (0,0) and (-1/2,1/2), respectively.



  §.§.§ Energy-conserving methods.

To develop mass and energy conserving schemes on the 10-point stencil, we approximate G_1 and the cubic term of 𝒬_3 on the most compact possible sub-stencils. This gives the one-parameter family of schemes found in <cit.>, namely (λ) with F_1=μ_mφ_-1,0, G_1=u_0,0, where

    φ_-1,0   =13(μ_n u_-1,0^2)(μ_n u_-1,0)+D_m^2μ_n u_-2,0+λ D_mD_nμ_mu_-2,0,

where λ=𝒪(Δ x^2,Δ t^2). The discrete local energy conservation law satisfied by each scheme in this family is

    _3=    D_mF_3+ D_nG_3=0,

where _3=φ_0,0 and

    F_3   =12{φ_-1,0φ_0,0+(D_mμ_nu_-1,0)D_nμ_mu_-1,0-(μ_mμ_nu_-1,0)D_mD_nu_-1,0+λ(D_nu_0,0)D_nu_-1,0},
    G_3   = 112u_0,0^4+12u_0,0D_m^2u_-1,0 .

This family of schemes is written as 
    EC_10(λ_3)=(λ_3Δ x^2).

Note that EC_10(0) is the AVF scheme AVF_EC that was derived in <ref>.



  §.§.§ Momentum-conserving methods.

The complexity of the symbolic computation for solving

    𝐄(_2)≡ 0,

is reduced by using the most compact second-order approximations of G_1 and _2. This gives a family depending on 27 parameters. For simplicity, we discuss here only members of this family having a compact approximation of the nonlinear term in F_1, ignoring parameters that are 𝒪(Δ x^4,Δ t^4). This yields a one-parameter family (λ) with

    F_1=13(μ_mμ_nu_-1,0) μ_m((μ_n u_-1,0)^2)+D_m^2μ_nμ_mu_-2,0+λ D_mD_nu_-1,0 ,   G_1=u_0,0,

where λ=𝒪(Δ x^2,Δ t^2). The discrete local conservation law for momentum is

    _2=    D_mF_2+ D_nG_2=0,

where 𝒬_2= μ_nμ_mu_-1,0 and

    F_2=    112(μ_nu_-1,0)(μ_nu_0,0){(μ_nu_-1,0)^2+(μ_nu_0,0)^2+(μ_nu_-1,0)(μ_nu_0,0)}
       +(μ_mμ_nu_-1,0)D_m^2μ_mμ_nu_-2,0-14(D_mμ_nu_-1,0)D_mμ_n(u_-2,0+u_0,0)
       +12 λ{(μ_mμ_nu_-1,0)D_mD_nu_-1,0-(D_mμ_nu_-1,0)D_nμ_mu_-1,0},
    G_2=    12u_0,0(u_0,0+λ D_m^2u_-1,0).

Again, all quantities that do not approximate any term in the corresponding continuous conservation law are identically zero when λ=0. 
We denote this family of schemes by 
    _10(λ_4)=(λ_4Δ x ^2).

The simplest scheme, MC_10(0), is AVF_MC.



§ NUMERICAL TESTS


In this section we consider the mKdV equation (<ref>) on a domain V=[a,b]×[0,T], with periodic boundary conditions. We use some benchmark tests to show the effectiveness of the schemes developed in <ref> compared with two well-known conservative methods (see <cit.>). The narrow box scheme is obtained by applying a standard finite volume discretization, giving

    D_n(μ_mu_-1,0)+D_m(13(μ_nu_-1,0)^3+D_m^2μ_nu_-2,0)=0.

The second method is multisymplectic and amounts to

    D_n(μ_m^3u_-2,0)+D_m(13 μ_m((μ_mμ_nu_-2,0)^3)+D_m^2μ_nu_-2,0)=0.

This compact one-step scheme is a more efficient version of the popular Preissmann scheme <cit.>.
Each of these schemes is in divergence form, preserving a discrete version of the mass conservation law (<ref>).


For every test in this section, the computational time is roughly the same for all of the numerical schemes. Therefore, the main difference between schemes is the solution error at the final time t=T, evaluated as

    .u-u_exact/u_exact|_t=T.

On a grid with M points in space and N points in time, we evaluate the error in the conservation laws by measuring the error in the global invariants in (<ref>) as follows:

    Err_ℓ=    Δ xmax_j=1,…,N|∑_i=1^M (G_ℓ(x_i,t_j)-G_ℓ(x_i,t_1))|,    ℓ=1,2,3.

Whenever G_2 or G_3 is not defined for one of the schemes considered here, the corresponding error is instead evaluated as

    Err_2   =Δ xmax_j=1,…,N|∑_i=1^M  12(v_i,j^2-v_i,0^2)|,
    Err_3   =Δ xmax_j=1,…,N|∑_i=1^M ( 112(v_i,j^4-v_i,0^4)+12v_i,jD_m^2v_i-1,j-12v_i,0D_m^2v_i-1,0)|,

where v_i,j=u(a+iΔ x,  jΔ t) for schemes defined on the 10-point stencil and

    v_i,j=12(u(a+(i-1)Δ x,  jΔ t)+u(a+iΔ x,  jΔ t))
 for 8-point schemes.

For each numerical test and family of schemes, we state the parameter value that minimizes the solution error (<ref>). None of the schemes preserve three conservation laws, so we state the parameter values λ_i that optimize the error in the unpreserved invariant given by (<ref>) or (<ref>). We also include the results for the simplest schemes, obtained by setting each λ_i to zero. As EC_10(0)≡ AVF_EC and MC_10(0)≡ AVF_MC, this enables comparison with the AVF schemes.

The first benchmark problem is the interaction of two solitons, with the exact solution

    u(x,t)=2√(6) κ(√(c_1)coshξ_2+√(c_2)coshξ_1)/(κ^2-1)+κ^2cosh(ξ_1-ξ_2)+cosh(ξ_1+ξ_2) ;

here

    κ=√(c_1)+√(c_2)/√(c_1)-√(c_2) ,    ξ_1=√(c_1)(x-c_1t+d_1),    ξ_2=√(c_2)(x-c_2t+d_2).

We set

    c_1=2.5,     c_2=0.5,     d_1=12,     d_2=2.5,

and solve this problem on Ω=[-20,20]×[0,10], using step lengths Δ x=0.1 and Δ t=0.025.

The schemes _8(λ_1), _8(λ_2), _10(λ_3) and _10(λ_4) give the minimal solution error when λ_1=1, λ_2=-0.077, λ_3=0.04 and λ_4=0.19. The values λ_1=-0.05, λ_2=-0.073, λ_3=0.20 and λ_4≪ 0 minimize the error in the invariant that is not preserved by the scheme. As the modulus of λ_4 is very large, the corresponding term in MC_10(λ_4) is not merely a perturbation; it undermines the correct behaviour of the solution. So for this problem, minimizing the error in energy turns out to be a poor criterion for choosing the parameter.














Table <ref> shows the errors at the final time T=10 in the conservation laws and the solution. We also show the error in the phase shift of the fastest and the slowest soliton, evaluated as, respectively,

    .Err_ϕ_1=(x_1-x̃_1)|_t=10,     .Err_ϕ_2=(x_2-x̃_2)|_t=10;

where x_i (resp. x̃_i) denotes the location of the soliton peak in the exact and (resp. numerical) solution, using piecewise cubic interpolation of approximations at the grid points.
The quantity 

    Err_ϕ= Err_ϕ_1-Err_ϕ_2,

measures the extent to which the numerical solution underestimates the phase shift produced by the interaction of the two solitons.


Table <ref> shows that the new conservative schemes each preserve two discrete invariants (up to rounding errors). Each family includes schemes that are highly accurate — more so than the multisymplectic and narrow box schemes. The small values of Err_ϕ_1, Err_ϕ_2 and Err_ϕ indicate that the best schemes also reproduce the correct phase shifts. Attempting to minimize the error in the unpreserved conservation law does not optimize the numerical solution; nevertheless, this approach gives MC_8 and EC_10 schemes that are more accurate than the narrow box and multisymplectic schemes. This is not true for EC_8 and MC_10 schemes.

The upper plot in Figure <ref> shows the initial condition and the numerical solution at T=10 given by the most accurate of our schemes, _10(0.04). The lower plot compares various numerical solutions with the exact solution at T=10, near to the top of the faster soliton. The approximations at the grid points are connected by piecewise cubic interpolation, except for the solution of _10(0.04) (for which the interpolation would cover the exact solution). The narrow box scheme is quite accurate for this problem, but it is outclassed by the AVF scheme _10(0) and more so by _10(0.04).

We now study the same problem on a coarser grid with Δ x=0.2 and Δ t=0.05. For each family, the parameter values minimizing the error in the solution or the unpreserved invariant are as given in Table <ref>; they are close to the optimal values for the finer grid. The solution error for the most accurate scheme in each family is around 4 times greater than on the finer grid, as expected. Figure <ref> compares the exact solution with various numerical solutions near to the top of the fastest soliton, again using piecewise cubic interpolation. 

 

As a second benchmark problem, we approximate the breather solution (see <cit.>),

    u(x,t)=∂/∂ x[-2√(6)arctan(√(3)sin(2x-64t-π/2)/cosh(2√(3)x))],

on the domain Ω=[-2,2]×[0,0.4], setting Δ x=0.02 and Δ t=0.002.

Table <ref> shows the error in the conservation laws and the solution at the final time: _8(2.22) is the most accurate scheme, while _8(-0.165), _10(0.92) and _10(1.15) are the best in their respective families at minimizing the solution error. The error in the unpreserved conservation law is minimized by choosing λ_1=0.49, λ_2=-0.128, λ_3=0.78 and λ_4≪ 0. As for the two-soliton problem, minimizing the error in the non-conserved invariant is a poor criterion for choosing the parameter in MC_10, so we do not show this result.
By contrast, choosing the values of the parameters in _8 and _10 that minimize the error in the unpreserved invariant yields fairly accurate approximations.




The upper part of Figure <ref> shows the numerical solution given by _8(2.22). The lower part compares the exact solution and the numerical solutions given by EC_8(2.22), the multisymplectic scheme and EC_10(0) (which is the most accurate AVF scheme for this problem) at the final time. The graph of the solution of the narrow box scheme is very close to the solution of EC_10(0), so is not shown in the figure. These schemes are more accurate than the multisymplectic scheme, but the frequency of the breather oscillations is best caught by EC_8(2.22), whose solution almost overlaps the exact profile.



§ CONCLUSIONS



The approach introduced in <cit.>, which uses a fast symbolic computation to find finite difference schemes that preserve two conservation laws, is not restricted to quadratic nonlinearity. By considering stencils with eight and ten nodes, we have introduced four new one-parameter families of schemes that preserve two conservation laws of the mKdV equation, which has a cubic nonlinearity.

The AVF schemes introduced by Quispel and co-workers are found by this approach; typically, each is the simplest member of its family. However, the need to obtain a skew-adjoint approximation of the skew-adjoint differential operator means that AVF methods are restricted to stencils with an odd number of points in space. 

Each of the families includes schemes that are very accurate. However, none of them preserve the first three conservation laws. Nevertheless, the value of the parameter can be chosen to minimize the error in the third invariant. Our numerical tests have shown that this criterion leads to reasonably accurate solutions for some, but not all, families.



§ ACKNOWLEDGMENTS
 The authors would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme 
Geometry, Compatibility and Structure Preservation in Computational Differential Equations,
when work on this paper was undertaken. This work was supported by EPSRC grant number EP/R014604/1.

99


Anco2 S. C. Anco, M. Mohiuddin, T. Wolf. Traveling waves and conservation laws for complex mKdV-type equations. Appl. Math. Comput. 219 (2012), 679–698.

AscherMcL U. M. Ascher, R. I.  McLachlan. Multisymplectic box schemes and the Korteweg–de Vries equation. Appl. Numer. Math. 48, (2004), 255–269.

Ascher U. M. Ascher, R. I.  McLachlan. On symplectic and multisymplectic schemes for the KdV equation. J. Sci. Comput. 25 (2005), 83–104.

Aydin A. Aydin, B. Karasözen. Multisymplectic box schemes for the complex modified Korteweg–de Vries equation. J. Math. Phys. 51, (2010), 083511.

BarlettiNLS L. Barletti, L. Brugnano, G. Frasca-Caccia, F. Iavernaro. Energy-conserving methods for the nonlinear Schrödinger equation. Appl. Math. Comput. 318 (2018), 3–18.

Bridges T. J. Bridges. Multisymplectic structures and wave propagation. Math. Proc. Cambridge Philos. Soc.
121 (1997), 147–190.

BHL T. J. Bridges, P. E. Hydon, J. K. Lawson. Multisymplectic structures and the variational bicomplex. Math. Proc. Camb. Phil. Soc. 148, (2010) 159–178.

BridgesReich2001 T. J. Bridges, S. Reich. Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity. Phys. Lett. A 284 (2001), 184–193.

BridgesReich T. J. Bridges, S. Reich. Numerical methods for Hamiltonian PDEs. J. Phys. A 39, (2006), 5287–5320.

BI L. Brugnano, F. Iavernaro. Line Integral Methods for Conservative Problems. Monograph and
Research Notes in Mathematics. Boca Raton, FL: CRC Press, 2016.

Cano B. Cano. Conserved quantities of some Hamiltonian wave equations after full discretization. Numer. Math. 103 (2006), 197–223.

Celledoni E. Celledoni, V. Grimm, R. I. McLachlan, D. I. McLaren, D. O'Neal, B. Owren, G. R. W.  Quispel. Preserving energy resp. dissipation in numerical PDEs using the “Average Vector Field” method.  J. Comput. Phys. 231 (2012), 6770–6789.

Cell2 E. Celledoni, R. I. McLachlan, B. Owren, G. R. W. Quispel. Energy-preserving integrators and the structure of B-series. Found. Comput. Math.
10 (2010), 673–693.

Cox D. Cox, J. Little, D. O’Shea. Ideals, Varieties and Algorithms. An Introduction to Computational
Algebraic Geometry and Commutative Algebra. New York, Springer, 1992.

Dahlby M. Dahlby, B. Owren. A general framework for deriving integral preserving numerical methods for
PDEs. SIAM J. Sci. Comput. 33 (2011), 2318–2340.

DuranLM A. Durán, M. A. López-Marcos. Conservative numerical methods for solitary wave interactions. J. Phys. A 36 (2003), 7761–7770.

DuranSS1 A. Durán, J. M. Sanz-Serna. The numerical integration of relative equilibrium solutions. Geometric theory. Nonlinearity 11 (1998), 1547–1567.

DuranSS2 A. Durán, J. M. Sanz-Serna. The numerical integration of relative equilibrium solutions. The nonlinear Schrödinger equation. IMA J. Numer. Anal. 20 (2000), 235–261.

FC18 G. Frasca-Caccia. Bespoke finite difference methods that preserve two local conservation laws of the modified KdV equation. AIP Conf. Proc. 2116 (2019), 140004.

FCHydon G. Frasca-Caccia, P. E. Hydon. Simple bespoke preservation of two conservation laws. IMA J. Numer Anal. (in press) doi:10.1093/imanum/dry087.

DeFrutos J. de Frutos, J. M. Sanz-Serna. Accuracy and conservation properties in numerical integration: the case of the Korteweg–de Vries equation. Numer. Math. 75 (1997), 421–445.

Furihata D. Furihata. Finite difference schemes for ∂ u/∂ t = (∂/∂ x)^αδ G/δ u that inherit energy conservation or dissipation property. J. Comput. Phys. 156 (1999), 181–205.

Furibook D. Furihata, T. Matsuo. Discrete variational derivative method: a structure-preserving numerical method for partial differential equations. CRC Press, Boca Raton, FL, 2011.

Gonz O. Gonzales. Time integration and discrete Hamiltonian systems. J. Nonlinear Sci. 6 (1996), 449–467.

Grant 
T. J. Grant.  
Bespoke finite difference schemes that preserve multiple conservation laws.
LMS J. Comput. Math. 18 (2015), 372–403.

GrantHydon 
T. J. Grant,  P. E.  Hydon.  
Characteristics of conservation laws for difference equations. Found. Comput. Math. 13 (2013), 667–692.

HLWbook E. Hairer, C. Lubich, G. Wanner. Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations.
2nd edn. Berlin, Springer.

Hydon P. E. Hydon. Conservation laws of partial difference equations with two independent variables.
J. Phys. A. 34 (2001), 10347–10355.

Hydonbook P. E. Hydon. Difference Equations by Differential Equation Methods. Cambridge University Press, Cambridge, 2014.

HydonMans  P. E. Hydon, E. L. Mansfield.  A variational complex for difference equations. Found. Comput. Math. 4 (2004), 187–217.

Kuper 
B. A. Kuperschmidt. Discrete Lax equations and differential-difference calculus, Astérisque No. 123, 1985.

Leimku B. Leimkuhler, S.  Reich. Simulating Hamiltonian dynamics. Cambridge University Press,
Cambridge, 2004.

Marsden J. E. Marsden, G. W. Patrick, and S. Shkoller. Multisymplectic geometry, variational integrators, and nonlinear PDEs. Commun. Math. Phys. 199 (1998), 351–395.

QuispMS F. McDonald, R. I.  McLachlan, B. E.  Moore, G. R. W.  Quispel. Travelling wave solutions of multisymplectic discretizations of semi-linear wave equations. J. Difference Equ. Appl. 22 (2016), 913–940.

McLQuisp R. I. McLachlan, G. R. W. Quispel. Discrete gradient methods have an energy conservation law. Discrete Contin. Dyn. Syst. 34 (2014), 1099–1104.

Robi R. I. McLachlan, G. R. W. Quispel, N. Robidoux. Geometric integration using discrete gradients. R. Soc. Lond. Philos. Trans. Ser. A Math.
Phys. Eng. Sci. 357 (1999), 1021–1045.

Miura R. Miura, C. S. Gardner, M. D. Kruskal. Korteweg-de Vries equation and generalizations. II. Existence of conservation laws and constant of motion. J. Math. Phys. 9 (1968), 1204.

oliver M. Oliver, M. West, C. Wulff. Approximate momentum conservation for spatial semidiscretization
of semilinear wave equations. Numer. Math. 97 (2004), 493–535.

Olverpaper P. J. Olver. Evolution equations possessing infinitely many symmetries. J. Math. Phys. 18 (1977), 1212–1215.

QMcLaren G. R. W. Quispel, D. I.  McLaren. A new class of energy-preserving numerical integration methods.
J. Phys. A 41 (2008), 045206.

QuispTurn G. R. W. Quispel, G. S. Turner. Discrete gradient methods for solving ODEs numerically while preserving a first integral. J. Phys. A 29 (1996), 341–349.

Zhang D. J. Zhang, S. L. Zhao, Y. Y. Sun, J.  Zhou. Solutions to the modified Korteweg-de Vries equation. Rev. Math. Phys. 26 (2014), 1430006. ]