Locally conservative finite difference schemes for the modified KdV equation

TL;DR

This paper develops new finite difference schemes for the modified KdV equation that preserve key conservation laws, using a symbolic approach extended to cubic nonlinearities, resulting in highly accurate, mass-preserving numerical methods.

Contribution

It extends a symbolic conservation-preserving scheme approach to cubic nonlinear PDEs, producing new second-order accurate, mass-preserving finite difference schemes for the modified KdV equation.

Findings

Several new families of schemes preserving mass and energy or momentum.

Some schemes are highly accurate compared to existing methods.

Includes Average Vector Field schemes adapted for the modified KdV.

Abstract

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 coworkers. 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.

Equations

$$\mathcal{A}(x,t,[u])=0,$$

$$\text{Div}\mathbf{F}=D_x\{F(x,t,[u])\}+D_t\{G(x,t,[u])\},$$

$$\text{Div}\mathbf{F}=0\text{when}[\mathcal{A}=0].$$

$$u_t=\mathcal{D}\frac{\delta }{\delta u}\mathcal{H},\mathcal{H}=\int H([u])\mathrm{d}x,$$

$$x_i=x(m+i)=x(m)+i\Delta x,t_j=t(n+j)=t(n)+j\Delta t.$$

$$D_n{\mathbf{u}}_0=\tilde{\mathcal{D}}\left({\int }_0^1\nabla \hat{\mathcal{H}}\left(\xi {\mathbf{u}}_1+(1-\xi ){\mathbf{u}}_0\right)\mathrm{d}\xi \right).$$

$$\mathcal{A}\equiv u_t+u^2{u}_x+u_{xxx}=0.$$

$$\mathcal{H}=-\int \frac{1}{12}{u}^4+\frac{1}{2}u{u}_{xx}\mathrm{d}x,\mathcal{D}=D_x,$$

$$\mathcal{H}=-\int \frac{1}{2}{u}^2\mathrm{d}x,\mathcal{D}=D_x^3+\frac{1}{2}(u^2{D}_x+u{u}_x).$$

$$D_x{F}_1+D_t{G}_1=D_x\left(\frac{1}{3}{u}^3+u_{xx}\right)+D_{t}u,$$

$$D_x{F}_2+D_t{G}_2=D_x\left(\frac{1}{4}{u}^4+u{u}_{xx}-\frac{1}{2}{u}_x^2\right)+D_t\left(\frac{1}{2}{u}^2\right),$$

$$D_x{F}_3+D_t{G}_3=D_x\left(\frac{1}{2}{\left(\frac{1}{3}{u}^3+u_{xx}\right)}^2+\frac{1}{2}{u}_x{u}_t-\frac{1}{2}u{u}_{xt}\right)+D_t\left(\frac{1}{12}{u}^4+\frac{1}{2}u{u}_{xx}\right).$$

$$\int G_1\mathrm{d}x=\int u\mathrm{d}x,\int G_2\mathrm{d}x=\frac{1}{2}\int u^2\mathrm{d}x,\int G_3\mathrm{d}x=\int \frac{1}{12}{u}^4+\frac{1}{2}u{u}_{xx}\mathrm{d}x.$$

$$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}).$$

$$D_m=\frac{1}{\Delta x}(S_m-I),D_n=\frac{1}{\Delta t}(S_n-I),{\mu }_m=\frac{1}{2}(S_m+I),{\mu }_n=\frac{1}{2}(S_n+I),$$

$$\hat{\mathcal{H}}=-\sum _j\frac{1}{12}{U}_j^4+\frac{1}{2}{U}_j{D}_m^2{U}_{j-1},\tilde{\mathcal{D}}=D_m{\mu }_m{S}_m^{-1}.$$

$${\text{AVF}}_{\text{EC}}=D_m{\mu }_m\left\{\frac{1}{3}({\mu }_n{u}_{-1,0}){\mu }_n\left(u_{-1,0}^2\right)+D_m^2{\mu }_n{u}_{-2,0}\right\}+D_n{u}_{0,0}.$$

