# Backward error analysis for multisymplectic discretizations of Hamiltonian PDEs

###### Abstract

Several recently developed multisymplectic schemes for Hamiltonian PDEs have been shown to preserve associated local conservation laws and constraints very well in long time numerical simulations. Backward error analysis for PDEs, or the method of modified equations, is a useful technique for studying the qualitative behavior of a discretization and provides insight into the preservation properties of the scheme. In this paper we initiate a backward error analysis for PDE discretizations, in particular of multisymplectic box schemes for the nonlinear Schrodinger equation. We show that the associated modified differential equations are also multisymplectic and derive the modified conservation laws which are satisfied to higher order by the numerical solution. Higher order preservation of the modified local conservation laws is verified numerically.

## 1 Introduction

When developing numerical integrators for Hamiltonian PDEs that possess a multisymplectic structure (i.e. symplectic in both space and time), it is natural to require the numerical scheme to preserve exactly a discrete version of the multisymplectic conservation law (MSCL) [2, 8]. However, this does not imply preservation of other dynamical invariants of the system such as the local energy and momentum conservation laws or global invariants which determine the phase space structure. A question that immediately arises then is, to what extent are the other invariants of the system preserved? Recent numerical experiments using multisymplectic integrators for nonlinear wave equations (e.g. the nonlinear Schrodinger (NLS), sine-Gordon, and Gross-Pitaevskii eqautions) show that the local conservations laws are preserved very well, although not exactly, over long times [5, 6, 7]. Further, the improved preservation of the local conservation laws is reflected in an improved preservation of complicated phase space structures [7]. This is reminiscent of the behavior of symplectic schemes for Hamiltonian ODEs. Symplectic integrators are designed to preserve the symplectic structure, not to preserve the energy. In fact, for general Hamiltonian systems, conservation of the symplectic structure and conservation of energy are conflicting requirements that, in general, are not solved simultaneously by a given scheme [4]. Even so, symplectic integrators preserve the Hamiltonian extremely well over very long times.

Backward error analysis (BEA), or the method of modified equations, is a particularly insightful technique for studying the qualitative behavior of a discretization as well as an alternative method for checking the accuracy of the numerical solution [12]. Since our main interest lies in the geometry preserving properties of multisymplectic schemes, the main question backward error analysis tries to answer (whether the distinguishing properties of the original equation carry over to the modified equation which the numerical solution satisfies to higher order) becomes relevant to our study. For a given scheme, the derivation of the associated modified equation is related to the calculation of the local truncation error and has, typically, been used to examine the dispersive, dissipative and diffusive properties of PDE discretizations. For example, in the numerical analysis of linear PDEs a backward error analysis of the Lax-Friedrichs method or the upwind method for the advection equation produces in both cases a modified equation that is an advection-diffusion equation. This helps one to understand the qualitative behavior of the methods and, from this perspective, explains why the numerical solution in both cases becomes smeared out as time evolves.

Likewise, BEA is an important tool in the study of geometric integrators [3, 4, 10, 9]. For Hamiltonian ODEs, symplectic methods lead to modified equations which are also Hamiltonian. In fact, the modified equation of a Hamiltonian ODE is also Hamiltonian if and only if the integrator is symplectic; this is then used to rigorously establish that a symplectic integrator almost preserves the total energy over an exponentially long period of time [4]. In striking contrast, nonsymplectic methods used to integrate Hamiltonian ODEs can introduce dissipation, a feature which is readily predicted by the dissipative form of the modified equations. Less has been established using BEA for Hamiltonian PDEs since there are a variety of ways to implement a BEA and the relevance of the analysis is open to interpretation. Spatial discretization of a PDE results in a system of ODES to which a standard BEA can be applied to derive a modified equation that is satisfied to higher order in one independent variable. Alternatively, a BEA can be used to derive modified equations for the PDE that are satisfied to higher order in both space and time [10, 9].

In this paper we implement a formal backward error analysis in both space and time of two multisymplectic box schemes, the Euler and the centered cell box schemes, as applied to the nonlinear Schrodinger equation. We find that the modified equations of these box schemes are also multisymplectic. The modified PDEs are used to derive modified conservation laws of energy and momentum that are approximated by the MS scheme to higher order in space and time. For the centered cell discretization of the NLS we numerically verify that the modified conservation laws are satisfied to higher order by the numerical solution. This provides a partial explanation of the superior resolution of the local conservation laws and global invariants by MS schemes (e.g. see the numerical experiments in section 5) and a deeper understanding of the local and global properties of MS integrators.

