Linearly implicit structure-preserving schemes for Hamiltonian systems

TL;DR

This paper investigates and compares linearly implicit, structure-preserving numerical schemes for Hamiltonian systems, demonstrating their effectiveness on specific PDEs like Korteweg-de Vries and Camassa-Holm equations.

Contribution

It introduces and analyzes linearly implicit schemes that preserve modified energy in Hamiltonian systems, with applications to important PDEs.

Findings

Schemes effectively preserve modified energy.

Numerical results demonstrate stability and accuracy.

Applicable to integrable PDEs like KdV and Camassa-Holm.

Abstract

Kahan's method and a two-step generalization of the discrete gradient method are both linearly implicit methods that can preserve a modified energy for Hamiltonian systems with a cubic Hamiltonian. These methods are here investigated and compared. The schemes are applied to the Korteweg-de Vries equation and the Camassa-Holm equation, and the numerical results are presented and analysed.

Figures

Equations

$$\begin{array}{rl}\dot{x}& =f(x)=S\nabla H(x),x\in {\mathbb{R}}^d,\\ x(0)& =x_0,\end{array}$$

$$\frac{d}{dt}H(x)=\nabla H(x{)}^T\frac{dx}{dt}=\nabla H(x{)}^{T}S\nabla H(x)=0,$$

$${\Psi }_{y_0}(t{)}^{T}S{\Psi }_{y_0}(t)=S,$$

$$H(y)-H(x)=(y-x{)}^T\overline{\nabla }H(x,y).$$

$$\frac{x^{n+1}-x^n}{\Delta t}=S\overline{\nabla }H(x^n,x^{n+1}),$$

$$\tilde{H}(x,x)$$

$$\tilde{H}(x,y)$$

$$\tilde{H}(y,z)-\tilde{H}(x,y)$$

$$\overline{\nabla }\tilde{H}(x,x,x)$$

$$\frac{x^{n+2}-x^n}{2\Delta t}=S\overline{\nabla }\tilde{H}(x^n,x^{n+1},x^{n+2}).$$

$$\tilde{H}(x^{n+1},x^{n+2})-\tilde{H}(x^n,x^{n+1})$$

$$=\Delta t\overline{\nabla }\tilde{H}(x^n,x^{n+1},x^{n+2}{)}^{T}S\overline{\nabla }\tilde{H}(x^n,x^{n+1},x^{n+2})$$

$$=0,$$

$${\overline{\nabla }}_{\text{AVF}}\tilde{H}(x,y,z)=2{\int }_0^1{\nabla }_x\tilde{H}(\xi x+(1-\xi )z,y)\mathrm{d}\xi ,$$

$$\overline{\nabla }\tilde{H}(x,y,z)=2{\nabla }_x\tilde{H}(\frac{x+z}{2},y).$$

$$\tilde{H}(x,y)=x^{T}A(y)x+b(y{)}^{T}x+c(y),$$

$${\nabla }_x\tilde{H}(x,y)=2A(y)x+b(y),$$

$${\nabla }_x\tilde{H}(\frac{x+z}{2},y{)}^T(z-x)$$

$$=z^{T}A(y)z+b(y{)}^{T}z-x^{T}A(y)x-b(y{)}^{T}x$$

$$=\tilde{H}(y,z)-\tilde{H}(x,y).$$

$$\overline{\nabla }\tilde{H}(x,x,x)=2{\nabla }_x\tilde{H}(x,x)=\nabla H(x).$$

$${\overline{\nabla }}_{\text{IA}}\tilde{H}(x,y,z{)}_i$$

$$\bar{\partial }\tilde{H}(x,y,z{)}_i=\frac{\tilde{H}((z_1,\dots ,z_i,x_{i+1},\dots ,x_d),y)-\tilde{H}((z_1,\dots ,z_{i-1},x_i,\dots ,x_d),y)}{z_i-x_i}.$$

$$\frac{x^{n+2}-x^n}{2\Delta t}=S\sum _{i,j=1}^3{\alpha }_{ij}(H{}^{\prime \prime }(x^{n-1+i})x^{n-1+j}+\beta (x^{n-1+i})),$$

