Symmetry-Preserving Numerical Schemes

Alexander Bihlo Francis Valiquette
Department of Mathematics and Statistics Department of Mathematics
Memorial University SUNY New Paltz
St. John’s, NL, Canada A1C 5S7 New Paltz, NY, USA  12561

Keywords: Differential equations, equivariant moving frames, finite difference equations, infinitesimal symmetry generators.

Mathematics subject classification: 39A99, 54H15, 65Q10

In these lectures we review two procedures for constructing finite difference numerical schemes that preserve symmetries of differential equations. The first approach is based on Lie’s infinitesimal symmetry generators, while the second method uses the novel theory of equivariant moving frames. The advantages of both techniques are discussed and illustrated with the Schwarzian differential equation, the Korteweg–de Vries equation and Burgers’ equation. Numerical simulations are presented and innovative techniques for obtaining better invariant numerical schemes are introduced. New research directions and open problems are indicated at the end of these notes.

1 Introduction

The aim of geometric numerical integration is to develop numerical integrators that preserve geometric properties of the system of differential equations under investigation. Classical examples include symplectic integrators, [38, 55], energy preserving methods, [77], and schemes that preserve a Lie–Poisson structure, [88]. The motivation behind geometric numerical integration is that, as a rule of thumb, such integrators will typically give better global or long term numerical results than standard methods since they incorporate qualitative properties of the system under consideration.

In mathematical physics, most fundamental differential equations are invariant under a certain collection of symmetry transformations. These symmetries can be point transformations, contact transformations, or generalized transformations, [68]. In all cases, the symmetries of a differential equation encapsulate important properties of the equation and its solutions. Furthermore, Lie group techniques are amongst the most effective methods for obtaining explicit solutions and conservation laws of nonlinear differential equations, [13, 68, 76].

When discretizing differential equations invariant under a certain symmetry group, there are different incentives for preserving the symmetries of these equations. From a physical standpoint, discrete spacetime models should preserve the symmetries of their continuous counterparts. Mathematically, Lie group techniques could then be used to find explicit solutions and compute conservation laws of the discrete models. From a more practical point of view, symmetry-preserving discretizations should share some exact solutions with the original differential equations, or at least provide better approximations than non-invariant numerical schemes.

In the last 30 years, the application of Lie group techniques to finite difference equations has become a very active field of research. To the best of our knowledge, Yanenko and Shokin were the first to use group theoretical methods to study finite difference schemes by introducing first differential approximations of difference equations, [82, 87]. The application of Lie group methods to finite difference equations, as we know it today, was first introduced by Dorodnitsyn in 1989, [25]. Early on, one of the main focuses in the field was to construct Lie point symmetry-preserving finite difference approximations of differential equations. Beside Dorodnitsyn, early contributors include Bakirova, Kozlov, Levi, and Winternitz who constructed symmetry-preserving schemes for heat transfer equations, [2, 3, 28], variable coefficient Korteweg–de Vries equations, [32], Burgers’ equation, [39], the nonlinear Schrödinger equation, [19], and second-order ordinary differential equations, [29]. Symmetry-preserving approximation of Euler–Lagrange equations and their corresponding Lagrangian have also been considered in [30, 31], and the application of Noether’s theorem to compute conservation laws has been extensively studied in the discrete setting, [26, 46, 48]. The applications of Lie point symmetries to finite difference equations have also been extended to generalized symmetries, [64, 65], -symmetries, [58, 59], and contact transformations, [60].

In recent years, more systematic efforts have been directed towards investigating the numerical performance of symmetry-preserving schemes. For ordinary differential equations, symmetry-preserving schemes have proven to be very promising. For solutions exhibiting sharp variations or singularities, symmetry-preserving schemes systematically appear to outperform standard numerical schemes, [15, 16, 20, 53]. For partial differential equations, the improvement of symmetry-preserving schemes versus traditional integrators is not as clear, [7, 52, 56, 78]. On one hand, it was shown in [85] that symmetry-preserving schemes do much better in tracking sharp interfaces in Hamilton–Jacobi equations. On the other hand, invariant numerical schemes for evolution equations generally require the use of time-evolving meshes which can lead to mesh tangling and thereby severely limit the use of symmetry-preserving schemes. In this case, special techniques have to be developed to avoid mesh singularities. For example, new ideas relying on r-adaptivity have been implemented to improve the performance of invariant integrators, [12]. Also, in [10, 11] an invariant evolution–projection strategy was introduced and invariant meshless discretization schemes were considered in [6].

The preceding references only provide a short bibliographical overview of the field. Many papers had to be omitted. More references on the subject can be found in the survey papers [63, 86], and the books [27, 46].

