Symmetry, Integrability and Geometry: Methods and Applications Invariant Discretization Schemes Using Evolution–Projection Techniques ⋆

Finite difference discretization schemes preserving a subgroup of the maximal Lie invariance group of the one-dimensional linear heat equation are determined. These invariant schemes are constructed using the invariantization procedure for non-invariant schemes of the heat equation in computational coordinates. We propose a new methodology for handling moving discretization grids which are generally indispensable for invariant numerical schemes. The idea is to use the invariant grid equation, which determines the locations of the grid point at the next time level only for a single integration step and then to project the obtained solution to the regular grid using invariant interpolation schemes. This guarantees that the scheme is invariant and allows one to work on the simpler stationary grids. The discretization errors of the invariant schemes are established and their convergence rates are estimated. Numerical tests are carried out to shed some light on the numerical properties of invariant discretization schemes using the proposed evolution-projection strategy.


Introduction
Discretization schemes for differential equations that are not solely constructed for the sake of reducing the local discretization error as much as possible, but rather to preserve some of the intrinsic properties of these differential equations have become increasingly popular over the last decades. While preserving one of these properties, namely conservation laws, led to the design of conservative discretization schemes which are quite popular in the scientific community [4,19,24] (and in particular in the geosciences, e.g. [18,31]), there are other geometric features of differential equations that can be attempted to be preserved as well that have received less attention from the side of numerical analysis so far. One of these features are symmetries of differential equations. While there have been theoretical advancements on the methodologies of finding numerical schemes that preserve the maximal Lie invariance groups of systems of differential equations over the past 20 years or so [13,22,30,35], little is known about the numerical properties of these invariant schemes. A part of the problem is that while conservation laws are always properties of the solutions of a differential equation, symmetries are by definition properties of differential equations. Therefore, it is a standing question whether a discretization This paper is a contribution to the Special Issue "Symmetries of Differential Equations: Frames, Invariants and Applications". The full collection is available at http://www.emis.de/journals/SIGMA/SDE2012.html scheme that preserves numerically a property of a differential equation also improves the quality of the numerical solution of that discretized differential equation.
The present paper is devoted to an investigation of this question and related problems exemplified with invariant discretization schemes for the linear heat equation. The heat equation is particularly suited for this investigation as it is a canonical example in the group analysis of differential and difference equations. Moreover, there are already several studies devoted to invariant numerical schemes for this equation [1,14,35]. At the same time, in none of these existing works a deeper background analysis of the numerical properties (e.g. order of approximation or stability) of the developed schemes was investigated. A first account on numerical properties of invariant numerical schemes for the heat equation was given in [12], in which a numerical comparison of invariant and non-invariant schemes for the heat equation regarding accuracy was presented.
There are several reasons why less attention has been paid so far on the numerical analysis of invariant schemes (with the exception of the works [10,12]). One of the reasons is that the field of invariant discretization schemes is still in its early stages, with new conceptual algorithms being developed only recently [2,5,6,13,21,22,25,30]. Another reason is that invariant finite difference schemes generally require the use of adaptive moving meshes, i.e. it is necessary to include a non-trivial mesh equation in the discretization problem. Moving meshes lead to nonuniform grid point distributions and, in the multi-dimensional case, to non-orthogonal grids. The analysis of schemes on such meshes is considerably more difficult than that for related difference schemes on fixed, uniform and orthogonal meshes. Due to this second reason, most invariant numerical schemes so far have been constructed only for (1+1)-dimensional single evolution equations, as in that case moving meshes can be handled still with limited effort. Although we will be concerned with a (1+1)-dimensional equation in the present paper too, the methods used subsequently can be employed in the multi-dimensional case without substantial modification.
The new approach we propose here is to use the invariant grid equation only for a single time step and then to interpolate the solution to the regular grid. The important observation is that this interpolation can be done in an invariant way, i.e. projecting the solution does not break the invariance of the scheme itself. At the same time, the possibility to project the solution of an invariant scheme to a regular grid is highly desirable as in the multi-dimensional case a freely evolving grid can cause severe numerical problems. Moreover, for realistic numerical models, as e.g. employed in weather and climate predictions, it is in general hard to use adaptive meshes as the discretization of the governing equations is only one part of such model. Other parts are related to the numerical data assimilation, i.e. the preparation of the initial conditions for the numerical model and this step usually involves the forecasting model itself. As the assimilation of the initial conditions cannot be done on an evolving mesh (because the data are given at fixed locations only) this at once renders invariant schemes on moving meshes impractical. Equally important, any realistic numerical model for a nonlinear system of partial differential equations has to contain subgrid-scale models, which mimic the effects of processes taking place at those scales that the numerical model is not capable of resolving explicitly [33,34]. The construction of subgrid-scale models for non-resolved physical processes involves in general ad-hoc arguments and these arguments rely on the particular scale on which the unresolved processes take place. As a moving mesh locally changes the resolution and thus impacts what processes are explicitly resolved by a numerical model, subgrid-scale schemes have to be designed that can operate on grids with varying resolution. For realistic processes (which are usually not self-similar), this might be difficult to achieve in practice.
All what was said above objects against invariant numerical schemes for multi-dimensional systems of differential equations using freely evolving meshes. Thus, whether mathematically feasible or not, such schemes would be of less practical interest. This is why other approaches should be sought that on the one hand allow one to retain the invariance group of a system of differential equations in the course of discretization and on the other hand yield schemes that are practical to avoid the above mentioned and related problems. The proposed invariant evolution-projection strategy we are going to introduce below may be considered as one such approach.
The further organization of this article is as follows. The subsequent Section 2 features a summary and some extensions on the various methods to construct invariant discretization schemes. In Section 3 the heat equation along with its maximal Lie invariance group G is presented. It is discussed which subgroup G 1 of G we aim to present when constructing invariant numerical discretization schemes. The selection of G 1 is based on preserving the class of periodic boundary value problems we are focussing on. Section 4 contains the construction of an equivariant moving frame for G 1 along with a presentation of some lower order differential invariants of G 1 . In Section 5 invariant discretization schemes for the heat equation in computational coordinates are found. The local discretization errors of these schemes are established in Section 6. In Section 7 we introduce the new idea of invariant interpolation schemes that can be used to project the numerical solution obtained from an invariant scheme on a moving mesh to the regular grid. The numerical analysis as well as some numerical tests for the schemes proposed in this paper are found in Section 8. The summary of this article is presented in the final Section 9.