$$\frac{x^{n+2}-x^n}{2\Delta t}=\frac{1}{4}S{H}^{\prime \prime }(x^{n+1})(x^n+x^{n+2}),$$

$$\frac{x^{n+1}-x^n}{\Delta t}=-\frac{1}{2}f(x^n)+2f(\frac{x^n+x^{n+1}}{2})-\frac{1}{2}f(x^{n+1}).$$

$$\frac{x^{n+2}-x^n}{2\Delta t}=-\frac{1}{4}f(x^n)-\frac{1}{2}f(x^{n+1})-\frac{1}{4}f(x^{n+2})+f(\frac{x^n+x^{n+1}}{2})+f(\frac{x^{n+1}+x^{n+2}}{2}).$$

$$(x^n{)}^T{H}^{\prime \prime }(x^n)x^{n+1}=(x^{n+1}{)}^T{H}^{\prime \prime }(x^{n+2})x^{n+2},$$

$$\tilde{H}(x,y)=\frac{1}{6}{x}^T{H}^{\prime \prime }(\frac{x+y}{2})y,$$

$$\frac{x^{n+2}-x^n}{2\Delta t}$$

$${\nabla }_x\tilde{H}(x,y)=\frac{1}{6}{H}^{\prime \prime }(\frac{x+y}{2})y+\frac{1}{6}{H}^{\prime \prime }(\frac{y}{2})x=\frac{1}{12}{H}^{\prime \prime }(2x+y)y,$$

$$\overline{\nabla }\tilde{H}(x,y,z)=2{\nabla }_x\tilde{H}(\frac{x+z}{2},y)=\frac{1}{6}{H}^{\prime \prime }(x+y+z)y=\frac{1}{6}{H}^{\prime \prime }(y)(x+y+z).$$

$$\dot{\bar{x}}$$

Full text

Keywords: Linearly implicit methods, Hamiltonian system, energy preservation, Camassa–Holm equation, Korteweg–de Vries equation.

1 Introduction

The field of geometric numerical integration has garnered increased attention over the last three decades. It considers the design and analysis of numerical methods that can capture geometric properties of the flow of the differential equation to be modelled. These geometric properties are mainly invariants over time; they are conserved quantities such as Hamiltonian energy, angular momentum, volume or symplecticity. Among them the conservation of energy is particularly important for proving the existence and uniqueness of solutions for partial differential equations (PDEs) [1]. Numerical schemes inheriting such properties from the continuous dynamical system have been shown in many cases to be advantageous, especially when integration over long time intervals is considered [2].

For general non-linear differential equations, one may use a standard fully implicit scheme to solve a problem numerically. Then a non-linear system of equations must be solved at each time step. Typically this is done by the use of an iterative solver where a linear system is to be solved at each iteration. This quickly becomes a computationally expensive procedure, especially since the number of iterations needed in general increases with the size of the system; see a numerical example comparing the computational cost for implicit and linearly implicit methods in [3]. A fully explicit method on the other hand, may over-simplify the problem and lead to the loss of important information, and will often have inferior stability properties. The golden middle way may be found in linearly implicit schemes, i.e. schemes where the non-linear terms are discretized such that the solution at the next time step is found from solving one linear system.

Our aim is to present and analyse linearly implicit schemes with preservation properties. We consider ordinary differential equations (ODEs) that can be written in the form

[TABLE]

where $S$ is a constant skew-symmetric matrix and $H$ is a cubic Hamiltonian function. The famous geometric characteristic for equations like (1.1) is that the exact flow is energy-preserving,

[TABLE]

and symplectic if $S$ is the canonical skew-symmetric matrix:

[TABLE]