Given a differential equation with symmetry group , the first step in constructing a symmetry-preserving numerical scheme is to compute difference invariants of the product action of on a chosen stencil. There are mainly two approaches for constructing those invariants. Most of the references cited above use the infinitesimal symmetry generators of the group action and Lie’s infinitesimal invariance criterion to construct difference invariants. Alternatively, the difference invariants can be constructed using the novel method of equivariant moving frames mainly developed by Olver, which was done in [6, 23, 53, 70, 78, 79]. Given sufficiently many difference invariants, an invariant numerical scheme is, in general, obtained by finding a suitable combination of these invariants that converges to the original differential equation in the continuous limit. When using Lie’s infinitesimal generator approach, a suitable combination is found by taking the Taylor expansion of the difference invariants and combining them in such a way to obtained the desired invariant scheme. With the method of moving frames, a suitable combination is found more systematically by invariantizing a non-invariant numerical scheme. Since the symmetry group of a differential equation will, in general, act on both the independent and dependent variables, a symmetry-preserving numerical scheme will usually not be defined on a uniform orthogonal mesh.

The application of Lie groups to finite difference equations is a vast and very dynamic field of study. While preparing these lecture notes we had to omit many interesting applications and important results. The focus of these notes will be on the theoretical construction of invariant numerical schemes and their numerical implementation. At the heart of all our considerations are differential equations, finite difference equations, symmetry groups, and invariants. These familiar notions are all reviewed in Sections 2, 3, and 4. As outlined above, there are two different approaches for computing invariants of a Lie group action. The infinitesimal approach based on Lie’s symmetry generators is introduced in Section 4.1, while the equivariant moving frame approach is explained in Section 4.2. Section 5 is devoted to weakly invariant equations, which can play an important role in the construction of symmetry-preserving schemes. The construction of symmetry-preserving numerical schemes is carefully explained in Section 6. To illustrate the implementation, we consider the Schwarzian differential equation and the Korteweg–de Vries (KdV) equation. In Section 7 we carry out numerical simulations for the Schwarzian equation, the KdV equation and Burgers’ equation. For partial differential equations, the invariance of a numerical scheme does not, in general, guarantee better numerical results. We will show that symmetry-preserving schemes can lead to mesh tangling, which limit their practical scope. To circumvent this shortcoming, we discuss new invariant numerical strategies. For the Korteweg–de Vries equation, we introduce invariant evolution–projection schemes and invariant adaptive numerical schemes. Unlike the KdV equation, solutions to Burgers’ equation can exhibit shocks. For these shock solutions we propose a new invariant finite volume type scheme. Finally, in Section 8 we identify some open problems and challenges in the field of symmetry-preserving numerical integrators.

2 Differential and difference equations

In this section we review the definitions of differential equations and finite difference equations. We take this opportunity to introduce some of the notation used throughout these notes.

2.1 Differential equations

Let be an -dimensional manifold. For , let denote the extended order jet space of dimensional submanifolds defined as the space of equivalence classes of submanifolds under the equivalence relation of order contact at a point, [68]. Local coordinates on are given by the -jet


where correspond to the independent variables and denotes the derivatives

In the above notation, is an ordered -tuple of non-negative integers, with entries indicating the number of derivatives taken in the variable . The order of the multi-index, denoted by , indicates how many derivatives are being taken.

Example 2.1.

In the case where and , we have two independent variables and one dependent variable . Then, the second order jet space is parametrized by

Definition 2.2.

A differential equation of order is the zero locus of a differential map . That is,


For later use, we introduce two regularity requirements on differential equations.

Definition 2.3.

A differential equation is said to be regular if the rank of its differential

is constant on the domain of definition of .

Example 2.4.

Any evolutionary partial differential equation

is regular since the rank of is one.

Definition 2.5.

A differential equation is locally solvable at a point if there exists a smooth solution , defined in the neighborhood of , such that . A differential equation which is both regular and locally solvable is said to be fully regular.

The above description assumes that a submanifold is locally represented as the graph of a function . Alternatively, when no distinction between independent and dependent variables is made, a submanifold is locally parameterized by variables such that

In numerical analysis, the independent variables are called computational variables, [41]. These are the variables that are discretized when finite difference equations are introduced in Section 2.2. We let denote the order jet space of submanifolds parametrized by computational variables. Local coordinates on are given by


with , , and .

Remark 2.6.

We hope that the jet notations and will not confuse the reader. The independent variable, that is and , respectively, indicates with respect to which variables the dependent variables (and in the second case) are differentiated in .

Example 2.7.

In the case where and , let denote the two computational variables and let be a local parametrization of . Then the second order jet space is parametrized by

The transition between the jet coordinates (2.1) and (2.3) is given by the chain rule. Provided


the implicit total derivative operators

are well-defined, and successive application of those operators to the dependent variables will give the coordinate expressions for the -derivatives of in terms of the -derivatives of and :


We note that the non-degeneracy constraint (2.4) implies that the change of variables is invertible.

