In recent years, much attention has been given to the design and analysis of numerical methods for differential equations that can capture geometric properties of the exact flow. The increased interest in this subject can mainly be attributed to the superior qualitative behaviour over long time integration of such structure-preserving methods, see [1, 2, 3]. A popular class of structure-preserving methods are energy-preserving methods. In particular, the energy preservation property has been found to be crucial in the proof of stability for several of these numerical methods, see e.g .
. It is also highly conceivable that the ideas behind the finite-dimensional setting can be extended to the infinite-dimensional Hamiltonian systems or Hamiltonian partial differential equations (PDEs)
. There are two popular ways to construct energy-preserving methods for Hamiltonian PDEs. One approach is to semi-discretize the PDE in space so that one obtains a system of Hamiltonian ordinary differential equations (ODEs), and then apply an energy-preserving method to this semi-discrete system, see for example. In this way, it is straightforward to generalise the energy-preserving methods for finite-dimensional Hamiltonian systems to Hamiltonian PDEs. However, such methods conserve only a global energy that relies on a proper boundary condition, such as a periodic boundary condition. If this is not present, the energy-preserving property will be destroyed. The other approach is based on a reformulation of the Hamiltonian PDE into a multi-symplectic form, which provides the PDE with three local conservation laws: the multi-symplectic conservation law, the energy conservation and the momentum conservation law [10, 11, 12]. Then one may consider methods that preserve the local conservation laws, see for example . These locally defined properties are not dependent on the choice of boundary conditions, giving the methods that preserve local energy an advantage over methods that preserve a global energy, especially since local conservation laws will always lead to global conservation laws whenever periodic boundary conditions are considered. The concept of a multi-symplectic structure for PDEs was introduced by Bridges in [10, 11], see also  for a framework based on a Lagrangian formulation of the Cartan form. Local energy-preserving methods were first studied in , and have garnered much interest recently, see for example [13, 16, 17].
Most of the local energy-preserving methods proposed so far are fully implicit methods, for which a non-linear system must be solved at each time step. This is normally done by using an iterative solver where a linear system is solved at each iteration, which can lead to computationally expensive procedures, especially since the number of iterations needed in general increases with the size of the system. A fully explicit method on the other hand, may over-simplify the problem and often has inferior stability properties, so that a strong restriction on the grid ratio is needed. A good alternative may therefore be to develop linearly implicit schemes, where the solution at the next time step is found by solving only one linear system.
One example of linearly implicit methods for Hamiltonian ODEs is Kahan’s method, which was designed for solving quadratic ODEs  and whose geometric properties have been studied in a series of papers by Celledoni et al. [19, 20, 21]. For Hamiltonian PDEs, Matsuo and Furihata proposed the idea of using multiple points to discretize the variational derivative and thus design linearly implicit energy-preserving schemes . Dahlby and Owren generalised this concept and developed a framework for deriving linearly implicit energy-preserving multi-step methods for Hamiltonian PDEs with polynomial invariants . A comparison of this approach and Kahan’s method applied to PDEs is given in 
. Recently, more work has been put into developing linearly implicit energy-preserving schemes for Hamiltonian PDEs, e.g. the partitioned averaged vector field (PAVF) method and schemes based on the invariant energy quadratization (IEQ) approach  or the multiple scalar auxiliary variables (MSAV) approach . However, little attention has been given to linearly implicit local energy-preserving methods. To the best of the authors’ knowledge, the only existing method is one based on the IEQ approach, specific for the sine-Gordon equation . In this paper, we use Kahan’s method to construct a linearly implicit method that preserves a discrete approximation to the local energy for multi-symplectic PDEs with a cubic energy function.
The rest of this paper is organized as follows. First, we give an overview of Kahan’s method and formulate it by using a polarised energy function. A brief introduction to multi-symplectic PDEs and their conservation laws are presented in Section 3. In Section 4, new linearly implicit local and global energy-preserving schemes are presented. Numerical examples for the Korteweg–de Vries (KdV) and Zakharov–Kuznetsov equations are given in Section 5, before we end the paper with some concluding remarks.
2 Kahan’s method
Consider an ODE system
where is an valued quadratic form, is a symmetric constant matrix, and is a constant vector. Kahan’s method is then given by
is the symmetric bilinear form obtained by polarisation of the quadratic form . Polarisation, which maps a homogeneous polynomial function to a symmetric multi-linear form in more variables, was used to generalise Kahan’s method to higher degree polynomial vector fields in .
Suppose we restrict the problem (2.1) to be a Hamiltonian system on a Poisson vector space with a constant Poisson structure:
is a constant skew-symmetric matrix, andis a cubic polynomial function. We first consider the Hamiltonian to be homogeneous. Then, following the result in Proposition of , Kahan’s method can be reformulated as
where is a symmetric
-tensor satisfying. Consider the -tensor , where , with being the Hessian of ; then we can rewrite Kahan’s method (2.3) as
where denotes the partial derivative with respect to the first argument of .
Consider then the cases where the Hamiltonian in problem (2.2) is non-homogeneous, i.e. of the general form
where is the linear part of and thus a symmetric matrix whose elements are homogeneous linear polynomials, is the constant part of and thus a symmetric constant matrix, is a constant vector and is a constant scalar. We follow the technique in , adding one variable to to get , extending to by adding a zero initial row and a zero initial column, considering a homogeneous function based on the non-homogeneous Hamiltonian such that , and finally solving instead of (2.2) the equivalent, homogeneous cubic Hamiltonian problem
with . In this way we can still get the reformulation of Kahan’s method as (2.4) with
The -valued function in (2.6) has the following properties:
is symmetric111Denote the elements in by , where , , are scalars and is the th element of . We have that satisfies since , which is unchanged under any permutation of . This provides the symmetry of . w.r.t. , and ,
is symmetric w.r.t. and .
In this paper, we will use the form of Kahan’s method in (2.4) to prove the energy preservation of the proposed methods.
3 Conservation laws for multi-symplectic PDEs
Many PDEs, including all one-dimensional Hamiltonian PDEs, can be written on the multi-symplectic form
where , are two constant skew-symmetric matrices and is a scalar-valued function. Following the results about multi-symplectic structure in , it can be shown that multi-symplectic PDEs satisfy the following local conservation laws : the multi-symplectic conservation law
the local energy conservation law (LECL)
and the local momentum conservation law (LMCL)
where and satisfy
Decomposition of the matrices is done to make deduction of the conservations laws for energy and momentum more efficient [12, Section 12.3.1].
The multi-symplectic form (3.1) can also be generalised to problems in higher dimensional spaces. Consider spatial dimensions; based on the work by Bridges , a multi-symplectic PDE can then be written as
where , are constant skew-symmetric matrices and is a smooth functional. Equation (3.3) has the following local energy conservation law:
where , , and are splittings of satisfying .
Say we have (3.3) defined on the spatial domain with periodic boundary conditions. Integrating over the domain on both sides of the equation (3.4) and using the periodic boundary condition then leads to the global energy conservation law for the multi-symplectic PDEs,
Korteweg–de Vries equation. Consider the KdV equation for modeling shallow water waves,
where . Introducing the potential , momenta and the variable by the covariant Legendre transform from the Lagrangian, we obtain
from which we find the multi-symplectic formulation (3.1) for the KdV equation with , the Hamiltonian , and
As for the conservation laws, there are many choices of and , for example , or and being the upper triangular parts of and , respectively.
Zakharov–Kuznetsov equation. Zakharov and Kuznetsov introduced in  a (2+1)-dimensional generalisation of the KdV equation which includes weak transverse variation,
A multi-symplectification of this leads to a system (3.3) for two spatial dimensions,
which is (3.9) with , the Hamiltonian , and the skew-symmetric matrices
4 New linearly implicit energy-preserving schemes
In , Gong, Cai and Wang present a scheme that preserves the local energy conservation law (3.4) of a one-dimensional multi-symplectic PDE, obtained by applying the midpoint rule in space and the averaged vetor field (AVF) method in time. They also present schemes that preserve the global energy, but not (3.4), obtained by considering spatial discretizations that preserve the skew-symmetric property of the difference operator . We build on their work by considering Kahan’s method for the discretization in time, ensuring linearly implicit schemes and also energy preservation.
To introduce our new schemes, we begin with some basic difference operators:
The operators satisfy the following properties :
All the operators commute with each other, e.g.
They satisfy the discrete Leibniz rule
One can obtain a series of similar commutative equations and discrete Leibniz rules that are not presented here, but which are also crucial in the proofs of the preservation properties of the schemes to be introduced in the remainder of this section.
4.1 A local energy-preserving scheme for multi-symplectic PDEs
In this section, we apply the midpoint rule in space and Kahan’s method in time to construct a local energy-preserving method for multi-symplectic PDEs. Introducing the concept by first considering the one-dimensional system (3.1), we apply the midpoint rule in space to get
Then applying Kahan’s method gives us the linearly implicit local energy-preserving (LILEP) scheme
The scheme (4.1) satisfies the discrete local energy conservation law
Taking the inner product with on both sides of (4.1) and using the skew-symmetry of matrix , we have
Taking the inner product with on both sides of (4.1), we get
Taking the inner product with on both sides of the scheme (4.1) for the next time step, we get
On the other hand, using the aforementioned commutative laws and discrete Leibniz rules for the operators, we can deduce
The polarised global energy may be considered as a function of the solution in time step only, similarly to the modified Hamiltonian defined in Proposition 3 of .