{'fraction_non_alphanumeric': 0.04544201750155034, 'fraction_numerical': 0.01633018672913939, 'mean_word_length': 4.374375115719311, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Linearly implicit local and global energy-preserving methods for Hamiltonian PDEs

We present linearly implicit methods that preserve discrete approximations to local and global energy conservation laws for multi-symplectic PDEs with cubic invariants. The methods are tested on the one-dimensional Korteweg-de Vries equation and the two-dimensional Zakharov-Kuznetsov equation; the numerical simulations confirm the conservative properties of the methods, and demonstrate their good stability properties and superior running speed when compared to fully implicit schemes.


Introduction
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 structurepreserving 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 [4].
Energy-preserving methods are well studied for finite-dimensional Hamiltonian systems [5,6,7,8]. 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) [9]. 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 [8]. 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 energypreserving 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 [13]. 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 [14] for a framework based on a Lagrangian formulation of the Cartan form. Local energy-preserving methods were first studied in [15], 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 [18] 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 [22]. 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 [23]. A comparison of this approach and Kahan's method applied to PDEs is given in [24]. 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 [25] and schemes based on the invariant energy quadratization (IEQ) approach [26] or the multiple scalar auxiliary variables (MSAV) approach [27]. 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 [28]. 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 Kortewegde Vries (KdV) and Zakharov-Kuznetsov equations are given in Section 5, before we end the paper with some concluding remarks.

Kahan's method
Consider an ODE systemẏ = f (y) =Q(y) +B y +ĉ, y ∈ R M , (2.1) whereQ(y) is an R M valued quadratic form,B ∈ R M ×M is a symmetric constant matrix, and c ∈ R M is a constant vector. Kahan's method is then given by is the symmetric bilinear form obtained by polarisation of the quadratic formQ [19]. 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 [29]. Suppose we restrict the problem (2.1) to be a Hamiltonian system on a Poisson vector space with a constant Poisson structure:ẏ = A∇H (y), where A is a constant skew-symmetric matrix, and H : R M → R is a cubic polynomial function. We first consider the Hamiltonian H to be homogeneous. Then, following the result in Proposition 2.1 of [29], Kahan's method can be reformulated as y n+1 − y n ∆t = 3AH (y n , y n+1 , ·), where ∂H ∂x denotes the partial derivative with respect to the first argument ofH. Consider then the cases where the Hamiltonian in problem (2.2) is non-homogeneous, i.e. of the general form H (y) = y T Q(y)y + y T B y + c T y + d , (2.5) where Q(y) is the linear part of ∇ 2 H (y) and thus a symmetric matrix whose elements are homogeneous linear polynomials, B is the constant part of ∇ 2 H (y) and thus a symmetric constant matrix, c is a constant vector and d is a constant scalar. We follow the technique in [19], adding one variable to y = (y 1 , . . . , y M ) T to getỹ = (y 0 , y 1 , . . . , y M ) T , extending A toÃ by adding a zero initial row and a zero initial column, considering a homogeneous functionH(ỹ) based on the non-homogeneous Hamiltonian H (y) such thatH(ỹ)| y 0 =1 = H (y), and finally solving instead of (2.2) the equivalent, homogeneous cubic Hamiltonian probleṁ y =Ã∇H (ỹ) with y 0 = 1. In this way we can still get the reformulation of Kahan's method as (2.4) with (2.6) Remark 1. The R-valued functionH(x, y, z) in (2.6) has the following properties: x, y and z,

3.
∂H (x,y,z) ∂x In this paper, we will use the form of Kahan's method in (2.4) to prove the energy preservation of the proposed methods.