Construction of invariant discretization schemes
The construction of invariant discretization schemes for differential equations can be seen as a part of the ongoing effort to turn group analysis into an efficient tool for the analysis of difference equations, see e.g. the review article [25]. As of now, there are three main methods that are used to construct invariant discretization schemes.

Dif ference invariant method
The first method was developed by V. Dorodnitsyn, see [1,13,14,25,35]. It uses the infinitesimal generators of one-parameter symmetry groups of the system of differential equations under consideration that span the maximal Lie invariance algebra g of this system. These generators are of the form where x = (x 1 , . . . , x p ) and u = (u 1 , . . . , u q ) are the tuples of independent and dependent variables, respectively. Here and in the following, the summation convention over repeated indices is used. Rather than prolonging v to higher order derivatives of u with respect to x, which is standard in the symmetry analysis of differential equations [3,26,29], in this method the vector fields are prolonged to all the points of the discretization stencil, i.e. the collection of grid points which are necessary to approximate a given system of differential equation up to a desired order. This prolongation is of the form . , x p i ) and u i = (u 1 i , . . . , u q i ), i.e. it is done by evaluating the vector field v at all m stencil points z i = (x i , u i ) and summing up the result. An example for such a prolongation is given in Remark 1 in Section 5.
As a next step, the invariants of the group action are found by invoking the infinitesimal invariance criterion [13,26], which in the present case is pr v(I) = 0. The functions I that fulfill this condition for all v ∈ g are termed difference invariants.
Once the difference invariants on the stencil space are found, one then has to assemble these invariants together to a finite difference approximation of the given system of differential equations. By construction, this procedure guarantees that the resulting numerical scheme is invariant under the symmetry group of the original system of differential equations.
The main drawback of this method is that it might be hard to find a combination of difference invariants that approximates a system of differential equations in the multi-dimensional case. The problem is, as discussed in the introduction, that invariant schemes generally require the use of moving and/or non-orthogonal grids. Formulating consistent discretization schemes using difference invariants as building blocks on moving meshes is rather challenging in higher dimensions and thus limited the application of this method to the case of single (1 + 1)-dimensional evolution equations. We stress though that this problem only enters at the stage of combining difference invariants to a discretization scheme. Computing difference invariants in the multi-dimensional case can be done as effectively as computing differential invariants for multidimensional problems using infinitesimal techniques.

Invariant moving mesh method
Retaining the invariance of finite difference schemes under the maximal Lie invariance groups of physically relevant time-dependent differential equations often requires the use of moving meshes. This is true both for the finite difference method discussed in the previous Section 2.1 and the moving frame method to be discussed in the next Section 2.3. This kind of mesh adaptation in which the number of grid points remains constant throughout the integration is referred to as r-adaptivity in the field of adaptive numerical schemes [7,20].
The standard strategy to handle r-adaptive meshes is to regard the grid adaptation as a timedependent mapping from a fixed reference space of computational coordinates to the physical space of the independent variables of the differential equation, i.e. x = x(ξ) for ξ = (ξ 1 , . . . , ξ p ) being the computational variables. Without loss of generality, we assume that ξ 1 = τ = t is the time variable. The dependent variables u are expressed in the computational space by settinḡ u(ξ) = u(x(ξ)). For the sake of simplicity we will omit the bars henceforth.
The significance of the computational coordinates is to provide a reference frame that remains stationary and orthogonal even in the presence of grid adaptation in the physical space of coordinates. In the course of discretization the variable ξ labels the position of the grid points in the mesh and this labeling stays unchanged during the mesh adaptation. Thus, the computational variables can be interpreted physically as Lagrangian coordinates and their invariance under the motion of the grid is equivalent to the identity of fluid particles in ideal hydrodynamics.
Because by construction the grid remains orthogonal in the ξ-coordinates, the usual finite difference approximations for derivatives can be used in the space of computational variables. This simplifies both the practical implementation of the discretization method and the numerical analysis of the resulting schemes.
The expression of the initial physical system of differential equations in terms of computational variables leads to a system of equations that explicitly includes the mesh velocity x τ , which is yet to be determined in order to close the resulting numerical scheme. A prominent strategy for determining the location of the grid points at the subsequent time level in the onedimensional case uses the equidistribution principle, which in its differential form is (ρx ξ ) ξ = 0, where ρ is a monitor function that determines the areas of grid convergence and divergence. For higher-dimensional problems, equidistribution has to be combined with heuristic arguments, see [20] for more details.
The invariance of the initial differential equations is brought into the scheme by adequately specifying the monitor function ρ. In [6] it was pointed out that using monitor functions that preserving the scale-invariance of a differential equation is particularly relevant in cases where the equation is capable of developing a blow-up solution in finite time, see also [7,8,20]. This finding is generalized upon requiring that the monitor function is chosen in a manner such that the equidistribution principle is invariant under the same symmetry group as the original differential equation. For a number of symmetry groups this appears to be possible, see [2] for an example.
The invariant moving mesh method was recently extended in [2]. The idea of this extension is to transform the initial system of differential equations to the space of computational coordinates and to determine the form of the symmetry transformations in the computational space. The equations in the computational space are then discretized such that the resulting scheme mimics the transformation behavior of the continuous case. The main advantage of this approach is that it allows one to retain an initial conserved form of the system of differential equations and thus to numerically preserve certain conservation laws in the invariant scheme. This is relevant as preserving conservation laws in the course of invariant numerical modeling is yet a pristine problem. An exception to this is the discretization of equations that follow from variational principles, which, if done in a proper way, can lead to the simultaneous preservation of both symmetries and associated conservation laws, owing to the discrete Noether theorem. See, e.g. [5] for an example of such an invariant Lagrangian discretization.
Another advantage of the extension proposed in [2] is that it allows one to find invariant numerical schemes without the detour of difference invariants. This is essential as it can happen that the single equations in a system of differential equations cannot be approximated directly in terms on differential invariants but only in combination with other equations of that system. If this is the case it is not natural to attempt to discretize the system using difference invariants as this would lead to rather cumbersome discretization schemes.