The paper is organized as follows. In the next section we recall the multisymplectic formulation of Hamiltonian PDEs and of the NLS equation. In section 3 we introduce the box schemes, establish multisymplecticity, and apply them to the NLS equation. We present a straightforward method for obtaining compact box schemes that is applicable to many multisymplectic PDEs. Section 4 contains the backward error analysis of the discretizations. In section 5 numerical experiments for the MS centered cell box scheme are discussed, illustrating the remarkable behavior of MS schemes. Higher order preservation of the modified local conservation laws is verified numerically, which supports the use of MS integrators in long time numerical simulations of Hamiltonian PDEs.

## 2 Multisymplectic Hamiltonian PDEs

A Hamiltonian PDE (in the “1+1” case) is said to be multisymplectic if it can be written as

(1) |

where are skew-symmetric matrices and is a smooth function of the state variable [11, 2]. The variational equation associated with (1) is given by

(2) |

The Hamiltonian system (1) is multisymplectic in the sense that associated with and are the 2-forms

(3) |

which define a symplectic space-time structure (symplectic with respect to more than one independent variable).

Any system of the form (1) satisfies conservation of symplecticity. Let be any solution of the variational equation (2). Then it can be shown that and , as defined in (3), satisfy the multisymplectic conservation law (MSCL):

(4) |

This result is obtained by noting that

since are skew-symmetric and is symmetric. The MSCL (4) is a local property and expresses the fact that symplecticity for Hamiltonian PDEs can vary locally over the spatial domain.

An important consequence of the MS structure is that when the Hamiltonian is independent of and , the PDE has local energy and momentum conservation laws [11, 2]

(5) | |||||

(6) |

For periodic boundary conditions, the local conservation laws can be integrated in to obtain global conservation of energy and momentum.

### 2.1 Multisymplectic formulation of the NLS equation

The focusing one dimensional nonlinear Schrödinger (NLS) equation,

(7) |

can be written in multisymplectic form by letting and introducing the new variables . Separating (7) into real and imaginary parts, we obtain the system [5]:

(8) |

which is equivalent to the multisymplectic form (1) for the NLS equation with

and Hamiltonian

Implementing (5)-(6) for the NLS equation yields the local energy conservation law (LECL)

(9) |

and the local momentum conservation law (LMCL)

(10) |

Additionally we have a norm conservation law for the NLS equation

(11) |

These three equations, when integrated with respect to , yield the classic global conservation of energy (Hamiltonian), momentum and norm .

## 3 Multisymplectic box schemes

Multisymplectic discretizations are numerical schemes for approximating (1) which preserve a discrete version of the multisymplectic conservation law (4). That is, if the discretization of the multisymplectic PDE and its conservation law are written schematically as

(12) |

and

(13) |

where , and and are discretizations of the corresponding derivatives and , then the numerical scheme (12) is said to be multisymplectic if (13) is a discrete conservation law of (12) [11, 2].

A standard method for constructing multisymplectic schemes is to apply a known symplectic discretization to each independent variable. For example, splitting the matrices and as

(14) |

and using the symplectic Euler forward-backward difference approximations on both space and time derivatives yields the Euler box scheme

(15) |

Similarly, applying the symplectic midpoint rule to both the time and space derivatives in (1) yields a “centered cell” box discretization

(16) |

where

(17) |

The local truncation error for the Euler box scheme is , while for the centered cell discretization it is .

Multisymplecticity of schemes (15) and (16) is easily established. For example, to do so for the centered cell scheme, we use the discrete variational equation associated with (16) given by

(18) |

Taking the wedge product of with (18), note that the right-hand side is zero, since is symmetric. The terms on the left-hand side can be simplified

whereas,

This implies that the numerical scheme (16) satisfies the discrete multisymplectic conservation law

### 3.1 Multisymplectic box schemes for the NLS equation