$$\begin{array}{rl}{\tilde{\mathcal{D}}}_{j}f=& D_m^3{\mu }_m{S}_m^{-2}f+\frac{1}{3}{S}_m^{-1}{\mu }_m\left\{\left(({\mu }_n{u}_{j,0}{)}^2+({\mu }_n{u}_{j,0})({\mu }_n{u}_{j+1,0})+({\mu }_n{u}_{j+1,0}{)}^2\right)D_{m}f\right\}\\ & +\frac{1}{3}\left({\mu }_n{u}_{j-1,0}+{\mu }_n{u}_{j,0}+{\mu }_n{u}_{j+1,0}\right)\left(D_m{\mu }_n{u}_{j-1,0}\right)f.\end{array}$$

$$\text{AVF}_{\text{MC}}=D_{m}\left\{\tfrac{1}{3}(\mu_{m}\mu_{n}u_{-1,0})\mu_{m}\big{(}(\mu_{n}u_{-1,0})^{2}\big{)}+D_{m}^{2}\mu_{m}\mu_{n}u_{-2,0}\right\}+D_{n}u_{0,0}.$$

$$\mathcal{A}(x,t,[u])=0,(x,t)\in V\subset {\mathbb{R}}^2,$$

$$\text{Div}\mathbf{F}=\mathcal{Q}\mathcal{A}.$$

$$\mathcal{E}=\sum _{i,j}(-D_x{)}^i(-D_t{)}^j\frac{\partial }{\partial u_{,ij}},\text{where}{u}_{,ij}=D_x^i{D}_t^j(u).$$

$$\mathcal{E}(\mathcal{Q}\mathcal{A})=0,$$

$$\tilde{\mathcal{A}}(m,n,[u_{0,0}])=0,$$

$$\text{mbox}Div\tilde{\mathbf{F}}=D_m\left\{\tilde{F}(m,n,[u])\right\}+D_n\left\{\tilde{G}(m,n,[u])\right\}=0,\text{when}[\tilde{\mathcal{A}}=0].$$

$$\text{mbox}Div\tilde{\mathbf{F}}=\tilde{\mathcal{Q}}(m,n,[u])\tilde{\mathcal{A}}(m,n,[u]);$$

$$\mathsf{E}=\sum _{i,j}{S}_m^{-i}{S}_n^{-j}\frac{\partial }{\partial u_{i,j}},$$

$$\mathsf{E}({\tilde{\mathcal{Q}}}_i\tilde{\mathcal{A}})=0,i=1,\dots p.$$

$$\mathsf{E}({\tilde{\mathcal{Q}}}_1\tilde{\mathcal{A}})=0$$

$${\mathcal{Q}}_1=1,{\mathcal{Q}}_2=u,{\mathcal{Q}}_3=\frac{1}{3}{u}^3+u_{xx}.$$

$$\tilde{\mathcal{A}}=D_m{\tilde{F}}_1+D_n{\tilde{G}}_1=0,$$

Full text

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.

1 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,

[TABLE]

is a divergence expression,

[TABLE]

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,

[TABLE]

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 $\int G\,\mathrm{d}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

[TABLE]

where $\mathcal{D}$ is a skew-adjoint differential operator, $\delta/\delta u$ denotes the variational derivative and $H$ is a given function. For every Hamiltonian PDE, the Hamiltonian functional $\mathcal{H}$ 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 [6, 7]. Multisymplectic finite difference schemes, which are extensions of symplectic methods for Hamiltonian ordinary differential equations (ODEs), preserve this conservation law [8, 9, 32, 33, 3, 2, 34].

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 [21, 16, 17, 18], 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 [8, 3, 11, 38], discrete line integral methods [5, 10] or discrete gradient methods [12, 15, 22, 23].