Moving frame method
The third method is the most recent one [10,11,21,22,23,30]. It relies on the notion of equivariant moving frames and their important property to provide a mapping that allows one to associate an invariant function to any given function. As we will mostly work with this method in the present paper, we describe it in greater detail here. We collect some important notions on moving frames below, a more comprehensive exposition can be found in the original references [9,15,16,27,28,30]. Definition 1. Let G be a finite-dimensional Lie group acting on a manifold M . A (right) moving frame ρ is a smooth map ρ : M → G satisfying the equivariance property for all z ∈ M and g ∈ G.
Theorem 1. A moving frame exists in the neighborhood of a point z ∈ M if and only if the group G acts freely and regularly near z.
Local freeness of a group action means thatz = g · z = z for all z from a sufficiently small neighborhood of each point on M only holds for g being the identity transformation, which implies that all the group orbits have the same dimension. Here and throughout the paper, a tilde over a variable denotes the corresponding transformed form of that variable. Regularity of a group action requires that there exists a neighborhood for each point z ∈ M , which is intersected by the orbits of G into a pathwise connected subset.
When a group G does not act freely on M , its action can be made free upon extending it to a suitably high-order jet space J n = J n (M, p) of M , 0 ≤ n ≤ ∞. Locally, the nth order jet space of a p-dimensional submanifold of M has coordinates z (n) = (x, u (n) ), where as in the previous subsections x = (x 1 , . . . , x p ) are considered as the independent variables, u = (u 1 , . . . , u q ), q = dim M − p, are the dependent variables and u (n) collects all the derivatives of u with respect to x of order not greater than n including u as the zeroth order derivatives. In practice, the prolongation of the group action of G on J n is implemented using the chain rule.
Moving frames are determined using a normalization procedure. The steps to find a moving frame for a group action G are the following: (i) Define a cross-section to the group orbits. A cross-section C is any submanifold C ⊂ M of complementary dimension to the dimension r of the group orbits, i.e. dim C = dim M −r, that intersects each group orbit once and transversally. Usually coordinate cross-sections are chosen in which some of the coordinates of M (or of J n if the group action is not free on M ) are set to constants, i.e. z 1 = c 1 , . . . , z r = c r . (ii) The algebraic systemz 1 = (g · z) 1 = c 1 , . . . ,z r = (g · z) r = c r is solved for the group parameters g = (g 1 , . . . , g r ). The resulting expression g = ρ(z) is the moving frame.
Moving frames can be used to map any given function to an invariant function by a procedure called invariantization.
That the function ι(f ) constructed in this way is indeed invariant follows from the equivariance property (1) of the moving frame ρ, which boils down to the definition of an invariant function I, i.e. I(g · z) = I(z). In practice, a function f (z) is invariantized by first transforming its argument using the transformations from G and then substituting the moving frame for the group parameters. By definition, an invariant that is defined on the jet space J n is a differential invariant.
Moving frames can also be constructed on a discrete space. In a finite difference approximation, derivatives of functions are approximated using a finite set of values of these functions, and all the points needed to approximate the derivatives arising in a system of differential equations are the points of the stencil introduced in Section 2.1. Because most of the interesting symmetries of differential equations that are broken in standard numerical schemes require the use of non-orthogonal discretization meshes, it is beneficial to both regard x and u as the dependent variables and the computational variables ξ as the independent variables as was discussed in the previous Section 2.2.
Sampling the tuples from the extended computational space , which can be identified as the joint product space of stencil variables. Because the identifier ξ i of the point w i is required to be unique, each element of M n ξ only includes distinct grid points in the physical space of equation variables. The dimension of the space M n ξ depends on the number of independent and dependent variables in the system of differential equations and the desired order of accuracy of the approximated derivatives.
It is possible to carry out the construction of the moving frame on M n ξ , i.e. to define the moving frame by an equivariant mapping ρ n ξ : M n ξ → G, where G acts on M n ξ by the product action, g · (w 1 , . . . , w n ) = (g · w 1 , . . . , g · w n ). Note that the extension of the group action to the computational variables ξ is trivial, i.e. they remain unaffected by G,ξ = g · ξ = ξ, see [2]. The compatibility between the moving frame ρ n ξ and the moving frame ρ on the space M (or an appropriate jet space J n ), i.e. that ρ n ξ → ρ in the continuous limit is assured provided that the cross-section defining the moving frame ρ n ξ in the continuous limit converges to the crosssection defining the moving frame ρ. Once the moving frame is constructed on the discrete space M n ξ of stencil variables, it can be used to invariantize any numerical scheme expressed in computational coordinates. This will be explicitly shown in Sections 5 and 6 where we will construct invariant schemes for the heat equation. It is essential that the construction of the moving frame on the grid point space is carried out in terms of computational coordinates rather than physical coordinates. This can be illustrated by the following simple example. Example 1. The Laplace equation u xx +u yy = 0 is, inter alia, invariant under the one-parameter group of rotations SO(2),x = x cos ε − y sin ε,ỹ = x sin ε + y cos ε. Let us obtain the moving frame ρ for this group action from the normalization condition u x = 0, i.e. we determine the moving frame on the first jet space J 1 (M, 2), ρ = ρ(x, y, u, u x , u y ). Prolonging the transformations from SO(2) to the derivative u x leads toũx = u x cos ε − u y sin ε and thus the moving frame is ε = arctan(u x /u y ).
Let us now find the product frame from the discrete normalization condition u d cannot be solved for the group parameter ε. On the other hand, setting u = u(x(ξ 1 , ξ 2 ), y(ξ 1 , ξ 2 )) and expressing u d x in terms of the computational variables ξ 1 , ξ 2 , the normalization u d This expression can be solved for ε and it gives ε = arctan which in the continuous limit goes to ε = arctan(u x /u y ) as required.