The multisymplectic centered cell box scheme was first developed for the NLS equation in [5] where an apparently ad hoc reduction provided a particularly compact form of the scheme. This reduction turns out to be generalizable as can be seen in McLachlan’s derivation of box schemes for the Korteweg de Vries equation [1]. Here we present a general approach for constructing compact box schemes which is applicable to many multisymplectic PDEs.

#### 3.1.1 Euler box scheme for the NLS equation

We begin by introducing the following finite difference operators

In terms of these operators the Euler box scheme (15) takes the form

(19) |

For the NLS, and are split using (14), where

Applying (19) to the NLS system (8) yields the system

After eliminating and the system reduces to

where we have set . When the second equation is shifted in time, the resulting six-point box scheme in stencil format is :

#### 3.1.2 Centered cell box scheme for the NLS equation

As before, we begin by introducing the appropriate finite difference operators

(20) |

In terms of these operators, the centered-cell discretization (16) becomes

(21) |

with discrete conservation law

The system which results upon applying (21) to (8) is

(22) |

Since the operators in (20) commute, by multiplying the first two equations in (22) by and back substituting and into the first two equations we obtain

Recombining these equations into a single complex equation (with ) yields the multisymplectic box scheme for the NLS equation

(23) |

or equivalently

(24) |

In finite difference stencil format the six-point box scheme is given by

The centered cell scheme naturally gives a two time level stencil for the NLS equation. If every term in (23) contained a common factor, e.g. or , further compactification would be possible. As it is, an additional reduction of (23) is not possible.

## 4 Backward Error Analysis

A useful method for analysing the qualitative behavior of symplectic methods for ODEs has been backward error analysis, where one interprets the numerical solution as the “nearly” exact solution of a modified Hamiltonian differential equation. In this section we implement a BEA in space and time for the multisymplectic box schemes. The modified differential equations are also multisymplectic and satisfy modified conservation laws.

### 4.1 BEA for the Euler box scheme

Let be a differentiable function that, when evaluated at the lattice points, satisfies the Euler box scheme (19). Using the Taylor series expansions in about

and equivalent expansions in , we obtain to first order the following modified equation

(25) |

If we introduce the new matrices

equation (25) can also be written in the multisymplectic form

where

and

Applying equation (25) to the NLS system and eliminating and yields the reduced system

or setting , the single equation

which is an perturbation of the NLS.

### 4.2 BEA for the centered cell box scheme

We now assume is a sufficiently smooth function that, when evaluated at the lattice points, is a solution to the centered cell scheme (21). Expanding in a Taylor series about the midpoints we obtain

where to simplify the notation and denote the grid points, denotes the midpoints, and . The symplectic midpoint rule approximation of the time derivative is given by

and, similarly, the space derivative is approximated by

Substituting these expansions into (16), one finds that, to order , satisfies the modified PDE

(26) |

where all quantities are evaluated at the midpoint .

When applying equation (26) to the NLS example, the resulting modified system of equations can be reduced to

which is an perturbation of NLS.

The modified local conservation laws can be obtained directly from equation (26) by multiplying the equation from the left by to obtain an energy conservation law and by to obtain a momentum conservation law. We prefer to show that the modified equation can be written in MS form and from this formulation obtain the associated local conservation laws via equations (5)-(6). Introducing the augmented variables

the modified equations (26) can be written in the MS form

(27) |

where are the skew-symmetric matrices given by

The modified multisymplectic PDE can be used to derive the modified LECL and LMCL. Substituting , and , into (5) and (6), the modified LECL and LMCL are found to be, respectively,

where and are given by equations (9)-(10). In the next section, we numerically verify that these modified local conservation laws are satisfied to higher order.

## 5 Numerical Results

For our numerical experiments we consider the NLS equation with periodic boundary conditions, . We use initial data for a multi-phase quasi-periodic (in time) solution, i.e., , . This initial data corresponds to a multi-phase solution, near the plane wave, characterized by one excited mode. We examine the performance of the centered cell box scheme (which we designate as MS-CC) for varying mesh sizes and time steps. The solution to equation (24) is found by writing it in matrix form and using an iteration technique to solve for .

The solution with and for is shown in Figure 1a.

The surface clearly exhibits the correct quasiperiodic behavior in time. In addition, we are interested in how well the local and global conservation laws are satisfied. To evaluate the local conservation laws, we use midpoint discretizations of the form