Example 2.8.

Combining Examples 2.1 and 2.7, assume that and are functions of the computational variables . Provided


the implicit derivative operators


are well-defined. It follows that


Relations for the higher order derivatives are obtained by applying (2.7) to (2.8).

Given a differential equation (2.2), the chain rule (2.5) can be used to re-express (2.2) in terms of , and their computational derivatives , :

Recall that for in (2.9a) while in . Equation (2.9a) can be supplemented by companion equations, [67],

The latter are introduced to impose restrictions on the change of variables . The system of differential equations (2.9) is called an extended system of the differential equation (2.2). For the extended system of differential equations (2.9) to share the same solution space as the original equation (2.2), the companion equations (2.9b) cannot introduce differential constraints on the derivatives .

Definition 2.5 is readily adapted to the computational variable framework.

Definition 2.9.

A differential equation (or system of differential equations) is regular if the rank of its differential

is constant on the domain of definition. The equation (or system of equations) is locally solvable at a point if there exists a smooth solution , , defined in the neighborhood of , such that and . The differential equation (or system of differential equations) is said to be fully regular if it is both regular and locally solvable.

Example 2.10.

As one of our main examples in these notes, we consider the Korteweg–de Vries (KdV) equation


We introduce the computational variables so that , . Then the implicit total derivative operators are given by (2.7). Before proceeding any further, we assume that


In other words,


where and are constants. The reasons for imposing the constraints (2.11) are explained in Example 3.10. The operators of implicit differentiation (2.7) then simplify to


and the KdV equation (2.10) becomes


together with the companion equations (2.11). The differential equation (2.13) is reminiscent of the equation one obtains when writing the KdV equation in Lagrangian form, [7]. In the classical Lagrangian framework, the differential constraint


is also imposed. The KdV equation then reduces to


together with the companion equations (2.11), (2.14). In particular, when in (2.12), we obtain the system of differential equations

2.2 Finite difference equations

We now move on to the discrete setting, which is the main focus of these lecture notes. In the previous section, we introduced two different jets spaces, namely and . The motivation for introducing computational variables and the corresponding jet space stems from the fact that the discrete framework is more closely related to than .

Let denote an integer-valued multi-index. Thinking of the multi-index as sampling the computational variables at integer values, the discrete notation should be understood as sampling the submanifold at the integer-valued points . In other words . To approximate the -jet at , we consider a finite collection of points


where . We require that the point is always included and that

so that no two discrete independent variables are the same. We refer to (2.16) as the order discrete jet at . In numerical analysis, a point in (2.16) is also called a stencil. For theoretical purposes, one can assume that the multi-index only takes non-negative values and that . The latter provides the minimal number of points required to approximate the -jet (or ) by first order forward differences. In applications, especially when constructing numerical schemes, it is generally preferable to consider points centered around and to include more than the minimum number of points in required to approximate for better numerical accuracy and stability. From now on, we will assume that a certain stencil (2.16) has been chosen. We denote by

the union over all the stencils and call the order discrete jet space as provides an approximation of . Since the jet coordinates of can be expressed in terms of the jet coordinates of using (2.5), it follows that the points in can be used to approximate .

Example 2.11.

Consider the case where and the dimension of the manifold is . Let be local coordinates on . In the continuous case, see Example 2.7, we introduced the computational variables . In the discrete case, let , which can be thought of as evaluating the computational variables at integer values. To make the multi-index notation more compact, we let

Working with forward differences, the simplest first order discrete jet is parametrized by

First order approximations of the first order derivatives and on a grid with unit spacing are then given by


Referring to (2.8) for the expressions of the and derivatives of in terms of the computational variable derivatives, and using (2.17) we have that


The latter expressions are first order forward approximations of the first order partial derivatives and on any mesh that satisfies

the latter being a discrete version of the non-degeneracy condition (2.6). The procedure can be repeated to obtain approximations of higher order derivatives on arbitrary meshes. For example, applying the implicit derivative operators (2.7) to the first order derivative expressions (2.8) one obtains formulas for the second order derivatives , , and expressed in terms of the computational derivatives. Substituting the approximations (2.17) and the second order derivative approximations

into the formulas obtained yields discrete approximations for , , and in the computational variables on an orthogonal grid with unit spacing.

Definition 2.12.

A finite difference equation is the zero locus of a discrete map . That is,

Definition 2.13.

A finite difference equation is said to be regular if the rank of the differential

is constant for all in the domain of definition of the equation.

Finite difference equations can be studied as mathematical objects of interest in their own, [33, 46, 66]. In the following we are interested in finite difference equations that approximate differential equations.

Definition 2.14.