where $\Psi_{y_{0}}(t):=\frac{\partial\varphi_{t}(y_{0})}{\partial y_{0}}$, with $\varphi_{t}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, $\varphi_{t}(y_{0})=y(t)$ the flow map of (1.1) [2]. A numerical one-step method is said to be energy-preserving if $H$ is constant along the numerical solution, and symplectic if the numerical flow map is symplectic. Both the energy-preserving methods and the symplectic methods, the latter of which has the ability to preserve a perturbation of the Hamiltonian $H$ of (1.1), have their own advantages. In particular, the energy-preserving property has been found to be crucial in the proof of stability for several such numerical methods, see e.g [4]. However, there is no numerical integration method that can be simultaneously symplectic and energy-preserving for general Hamiltonian systems [5]. In this paper we will focus on energy-preserving numerical integration.

We wish to study and compare two types of existing methods with geometric properties. The first one is Kahan’s method for quadratic ODE vector fields [6], which by construction is linearly implicit, and for which the geometric properties have been studied in [7]. Kahan’s method has not been extensively studied for solving PDEs so far, with the notable exception [8]. This is a one-step method, but we will also give its formulation as a two-step method in this paper, for easier comparison to the other method to be studied. That method, which we call the polarised discrete gradient (PDG) method, is based on the multiple points discrete variational derivative method for PDEs presented by Furihata, Matsuo and coauthors in the papers [9, 10, 11] and the monograph [12]. A more general framework for such schemes is given by Dahlby and Owren in [3]. With the aim of easing the comparison to Kahan’s method, we present here the two-step method of [12, 3] as it looks for ODEs of the form (1.1). When Hamiltonian PDEs are considered, by semi-discretizing in space to obtain a system of Hamiltonian ODEs and then applying the PDG method, one may obtain the schemes of the aforementioned references; a specific scheme will depend on the choice of spatial discretization as well as the choices of some functions to be explained in the next section: the polarised energy and the polarised discrete gradient.

This paper is divided into two main parts. In the next chapter, we present the methods in consideration, and give some theoretical results on their geometric properties. In Chapter 3, we present numerical results for the Camassa–Holm equation and the Korteweg–de Vries equation, including analysis of stability and dispersion, comparing the methods.

2 Linearly implicit schemes

We will present an ODE formulation of the linearly implicit schemes presented by Furihata, Matsuo and coauthors in [9, 10, 11, 12] and by Dahlby and Owren in [3]. Inspired by the nomenclature of the latter reference, we call these schemes polarised discrete gradient methods. Then we present a special case of this polarisation method in the same framework as Kahan’s method, with the goal of obtaining more clarity in comparison of the methods.

2.1 Polarised discrete gradient methods

The idea behind the PDG methods is to generalise the discrete gradient method in such a way that a relaxed variant of the preservation property is intact, while nonlinear terms are discretized over consecutive time steps to ensure linearity in the scheme. Let us first recall the concept of discrete gradient methods. A discrete gradient is a continuous map $\overline{\nabla}H:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ such that for any $x,y\in\mathbb{R}^{d}$

[TABLE]

The discrete gradient method for (1.1) is then given by

[TABLE]

which will preserve the energy of the system (1.1) at any time step. Here and in what follows, $x^{n}$ is the numerical approximation of $x$ at $t=t_{n}$ and $x_{k}^{n}$ is the numerical approximation of the $k$th component of $x$ at $t=t_{n}$.

Restricting ourselves to two-step methods, we define the PDG methods as follows.

Definition 1**.**

For the energy $H$ of (1.1), consider the polarised energy as a function $\tilde{H}:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ satisfying the properties

[TABLE]

A polarised discrete gradient (PDG) for $\tilde{H}$ is a function $\overline{\nabla}\tilde{H}:\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ satisfying

[TABLE]

and the corresponding polarised discrete gradient scheme is given by

[TABLE]

Proposition 1**.**

The numerical scheme (2.2) preserves the polarised invariant $\tilde{H}$ in the sense that $\tilde{H}(x^{n},x^{n+1})=\tilde{H}(x^{0},x^{1})$ for all $n\geq 0$.

Proof.

[TABLE]

where the last equality follows from the skew-symmetry of $S$. ∎