Lie symmetries of the heat equation
The one-dimensional linear heat transport equation is where we scaled the thermal diffusivity ν to 1, i.e. equation (2) is in non-dimensional form. The heat equation (2) admits the following infinitesimal generators of one-parameter groups, which generate the maximal Lie invariance algebra g of equation (2): where α runs through the set of solutions of equation (2), see e.g. [26]. These vector fields generate (i) time-translations, (ii) space translations, (iii) scalings in u, (iv) simultaneous scalings in t and x, (v) Galilean boosts, (vi) inversions and (vii) the superposition principle symmetry.
In this paper, we will construct invariant numerical schemes for a class of initial value problems of the heat equation using periodic boundary conditions. This class of initial-boundary value problems only admits a subgroup of the symmetry group of the heat equation as inversions are no longer admitted; inversions do not send an initial value problem from the considered class to another initial value problem. The symmetries associated with the first five vector fields in (3) are compatible with the class of initial-boundary value problems we are interested in, i.e. they map the class of initial-boundary value problems for the heat equation under consideration to itself. This re-interpretation of symmetries of differential equations without initial and boundary conditions as equivalence transformations for a class of initial-boundary value problems was recently pointed out in [2].
In what follows we will thus focus our attention on constructing numerical schemes that preserve the symmetries generated by the first five operators of (3). The associated subalgebra of g will be denoted by g 1 . We do not require to preserve the linearity operator here by construction. At the same time, as will be shown in Section 8 the numerical schemes we propose in this paper preserve the linearity property up to the discretization error expected.

Moving frame and dif ferential invariants for the heat equation
We determine the moving frame for the subgroup G 1 of transformations associated with the subalgebra g 1 . Transformations of G 1 are of the form Because the group action of G 1 is not free on M = {(t, x, u)} we construct the moving frame on a suitably high-order jet space. In the present case, the group action of G 1 becomes free on J 1 = J 1 (M, 2). Thus, it is necessary to extend the transformations (4) to derivatives of u with respect to t and x.
Using the chain rule we can compute the transformed operators of total differentiation, which read as where D x = ∂ x + u x ∂ u + u tx ∂ ut + u xx ∂ ux + · · · and D t = ∂ t + u t ∂ u + u tt ∂ ut + u tx ∂ ux + · · · denote the usual operators of total differentiation. With the transformed total differentiation operators at hand it is possible to compute the transformed partial derivatives of u with respect to t and x. The transformation rules for the lowest order derivatives arẽ In fact, for the construction of the moving frame already the knowledge of the first order derivatives is sufficient. We compute the moving frame for the five-parameter group of transformations of the form (4) using the following normalization conditions which determine a valid cross-section to the group orbits of the prolonged action of G 1 on J 1 , The moving frame is computed by taking the transformed form of the normalization conditions, t = 0,x = 0,ũ = 1,ũt = 1 andũx = 0 and by solving the resulting algebraic system for the group parameters ε 1 , . . . , ε 5 . The result of this computation is the following moving frame g = ρ(z (1) ), With the moving frame at hand, we can invariantize any of the partial derivatives of u with respect to t and x and thus obtain a complete set of differential invariants for the subgroup G 1 of the maximal Lie invariance group of the heat equation. As an example, invariantizing the derivative u xx , i.e. computing ι(u xx ) as (g · u xx )| g=ρ(z (1) ) we produce the differential invariant Invariantizing the heat equation, i.e. computing ι(u t − u xx ) = 0 and recalling that ι(u t ) = 1, we obtain which yields the original heat equation expressed in terms of differential invariants. This reexpression of a differential equation using the differential invariants of its symmetry group is known as the replacement theorem [9].