A finite difference equation is said to be consistent with the differential equation (or if in the continuous limit ,

Remark 2.15.

The process of taking continuous limits is discussed in more details in Section 6.

Definition 2.16.

Let be a differential equation with extended system . A numerical scheme is a system of finite difference equations

where approximates the differential equation

and the equations provide an approximation of the companion equations

Intuitively, the difference equations provide constraints on the mesh used to approximate the differential equation . The latter should not yield any restrictions on the discrete dependent variables .

Example 2.17.

To illustrate Definition 2.16, let us consider the KdV equation (2.10). Assume the equation is to be discretized on the orthogonal mesh


where , , and , are arbitrary constants. The mesh (2.19) can be encapsulated in a system of finite difference equations in different ways. For example, it is not difficult to see that (2.19) is the solution to the system of equations


The mesh (2.19) is also a solution to


The difference between the two systems of mesh equations is that in (2.20) the time step and the spatial step are fixed by the system whereas in (2.21) those steps corresponds to constants of integration. In both cases, the KdV equation can be approximated by


The systems of equations (2.20)–(2.22) or (2.21)–(2.22) provide two examples of Definition 2.16. They also illustrate the fact that, in general, the equations specifying the mesh are not unique.

3 Lie symmetries

Let be an -dimensional Lie group, and let be a -dimensional manifold with local coordinates . In the following, the manifold can represent the submanifold jet spaces or or the discrete jet space . In the latter case, should in fact be called a lattifold or lattice variety, that is a manifold-like structure modeled on , [5, 47].

Definition 3.1.

A transformation group acting on a manifold is given by a Lie group and a smooth map , such that , which satisfies the following two properties


and where denotes the identity element.

It follows from (3.1) that the inverse of the transformation defined by the group element is given by the inverse group element . Therefore induces a diffeomorphism from to itself.

Remark 3.2.

Definition 3.1 assumes that the group action is global, meaning that is defined for every and every . In practice, group actions may only be defined locally, meaning that for a given , the transformation is only defined for group elements sufficiently near the identity. For a local transformation group, the map is defined on an open subset with , and the conditions (3.1) of Definition 3.1 are imposed wherever they are defined.

In the following, we use capital letters to denote the image of a point under a group transformation. For example,

At the infinitesimal level, let denote the Lie algebra of vector fields corresponding to the infinitesimal generators of the group action. A vector field

will be in if it is tangent to the orbits of a one-parameter subgroup of transformations of . The flow through the point generated by a vector field , is found by solving the initial value problem

The maximal integral curve is denoted , and is called the exponentiation of the infinitesimal generator .

Definition 3.3.

Let be a local Lie group of transformations acting on . The Lie group is a (local) symmetry group of the (fully) regular equation111Depending whether represents , , or , we refer to Definitions 2.5, 2.9, or 2.13 for the notion of (fully) regular equation. if and only if

for all such that the local action is defined. Infinitesimally, a connected Lie group of transformations acting on is a local symmetry group of if and only if

Remark 3.4.

Definition 3.3 extends to systems of equations and more general local groups of transformations by including discrete transformations as well, [9, 42, 43, 45]. In the following we restrict all our considerations to Lie point symmetries and omit the interesting case of discrete symmetries.

3.1 Symmetries of differential equations

Symmetries of differential equations are covered extensively in many excellent textbooks such as [13, 14, 43, 68, 71, 76]. We refer to these references for a more detailed exposition.

If , then the local group action is given by the prolonged action on the submanifold -jet . Let


denote the local group action of on the manifold locally coordinatized by . To compute the prolonged action, we introduce the implicit differentiation operators, [35],


denotes the entries of the inverse Jacobian matrix and

is the total differentiation operator with respect to the independent variable . In the above formula, denotes the unit vector with zeros everywhere except in the component. We note that the operators (3.4) mutually commute

Successively applying the implicit differentiation operators (3.4) to yields the expressions for the prolonged action

At the infinitesimal level, let


denote an infinitesimal generator of the group action (3.3). The prolongation of (3.5) to is given by

where the prolonged vector field coefficients are defined recursively by the standard prolongation formula

Given a differential equation , the Lie point symmetries of the equation are found from the infinitesimal invariance criterion


The latter yields a differential equation in , and the derivatives of with respect to , as well as and and their partial derivatives with respect to and . After eliminating any dependencies among the derivatives of the ’s due to the equation , one can equate the coefficients of the remaining unconstrained partial derivatives of to zero. This yields a system of linear partial differential equations for the coefficients and , called the determining equations of the (maximal) Lie symmetry algebra. The procedure for obtaining and solving the determining equations has been implemented in all major computer algebra systems such as Macsyma, Maple, Mathematica, MuMath and Reduce. An extensive list of references on the subject can be found in [21].

Example 3.5.

To illustrate the algorithm outlined above, we compute the infinitesimal generators of the KdV equation (2.10). Let denote a general vector field on . The third order prolongation of is