We remark here that in the cases where we seek a time-stepping scheme for the system of Hamiltonian ODEs resulting from discretizing a Hamiltonian PDE in space in an appropriate manner, e.g. as described in [13], $H$ will be a discrete approximation to an integral $\mathcal{H}$. Thus a two-step PDG method and a standard one-step discrete gradient method, the latter in general fully implicit, will preserve two different discrete approximations separately to the same $\mathcal{H}$.

The task of finding a PDG satisfying (2.1) is approached differently in our two main references, [9, 10, 11, 12] and [3]. Furihata, Matsuo and coauthors apply a generalisation of the approach introduced by Furihata in [14] for finding discrete variational derivatives, while Dahlby and Owren suggest a generalisation of the average vector field (AVF) discrete gradient [15], given by

[TABLE]

where $\nabla_{x}\tilde{H}(x,y)$ is the gradient of $\tilde{H}(x,y)$ with respect to its first argument. Provided that the spatial discretization is performed in the same way, these two approaches lead to the same scheme for an $\tilde{H}$ quadratic in each of its arguments, as does a generalisation of the midpoint discrete gradient of Gonzalez [16]. Based on this, we present the most straightforward approach for finding this specific PDG for the cases we are studying in this paper:

Proposition 2**.**

Given an $\tilde{H}(x,y)$ that is at most quadratic in each of its arguments, define $\nabla_{x}\tilde{H}(x,y)$ as the gradient of $\tilde{H}$ with respect to its first argument. Then a PDG for $\tilde{H}$ is given by

[TABLE]

Proof.

We may write

[TABLE]

for some symmetric $A:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\times\mathbb{R}^{d}$, $b:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and $c:\mathbb{R}^{d}\rightarrow\mathbb{R}$. Then

[TABLE]

and

[TABLE]

Furthermore,

[TABLE]

∎

As remarked in Theorem $4.5$ of [3]: if the polarised energy $\tilde{H}(x,y)$ is at most quadratic in each of its arguments, the scheme (2.2) with the PDG (2.3) is linearly implicit.

An alternative to (2.3) could be a generalisation of the Itoh–Abe discrete gradient [17], defined by its $i$-th component

[TABLE]

where

[TABLE]

A symmetrized variant of this, given by $\overline{\nabla}_{\text{SIA}}\tilde{H}(x,y,z):=\frac{1}{2}(\overline{\nabla}_{\text{IA}}\tilde{H}(x,y,z)+\overline{\nabla}_{\text{IA}}\tilde{H}(z,y,x))$ is again identical to (2.3), whenever $\tilde{H}$ is quadratic in each of its arguments.

2.2 A general framework and Kahan’s method

For ODEs of the form (1.1), consider the two-step schemes of the form

[TABLE]