Invariant discretization of the heat equation
The invariant discretization of equation (2) cannot be done on a fixed, uniform grid. To see this, let us check the transformation behavior of the grid equation x n+1 i − x n i = 0, which is the definition of a stationary grid, under the transformations (4). This yields which is only zero in the case when ε 5 = 0. Stated in another way, a discretization on a fixed grid can at most preserve the symmetry subgroup of G, which is generated by the first four elements of the maximal Lie invariance algebra g of the heat equation (3). Thus, the discretization of (2) preserving G 1 will require the use of moving grids. For this reason it is convenient to express (2) in terms of computational coordinates initially, i.e. we set u(τ, ξ) = u(τ, x(τ, ξ)), where ξ is the single spatial computational variable and τ = t. The heat equation in this set of coordinates reads So as to find the invariant discretization of the heat equation in the form (7), we determine the moving frame in the space of stencil variables M 4 ξ using the discrete analogs of the normalization conditions (5) expressed in terms of computational coordinates.
For the sake of convenience we introduce the notation h The discretization stencil we use is depicted in Fig. 1.
The appropriate normalization conditions for a compatible moving frame ρ 4 ξ are where  in the above normalization conditions by their respective transformed expressions and solving the resulting algebraic system for the group parameters we obtain the following moving frame on the space of stencil variables M 4 ξ , where we introduced (ln u) d x = (ln u n i+1 − ln u n i−1 )/(h + + h − ). This moving frame is compatible with the moving frame (6) in that it converges to (6) in the continuous limit ∆ξ → 0 and ∆τ → 0 upon using The moving frame (9) can now be used to invariantize any non-invariant finite difference discretization of (7) on M 4 ξ . To illustrate this, we invariantize the standard FTCS (forward in time centered in space) scheme This is done by first replacing all terms by their respective transformed expressions and substituting the moving frame (9) for the arising group parameters. The result of this procedure is the invariant scheme Again, it can be checked that the above scheme (11) indeed converges to (7) in the limit of ∆ξ → 0 and ∆τ → 0. This will be shown explicitly in Section 6, where we establish the order of approximation of (11).
So as to complete the scheme (11) it is necessary to determine x n+1 i , which is the ingredient missing in (11). There are different ways to determine a grid equation, such as using the equidistribution principle as outlined in Section 2.2. The problem with this strategy in the present case is that while it might be beneficial from the numerical point of view, it might not be easy to obtain an invariant discretization of this principle which does not lead to a fully coupled equation-grid system. In other words, it can happen that the grid equation includes values of u at both t n and t n+1 . While this coupling is not a problem in the one-dimensional case, it can lead to a severe restriction of the applicability for multi-dimensional equations as solving the coupled equation-grid system might then be too expensive.
A G 1 -invariant grid equation that circumvents the aforementioned coupling problem can be derived from the invariantization of x n+1 i . This invariantization yields where we did not explicitly substitute the frame value for ε 4 . An appropriate grid is then given This grid equation is quite similar to the grid equation which was found in [1,14] using the method of difference invariants. This last grid (13) is not only invariant under the subgroup G 1 but under the whole maximal Lie invariance group G of the heat equation. In the continuous limit, both equation (12) and (13) converge to We have tested all our numerical schemes with both (12) and (13) and found that the resulting schemes give asymptotically the same numerical results. In fact, as in the evolution-projection strategy that will be introduced in Section 7 we have h + = h − = h, equation (12) and (13) coincide.

Remark 1.
While the invariantization algorithm guarantees that the scheme (11) is indeed invariant under the subgroup G 1 of the maximal Lie invariance group of the heat equation, the invariance can be checked in a straightforward fashion using the infinitesimal invariance criterion as invoked in the Dorodnitsyn method discussed in Section 2.1. Let us recall that this criterion states that an invariant I of a group action is annihilated by the associated infinitesimal generators, i.e. v(I) = 0 for all v ∈ g. Because in the present case, the invariants are defined on the stencil space with coordinates τ n , ∆τ , x n i , x n i+1 , x n i−1 , x n+1 i , u n i , u n i+1 , u n i−1 and u n+1 i , we have to prolong the operators of g accordingly. The prolongations of the first five operators of (3) to the variables of the stencil are i , see [13,25] for more details. It can be checked that pr v(S) = 0 and pr v(M ) = 0 hold on S = 0 and M = 0 for all the prolonged infinitesimal generators and thus S = 0 is a proper invariant numerical scheme and M = 0 an invariant grid equation.

Remark 2.
The heat equation is a linear partial differential equation in two independent variables. One might thus consider to set up a grid equation not only for spatial but also for temporal adaptation. The reason why we refrain from spatial-temporal adaptation here is that the symmetry group G of the heat equation is compatible with flat time layers, i.e. t n i+1 − t n i = 0 is a G-invariant equation. Improving the invariant numerical scheme constructed above using temporal adaptation would thus not allow one a fair comparison against the original noninvariant FTCS scheme for the heat equation. Moreover, flat time layers are well-agreed with the physics of the heat transfer problem, which affects all points of the domain simultaneously.