In general, these residuals are not zero (Figures 1b-c). The errors in the local conservation laws, the LECL and LMCL (9)-(10), are and , respectively, and are concentrated around the regions where there are steep gradients in the solution. If were a quadratic functional of , , with a symmetric matrix, then the local conservation laws would be conserved exactly [2]. In general, as in the present case, the PDE is nonlinear and is not a quadratic functional. Therefore, the local energy and momentum conservation laws will not be preserved exactly. However the numerical experiments show that the local conservation laws are preserved very well over long times. In addition to resolving the LECL and LMCL very well, the MS scheme preserves the global errors extremely well. The error in the global energy oscillates in a bounded fashion, as expected of a symplectic integrator (Figure 1d) while the error in the global momentum (Figure 1e) and the norm (not shown) are conserved exactly (up to the error criterion of in the solver) since they are quadratic invariants.

The maximum error in the LECL and LMCL and in the global energy and momentum for the MS scheme are provided in Table 1 for varying mesh sizes and time steps. The error in the LECL depends only on the timestep and is second order, while the error in the LMCL depends only on the spatial mesh size and is second order.

N | 32 | 32 | 32 | 64 | 64 | 64 |
---|---|---|---|---|---|---|

2.0E-02 | 1.0E-02 | 5.0E-03 | 2.0E-02 | 1.0E-02 | 5.0E-03 | |

LE | 6.0E-05 | 1.5E-05 | 4.0E-06 | 8.0E-05 | 2.0E-05 | 5.0E-06 |

LM | 1.7E-02 | 1.7E-02 | 1.7E-02 | 4.8E-03 | 4.8E-03 | 4.8E-03 |

GE | 7.3E-05 | 2.0E-05 | 5.0E-06 | 7.6E-05 | 2.2E-05 | 5.0E-06 |

GM | 1.2E-13 | 2.5E-14 | 2.0E-13 | 1.3E-13 | 1.0E-13 | 4.5E-13 |

We next examine whether the modified local conservation laws obtained using the MS-CC discretization of the NLS are preserved to a higher order than the original local conservation laws. Since our solution is quasiperiodic, we compute the solution for , where is chosen to include a characteristic cycle. From Figure 1, is sufficient. Since the ECL is independent of (see Table 1), for a fixed , we let . That is, start with and let .

We compute the LECL and the modified LECL at each time step using centered approximations of the derivatives of sufficiently high order so as not to affect the order of the MS-CC discretization of the residuals. Figure (2) shows the loglog plot of the maximum error as a function of the timestep for the LECL and the modified LECL.

Clearly we can see that while the LECL is satisfied to 2nd order, the modified LECL is satisfied to 4th order. Verification of higher preservation of the LMCM as a function of the mesh size is similar.

## Acknowledgements

This work was partially supported by the NSF, grant number DMS-0204714.

## References

- [1] U. Ascher and R. McLachlan, Multisymplectic box schemes and the Korteweg-de Vries equation, preprint 2003.
- [2] T.J. Bridges and S. Reich, Physics Letters A, 284, 184-193 (2001).
- [3] E. Hairer and Ch. Lubich, Numer. Math., 76, 441 (1997).
- [4] E. Hairer, Ch. Lubich and G. Wanner, Geometric Numerical Integration, Springer Verlag, Berlin, 2002.
- [5] A.L. Islas, D.A. Karpeev and C.M. Schober, J. of Comp. Phys. 173, 116–148 (2001).
- [6] A.L. Islas and C.M. Schober, Fut. Gen. Comp. Sys, 19, 403 (2003).
- [7] A.L. Islas and C.M. Schober, On the preservation of phase space structure under multisymplectic discretization, accepted J. of Comp. Phys. 2003.
- [8] J.E. Marsden and S. Shkoller, Math. Proc. Camb. Phil. Soc. 125, 553–575 (1999).
- [9] B. Moore and S. Reich, Fut. Gen. Comp. Sys., 19, 395 (2003).
- [10] B. Moore and S. Reich, Num. Mathematik, 95, 625 (2003).
- [11] S. Reich, J. of Comp. Phys. 157, 473–499 (2000).
- [12] J.W. Thomas, Numerical Partial Differential Equations, Springer Verlag, New York, 1995.