where $H{{}^{\prime\prime}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\times\mathbb{R}^{d}$ is the Hessian matrix of $H$ and $\beta(x):=2\nabla H(x)-H{{}^{\prime\prime}}(x)x$. For cubic $H$, this scheme is linearly implicit if and only if $\alpha_{33}=0$.

In this section, we first consider the case when the Hamiltonian is a cubic homogeneous polynomial, in which case the term $\beta(x)$ in (2.4) will disappear, and then generalise the results to the non-homogeneous case.

Theorem 1**.**

The scheme (2.4) with $\alpha_{21}=\alpha_{23}=\frac{1}{4}$, $\alpha_{ij}=0$ otherwise, i.e.

[TABLE]

where $x^{1}$ is found from $x^{0}$ by Kahan’s method, is equivalent to Kahan’s method over two consecutive steps, when applied to ODEs of the form (1.1) with homogeneous cubic $H$.

Proof.

As shown in [7], Kahan’s method can be written into a Runge–Kutta form

[TABLE]

Two steps of this can be written as

[TABLE]

Using that for a homogeneous cubic $H$ we have $\nabla H(x)=\frac{1}{2}H^{\prime\prime}(x)x$, $H^{\prime\prime}(x)y=H^{\prime\prime}(y)x$ and $H^{\prime\prime}(x+y)=H^{\prime\prime}(x)+H^{\prime\prime}(y)$, and inserting $f(x)=S\nabla H(x)$ in (2.7), we get (2.5). On the other hand, if we have found $x^{n+1}$ by Kahan’s method and $x^{n+2}$ by (2.7), we see that by subtracting (2.6) from (2.7) we get (2.6) with $n$ replaced by $n+1$. ∎

Remark 1**.**

The scheme (2.5) with the first step computed by Kahan’s method preserves the polarised invariant $\tilde{H}(x^{n},x^{n+1})=\frac{1}{6}(x^{n})^{T}H^{\prime\prime}(\frac{x^{n}+x^{n+1}}{2})x^{n+1}$, since Kahan’s method preserves this polarised invariant [7]. We note that the scheme (2.5) satisfies

[TABLE]

independent of how $x^{1}$ is found, following from the skew symmetry of the matrix $S$. However, it preserves the polarised invariant $\frac{1}{6}(x^{n})^{T}H^{\prime\prime}(\frac{x^{n}+x^{n+1}}{2})x^{n+1}$ only if Kahan’s method or an equivalent scheme is used to calculate $x^{1}$ from $x^{0}$.

A special case of the PDG method which preserves the same polarised Hamiltonian as Kahan’s method, can also be written on the form (2.4):

Theorem 2**.**

For a homogeneous cubic $H$ and the polarised energy given by

[TABLE]

the scheme (2.2) with the PDG (2.3) applied to (1.1) is equivalent to (2.4) with $\alpha_{21}=\alpha_{22}=\alpha_{23}=\frac{1}{6}$, $\alpha_{ij}=0$ otherwise, i.e.

[TABLE]

Proof.

[TABLE]

and thus

[TABLE]

∎

It can be shown that many well known Runge–Kutta methods performed over two consecutive steps are methods in the class (2.4) when applied to (1.1) with $H$ cubic. As two examples, the implicit midpoint method over two steps is (2.4) with $\alpha_{11}=\alpha_{33}=\frac{1}{16},\alpha_{21}=\alpha_{22}=\alpha_{23}=\frac{1}{8}$, $\alpha_{ij}=0$ otherwise, while the trapezoidal rule is (2.4) with $\alpha_{11}=\alpha_{33}=\frac{1}{8},\alpha_{22}=\frac{1}{4}$, $\alpha_{ij}=0$ otherwise. The integral-preserving average vector field method [18] over two steps is (2.4) with $\alpha_{11}=\alpha_{21}=\alpha_{23}=\alpha_{33}=\frac{1}{12},\alpha_{22}=\frac{1}{6}$, $\alpha_{ij}=0$ otherwise.

Now, in the cases where $H$ is non-homogeneous, one can use the technique employed in [7], i.e. adding one variable $x_{0}$ to generate an equivalent problem to the original one, for a homogeneous Hamiltonian $\bar{H}:\mathbb{R}^{d+1}\rightarrow\mathbb{R}$ defined such that $\bar{H}(1,x_{1},\ldots,x_{d})=H(x_{1},\ldots,x_{d})$. Also constructing the $(d+1)\times(d+1)$ skew-symmetric matrix $\bar{S}$ by adding a zero initial row and a zero initial column to $S$, we get that solving the system

[TABLE]

is equivalent to solving (1.1). Following the above results for the homogeneous $\bar{H}$ and (2.9), we can generalise Theorem 1 and Theorem 2 for all cubic $H$. Generalisations of the preservation properties follow directly; e.g., Kahan’s method and the PDG method can preserve the perturbed energy $\tilde{H}(x^{n},x^{n+1}):=\frac{1}{6}(\bar{x}^{n})^{T}\bar{H}^{\prime\prime}(\frac{\bar{x}^{n}+\bar{x}^{n+1}}{2})\bar{x}^{n+1}$ also for non-homogeneous cubic $H$.

3 Numerical experiments

To have a better understanding of the above methods, we will apply them to systems of two different PDEs: the Korteweg–de Vries (KdV) equation and the Camassa–Holm equation. We will compare our methods to the midpoint method, which is a symplectic, fully implicit method. We solve the two PDEs by discretizing in space to obtain a Hamiltonian ODE system of the type (1.1) and then applying the PDG method (denoted by PDGM), Kahan’s method (Kahan) and the midpoint method (MP) to this.

3.1 Camassa–Holm equation

In this section, we consider the Camassa–Holm equation

[TABLE]

defined on the periodic domain $\mathbb{S}:=\mathbb{R}/L\mathbb{Z}$. It has the conserved quantities

[TABLE]

Here we consider the variational form of the Hamiltonian $\mathcal{H}_{2}$:

[TABLE]

We follow the approach presented in [13] and semi-discretize the energy $\mathcal{H}_{2}$ of (3.1) as

[TABLE]

where the difference operators $\delta^{+}_{x}$ and $\delta^{-}_{x}$ are defined by

[TABLE]

For later use, we here also introduce the notation

[TABLE]

[TABLE]

and the matrices corresponding to the difference operators $\delta^{+}_{x}$, $\delta^{-}_{x}$, $\delta^{\langle 1\rangle}_{x}$, $\delta^{\langle 2\rangle}_{x}$, $\mu^{+}_{x}$ and $\mu^{-}_{x}$, which are denoted by $D^{+}$, $D^{-}$, $D^{\langle 1\rangle}$, $D^{\langle 2\rangle}$, $M^{+}$ and $M^{-}$. Denoting the numerical solution $U=[u_{1},\ldots,u_{K}]^{T}$, and by using the properties of the above difference operators, we thus get

[TABLE]

where $U_{\boldsymbol{\cdot}}^{2}$ is the elementwise square of $U$. Then the semi-discretized system for the Camassa–Holm equation becomes

[TABLE]

The above-mentioned schemes applied to (3.3) give us

[TABLE]

where $H_{2}^{{}^{\prime\prime}}(U)=3\,\text{diag}(U)+M^{-}\,\text{diag}(D^{+}U)D^{+}-D^{\langle 2\rangle}\,\text{diag}(U)$ is the Hessian of $H_{2}(U)$ and

$\overline{\nabla}\tilde{H}_{2}(U^{n},U^{n+1},U^{n+2})$ is the PDG of Proposition 2 with polarised discrete energy

[TABLE]

for some $a\in\mathbb{R}$, typically between $-1$ and $2$.

Remark 2**.**

We performed numerical experiments for finding a good choice of the parameter $a$ in PDGM (3.6) and based on these set $a=\frac{1}{2}$ in the following.

3.1.1 Numerical tests for the Camassa–Holm equation

Example 1 (Single peakon solution): In this numerical test, we consider the same experiment as in [19], where multisymplectic schemes are considered for the Camassa–Holm equation with

[TABLE]

$x\in[0,L]$, $L=40$, $t\in[0,T]$, $T=5$, spatial step size $\Delta x=0.04$ and time step size $\Delta t=0.0002$. All our methods keep a shape close to the exact solution except some small oscillatory tails, also observed in [19], resulting from the semi-discretization, see Figure 2 (the right two plots). The numerical simulations show that the global error is mainly due to the shape error111 Shape error is defined by $\epsilon_{\text{shape}}:=\underset{\tau}{\text{min}}\parallel U^{n}-u(\cdot-\tau)\parallel_{2}^{2}$, and phase error is defined by $\epsilon_{\text{phase}}:=\lvert\underset{\tau}{\text{argmin}}\parallel U^{n}-u(\cdot-\tau)\parallel_{2}^{2}-ct_{n}\rvert$, [3]., see Figure 1. In Figure 2 (the left plot), we can see that the numerical energy for all the methods oscillate, but it appears to be bounded. Here we consider also coarser grids. We observe that there appear some small wiggles for both PDGM and Kahan’s method for $\Delta t=0.02$ and long time integration $T=100$. However, the wiggles in the solution by PDGM are much more evident than those in the solution of Kahan’s method, see Figure 3 (the left two plots). We keep on increasing $\Delta t$ to $0.15$ and $0.2$; we observe that the numerical solution obtained with the PDG method with $\Delta t=0.15$ suffers from evident numerical dispersion, while Kahan’s method seems to keep the shape well when comparing to the exact wave. Spurious oscillations appear also in Kahan’s method when the time-step is increased to the value $\Delta t=0.2$, see Figure 3 (right).

Example 2 (Two peakons solution): Now we consider the initial condition

[TABLE]

where $x\in[0,L]$, $L=40$, $t\in[0,T]$, $T=5$, and $\Delta x=0.04$, $\Delta t=0.0002$. We observe that all the methods keep the shape of the exact solution very well and the numerical energy appears bounded, see Figure 5. The numerical simulation shows that the global error is mainly due to the shape error, see Figure 4. When a coarser time grid and longer time integration is considered, $\Delta t=0.02$ and $T=100$, small wiggles appear in the solution of PDGM and Kahan’s method, see Figure 6 (the left two figures). We increase $\Delta t$ to $0.2$, and observe that PDGM fails to preserve the shape of the solution, while Kahan’s method can still keep a shape close to the exact solution even though also for this method the numerical dispersion increases, see Figure 6 (right).

3.2 Korteweg–de Vries equation

For the Camassa–Holm equation, the vector field of the semi-discretized system is a homogeneous quadratic polynomial. In this section, we deal with the KdV equation, for which the vector field of the semi-discretized equation is a non-homogeneous quadratic polynomial. Kahan’s method has also previously been used for the temporal discretization of this equation, see [8].

The KdV equation

[TABLE]

on the periodic domain $\mathbb{S}:=\mathbb{R}/L\mathbb{Z}$ has the conserved Hamiltonians

[TABLE]

In the following we consider the variational form based on the Hamiltonian $\mathcal{H}_{2}$:

[TABLE]

3.2.1 Numerical schemes for the KdV equation

We discretize the energy $\mathcal{H}_{2}$ for the KdV equation (3.8) as

[TABLE]

From simple calculations, the corresponding gradient is given by

[TABLE]

and thus we have the semi-discretized form for (3.8):

[TABLE]

Applying the schemes under consideration to (3.9) gives

[TABLE]

where $H_{2}^{{}^{\prime\prime}}(U)=-6\,\text{diag}(U)-D^{\langle 2\rangle}$ is the Hessian of $H_{2}(U)$ and $\overline{\nabla}\tilde{H}_{2}(U^{n},U^{n+1},U^{n+2})$ is found as in Proposition 2, with polarised discrete energy

[TABLE]

Remark 3**.**

We perform several numerical simulations to find a good choice of the parameter $a$, and we take $a=-\frac{1}{2}$ for PDGM in the following numerical examples for the KdV equation.

3.2.2 Stability analysis of the schemes

To analyse the stability of the above methods, we perform the von Neumann stability analysis for the Kahan and PDGM schemes applied to the linearized form of the KdV equation (3.7)

[TABLE]

The equation for the amplification factor for Kahan’s method is

[TABLE]

and its root is

[TABLE]

where $\lambda:=\frac{\Delta t}{{\Delta x}^{3}}$. Since $g$ is a simple root on the unit circle, Kahan’s method is unconditionally stable for the linearized KdV equation.

The equation for the amplification factor for PDGM is

[TABLE]

The two roots of the above equation are thus

[TABLE]

where $b=\lambda(1-\cos\theta)\sin\theta$. We observe that $\lvert g_{1}\rvert=\lvert g_{2}\rvert=1$, and $g_{1}\neq g_{2}$, therefore PDGM is unconditionally stable for the linearized KdV equation.

3.2.3 Numerical tests for the KdV equation

Example 1 (One soliton solution): Consider the initial value

[TABLE]

where $x\in[0,L]$, $L=40$. We apply our schemes over the time interval $[0,T]$, $T=100$, with step sizes $\Delta x=0.05$, $\Delta t=0.0125$. From our observations, all the methods behave well. The shape of the wave is well kept by all the methods, also for long time integration, see Figure 7. The energy errors of all the methods are rather small and do not increase over long time integration, see Figure 8 (left). We then use a coarser time grid, $\Delta t=0.035$, and both methods are still stable, see Figure 9 (left two). However we observe that the global error of PDGM becomes much bigger than that of Kahan’s method. When an even larger time step-size, $\Delta t=0.04$, is considered, the solution for PDGM blows up while the solution for Kahan’s method is rather stable. In this case, the PDG method applied to the nonlinear KdV equation is unstable and the numerical solution blows up at around $t=8$. Even if we increase the time step-size to $\Delta t=0.1$, Kahan’s method still works well, see Figure 9 (middle). When $\Delta t=0.15$ is considered, we observe evident signs of instability in the solution of Kahan’s method. The solution will blow up rapidly when $\Delta t=0.2\gg\Delta x$.

Example 2 (Two solitons solution): We choose initial value

[TABLE]

and consider periodic boundary conditions $u(0,t)=u(L,t)$, where $x\in[0,L]$, $L=40$. We set the space step size $\Delta x=0.05$ and apply the aforementioned schemes on time interval $[0,T]$ with $T=100$, $\Delta t=0.001$. All the methods behave stably. The profiles of Kahan’s method and the midpoint method are almost indistinguishable, and the profiles for the midpoint method are thus not presented here. Kahan’s method and PDGM preserve the modified energy, and accordingly the energy error of all the methods are rather small over long time integration, see Figure 10 (left). After a short while the solution has two solitons; one is tall and the other is shorter, see Figure 10 (the right two plots).

When we consider a coarser time grid, $\Delta t=0.00375$, both methods are still stable, see Figure 11 (the left two plots). However, there appear more small wiggles in the solution by PDGM and we observe that the solution of PDGM will blow up rather soon, around $t=1$, for an even coarser time grid $\Delta t=0.005$. When we increase the time step size to $\Delta t=0.0125$ and consider $T=100$, the shape of the exact solution is still well preserved by Kahan’s method, even though there appear some small wiggles in the solution at around $t=100$. We observe that the solution of Kahan’s method will blow up when $\Delta t=0.05$ is considered. Similar experiments as in this subsection, but for the multisymplectic box schemes, can be found in a paper by Ascher and McLachlan [20]. However, here we consider even coarser time grid than there, and the numerical results show that Kahan’s method is quite stable, even though it is linearly implicit.

3.2.4 Dispersion analysis

We consider the traditional linear analysis of numerical dispersion relations for the numerical schemes applied to the KdV equation, getting the dispersion relation of frequency $\omega$ and wave number $\xi$ to be

[TABLE]

where $\lambda=\frac{\Delta t}{{\Delta x}^{3}}$. The dispersion curve is displayed in Figure (9) (right). We observe that Kahan’s method is better than PDGM at preserving the exact dispersion relation. This coincides with the behaviour of the methods applied to the nonlinear KdV equation shown in Section 3.2.3.

4 Conclusion

In this paper we perform a comparative study of Kahan’s method and what we call the polarised discrete gradient (PDG) method. To that end, we present a general form encompassing a class of two-step methods that includes both a specific case of the PDG method and Kahan’s method over two steps. We also compare the methods for two Hamiltonian PDEs: the KdV equation and the Camassa–Holm equation. Both Kahan’s method and the PDG method are linearly implicit methods, which will save computational cost. A series of numerical experiments has been performed here, for the KdV equation with one and two solitons, and the Camassa–Holm equation with one and two peakons. These experiments show that Kahan’s method is more stable than the PDG method. They also indicate that Kahan’s method yields more accurate results, as we have witnessed in the energy error and the shape and phase error when comparing to analytical solutions. Based on our results, we would recommend the use of Kahan’s method if one seeks a linearly implicit scheme for a Hamiltonian system with $H$ cubic.

Acknowledgements

The authors would like to thank Elena Celledoni, Takayasu Matsuo and Brynjulf Owren for initiating the work that led to this paper, and for their very helpful comments on the manuscript. The first and second authors acknowledge support from the European Union Horizon 2020 research and innovation programmes, under the Marie Skłodowska-Curie grant agreement no. 691070.