{'fraction_non_alphanumeric': 0.036913749034060514, 'fraction_numerical': 0.013225940676454854, 'mean_word_length': 4.248323194509437, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A general framework for deriving integral preserving numerical methods for PDEs

A general procedure for constructing conservative numerical integrators for time dependent partial differential equations is presented. In particular, linearly implicit methods preserving a time discretised version of the invariant is developed for systems of partial differential equations with polynomial nonlinearities. The framework is rather general and allows for an arbitrary number of dependent and independent variables with derivatives of any order. It is proved formally that second order convergence is obtained. The procedure is applied to a test case and numerical experiments are provided.


Introduction
Schemes that conserve geometric structure have been shown to be useful when studying the long time behaviour of dynamical systems. Such schemes are sometimes called geometric or structure preserving integrators [18,19]. In this paper we shall mostly be concerned with the conservation of first integrals.
Even if a presumption in this work is that the development of new and better integral preserving schemes is useful, we would still like to mention some situations where schemes with such properties are of importance. In the literature one finds several examples where stability of a numerical method is proved by directly using its conservative property, one example is the scheme developed for the cubic Schrödinger equation in [11]. Another application where the exact preservation of first integrals plays an important role is in the study of orbital stability of soliton solutions to certain Hamiltonian partial differential equations (PDEs) as discussed by Benjamin and coauthors [1,2].
For ordinary differential equations (ODEs) it is common to devise relatively general frameworks for structure preservation. This is somewhat to the contrary of the usual practice with partial differential equations where each equation under consideration normally requires a dedicated scheme. But there exist certain fairly general methodologies that can be used for developing geometric schemes also for PDEs. For example, through space discretisation of a Hamiltonian PDE one may obtain a system of Hamiltonian ODEs to which a geometric integrator may be applied. Another approach is to formulate the PDE in multi-symplectic form, and then apply a scheme which preserves a discrete version of this form, see [4] for a review of this approach.
In this paper we consider methods for PDEs that are based on the discrete gradient method for ODEs. The discrete gradient method was perhaps first treated in a systematic way by Gonzalez [16], see also [18,25]. For PDEs one may derive discrete gradients either for the abstract Cauchy problem, where the solution at any time is considered as an element of some infinite dimensional space, or one may semidiscretise the equations in space and then derive the corresponding discrete gradient for the resulting ODE system. This last procedure has been elegantly presented in several articles by Furihata, Matsuo and collaborators, see e.g. [12,13,14,21,22,23,24], using the concept of discrete variational derivatives. See also the monograph [15]. The first part of this paper develops a similar framework that is rather general and allows for an arbitrary number of dependent and independent variables with derivatives of any order. The suggested approach does not require the equations to be discretised in space.
We consider a class of conservative schemes which are linearly implicit. By linearly implicit we mean schemes which require the solution of precisely one linear system of equations in each time step. This is opposed to fully implicit schemes for which one typically applies an iterative solver that may require a linear system to be solved in every iteration. For standard fully implicit schemes one would typically balance the iteration error in solving the nonlinear system with the local truncation error. However, for conservative schemes the situation is different since exact conservation of the invariant requires that the nonlinear system is solved to machine precision. This work can be seen as a generalisation of ideas introduced in [22, Section 6].
It may not, in general, be an easy task to quantify exactly what can be expected of gain in computational cost, if any, when replacing a fully implicit scheme with a linearly implicit one. For illustration we present an example where the KdV equation is solved on a periodic domain using a fully implicit scheme and a linearly implicit scheme These schemes are derived in Section 6. For space discretisation centered differences are used for both schemes. Note that the linearly implicit scheme (1.3) has a multistep nature, but should not be confused with standard linear multistep methods. Furihata et al. sometimes use the term "multiple points linearly implicit schemes" to emphasise this fact.
The schemes are both second order in time, but in our example the linearly implicit multistep scheme has an error constant which is about 3-4 times larger than the fully implicit one-step scheme. In Figure 1.1 we plot the global error versus the number of linear solves for the two schemes. The linearly implicit scheme solves one linear system in each time step. The fully implicit scheme, on the other hand, solves a linear system for each Newton iteration which is repeated to machine precision in each time step. For the largest time step in our experiment this amounts to 561 linear solves per time step. The linear systems in each of the two cases have the same matrix structure, they are both penta-diagonal, and we therefore assume that the cost of solving the linear system is approximately the same for both methods. The x-axis in Figure 1.1 can thus be interpreted as a measure of the computational cost in each scheme. The plot shows that for a given global error the linearly implicit scheme is computationally cheaper than the fully implicit scheme.
There are situations in which the results from this example may be less relevant. For instance, the iteration method used in a fully implicit schemes may use approximate versions of the Jacobian for which faster solvers can be applied, therefore the cost of a linear solve may not be the same for the two types of schemes. For large time steps, both types of schemes are likely to encounter difficulties, but for slightly different reasons. The fully implicit scheme may experience slow or no convergence at all of the iteration scheme, whereas the linearly implicit scheme may become unstable for time steps over a certain threshold [9]. For instance, in the case of the stiff ODE considered by Gonzales and Simo [17] one will observe that the stability properties of the linearly implicit schemes will be completely lost whereas fully implicit conservative schemes behave remarkably well. For large-scale problems, one may have situations in which iterative linear solvers are required and where one cannot afford to solve these systems to machine accuracy, in such cases the linearly implicit schemes are less useful. In conlusion, we believe that which of the two types of schemes that is preferable depends on the PDE and the circumstances under which it is to be solved.
The two schemes used in the example above have slightly different conservation properties. The first one (1.2) conserves the exact Hamiltonian whereas the second scheme (1. 3) conserves what we will define as the polarised Hamiltonian Both of these functions are approximations to the true Hamiltonian, the first is a spatial approximation for a fixed time, and the second also includes an averaging over time. The intention is that in both cases one can see the methods as exactly preserving a slightly perturbed first integral over very long times. That this seems to work for the chosen example is clearly seen in the first plot in Figure  1.2 where we plot the error in H as a function of the solution obtained by the linearly implicit scheme (1.3). We integrate to t = 1000 and the error is plotted from t = 980 to t = 1000. Notice that in this example there is no drift in the energy error. The corresponding error plots for H as a function of the solution of (1.2) and H as a function of the solution of (1.3) are omitted since they are both preserved up to round-off error, and thus not that interesting. The second plot in Figure 1.2 shows how the error in H at the endpoint depends on ∆t. Empirically, we have the relation where C is a constant that depends on the solution, but not on n. See Section 6 for another example that tests the long time structure preserving properties of these schemes. This and similar examples show that there are situations where linearly implicit schemes can be a better choice than their fully implicit counterparts.  [11] and by Besse [3]. Examples of methods for other PDEs can be found in the monograph [15] and the papers [20] and [29].
In the next section we define the PDE framework that we use. Then, in Section 3 we consider discrete gradient methods and how they can be applied to PDEs. We study in particular the average vector field method by Quispel and McLaren [27] and the discrete variational derivative method by Furihata, Matsuo and coauthors [12,13,21,22,23,24]. We develop a framework that works for a rather general class of equations.
The key tools for developing linearly implicit methods for polynomial Hamiltonians are treated in Section 4, introducing the concept of polarisation. There is some freedom in this procedure, and we show through a rather general example term how the choice may significantly affect the stability of the scheme. We defer the introduction of spatial discretisation until Section 5. This is done mostly in order to keep a simpler notation, but also because our approach concerns conservative time discretisations and is essentially independent of the choice of spatial discretisation. The last section offers some more details on the procedure for constructing schemes and we give some indication through numerical tests on the long term behaviour of the schemes.

Notation and preliminaries
We consider integral preserving PDEs written in the form is the preserved quantity and D is a skew-symmetric operator that may depend on u. We write dx = dx 1 · · · dx d . We remark in passing that the class of PDEs which can be written in the form (2.1) contains the class of Hamiltonian PDEs, however we do not make the additional assumption that D satisfies the Jacobi identity [26]. By (u α J ) we mean u itself, which may be a vector u = (u α ) ∈ R m , and all its partial derivatives with respect to all independent variables, (x 1 , . . . , x d ), up to and including some degree ν. Thus, J is a multi-index, we let J = (j 1 , . . . , j r ), where r = |J| the number of components in J, and As in [26], the square brackets in (2.2) are used to indicate that a function depends also on the derivatives of its arguments with respect to the independent variables. In one dimension d = 1 and m = 1, for example, we have The variational derivative δH δu is an m-vector depending on u α J for |J| ≤ ν where ν ≥ ν. It may be defined through the relation [26, p. 245] for any sufficiently smooth m-vector of functions ϕ(x). One may calculate δH δu by applying the Euler operator to G[u], the α-component is given as so that the sum ranges over all J corresponding to derivatives u α J featuring in G. We have used total derivative operators, In parts of the paper we refer to Hamiltonians as polynomial, or specifically quadratic. By this we mean that H is of a form such that G is a multivariate polynomial in the indeterminates u α J , which in the quadratic case is of degree at most two. For example, the KdV equation (1.1) has a polynomial Hamiltonian of degree 3 In this case G = G(u, u x ) and thus m = d = ν = 1, and we get We always assume sufficient regularity in the solution and that the boundary conditions on Ω are such that the boundary terms vanish when doing integration by parts, for example periodic boundary conditions. The operator D should be skew-symmetric with respect to the L 2 inner product For the KdV case we simply have D = ∂ ∂x .
Furthermore, to be a true Hamiltonian system it should induce a Poisson bracket on the space of functionals as described e.g. in [26,Ch 7.1], meaning that the Jacobi identity must be satisfied. However, the approach presented here only requires D to be skew-symmetric so that the functional H is a conserved quantity. In the case that the PDE has more than one Hamiltonian formulation, we may make a choice of which of the integrals to preserve. Our approach does not in general allow for the preservation of more than one Hamiltonian at the same time, for this see the upcoming paper [10].
PDEs such as the wave equation are typically written with u tt appearing on the left hand side, in such cases we double the dimension of u in order to apply the stated framework. For complex equations one may do something similar, splitting either into a real and an imaginary part, or adding in the complex conjugate as a separate variable.

Discrete gradient and variational derivative methods
Discrete gradient methods for ODEs were introduced by Gonzalez [16]. See also [6], [7], [25], and [18,Chapter V.5]. Recently this idea has been applied to PDEs in the form of the average vector field (AVF) method [5] and in a somewhat more general setting, the discrete variational derivative (DVD) method. We recall the definition of a discrete gradient as presented for ODEs. If Since an ODE system preserving H can be written in the form for some skew-symmetric matrix S(y), one obtains a conservative method simply by defining approximations y n ≈ y(t n ) = y(t 0 + n∆t) through the formula whereS, typically allowed to depend on y n and y n+1 , is some skew-symmetric matrix approximating the original S. There are many possible choices of discrete gradients for a function H, see for instance [18,25]. A particular example is the one used in the AVF method defined as When applying this approach to PDEs the obvious strategy is to discretise the Hamiltonian H[u] in space, replacing each derivative by a suitable approximation like e.g. finite differences, to obtain a Hamiltonian H d (u) as for ODEs. Similarly, the skew-symmetric operator D is replaced by a skew-symmetric M × M -matrix D d to yield the scheme for advancing the numerical solution u n at time t n to u n+1 at time t n+1 . Examples are worked out for several PDEs in [5]. Furihata, Matsuo and coauthors present a whole framework for discretising PDEs in the variational setting in a series of papers, providing a discrete analogue of the continuous calculus, see for instance [12]. They discretise G to obtain G d using difference operators, and then the integral in H is approximated by a sum to yield H d . Then they derive a discrete counterpart to the variational derivative, and finally state the difference scheme in a form which is a perfect analogue to the Hamiltonian PDE system (2.1), letting The use of integration by parts in deriving the Euler operator is mimicked by similar summation by part formulas for the discrete case. Their discrete variational derivative is in fact rather similar to a discrete gradient, as it satisfies the relation for the discrete L 2 inner product.
In the present paper, we focus on the time dimension in most of what follows, thus we shall defer the steps in which H and thereby G are discretised in space. But (3.2) makes perfect sense after removing the subscript d, replacing u and v by functions u and v, and the discrete L 2 inner product by the continuous one. A discrete variational derivative (DVD) is here defined to be any continuous function δH δ(v,u) of (u (ν) , v (ν) ) satisfying The integrator yields a continuous function .
By combining (3.3) and (3.5) we see that the method preserves H. The AVF scheme can of course also be interpreted as a discrete variational derivative method where The fact that (3.6) verifies the condition (3.3) is seen from the elementary identity The derivative under the integral is written Now substitute this into (3.7) and interchange the integrals to obtain (3.3).
In most of the cited papers by Furihata, Matsuo and coauthors, the notion of a DVD method is less general than what we just presented, in the sense that the relation (3.2) is not actually used as the defining equation for a discrete variational derivative. Instead the authors present a relatively general format that can be used for discretising H, this format is depending on the class of PDEs under consideration, and they work out the explicit expression for a particular discrete variational derivative. To give an idea of how the format may look like, we briefly review some points from [12] where PDEs of the form (2.2) are considered with d = m = ν = 1 such that G = G(u, u x ). G is assumed to be written as a finite sum where f and g are differentiable functions. 1 The form (3.3) is then derived through followed by an integration by part on the second term. This technique can be extended in any number of ways to allow for more general classes of PDEs. For instance, one may allow for more factors in (3.8), like and repeated application of the formula ab − cd = a+c 2 (b − d) + b+d 2 (a − c) to this equation combined with integration by parts will result in a discrete variational derivative.
Schemes which are built on this particular type of discrete variational derivative will be called the Furihata methods in the sequel since it was first introduced in [12]. Matsuo et al. extend the method to complex equations in [22], while [13,21] derive methods for equations with second order time derivatives. Other papers using the discrete variational derivative approach include [14], [24], and [30].
The lack of a general formalism in the papers just mentioned, makes it somewhat difficult to compare the approach to the AVF method and characterise 1 In [12] the expression is discretised in space and g (ux) is replaced by a product g + (δ + k U k )g − (δ − k U k ) where δ + k and δ − k are forward and backward divided differences respectively.
in which cases they lead to the same scheme. Taking for instance the KdV equation (1.1) one easily finds that both approaches lead to the scheme (1.2), however, considering for instance the Hamiltonian one would obtain two different types of discrete variational derivative in the Furihata method and the AVF method, that is δH F δ(v,u) = δHAVF δ(v,u) . In some important cases, the Furihata method and the AVF method lead to the same scheme. . Ω ∂ J u · ∂ K u dx for multi-indices J and K, or Then the AVF and the Furihata methods yield the same scheme Proof. It suffices to check one general term of each type.
1. We find the variational derivative using (2.3) Inserting the variational derivative into (3.6) gives To find the discrete variational derivative of the Furihata method we compute After integration by parts we get and we see that

In this case we get
For the Furihata method one would here just compute and integration by parts yields 4 Linearly Implicit Difference Schemes

Polarisation
The key to constructing conservative linearly implicit schemes will be to portion out the nonlinearity over consecutive time steps. In effect, this means that we replace the original Hamiltonian H with an approximate one H. We shall call H a polarisation of H since its definition resembles the way an inner product is derived from a quadratic form. We shall see that the difference scheme resulting from such a polarised Hamiltonian will be a multistep method. This method will now preserve exactly H, as opposed to H for the methods in the previous section. The requirements on H are given in the following definition. In Section 4.2 we will discuss local order of consistency, it will then be convenient to make the stronger assumption that H and H are at least twice Fréchet differentiable. To distinguish from the weaker notion of variational (Gâteaux) derivative, we replace δ by ∂, noting that the first derivative in the two definitions are the same when they both exist. We then find from (4. which is used to compute

Polynomial Hamiltonians
The polarisation of polynomial Hamiltonians will be key to constructing linearly implicit schemes. We will now explain in detail how to do this, and we begin with an example term in the integrand G[u] = G(∂ α J u) depending on just one scalar indeterminate, namely G(z) = z p where z = ∂ J u α for some (J, α) and where p ≤ 4. This example is important not only as a simple illustration of the procedure, but also because terms of this type are common in many of the Hamiltonians found in physics. As we will see in the next section, it will be natural to use two arguments, k = 2, in the polarised Hamiltonian. In fact, we need to restrict ourself to cases with polynomial Hamiltonians for our technique to yield linearly implicit schemes. Then, by using k ≥ p/2 , we can obtain polarised Hamiltonians which are at most quadratic in each argument. We call these quadratic polarisations. We see that if k = 2 then cyclic is the same as symmetric G(u, v) = G(v, u), and the possible quadratic polarisations for p = 2, 3, 4 are respectively, Note that for these monomials both the third and fourth degree case are uniquely given, but the second degree case is not. In Section 4.3 we will consider how the choice of θ influences the stability of the scheme for a term which appears frequently in PDEs. We now consider the general case when G[u] is a multivariate polynomial in N ν variables of degree p. It suffices in fact to let G((u α J )) be a monomial since each term can be treated separately, for u ∈ R Nν . For a convenient notation, we rename the vector of indeterminates (u α J ) by using a single index i.e. u = (u 1 , . . . , u Nν ) and write (4.9) One may use the following procedure for obtaining a quadratic polarisation from (4.9) 1. Group the factors of the right hand side of (4.9) into pairs z r = u i2r−1 u i2r and if p is odd z k = u ip . Set Note that there are potentially many ways of ordering the factors in (4.9) which give rise to different polarisations.

Linearly Implicit Methods
We may now define the discrete variational derivative for this polarised Hamiltonian as a generalisation of (3.3) and (3.4). We let Our standard example will be a generalisation of the AVF discrete variational derivative, which we define as Here the variational derivative on the right hand side, δH δw1 is defined as before, considering H as a function of its first argument only, leaving the others fixed. Similar discrete variational derivatives could be derived in a number of different ways. In particular one finds that when the function H is quadratic in all its arguments, the approach used in deriving the Furihata methods would lead to a discrete variational derivative which is identical to that of the AVF-method. Now we define the polarised discrete variational derivative scheme and prove that, under some assumptions, this scheme is conservative, linearly implicit and has formal order of consistency two.  1) and (4.2), and suppose that approximations U j to u(j∆t, ·) are given for j = 0, . . . , k − 1.
• The polarised DVD (PDVD) scheme is given as , n ≥ 0. (4.13) • D is a skew-symmetric operator approximating D. In • If the discrete variational derivative is given by (4.12), then the scheme is called the polarised AVF (PAVF) scheme. for any polarised Hamiltonian function H.
Proof. By induction, this is an immediate consequence of (4.10) In a framework as general as the one presented here, it is not possible to present a general analysis for convergence or the order of the truncation error. However, it seems plausible that a necessary condition to obtain a prescribed order of convergence can be derived through a formal Taylor expansion of the local truncation error, we denote this the formal order of consistency. • The PAVF scheme has formal order of consistency one for any polarised Hamiltonian, and skew-symmetric operator D satisfying (4.14).
• If in addition (4.15) is satisfied, the scheme has formal order of consistency two.
Proof. We show that when the exact solution is substituted into (4.13) where the discrete variational derivative is given by (4.12), then the residual is O(∆t 2 ). Throughout the proof we assume the existence of Fréchet derivatives. Writing, for any j, u j = u(·, t j ) for the exact local solution at t = t j , we get for the left hand side Next we expand (4.12) to get δH AVF δ(u n , . . . , u n+k ) = 1 0 δH δw 1 (ξu n+k + (1 − ξ)u n , u n+1 , . . . , u n+k−1 ) dξ where u = (u n , . . . , u n ). Using first (4.4) and then (4.3), (4.5) we find Expanding D we get If the cyclicity condition (4.15) holds for D, we can simplify (4.18) to obtain By substituting into (4.13) the expressions (4.16), (4.17) and (4.19), all terms of zeroth and first order cancel and we are left with O (∆t 2 ) .
Theorem 4.5. Suppose that the polarised Hamiltonian H is a quadratic polynomial in each of its arguments, then the PAVF scheme is linearly implicit Proof. Since H is at most quadratic in the first argument, it follows from (2.5) that δH δw1 is of degree at most 1 in its first argument, and so we see from (4.12) that δH δ(U n , . . . , U n+k ) is linear in U n+k . Since D does not depend on U n+k we conclude that the scheme (4.13) is linearly implicit.
In some cases one wishes to have time-symmetric numerical schemes, see for example [18]. The numerical scheme (4.13) will in general not be symmetric, however it is not hard to modify the procedure to yield symmetric schemes. One needs to polarise H such that H is invariant also when the order of its arguments is reversed, it turns out that this can be achieved by symmetrising over the dihedral group rather than just the cyclic one. A similar adjustment must be made for D.
We remark that one can construct explicit schemes by using p time steps (as opposed to k) in H such that H becomes p-linear (as opposed to k-quadratic). The rest of the procedure for the explicit case is the same as for the linearly implicit case. Clearly, one expects that explicit schemes will have more severe stability restrictions than the linearly implicit ones.
Since these multistep schemes need the k previous values, it is not selfstarting. We have to provide the starting-values U 1 , . . . , U k in addition to the initial value U 0 . Usually these are computed using another sufficiently accurate conservative scheme, such as for example the AVF scheme. Another possibility is to use any integrator and integrate to machine precision.

Stability
In [9] we studied linearly implicit schemes for the cubic Schrödinger equation, and found that two-step schemes can develop a two-periodic instability in time.
We also saw that this can be remedied by choosing a different polarisation of the Hamiltonian.
As it turns out, a common case is when the Hamiltonian is a univariate polynomial of degree 4 or less. If we polarise this Hamiltonian using two time-steps, we get three linearly independent H, corresponding to (4.6), (4.7), and (4.8). The third and fourth degree Hamiltonians are uniquely given. However, in the second degree case we can choose θ ∈ [0, 1] such that the scheme becomes unstable. Since Hamiltonians of the type (4.6) appear in many important PDEs it may be useful to determine which θ ∈ [0, 1] lead to unstable schemes.
We choose to study the test equation with Hamiltonian and a skew-symmetric operator D which satisfies the eigenvalue equation for all integers k. The Airy equation Rewriting (4.6) gives And the numerical scheme is (4.20) Since this is a linear equation we can use von Neumann stability analysis [8].
We insert the ansatz U n (x) = ζ n e ikx to obtain the quadratic equation A necessary condition for stability is |ζ| ≤ 1 which implies Assuming that {λ k k 2 } k∈Z is unbounded, we must require that θ is chosen greater than or equal 1 2 . This is exactly the condition found in [9] for the cubic Schrödinger equation. When θ ≥ 1 2 the roots of (4.21) satisfy |ζ 1 | = |ζ 2 | = 1. In Figure 4.1 we solve the Airy equation with the scheme (4.20) using θ = 0.5 and θ = 0.49. We use the initial value u(x, 0) = sin(x), which, in the exact case, yields the traveling wave solution u(x, t) = sin(x + t). The θ = 0.49 solution blows up in few time steps, while the θ = 0.50 solution shows no signs of instability. Doing a discrete Fourier transform of the unstable solution we see that the instability starts at high frequencies, that is large k, which corresponds to the results shown above.
There might be cases where the scheme develop instabilities due to spurious modes no matter what polarisation one chooses. A full stability analysis of either the fully or linearly implicit schemes is to our knowledge not been done. Standard linearisation techniques will usually lead to the conclusion that the schemes are neutrally stable. The nonlinear effects, however small, may still cause the scheme to be unstable. The tests we have done on a wide range of PDEs seem to indicate that the stability usually is very good, future work may shed a light on this issue.

Space discretisation
Until now we have mostly considered the situation where the PDE is discretised in time while remaining continuous in space. The methodology developed in the previous sections apply equally well to systems of ODEs. Arguably, the most straightforward approach is simply to discretise the space derivatives in the Hamiltonian, for instance by finite differences. This leads to One also needs to replace the skew-symmetric operator D by a skew-symmetric matrix D d . The fully implicit method reviewed in Section 3 is then just the discrete gradient method (3.1), which conserves the discretised Hamiltonian H d (u) in every time step.
We consider now finite difference approximations. The function space to which the solution u belongs, is replaced by a finite dimensional space with functions on a grid indexed by I g ⊂ Z d . We use boldface symbols for these functions. Let there be N r grid points in the space direction r so that N = N 1 · · · N d is the total number of grid points. We denote by u α the approximation to u α on such a grid, and by u the vector consisting of (u 1 , . . . , u m ). We will replace each derivative u α J by a finite difference approximation δ J u α , and replace the integral by a quadrature rule. We then let Here ∆x is the volume (length, area) of a grid cell and b = (b i ) i∈Ig are the weights in the quadrature rule. The discretised G d has the same number of arguments as G, and each input argument as well as the output are vectors in R N . We have here approximated the function u α J by a difference approximation δ J u α , where δ J : R N → R N is a linear map. As in the continuous case, we use square brackets, say F [u], as shorthand for a list of arguments involving difference operators F [u] = F (u, δ J1 u, . . . , δ Jq u). We compute B is the diagonal linear map B = diag(b i ), i ∈ I g , and the discrete inner product used in (5.2) is Notice the resemblance between the operator acting on G d in (5.2) and the continuous Euler operator in (2.5). Alternatively, suppose that

Examples
To illustrate the procedures for constructing conservative schemes presented in this paper we consider as an example the generalised Korteweg-de Vries (gKdV) equation for an integer p ≥ 3, see for example [28]. The case p = 3 is the KdV equation (1.1), the case p = 4 is known as the modified KdV equation, and p = 6 is sometimes referred to as the mass critical generalised KdV equation. The gKdV can be written as (2.1) with The AVF discrete variational derivative (3.6) gives rise to the fully implicit scheme (3.5) After applying the polarising procedure of Section 4.1 to H = Ω G dx we get H = Ω G dx which depends on k = p/2 arguments wi , p odd, Notice that U n+k is indeed only appearing as linear terms in this scheme. The schemes (1.2) and (1.3) are found by setting p = 3 (k = 2) in (6.1) and (6.2), respectively. Following the procedure of Section 5 one can get a fully discretised scheme by replacing U by U and the first and third derivative operators by δ 1 and δ 3 respectively. In the Figures 6.1  We test the four methods on a traveling wave solution Φ(x − ct) = 3c 2 sech 2 3 √ c 2 (x − ct) , c > 0 using the parameters c = 1, x = (−5, 5), ∆x = 10 32 and ∆t = 0.1. As an indication of the long time behaviour of the presented schemes we consider to which extent the methods are able to preserve the shape and propagation speed of a traveling wave solution. We define the two quantities ε shape = min τ U n − Φ(· − τ ) 2 2 (6.5) and ε distance = |argmin τ U n − Φ(· − τ ) 2 2 − ct n |. (6.6) Thus ε shape measures the shape error of the numerical solution, and ε distance measures the error in the travelled distance of the numerical solution. We see in Figure 6.1 the fully implicit schemes preserves the shape better than the linearly implicit ones, and that the conservative schemes perform better than the non-conservative ones. In Figure 6.2 we see that the linearly implicit schemes have a more accurate phase speed than the fully implicit ones. Figure  6.3 shows the global error as a function of the time step. As expected the plot shows that the four methods are second order and that the linearly implicit schemes are inaccurate for large ∆t. In conclusion we see that the linearly implicit conservative scheme performs comparably to the other methods while being more efficient (the latter is shown in Figure 1