Numerical properties of the invariant scheme
In this section we investigate the numerical properties of the scheme (11) and related schemes. We start our consideration with the estimation of the local truncation error of the scheme. The study of this question is relevant because so far little is known about the relation between the order of a non-invariant scheme and its invariantized counterpart. The discretization of the heat equation in computational coordinates (7) can be formally represented as where in the present case we assume that derivatives are approximated with the aid of a standard FTCS scheme. More general schemes will be considered after the order of the invariantized FTCS scheme is established.
Theorem 2. The order of the invariant scheme (11) is the same as the order of the scheme (14), namely first order in time and second order in space, provided that an Euler forward step and second order centered differences are used to approximate the time and space derivatives arising in both the differential equation (14) and the normalization conditions (8).
Proof . Invariantizing the scheme (14) using the normalization conditions (8) leads to where ι(f )(z) denotes the invariantization of the function f (z). By definition, invariantization of a function f (z) means to transform the argument z and plug in the moving frame for the group parameters. In the present case, the transformed form of (15) can be written as Using the normalization condition u n i = 1 we obtain that u n i = 1 = e ε 3 −ε 5 x−ε 2 5 t u n i and thus the last expression can be recast as Let us now determine the local discretization error in the parameter ε 5 . The respective moving frame component is which upon using (10) expands to Substituting ε 5 into the second term of equation (16) and expanding the exponential functions in the same term into Taylor series, we obtain after some rearranging This can be further simplified to It now remains to expand the first term in equation (18). The moving frame component for ε 4 in (9) can be recast as Using x n+1 i = x n i + x τ ∆τ + O(∆τ 2 ) and u n+1 i = u n i + u τ ∆τ + O(∆τ 2 ) and again expanding the exponential function into a Taylor series, we derive Plugging this into equation (18) we arrive at which completes the proof of the theorem.
A more general statement is the following one: Theorem 3. The order of spatial discretization of an invariant finite difference scheme for the heat equation in computational variables equals the order p ∈ N of the spatial discretization of the associated non-invariant finite difference scheme provided that centered differences of order p are used to approximate both the derivatives in the heat equation and in the normalization conditions (8).
Proof . In view of the general form (15) of the invariantization of scheme (14), we study the invariantization of the terms x ξ and u ξξ . The invariantization of (x d ξ ) 2 = (x ξ + O(∆ξ p )) 2 is ι((x d ξ ) 2 ) = e 2ε 4 (x 2 ξ + O(∆ξ p )) and is of the same order p if, as required, the moving frame component ε 4 stems from the approximation of u d τ = 1 using pth order accuracy and thus only includes approximations of derivatives with that accuracy.
See [17] for a discussion of the algorithm for finding the weights c k p,j in higher-order finite difference approximations of the kth derivative of u. The invariantization of u ξξ is ι(u ξξ ) = 1 ∆ξ 2 p/2 j=−p/2 c 2 p,j exp ε 3 − ε 5 x n j − ε 2 5 τ n u n j , or, upon using the normalization condition u n i = 1, where we expand Using the expressions for ∆x j and u n i , the expression (19) can be expanded and rearranged in powers of j∆ξ in the form The expressions for A k , k = 2, are not required subsequently. The proof is completed upon substituting for ε 5 the corresponding moving frame component (which is of order p if the normalization u d x = 0 is approximated with pth order accuracy) and by noting that where the precise values of the constants c k follow from evaluating the respective sums.
The scheme (11) is only of first order in time τ = t. To construct a scheme that is second order in time, we can start with a non-invariant scheme (14) and discretize the time derivative u d τ with second order accuracy, i.e. we set u d is the value of u at the previous time step τ n−1 = τ n − ∆τ . It is now necessary to check whether invariantizing this leapfrog discretization leads to an invariant scheme that is also second order in time.
Theorem 4. Invariantization of the scheme (14) in which a leapfrog step and second order centered differences are used to approximate the time and space derivatives, leads to an invariant scheme that is both second order in time and space provided that the normalization conditions (8) are approximated using discretizations that are of second order.
Proof . To prove this theorem it is sufficient to establish the order of the first term in equation (18). We proceed in an analog manner as in the proof of Theorem 2, i.e. we discretize the normalization condition u d τ = 1 but now with second order accuracy. This yields Using the normalization condition u n i = 1 as before and expanding the exponential functions we derive where we have substituted the expression (17) for ε 5 . Plugging this result into equation (18) completes the proof of the theorem.
The actual form of the resulting invariant leapfrog scheme is )/∆τ . Higher order in time schemes can be constructed upon invariantizing multi-stage schemes. Combining this result with the result established in Theorem 3 we have found the following: Corollary 1. Invariantizing a non-invariant finite difference scheme for the heat equation in computational coordinates preserves the spatial and temporal order of the initial non-invariant finite difference scheme provided that centered differences are used and the normalization conditions for the moving frame are discretized with the same order as the respective derivatives in the non-invariant finite difference scheme.