Among the most useful discrete gradient methods for approximating Hamiltonian PDEs is the Average Vector Field (AVF) method introduced in [12], which generalizes Quispel and McLaren’s energy-preserving approach for Hamiltonian ODEs [40]. According to [12], 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 (1).

Relative to a generic lattice point $\mathbf{n}=(m,n)$, the uniform grid points are

[TABLE]

In the following, $u$ is approximated by either the semidiscretization $U_{i}\approx u(x_{i},t)$ or the full discretization $u_{i,j}\approx u(x_{i},t_{j})$. The vectors whose $i$-th components are $U_{i}$ and $u_{i,j}$ are denoted by $\mathbf{U}$ and $\mathbf{u}_{j}$ respectively, and $\nabla$ is the differential operator whose $i$-th component is $\partial/\partial U_{i}$. Let $\widehat{\mathcal{H}}(\mathbf{U})\Delta x$ be a semidiscretization of $\mathcal{H}$ and let $\widetilde{\mathcal{D}}$ be a discrete skew-adjoint operator (typically depending on $\mathbf{u}_{0}$ and $\mathbf{u}_{1}$) that approximates $\mathcal{D}$. The AVF method approximates (1) by

[TABLE]

In [35], 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 [24, 13, 36, 27, 41] 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,

[TABLE]

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

[TABLE]

and with

[TABLE]

The mKdV equation has infinitely many independent conservation laws [39]. Those for mass, momentum and energy, respectively, are [1, 37]:

[TABLE]

For suitable (e.g. zero) boundary conditions, integrating (7)–(9) over the spatial domain yields the invariants

[TABLE]

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

[TABLE]

The forward difference operators $D_{m},D_{n}$ and forward average operators $\mu_{m},\mu_{n}$ are

[TABLE]

where $I$ is the identity operator.

The AVF method (3) for (5) uses

[TABLE]

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

[TABLE]

To apply the AVF method (3) to (6), let $\widehat{\mathcal{H}}=-\sum_{j}\tfrac{1}{2}U_{j}^{2}$ and let $\widetilde{\mathcal{D}}$ be the diagonal matrix whose entry $\widetilde{\mathcal{D}}_{j}$ in the row and column indexed by $j$ acts on functions $f(m,n,[u])$ as

[TABLE]

This yields the 10-point momentum-conserving scheme

[TABLE]