Conservation laws for multi-symplectic PDEs
Many PDEs, including all one-dimensional Hamiltonian PDEs, can be written on the multisymplectic form where K , L ∈ R l ×l are two constant skew-symmetric matrices and S : R l → R is a scalar-valued function. Following the results about multi-symplectic structure in [10], it can be shown that multi-symplectic PDEs satisfy the following local conservation laws [30]: the multi-symplectic conservation law the local energy conservation law (LECL) and the local momentum conservation law (LMCL) 1 Denote the elements in Q(y) by q i j y = k q k i j y k , where q k i j , i , j , k = 1, · · · , M , are scalars and y k is the kth element of y. We have that q k i j satisfies q k , which is unchanged under any permutation of i , j , k. This provides the symmetry ofH(x, y, z).
where K + and L + 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 d spatial dimensions; based on the work by Bridges [10], a multi-symplectic PDE can then be written as where K , L α ∈ R l ×l (α = 1, . . . , d ) are constant skew-symmetric matrices and S : R l → R is a smooth functional. Equation (3.3) has the following local energy conservation law: 3) defined on the spatial domain Ω ∈ R d 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, where η, γ ∈ R. Introducing the potential φ x = u, momenta v = γu x and the variable w = γv x φ t + γ 2 u 2 2 by the covariant Legendre transform from the Lagrangian, we obtain

7)
from which we find the multi-symplectic formulation As for the conservation laws, there are many choices of K + and L + , for example K + = K 2 , L + = L 2 , or K + and L + being the upper triangular parts of K and L, respectively.  A multi-symplectification of this leads to a system (3.3) for two spatial dimensions, (3.9) Following [9], we have that (3.8) is equivalent to a system of first-order PDEs,

New linearly implicit energy-preserving schemes
In [16], Gong, Cai and Wang present a scheme that preserves the local energy conservation law (3.2) 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.2), obtained by considering spatial discretizations that preserve the skew-symmetric property of the difference operator ∂ x . 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 [13]: 1. All the operators commute with each other, e.g.
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.

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 . (4.1) Here we consider S of the form S(y) = y T Q(y)y + y T B y + c T y + d , as in (2.5), and accordinglȳ S(x, y, z) of the form (2.6).
Theorem 1. The scheme (4.1) satisfies the discrete local energy conservation law Proof. Taking the inner product with 1 3 δ t µ x z n j on both sides of (4.1) and using the skew-symmetry of matrix K , we have .

(4.4)
Taking the inner product with 1 j on both sides of (4.1), we get . (4.5) Taking the inner product with 1 3 δ t µ x z n j on both sides of the scheme (4.1) for the next time step, we get . (4.6) Adding equations (4.4), (4.5) and (4.6) and using the skew-symmetry of matrix K , we obtain .
On the other hand, using the aforementioned commutative laws and discrete Leibniz rules for the operators, we can deduce (4.8) Using the above relations (4.8), the fact that L = L + − L T + and the result (4.7), we obtain where ∆x = P /M and (Ē L ) n j is given by (4.3).
The polarised global energyĒ n L may be considered as a function of the solution in time step n only, similarly to the modified Hamiltonian defined in Proposition 3 of [19]. Proposition 1. With the solution z n+1 found from z n by (4.1), the discrete global energyĒ n L of (4.9) satisfiesĒ Proof. Note that Inserting this in (4.3), we get (4.10) from (4.9). Furthermore, observing that we may rewrite (4.1) as which is (4.12).
Note that (4.11) is the discrete energy preserved by the fully implicit local energy-preserving method of [16]. Also, for methods based on the multi-symplectic structure, instead of solving for z directly, the normal procedure is to eliminate the auxiliary variables from the scheme and get an equation for one variable u. Therefore we do not give an explicit expression for the modified energy in z n . However, in Section 5, we present an explicit expression for the modified energy in u n when our scheme is applied to the KdV equation.
The results about the energy conservation for the LILEP method applied to one-dimensional multi-symplectic PDEs can be generalised to problems in spatial dimensions of any finite degree. Consider for example a 2-dimensional multi-symplectic PDE for which we have the following corollary. This is presented without its proof, which is rather technical but similar to the proof of Theorem 1.

Global energy-preserving methods for multi-symplectic PDEs
As shown in Section 3, Hamiltonian PDEs of the form (3.1) with periodic boundary conditions have global energy conservation which can be deduced from the local conservation law. On the other hand, the local conservation law is not inherent in the global conservation law. In this section, we will focus on giving a systematic method that preserves the global energy conservation law directly. We discretize ∂ x with an antisymmetric differential matrix D and get the semidiscretized variant of (3.1), where z := (z 0 , z 1 , . . . , z M −1 ) T ∈ R M ×l and (D z) j = M −1 k=0 D j ,k z k . We then apply Kahan's method to (4.15) and obtain the linearly implicit global energy-preserving (LIGEP) scheme . (4.16) Define the polarised energy density bȳ and we get the following result.  Proof. Taking the inner product with 1 3 δ t z n j on both sides of equation (4.16) and using the skew-symmetry of the matrix K , we get .