Invariant interpolation schemes
A common property of invariant numerical schemes for evolution equations possessing a nontrivial maximal Lie invariance group is that it is not possible to use a fixed, orthogonal discretization mesh. The continuous evolution of the mesh, if not handled properly, can lead to several undesirable properties, such as an overly strong concentration of grid points in certain regions and therefore too poor a resolution in other parts of the integration domain. The problem gets worse in the multi-dimensional case where mesh tangling or strongly skewed meshes can occur. But even if the mesh movement can be managed in an optimal way there are various physical problems for which continuously adapting grids pose a severe challenge. An example for this are practically all models that are in operational use in weather and climate prediction. These models employ sophisticated data assimilation strategies and are coupled to subgrid-scale parameterizations that aim to mimic the effects of unresolved processes on the grid scale variables. Attempting to make use of data assimilation or parameterization schemes on moving meshes is not only a technical problem that would cause a significant computational overhead compared to standard schemes but also a conceptual challenge for it is unclear on how to design parameterization schemes that can operate on grids with varying resolution. In order to promote the ideas of invariant numerical discretization schemes beyond their application to simple evolution equations it is thus instructive to study possible ways of overcoming the limitations imposed by the requirement of using moving meshes.
One straightforward idea is to use invariant schemes on fixed (i.e. non-invariant grids). As was shown in [30] this can lead to improved numerical solutions compared to non-invariant integrators, while still being excelled by the results that can be obtained using completely invariant schemes. On the other hand, if moving (invariant) meshes are not tractable for a particular class of problems, preserving the invariance of a system of differential equations at least for the discretization of the system itself might be a possible trade-off to take.
Another idea is to use an evolution-projection strategy, which will be proposed in the following. This concept relies on using the invariant scheme for the system of differential equations together with the invariant mesh equations for a single time step followed by the projection of the numerical solution back to the regular mesh. A similar strategy has proven successful in semi-Lagrangian time integration schemes [32].
In the present case, the projection step can be practically realized by using interpolation. Obviously, any standard interpolation scheme can be used to map the numerical solution u n+1 i defined at x n+1 i to the uniformly spaced ξ-grid. This, however, can break the invariance of the numerical scheme as a whole and so the question arises whether it is possible to accomplish the interpolation step in a symmetry-preserving fashion.
In the following we discuss two possible ways of formulating invariant interpolation schemes, both of which can be used for finding interpolations that allow the re-mapping of the numerical solution on a moving mesh to a fixed, Cartesian, equally-spaced grid. These ways are the invariantization of non-invariant interpolation schemes with moving frames and the construction of interpolations using difference invariants.
Invariantization of interpolation schemes. The moving frame constructed in the course of invariantizing a finite difference scheme can also be used to invariantize a certain interpolation method. We exemplify this idea by invariantizing the formula for linear interpolation, ]. The invariantization of this expression using a moving frame associated with G 1 yields the invariant interpolation formula Note that we have used a slightly different moving frame for the invarianization as we have used for invariantizing the finite difference scheme for the heat equation. Specifically, this moving frame is constructed by replacing the normalization condition u d x = 0 withû d x = 0. The reason for this is that the moving frame used earlier yielded ε 5 = (ln u x ) d , i.e. it involves the solutions of u i at the time step τ n rather than at τ n + ∆τ . Irrespectively of what normalization is used, both interpolations are invariant. Setting y = ξ i in the above interpolation formula yields u n+1 i on the regular computational grid. Note that the interpolation (20) In a similar manner more sophisticated interpolation schemes can be invariantized. In the following, we will use the invariantization of quadratic interpolation. Usual quadratic interpolation is based on the expression where y ∈ [x n+1 i−1 , x n+1 i+1 ] and L j (y) are the Lagrangian interpolation polynomials. Invariantizing this formula using the same moving frame as above we get where , as in the case of the invariant linear interpolation (20). Numerical examples using the invariant quadratic interpolation will be given in Section 8.
Interpolation using dif ference invariants. The product frame on the grid point space allows invariantizing the elementary variables x n i and u n i , which yields the system of joint invariants. In the continuous limit these invariants take the normalization values chosen for x and u to construct the usual moving frame ρ [27]. On the other hand, on the discrete space M n ξ we only normalize one x n i (i, n fixed) among all the grid points x k l and the analog statement is true for the associated values u k l . This means that the joint invariants ι(x k l ) and ι(u k l ), l = i, k = n, are nontrivial and can be used to assemble invariant interpolation schemes.
In the present case, while we have normalized u n i = 1 in the course of constructing the moving frame ρ 4 , we are free to use the moving frame to invariantized any u k l where l = i, k = n and this will yield a proper (nontrivial) invariant on the discrete space M 4 ξ . As above, we again recompute the moving frame for G 1 by replacing the normalization conditions u n i = 1 and u d x = 0 with u n+1 i = 0 andû d x = 0, respectively, which yields new expressions for the moving frame components of ε 3 and ε 5 given by Using this modified moving frame, we then invarianize the variable u n+1 (ξ i ), which is the sought value of u at the point (τ n+1 , ξ i ) of the computational domain. This invariantization yields Because in the continuous limit the invariantization of u n+1 (ξ i ) must reproduce the normalization condition u = 1, we restrict the difference invariant to the manifold ι(u n+1 (ξ i )) = 1. The invariant interpolation is thus and it is again consistent as u n+1 (x n+1 i ) = u n+1 i . More accurate interpolations could be constructed by combining the invariants ι(x k l ) and ι(u k l ) in a suitable way. The advantage of the interpolation methods introduced in this section is that they are invariant under the group G 1 , i.e. using these interpolation formulas to map u n+1 i back to the ξ-grid does not break the invariance of the numerical schemes for the heat equation, while still allowing one to use a regular grid. Invariant interpolations thus allow avoiding the complications that moving meshes impose on the applicability of symmetry-preserving finite difference discretizations.