A new direct approach that enables multiple conservation laws to be preserved was introduced recently in [25, 26] and greatly simplified in [20]. 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 [20] for a non-Hamiltonian, non-integrable example. In the next sections we show that the two AVF schemes $\text{AVF}_{\text{EC}}$ and $\text{AVF}_{\text{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 2 describes the simplified direct approach to finding finite difference schemes that preserve two local conservation laws, for mass and either momentum or energy. In §3, we apply this strategy to the mKdV equation and obtain several families of conservative schemes. In §4, 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.

2 A strategy for preserving conservation laws

Given a (not necessarily Hamiltonian) PDE,

[TABLE]

a conservation law is in characteristic form if it satisfies

[TABLE]

The function $\mathcal{Q}$ is called the characteristic of the conservation law. For simplicity, we assume throughout that the domain $V$ is contractible.

Remark 1

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

[TABLE]

From Remark 1, if

[TABLE]

there exists $\mathbf{F}$ such that $\text{Div}\,\mathbf{F}=\mathcal{Q}\mathcal{A}$ is a conservation law.

A finite difference approximation of (11) on the grid is a partial difference equation,

[TABLE]

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,

[TABLE]

This difference conservation law is in characteristic form if

[TABLE]

just as for the continuous case, the multiplier function $\widetilde{\mathcal{Q}}$ is called the characteristic (see [29] for more details). The following result, due to Kuperschmidt [31] and generalized in [30], is a discrete version of Remark 1 and plays a pivotal role in our approach.

Remark 2

The kernel of the difference Euler operator

[TABLE]

is the vector space of all difference divergences (14).

Given a PDE (11) with $p$ conservation laws in characteristic form (12) that one wishes to preserve, our strategy is to seek approximations $\widetilde{\mathcal{Q}}_{i}$ and $\widetilde{\mathcal{A}}$ such that

[TABLE]

From Remark 2, there exist $\widetilde{\mathbf{F}}_{i}$ such that each $\widetilde{\mathcal{Q}}_{i}\widetilde{\mathcal{A}}=\text{Div}\,\widetilde{\mathbf{F}}_{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 $\mathcal{A}$ and all $\mathcal{Q}_{i}$. 2. 2.

On the given stencil, the most general finite difference approximations $\widetilde{\mathcal{A}}$ and $\widetilde{\mathcal{Q}}_{1}$ will depend on a large number of coefficients. 3. 3.

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

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. 5. 5.

Some of the remaining parameters are determined by solving

[TABLE]

symbolically, whenever a solution exists. The discrete flux $\widetilde{F}_{1}$ and density $\widetilde{G}_{1}$, which satisfy $\widetilde{\mathcal{Q}}_{1}\widetilde{\mathcal{A}}=D_{m}\widetilde{F}_{1}+D_{n}\widetilde{G}_{1}$, can be reconstructed from the characteristics [28]. 6. 6.

Steps 2 onward are iterated to preserve more conservation laws. If $\mathsf{E}(\widetilde{\mathcal{Q}}_{i}\widetilde{\mathcal{A}})\!=\!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 (15) 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 [14]) 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 [20] 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 [19] 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.

3 Conservative methods for the mKdV equation

Each of the conservation laws (7)–(9) is in characteristic form; their characteristics are

[TABLE]

For simplicity, we consider only one-step schemes for the mKdV equation (4). Therefore, our stencils are as shown in Figure 1, with $B\!-\!A\geq 4$.

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

[TABLE]

so that the mass conservation law (7) is preserved.

From this starting-point, the strategy described in §2 is used to preserve either (8) or (9). Linear terms in $F_{1}$, $G_{1}$ and $\mathcal{Q}_{2}$ or $\mathcal{Q}_{3}$ are approximated by linear combinations of the values of $u$ at the stencil points, with undetermined coefficients:

[TABLE]

Similarly, the cubic terms in $\mathcal{Q}_{3}$ and $F_{1}$ are approximated by cubics:

[TABLE]

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

Remark 3

An approximation of a conservation law in the form (14) is second-order accurate at the centre $(x,t)$ if the approximations of $\widetilde{F}$ and $\widetilde{G}$ are second-order accurate at the points $(x-\Delta x/2,t)$ and $(x,t-\Delta t/2)$, respectively.

We now use the simplified strategy to obtain families of schemes, all of which depend on parameters that are $\mathcal{O}(\Delta x^{2},\Delta 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 1, 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:

[TABLE]

where $\lambda=\mathcal{O}(\Delta x^{2},\Delta t^{2})$. All other terms in $\widetilde{\mathcal{A}}$ and $\widetilde{Q}_{3}$ are of the form (19) or (20), subject only to second-order consistency and

[TABLE]

In this way, we find a family of schemes that depends on two parameters. One of these parameters is $\mathcal{O}(\Delta x^{4},\Delta t^{4})$; as this is negligible on fine grids, we set this parameter to zero, obtaining a family of schemes $\widetilde{\mathcal{A}}(\lambda)$ with

[TABLE]

Each of these schemes preserves the following discrete version of the conservation law (9):

[TABLE]

with $\widetilde{\mathcal{Q}}_{3}=\mu_{m}\widetilde{F}_{1}$ and

[TABLE]

In both $\widetilde{G}_{3}$ and $\widetilde{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 $\Delta t<\Delta x$, the leading error terms are $\mathcal{O}(\Delta x^{2})$. We use the notation

[TABLE]

the parameter $\lambda_{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

[TABLE]

and choosing the following compact approximations of $G_{1}$ and $Q_{2}\,$:

[TABLE]

Then $F_{1}$ is determined from the general forms (19) and (20) by requiring consistency and

[TABLE]

This gives a three-parameter family of schemes that preserve mass and momentum.111Without the constraint on $\widetilde{G}_{1}$ in (21), there is another parameter. Removing (instead) the constraint on $\widetilde{Q}_{2}$ yields two other families of schemes. Two of the parameters are $\mathcal{O}(\Delta x^{4},\Delta t^{4})$, so we set these to zero and obtain the one-parameter family $\mathcal{A}(\lambda)$ with $\widetilde{G}_{1}$ given in (21) and

[TABLE]

where $\lambda=\mathcal{O}(\Delta x^{2},\Delta t^{2})$. For any value of $\lambda$, the discrete momentum conservation law

[TABLE]

is preserved, with $\widetilde{\mathcal{Q}}_{2}$ given in (21) and

[TABLE]

All terms in $\widetilde{F}_{1}$, $\widetilde{F}_{2}$ and $\widetilde{G}_{2}$ that do not approximate any quantity in the continuous expressions vanish as the stepsizes tend to zero. They are identically zero if $\lambda=0$.

If $\Delta t\!<\!\Delta x$, the leading term in the local truncation error is $\mathcal{O}(\Delta x^{2})$; by choosing $\lambda=\lambda_{2}\Delta x^{2}$ optimally, one may be able to remove at least part of this error. In the numerical tests section, we use the notation

[TABLE]

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 1. 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 $\widetilde{G}_{1}$ and the cubic term of $\widetilde{\mathcal{Q}}_{3}$ on the most compact possible sub-stencils. This gives the one-parameter family of schemes found in [19], namely $\widetilde{\mathcal{A}}(\lambda)$ with $\widetilde{F}_{1}=\mu_{m}\varphi_{-1,0}$, $\widetilde{G}_{1}=u_{0,0}$, where

[TABLE]

where $\lambda=\mathcal{O}(\Delta x^{2},\Delta t^{2})$. The discrete local energy conservation law satisfied by each scheme in this family is

[TABLE]

where $\widetilde{\mathcal{Q}}_{3}=\varphi_{0,0}$ and

[TABLE]

This family of schemes is written as

[TABLE]

Note that EC${}_{10}(0)$ is the AVF scheme AVF${}_{\text{EC}}$ that was derived in §1.

Momentum-conserving methods.

The complexity of the symbolic computation for solving

[TABLE]

is reduced by using the most compact second-order approximations of $\widetilde{G}_{1}$ and $\widetilde{\mathcal{Q}}_{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 $\widetilde{F}_{1}$, ignoring parameters that are $\mathcal{O}(\Delta x^{4},\Delta t^{4})$. This yields a one-parameter family $\widetilde{\mathcal{A}}(\lambda)$ with

[TABLE]

where $\lambda=\mathcal{O}(\Delta x^{2},\Delta t^{2})$. The discrete local conservation law for momentum is

[TABLE]

where $\widetilde{\mathcal{Q}_{2}}=\,\mu_{n}\mu_{m}u_{-1,0}$ and

[TABLE]

Again, all quantities that do not approximate any term in the corresponding continuous conservation law are identically zero when $\lambda=0$. We denote this family of schemes by

[TABLE]

The simplest scheme, MC${}_{10}(0)$, is AVF${}_{\text{MC}}$.

4 Numerical tests

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

[TABLE]

The second method is multisymplectic and amounts to

[TABLE]

This compact one-step scheme is a more efficient version of the popular Preissmann scheme [4, 2]. Each of these schemes is in divergence form, preserving a discrete version of the mass conservation law (7).

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

[TABLE]

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 (10) as follows:

[TABLE]

Whenever $\widetilde{G}_{2}$ or $\widetilde{G}_{3}$ is not defined for one of the schemes considered here, the corresponding error is instead evaluated as

[TABLE]

where $v_{i,j}=u(a+i\Delta x,\,j\Delta t)$ for schemes defined on the 10-point stencil and

[TABLE]

for 8-point schemes.

For each numerical test and family of schemes, we state the parameter value that minimizes the solution error (22). None of the schemes preserve three conservation laws, so we state the parameter values $\lambda_{i}$ that optimize the error in the unpreserved invariant given by (24) or (25). We also include the results for the simplest schemes, obtained by setting each $\lambda_{i}$ to zero. As EC${}_{10}(0)\equiv$ AVF${}_{\text{EC}}$ and MC${}_{10}(0)\equiv$ AVF${}_{\text{MC}}$, this enables comparison with the AVF schemes.

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

[TABLE]

here

[TABLE]

We set

[TABLE]

and solve this problem on $\Omega=[-20,20]\times[0,10]$, using step lengths $\Delta x=0.1$ and $\Delta t=0.025$.

The schemes $\mbox{EC}_{8}(\lambda_{1})$, $\mbox{MC}_{8}(\lambda_{2})$, $\mbox{EC}_{10}(\lambda_{3})$ and $\mbox{MC}_{10}(\lambda_{4})$ give the minimal solution error when $\lambda_{1}=1$, $\lambda_{2}=-0.077$, $\lambda_{3}=0.04$ and $\lambda_{4}=0.19$. The values $\lambda_{1}=-0.05$, $\lambda_{2}=-0.073$, $\lambda_{3}=0.20$ and $\lambda_{4}\ll 0$ minimize the error in the invariant that is not preserved by the scheme. As the modulus of $\lambda_{4}$ is very large, the corresponding term in MC${}_{10}(\lambda_{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 1 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,

[TABLE]

where $x_{i}$ (resp. $\tilde{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

[TABLE]

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

Table 1 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${}_{\phi_{1}}$, Err${}_{\phi_{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 MC8 and EC10 schemes that are more accurate than the narrow box and multisymplectic schemes. This is not true for EC8 and MC10 schemes.

The upper plot in Figure 2 shows the initial condition and the numerical solution at $T=10$ given by the most accurate of our schemes, $\mbox{EC}_{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 $\mbox{EC}_{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 $\mbox{EC}_{10}(0)$ and more so by $\mbox{EC}_{10}(0.04)$.

We now study the same problem on a coarser grid with $\Delta x=0.2$ and $\Delta t=0.05$. For each family, the parameter values minimizing the error in the solution or the unpreserved invariant are as given in Table 2; 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 3 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 [42]),

[TABLE]

on the domain $\Omega=[-2,2]\times[0,0.4]$, setting $\Delta x=0.02$ and $\Delta t=0.002$.

Table 3 shows the error in the conservation laws and the solution at the final time: $\mbox{EC}_{8}(2.22)$ is the most accurate scheme, while $\mbox{MC}_{8}(-0.165)$, $\mbox{EC}_{10}(0.92)$ and $\mbox{MC}_{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 $\lambda_{1}=0.49$, $\lambda_{2}=-0.128$, $\lambda_{3}=0.78$ and $\lambda_{4}\ll 0$. As for the two-soliton problem, minimizing the error in the non-conserved invariant is a poor criterion for choosing the parameter in MC10, so we do not show this result. By contrast, choosing the values of the parameters in $\mbox{MC}_{8}$ and $\mbox{EC}_{10}$ that minimize the error in the unpreserved invariant yields fairly accurate approximations.

The upper part of Figure 4 shows the numerical solution given by $\mbox{EC}_{8}(2.22)$. The lower part compares the exact solution and the numerical solutions given by EC8(2.22), the multisymplectic scheme and EC10(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.

5 Conclusions

The approach introduced in [20], 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.