(4.22)
By using the commutative laws and discrete Leibniz rules, (4.23) Based on the above equations (4.22) and (4.23), we obtain δ t E n j = δ tS (z n j , z n j , z n+1 Since D is skew-symmetric and G j ,k = G k, j , we get which implies that the discrete global energy conservation lawĒ n+1 =Ē n is satisfied. The polarised energyĒ preserved by (4.16) may also be expressed as a modification of the discrete energy which is preserved by the fully implicit global energy-preserving scheme of [16]. The proof of the following proposition is similar to the proof of Proposition 1, and hence omitted.

Proposition 2.
If the solution z n+1 is found from z n by (4.16), the discrete global energyĒ n of (4.18) satisfiesĒ where g (z n j ) = ∇S(z n j ) − L(D z) n j and R(z n j ) = K + ∆t 2 ∇g (z n j ).
The above global conservation results can be generalised to multi-symplectic formulations in higher spatial dimensions, as demonstrated for the two-dimensional case by the following corollary, whose omitted proof is in the same vein as the proof of Theorem 2.

Corollary 3.
Discretizing ∂ x and ∂ y by skew-symmetric differential matrices D x and D y in equation (4.13) and then applying Kahan's method to the semi-discrete system, one obtains the linearly implicit global energy-preserving (LIGEP) scheme where j = 0, . . . , M x − 1 and k = 0, . . . , M y − 1. For periodic boundary conditions z(x + P x , y, t ) = z(x, y, t ), z(x, y + P y , t ) = z(x, y, t ), the scheme (4.25) satisfies the discrete global energy conservation lawĒ

Numerical examples
In this section, we apply our proposed new linearly implicit energy-preserving schemes to the KdV equation and Zakharov-Kuznetsov equation, and compare them with fully implicit schemes. Among our reference methods are the methods introduced in [16], for which the local energypreserving method is denoted by LEP, and the global energy-preserving method by GEP. For the GEP and LIGEP schemes, two different choices are considered for approximating the spatial derivative: the central difference operator δ c x defined by δ c x v n j := 1 2 (δ x v n j −1 + δ x v n j ) and the first order Fourier pseudospectral operator [9]. The latter results in the M × M matrix D, given explicitly by its elements evaluated on the domain [0, P ], where we assume M even and periodic boundary conditions [32]. If M is odd, we have instead

Korteweg-de Vries equation
Consider the multi-symplectic structure of the KdV equation as presented in Example 1. Applying the LILEP scheme (4.1) to (3.7), we obtain By eliminating the auxiliary varibles φ, v and w, we see that this is equivalent to Omitting the average operator µ t gives us The polarised discrete energy preserved by this scheme is On the other hand, the discrete energy preserved by the LEP method of [16] is

(5.3)
By Proposition 1 and elimination of the variables φ, v and w, (5.2) can be expressed as a modification of (5.3): we may rewrite (5.1) as where (µ x u n ) 2 denotes the element-wise square of µ x u n . Inserting this in (5.2), we get where ∇E n L means the gradient of E n L with respect to u n , and ζ L (u n ) denotes the Jacobian matrix of ζ L (u n ).
Similarly for the LIGEP method (4.16); applying it to the the multi-symplectic KdV equations (3.7) and eliminating the auxiliary varibles φ, v and w, we obtain where u n u n+1 denotes element-wise multiplication of the vectors. Omitting the average operator µ t , we get δ t u n j + η 2 (D(u n u n+1 )) j + γ 2 µ t (D 3 u n ) j = 0.

(5.4)
The discrete global energy preserved by the GEP method is while the polarised discrete energy preserved by (5.4) is In the first numerical experiment, we consider the problem introduced in [33] and then used by Zhao and Qin [34] and Ascher and McLachlan [35] to test various symplectic and multisymplectic schemes: the KdV equation with γ = 0.022, η = 1, and initial value This problem is also considered in Example 3 of [16], where it is solved by implicit schemes that preserve local and/or global energy. As observed by Gong et al., the global energy-preserving scheme (GEP) with the central difference operator used to approximate ∂ x gives unsatisfactory results for this problem; we observed that the same is true for the LIGEP scheme. Therefore, the Fourier pseudospectral operator is used to approximate the spatial derivatives in the GEP and LIGEP schemes. This seems to result in more accurate solutions than the LEP and LILEP schemes for the same number of discretization points, but at a considerably higher computational cost, as seen from Table 1. From Figure 1, we can conclude that our linearly implicit schemes give results close to their fully implicit counterparts introduced in [16], and that the different schemes converge to the same solution. Here and in the following test problem, we have solved the fully implicit schemes in each step by Newton's method until F (u n ) 2 < 10 −10 . Compared to the schemes tested in [34,35], our schemes do also perform well; see Figure  2, where we have plotted solutions by our schemes for the same discretization parameters used in Example 5.3 of [35]. The reference solution is found by the implicit midpoint scheme of [35] with θ = 1 and very fine discretization in space and time: M = 2000 and ∆t = 0.0001. We observe that the LILEP scheme behaves similarly to the multi-symplectic box scheme of Arscher and McLachlan (see figures 3 and 4 in [35]), seemingly with the same superior stability for rough discretization in space and time. The LIGEP scheme, on the other hand, starts to blow up at around t = 1 when M = 60, ∆t = 1/150, but produces for M = 100, ∆t = 0.004 a solution that is much closer to the correct solution than any of the schemes tested in [35] (see Figure 3 in that paper for comparison).

Test problem 2
To get quantitative results on the performance of our methods, we wish to study a problem with a known solution. For the KdV equation with γ = 1, η = 6, initial value u 0 (x) = 1 and periodic boundary conditions u(x + P, t ) = u(x, t ), the exact solution is a soliton moving with a constant speed c in the positive x-direction while keeping its initial shape. That is, In our numerical experiments, c = 4 and P = 20. For this problem, we have used the central difference operator to approximate ∂ x in the GEP and LIGEP schemes, since it gives good results and yields considerably shorter computational time than if the pseudospectral operator is used. The proposed methods all show very good stability conditions when applied to this problem, as expected by methods conserving some invariant. The initial shape of the wave is well kept for long integration times, even when quite large step sizes in space and time are used; Figure 3 gives a good illustration of this. As in the previous example, we again observe that little is lost in accuracy by choosing linearly implicit over fully implicit time integration. A close inspection of Figure 3 also indicates that the local energy-preserving schemes preserve the shape of the wave better than the global energy-preserving schemes, while on the other hand, the GEP and LIGEP schemes are better than the LEP and LILEP schemes at preserving the speed of the wave. This is confirmed in Table 2 by measuring the shape error and phase error where U N is the numerical solution at end time t .    In Figure 4, we have plotted the computational time required to reach a certain accuracy in the global error for the different methods, both at time t = 0.5 and at time t = 10. We compare our methods to the fully implicit LEP and GEP schemes of [16], and also to two of the schemes studied in [35]: the multi-symplectic box scheme (MSB) and the implicit midpoint scheme (IMP). Most notably we see from both plots in Figure 4 that the linearly implicit schemes perform better than the fully implicit schemes. Also, we see that at time t = 0.5 the global error is lowest for the LILEP scheme, while at t = 10 it is lowest for the LIGEP scheme. This is in accordance with the schemes' phase and shape errors, which can be observed from Figure 3 and Table 2; with increasing time, the phase error becomes more dominant, and thus the scheme with the smallest phase error becomes increasingly advantageous.  Figure 5 illustrates how the different schemes preserve a discrete approximation to the energy to machine precision. That is, the linearly implicit schemes LILEP and LIGEP preserve exactly the discrete energies (5.2) and (5.6), respectively, while keeping the discrete energies (5.3) and (5.5), respectively, within some bound which depends on the discretization parameters. Likewise, the reverse is true for the fully implicit schemes. These observations fit well with our above results about the different discrete approximations to the energy: that for both the local energy preserving and the global energy preserving schemes, either discrete energy given can be seen as a modification of the other approximation. Finally, we have included plots in Figure 6 which confirm that our schemes are of second order in space and time. Step size in time

Zakharov-Kuznetsov equation
Kahan's method is previously shown to have nice properties when applied to integrable ODE systems [19,20], and to perform well compared to other linearly implicit methods when applied to the KdV and Camassa-Holm equations [24], which are completely integrable PDEs. We wish to test our methods also on non-integrable systems, as well as on higher-dimensional problems. Therefore we consider the Zakharov-Kuznetsov equation, which is a non-integrable PDE [36,37]. This two-dimensional generalisation of the KdV equation has a variety of applications, see e.g. [38] for a brief summary. Applying the (2+1)-dimensional LILEP method (4.14) to the Zakharov-Kuznetsov equation (3.8) multi-symplectified as described in Example 2, we find Upon eliminating all variables except u, we are left with The operator µ t is again superfluous. Hence we get the scheme This scheme preserves which is a two-step discrete approximation of the energy Similarly, applying the linearly implicit global energy-preserving method (4.25) to (3.10), we get the scheme which preserves the two-step discrete energy approximation

Test problem
Taking a note from a numerical experiment performed in [9], we study the formation of cylindrical soliton pulses on the domain [0, P ] × [0, P ], P = 30, following the initial condition where ξ(y) is a random perturbation. Upon trying the different schemes we can immediately conclude that the local energypreserving schemes are superior for this problem when compared to the global energy-preserving schemes. The GEP and LIGEP schemes are too costly when the pseudospectral operator is used, and gives oscillatory behaviour in the y-direction when the central difference operator is used, unless the discretization in this direction is very fine. Although the global energy-preserving schemes with the central difference operator are slightly faster then the local energy-preserving schemes, as can be seen in Table 3, this is undermined by the cost of the extra discretization points needed to avoid oscillations in the former case. As was the case for the KdV problem, we see little difference between the linearly implicit schemes and their fully implicit counterparts. This can be seen in Figure 7, as can the oscillations in y-direction of the solution found by the GEP and LIGEP methods. The plots in Figure 7 can be compared to the plot in Figure 8, where the same problem is solved by the LILEP method using finer discretization in space and time. The initial random perturbation in y-direction over 75 points is then transferred over to 225 points using linear interpolation.  Table 3: Running time, in seconds, for computing 100 steps in time by the various schemes and various number of discretization points M = M x = M y in each spatial direction, solving our test problem for the Zakharov-Kuznetsov equation. As for the KdV equation test problems, a tolerance of 10 −10 is used when solving the fully implicit schemes by Newton's method.

Concluding remarks
In this paper, we propose two types of linearly implicit methods with conservation properties for cubic invariants of multi-symplectic PDEs. The linearly implicit local energy-preserving (LILEP) method preserves a discrete approximation to the local energy conservation law, and by extension, the global energy whenever periodic boundary conditions are considered. The linearly implicit global energy-preserving (LIGEP) method preserves the global energy without inheriting the local preservation from the continuous system. We test our methods on two PDEs: the one-dimensional, integrable Korteweg-de Vries (KdV) equation and the two-dimensional, non-integrable Zakharov-Kuznetsov equation. The numerical experiments confirm that the proposed methods are of second order both in space and time and that they preserve the expected local and global energy conservation laws. We have observed excellent stability properties for the LILEP scheme in particular, and very high accuracy in the LIGEP scheme even for quite coarse discretization when a Fourier pseudospectral operator is used to approximate the spatial derivative. Compared to the fully implicit methods of Gong et al. in [16], which was an inspiration for this paper, our methods show comparable wave profiles, global errors and energy errors, at a significantly lower computational cost. For twodimensional problems, where fully implicit schemes quickly become very expensive to compute, the combination of local energy-preservation and a linearly implicit method seems to provide for a very competitive method.
Although we have only considered the preservation of cubic invariants in this paper, our schemes can be extended to preserve higher order polynomials by the polarisation techniques for generalising Kahan's method suggested in [29]. This would result in (p − 2)-step methods for preservation of a discrete p-order polynomial invariant.