Numerical verif ication
In order to verify the accuracy predicted above for the various schemes proposed, we set-up the following problem. On a periodic domain x ∈ [0, 2π[, consider On a sequence of grids with N ∈ {2, 4, 8, . . . , 256}, the number of grid points, we compute the error in the maximum norm between the numerical solution and the exact solution at t = 1. The time step ∆τ is taken as proportional to h 2 , h = h + = h − , in all simulations.
In each of the following figures, we plot the reference line corresponding to O h 2 dashdotted, and the L ∞ error as black line where, Note that for all approaches described below we expect a second order convergence since the numerical scheme being used is of second order, its invariantization was shown to preserve this order and the quadratic interpolation and its invariantization is also of second order.

Invariant scheme without projection
In this test run we use the scheme (11) without projection. As a result, the solution is evolving along the trajectories of the grid equation (12). Since the spacing between trajectories is not constant, i.e. h + and h − are changing in time, we choose to plot the error versus 1/N . In Fig. 2, we observe the second order convergence expected.

Invariant scheme with non-invariant quadratic interpolation
In this scheme we interpolate the solution of the invariant scheme (11) at every step back onto the regular grid using standard Lagrange quadratic interpolation (21). As a result, we have the solution on a regular grid with step-size h = 1/N . In Fig. 3, we observe the second order convergence expected.

Invariant scheme with invariant quadratic interpolation
In this scheme, we interpolate the solution of the invariant scheme (11) at every step back onto the regular grid using the invariant Lagrange quadratic interpolation (22) described in the previous section. As for the case above, at each time step we have the solution on a regular grid with step-size h = 1/N . In Fig. 4, we observe the second order convergence expected.

Linearity preservation in the invariant numerical scheme
Linearity is not preserved by construction in the schemes proposed. For this reason it is instructive to check numerically whether or not linearity is preserved in the fully invariant scheme.
Consider the initial value problem, u t = u xx u(x, t = 0) = (sin(x − 1) + 2) + (cos x + 2), the solution of which we call u exact . We then solve numerically the following two equations u a t = u a xx with u a (x, 0) = sin(x − 1) + 2, u b t = u b xx with u b (x, 0) = cos x + 2, and define u s = u a + u b . Fig. 5 depicts the L ∞ error between u exact and u s . We observe a convergence rate of second order. In other words, despite the fully invariant numerical scheme does not explicitly preserve the symmetry associated with the linear superposition principle, we observe that the linearity property is preserved approximately to the order of the method.

Conclusions
In this paper we construct invariant discretization schemes using the method of invariantization via equivariant moving frames. The advantage of this technique is that it allows one to start with a given non-invariant scheme and convert this initial scheme into a finite difference approximation of a system of differential equations L that is invariant under the same maximal Lie invariance group G (or a suitably chosen subgroup of G) as admitted by L.
The possibility of converting non-invariant numerical schemes into invariant discretizations may lead to the overly optimistic speculation that the schemes constructed by invariantization could be easily included in existing numerical models using the original scheme. The hurdle preventing this in practice is that preserving symmetry groups of systems of evolution equations more complicated than scalings or translations requires the use of moving grids. Converting numerical models that use standard discretization schemes based on fixed lattices to (invariant) discretization schemes on moving meshes is not an easy task. At the same time, rewriting numerical models from scratch for the simulation of involved physical processes using symmetrypreserving schemes might be a time-consuming and costly task too and it not certain that this is feasible at all. Moreover, it is as of now unclear whether preserving symmetries in numerical schemes for multi-dimensional systems of partial differential equations gives enough added value compared to standard schemes that one might justify such an undertaking in practice. This is why one relies on finding methods allowing one to efficiently include invariant discretization schemes into existing numerical models without the need to rewrite new models from scratch that incorporate the invariance methodology. The method proposed in this article solves this problem by breaking the integration procedure into two steps, the time-stepping using the invariant numerical scheme with an invariant numerical grid equation and the projection (interpolation) of the results obtained at intermediate grid points to the regular mesh. This interpolation can be done in an invariant way by applying the moving frame map used to invariantize the initial discretization scheme also to a particular interpolation method. An alternative is to assemble the invariant interpolation method using joint invariants. Either way, it is worthwhile pointing out that interpolations requiring only data given at a single time level are already invariant under most symmetry groups as admitted by physical systems of differential equations. Thus, invariantization of interpolation formulas will often only lead to minor modifications of the initial interpolation method chosen and the influence on the numerical solution might be rather small. In the numerical tests carried out above for the heat equation, the difference in the convergence properties we found when using invariant or non-invariant interpolation methods is indeed small although using the invariant interpolation gave slightly better numerical results. This is encouraging and the reason why we plan to further investigate invariant numerical schemes using the projection procedure.
We illustrate the evolution-projection strategy by integrating the one-dimensional linear heat equation with an invariant numerical scheme. The heat equation has been studied quite extensively in light of its invariance properties and in particular it is a standard model for the construction of invariant numerical schemes [1,14,35]. At the same time, a comprehensive numerical analysis of such schemes was not given before and thus seems relevant to be reported. This is another aim of the present paper. Again, the analysis of numerical properties of discretization schemes is considerably easier if one can use non-evolving meshes.
Further work we intend to do is to employ the evolution-projection strategy to multidimensional systems of differential equations using both higher-order discretization and interpolation schemes.