On the Rate of Error Growth in Time for Numerical Solutions of Nonlinear Dispersive Wave Equations

We study the numerical error in solitary wave solutions of nonlinear dispersive wave equations. A number of existing results for discretizations of solitary wave solutions of particular equations indicate that the error grows quadratically in time for numerical methods that do not conserve energy, but grows only linearly for conservative methods. We provide numerical experiments suggesting that this result extends to a very broad class of equations and numerical methods.


Introduction
In this work we study the behavior of numerical approximations of solitary wave solutions of nonlinear wave partial differential equations (PDEs). These solitary waves evolve in time only by translation; this translation symmetry is, in accordance with Noether's Theorem [50], associated with a particular conserved functional of the PDE. We investigate in a broad setting the phenomenon, already studied for certain equations, that numerical methods that conserve (up to rounding errors) this functional (or the Hamiltonian of the system) exhibit much smaller errors over long times. This behavior has been demonstrated for numerical approximations of certain solutions of Hamiltonian systems. These solutions are referred to as relative equilibria, and are those whose trajectory lies on a manifold defined by a symmetry group. For such solutions, numerical methods that exactly conserve certain invariants give rise to a leading term in the global error that grows linearly in time [24]. In contrast, for general (nonconservative) numerical solutions approximating relative equilibria, the error growth is typically quadratic (this should be contrasted with more general estimates of numerical error growth, which are typically exponential in time). This behavior has been shown for finite-dimensional Hamiltonian systems [12,14,15,24]. For dispersive nonlinear wave equations, it was first observed in the context of solitary wave solutions of the Kortewegde Vries (KdV) equation, by Sanz-Serna [72] and later by Bona, Dougalis, & Karakashian [10]. This behavior was eventually proved for discretizations of the KdV equation by de Frutos & Sanz-Serna [18].
Similar results were later obtained for solitary wave solutions of the nonlinear Schrödinger equation [25], Peregrine's regularized long wave (RLW) equation (also known as the Benjamin-Bona-Mahoney (BBM) equation) [7,49], and the generalized BBM equation [3,23]. An extension to certain other dispersive wave equations, including the generalized KdV, generalized BBM, and generalized Benjamin-Ono equations, is presented in [3]. The results in [49] are purely experimental, while the remaining results just cited include proofs along with numerical demonstrations.
Most of the PDE results in this direction are specific to solitary wave solutions. The conservative discretizations that yield linear error growth also exhibit modified solitary waves [7,23]. These are discrete solitary waves that are close to but different from the exact solitary wave solutions. The ability of these modified solitary waves to propagate discretely without changing shape or amplitude means that conservative methods are appealing for long-time simulations of solitary waves and their interactions [22]. The concept of a modified solitary wave also provides, for some systems, an intuitive rough explanation of the observed error growth rates. In a non-conservative method, the solitary wave is approximated by a wave whose amplitude grows or diminishes linearly in time.
Since the speed of the solitary wave is (at least to first approximation) a linear function of the amplitude, the error in the speed of the wave also grows linearly in time. This leads to a quadratic error in the location (phase) of the wave. Meanwhile, it was shown in [3,7,23] that a conservative scheme yields a modified solitary wave whose amplitude and speed are constant in time (though different from those of the true solitary wave). Hence the phase error grows only linearly.
The foregoing results serve as the motivation for the present study. Here we take an experimental approach, allowing the consideration of several other dispersive wave equations as well as non-dispersive hyperbolic systems that possess similar kinds of solutions. In this work we provide strong evidence that the phenomenon of linear error growth in time for conservative numerical solutions extends to a very large class of dispersive wave equations and nonlinear invariants, including invariants that are not quadratic, and even certain solutions of first-order hyperbolic systems that have no dispersive terms.
We briefly introduce the spatial and temporal discretization methods we will use in Section 3. In Section 4, we analyze the temporal error growth of solitary wave solutions to several nonlinear dispersive wave equations. Our experimental results include the Fornberg-Whitham [80], Camassa-Holm [13], Degasperis-Procesi [19], Holm-Hone [36], and BBM-BBM equations [9,11]. For each of these models, we apply recently-constructed conservative schemes and show experimentally that they yield improved (linear) error growth when compared to non-conservative schemes for the same model. To demonstrate the importance of nonlinearity for the different error growth behaviors of conservative and non-conservative methods, we investigate a linear dispersive wave equation in Section 5.2. There, both kinds of methods result in the same asymptotic error growth rate, although the absolute error of the conservative method is smaller. We then go further outside the realm of available results and consider two first-order nonlinear hyperbolic systems that have been shown to possess solitary-wave-like solutions. In Section 6 we study the variablecoefficient -system in one dimension, with periodically varying coefficients. Previous work on this system has noted similarities with dispersive nonlinear wave equations; see [40,46]. In Section 7, we study the two-dimensional shallow water equations with varying bathymetry. Finally, we summarize and discuss our results in Section 8.
The main contributions of this work are as follows: 1. We illustrate via computational experiments that the favorable behavior proved for conservative numerical solutions of certain nonlinear dispersive PDEs [3,7,18,25] holds for several other such PDEs.
2. We show through numerous examples that not only do conservative schemes exhibit a better asymptotic error behavior for large times, but also a much smaller size of the global error even for short times.
3. We show through a computational example that the condition of mass conservation, assumed in the main theorems of [7,18], is necessary in practice (see Section 5.1).
4. We demonstrate for the first time that the aforementioned favorable behavior can be observed also in the numerical solution of a system that does not possess solitary wave solutions with simple translation or rotation symmetries (see Section 6).
5. We demonstrate for the first time that dramatically better error growth for conservative methods can be observed also in systems with two spatial dimensions (see Section 7).
The source code used to generate the results of this study is available online [67]. The numerical methods studied in this article are implemented in Julia [8] and make use of the open source library SummationByPartsOperators.jl [63], which utilizes FFTW [29] for Fourier collocation methods. We use time integration methods from [58]. The plots are generated using Matplotlib [37].

Relative equilibrium solutions
One framework for understanding the rates of error growth observed in the present work is that of relative equilibrium theory. Here we review very briefly this theory in the context of finite-dimensional dynamical systems (where it is most completely developed) and Hamiltonian PDEs.
Relative equilibria are equilibria of a reduced dynamical system obtained by reduction modulo a set of symmetry groups of the original dynamical system. Relative equilibria evolve only along the manifold corresponding to the symmetry groups. In the case of finite-dimensional Hamiltonian systems, each symmetry group corresponds to an invariant quantity. A relative equilibrium 0 of a Hamiltonian system with Hamiltonian ℋ ( ) and invariant quantities ℐ ( ) satisfies [24,Eqn. (36)] where the values are Lagrange multipliers. For a relative equilibrium solution of a Hamiltonian system, perturbations in directions that lie on the manifold corresponding to the symmetry group grow linearly in time, while perturbations in other directions grow quadratically [24,Lemma 2.4]. Therefore, if a numerical method conserves the invariants corresponding to the symmetry group, the numerical error will grow only linearly in time [24,Theorem 3.2].
The situation is analogous but more involved for Hamiltonian PDEs like the dispersive wave equations considered in Section 4. Each of these equations can be written as a Hamiltonian system where ℋ is the Hamiltonian, is the variational derivative, and is a skew-symmetric operator. Each system possesses an additional invariant ℐ( ), such that This implies that the symmetry group of translations is associated with the invariant ℐ; the relative equilibrium condition (analogous to (2.1)) is It is possible to show also for these systems that numerical methods whose errors lie entirely within the appropriate manifold yield an error that grows only linearly in time; see e.g. [3,18,25] for examples of such analysis. Due to the condition (2.4), it is sufficient to conserve either ℋ or ℐ, in order to obtain linear growth in time of the leading term in the global error [3,18]. Rigorous results in this vein require a detailed perturbation analysis specific to the PDE in question. We do not pursue such analysis here, and focus instead on an experimental study.
We remark here that the results proving linear error growth for KdV and RLW solitary waves [7,18] also require the assumption of mass conservation. The necessity of this assumption has not been investigated experimentally, perhaps since "sensible methods exactly conserve the mass" [18]. However, the well-known orthogonal projection schemes can be used to conserve a nonlinear invariant at the cost of not conserving mass; in Section 5.1 we investigate the asymptotic error growth for such schemes.

Discretization methods
Next, we introduce the type of discretization methods employed in this article.

Spatial semidiscretizations
Proving the conservation of invariants of partial differential equations usually requires the product/chain rule and integration by parts. It is often useful to mimic the same procedure at the discrete level by using summation by parts (SBP) operators, which are constructed specifically to satisfy a discrete analogue of integration by parts. A review of the relevant theory can be found in [16,26,76]. Many classes of numerical methods can be formulated within the SBP framework, including finite difference [75], finite volume [52,53], continuous Galerkin [1,34,35], discontinuous Galerkin [30], and flux reconstruction methods [69]. A brief review of how to formulate these methods in the SBP framework with application to structure-preserving numerical methods used in this article is given in [68].
Next, we will briefly introduce SBP operators in periodic domains in one space dimension. Extensions to multiple space dimensions can be achieved via tensor products. We consider a grid = ( 1 , . . . , ) where 1 = min ≤ 2 · · · ≤ = max . All nonlinear operations will be performed pointwise, i.e. we use a collocation approach. with the convention 0 = 1 1 1 and 0 = 0 0 0. We say is consistent if ≥ 0. ⊳ For periodic boundary conditions (under whichsolutions at min and max are identical), integration by parts is basically a statement about the symmetry of a derivative operator with respect to the 2 scalar product. Discretely, an approximation of this scalar product is represented by a so-called mass matrix . 1

Definition 3.2.
A periodic first-derivative SBP operator consists of a grid , a consistent firstderivative matrix 1 , and a symmetric and positive-definite matrix such that ⊳ We will often refer to an operator as a (periodic) SBP operator if the other operators (such as the mass matrix ) are clear from the context. We always assume derivative operators are consistent, but we will usually omit this term. In this work, we use Fourier (pseudospectral) collocation methods [28,44] because of their efficiency and accuracy for smooth problems. This allows us to ensure, with reasonable grids, that the temporal discretization error dominates the spatial error for all of our 1D numerical tests. Nevertheless, all classes of methods within the SBP framework can be used and analyzed interchangeably to construct conservative semidiscretizations.
Semidiscretizations reduce an initial boundary value PDE to an initial value ODE Throughout this work we use superscripts to denote the time step, and here 0 is the initial condition. To conserve a nonlinear invariant functional ( ) discretely in a one-step method, we require ( ) = ( −1 ) = ( 0 ). For a Runge-Kutta method we define the update direction where we use the abbreviation := ( + Δ , ). Since the new solution will not be conservative in general, we modify the update formula (3.6b) by introducing a relaxation parameter and use To guarantee conservation of , is computed as root of This scalar nonlinear equation can be solved efficiently using algorithms such as the one of [2]. By the general theory on relaxation methods, there is exactly one root = 1 + (Δ −1 ) of (3.9) [66, Theorem 2.14]. By construction of the relaxation parameter , the resulting solution +1 conserves the invariant . Moreover, the relaxation approach also automatically conserves linear invariants (as long as the semi-discretization conserves them), which is an important prerequisite of analytical results on linear vs. quadratic error growth in time. Finally, the solution (3.8) has the same local order of accuracy as that given by the original Runge-Kutta method (3.6).

Nonlinear dispersive wave equations
In this section, we consider several nonlinear dispersive wave equations. Each possesses an additional invariant ℐ( ), and possesses stable solitary wave solutions that satisfy the relative equilibrium condition (2.4). Each equation also possesses one or more linear invariants, which are preserved by the discretizations we use. Since these invariants do not require any special numerical treatment, we do not discuss them explicitly.
For each equation we give a numerical method that conserves ℐ( ), based on SBP operators applied to a split form of the equation. In each case, always denotes a periodic th-derivative SBP operator. The matrices used for discretization of a given equation are also assumed to have in common a corresponding diagonal mass matrix . We will sometimes require an assumption that the derivative matrices commute; this will be stated explicitly when it is required.
Conservative numerical methods for these equations have been developed and analyzed in [68] using various classes of SBP operators. We will use split forms to preserve local conservation laws, since the chain and product rules cannot hold discretely for many high-order discretizations [61]. These are related to entropy-conservative methods in the sense of Tadmor [77]. While certain split forms have been known for some time [71, eq. (6.40)], they are still state of the art and enable the construction of numerical methods with desirable properties [32,60,74,82].
For each equation, after reviewing the Hamiltonian structure and split-form conservative spatial discretization, we present a numerical test of error growth for a solitary wave solution. The solitary wave solutions are obtained via the Petviashvili method [56] using a Fourier collocation method with 2 16 nodes. The resulting solitary wave profile is interpolated to a grid using fewer nodes and used as initial condition. A conservative semidiscretization is obtained by using Fourier collocation methods in space. This semidiscretization is integrated in time using the fifth-order accurate Runge-Kutta method of Tsitouras [78] with adaptive time stepping based on local error estimates. We show results for both this method without relaxation (non-conservative) and with relaxation in time (conservative). The errors plotted are discrete 2 errors computed using the discrete norm induced by the mass matrix , which is the identity matrix times the grid spacing for Fourier collocation methods.

Fornberg-Whitham equation
For convenience we define = I − 2 . Then the Fornberg-Whitham equation [80] ( , ) + ( ( , )) + ( , ) = 0, with periodic boundary conditions can also be written as where −1 is the inverse of the elliptic operator I − 2 with periodic boundary conditions. This equation possesses the Hamiltonian and additional nonlinear invariant It can be written in the form (2.2) by taking = − . The discrete equivalent of ℐ( ) is conserved by semidiscretizations of the form where 1 and 2 commute [68].

Camassa-Holm equation
The Camassa-Holm equation [13] ( , ) with periodic boundary conditions can also be written as where ∈ R is a parameter determining the split form [68]. This equation possesses the Hamiltonian and an additional nonlinear invariant

BBM-BBM system
The BBM-BBM system [5,6,9,11] ( , ) with periodic boundary conditions can also be written as This equation admits the nonlinear invariants and can be written in the form (2.2) by taking As mentioned in Section 2, we expect to observe linear error growth whenever the numerical method conserves either ℋ ( , ) or ℐ( , ). For this example, we demonstrate this by testing methods that conserve each of these quantities separately. The discrete equivalent − 1 2 ( + (1 1 1 + ) ) of ℋ ( , ) is conserved by semidiscretizations of the form + (I − 2 ) −1 1 ( + ) = 0 0 0, where 1 , 2 commute [68]. Conservation of the discrete equivalent (I − 2 ) of ℐ( , ) can be achieved with the standard Galerkin method [48] or SBP methods using the split form Applying the SBP property (3.2) additionally and using that is diagonal yields Finally, similar arguments based on the symmetry of I − 2 with respect to the mass matrix lead to (4.20)

Holm-Hone equation
As our last example in this section, we consider the Holm-Hone equation [36]: with periodic boundary conditions. It can also be written as

Numerical results
For each of the five equations just studied, we perform a similar numerical experiment. We use a numerically generated smooth solitary wave solution as initial condition in a periodic domain, and apply the corresponding SBP spatial discretization with Fourier collocation methods in space and Tsitouras' fifth-order Runge-Kutta method [78] in time, with or without using relaxation to enforce conservation of a nonlinear invariant ℐ. Specific parameters for each simulation are shown in Table 1.
In each case, we observe approximately linear error growth for the conservative method and approximately quadratic growth for the non-conservative method, until eventually the error for the non-conservative method saturates when the exact and numerical solutions no longer overlap.
An example of the numerical solutions themselves is shown in Figure 3, using BBM-BBM. The numerical solutions obtained using the non-conservative method have a visible amplitude error and a significant phase error. In contrast, the conservative method yields numerical solutions that are visually indistinguishable from the reference solution. Results for the other equations are similar.
For the BBM-BBM system, we also conduct numerical experiments using the discretization (4.17) and applying relaxation to conserve the corresponding functional. The error growth in time using this method is visually indistinguishable from that shown for the method (4.16). As expected, the error of the associated conservative method grows linearly in time while the error of the non-conservative method grows quadratically.
Additionally, we perform numerical experiments with non-smooth traveling waves. While these are available in closed form for some of the nonlinear dispersive wave equations studied in this article, their reduced regularity makes it more difficult to approximate them numerically well enough in space. Note that existing analytical results about the error growth in time described in Section 2 are based on discretization only in time and assume exactness in space. Thus, the spatial semidiscretization must be well-resolved in order for these results to apply to fully discrete schemes.
Here, we use eighth-order accurate finite difference methods with 2 13 nodes for peakon solutions of the form ( , ) = exp(−| − |) with = 1.2 of the Camassa-Holm equation [13,45] in the periodic domain [− 35,35]. The fifth-order Runge-Kutta method uses a local error tolerance of 10 −7 . The results visualized in Figure 4 are in accordance with the expectations: We observe quadratic error growth for non-conservative methods and linear error growth for conservative methods. ( 1 )

Significance of linear invariants
In this section we investigate the importance of conservation of linear invariants. First, we study the behavior of a method that preserves energy but not mass for the KdV equation. Then, we consider a linear PDE in which nonlinear invariants play no role.

Projection methods: non-conservation of mass
Here we consider the propagation of a soliton solution of the KdV equation, using three methods. All three methods use the SBP spatial discretization of [65] with eighth-order accurate finite differences that is mass-and energy-conserving. The time discretizations are: • An IMEX RK method of [38] that conserves mass but not energy • The corresponding relaxation method that conserves both mass and energy ( 1 )

Error of (non-conservative)
Error of (conservative) Error of (non-conservative) Error of (conservative) • The corresponding orthogonal projection method that conserves energy but not mass For details on orthogonal projection, we refer to [33,Section IV.4]. A similar comparison was conducted in [65]; here we conduct a more detailed study of the error growth for the projection scheme.
Results are shown in Figure 5. The first two methods behave as expected, in accordance with theoretical results in [18]. The projection method initially performs similar to the relaxation method. However, its error grows quadratically after some time and eventually leads to a much more inaccurate solution, in accordance with the theory of [18]. Interestingly, the large global error in total mass at late times is exhibited not primarily in a change in mass of the soliton, but in a gradual downward drift of the ambient value of , as shown in Figure 6. We observed a similar behavior of orthogonal projection methods also for other nonlinear dispersive wave equations such as the BBM equation.

A linear dispersive equation
As alluded to in the introduction, nonlinearity is an important condition to observe different behaviors of the error growth in time. To demonstrate that, we consider the linear dispersive equation   Figure 7. The spatial semidiscretization using a Fourier pseudospectral method with = 2 6 nodes is integrated in time with the fifth-order accurate Runge-Kutta method of [78] and a tolerance of 10 −5 . While the conservative method has a significantly reduced amplitude error, the phase error of the conservative and non-conservative methods are visually indistinguishable. Consequently, the error growth rate in time for both methods is linear. This is in accordance with analytical results available for other linear equations such as hyperbolic systems with constant coefficients in periodic domains; for these equations, the error can also be bounded in time if non-periodic domains are used and the boundary conditions are imposed appropriately [43,51,54,55].

The variable-coefficient -system
The examples in Section 4 all involve dispersive nonlinear wave equations with constant coefficients. The results there are in line with the results available in the literature, all of which involve similar dispersive nonlinear wave equations. We now turn to a very different kind of example that demonstrates the generality of the behavior in question. Specifically, we study a system of nonlinear first-order hyperbolic conservation laws with spatially-varying coefficients. This system has no dispersive terms, and in a strict sense has no traveling wave solutions. It has been observed to possess solutions (known as "stegotons") that are similar to solitary waves, although they change shape periodically in time; see [46]. Specifically, we investigate the behavior of numerical solutions to the first-order system where is the velocity, the prescribed density, the strain, and the stress [46]. Here we have used notation corresponding to elasticity, but the same system arises in Lagrangian gas dynamics and is referred to as the -system. The energy is conserved for strong solutions of (6.1). Generically, solutions of this system give rise to discontinuities (shocks) and a theory of weak solutions must be invoked. If ( ) and ( , ) do not depend explicitly on then the system has no traveling wave solutions.   However, with appropriate initial data and periodic variation of the PDE coefficients, solutions appear to be regular for all time and the energy (6.3) is conserved [40]. A straightforward application of an SBP operator 1 in a periodic domain results in the semidiscretization = 1 , where and the division by are evaluated pointwise.   Furthermore, where we have used that the mass matrix is diagonal.
We use smooth coefficients given by in the domain [0, 20]. To study the error growth in time, we construct a stegoton solution numerically using Clawpack [17,41,42,47] as follows. We start with a zero initial condition and a left boundary condition given by where 0 = −2.5 2.5 . This corresponds to a moving wall that generates a pulse that eventually becomes a train of stegotons. We solve this problem to a late time (so that the stegotons are well separated) on a highly-refined grid, and then isolate the first resulting stegoton. For more details see [67].
Next, we use the numerically constructed single-stegoton solution as initial condition. We apply an energy-conservative Fourier pseudospectral semidiscretization (6.4) with 2 8 nodes. We solve the resulting ODE system using the fifth-order Runge-Kutta method of [78] with adaptive time stepping and a tolerance of 10 −6 . The error growth is shown in Figure 8 for solutions obtained with relaxation (conservative) and without relaxation (nonconservative). After an initial transient period, the error of the non-conservative scheme grows quadratically while the corresponding conservative method results in approximately linear error growth. The error of the non-conservative method starts to saturate at ≈ 10 4 since the numerical and reference waves do not overlap anymore, as can be seen in Figure 9. The numerical solutions obtained using the non-conservative method have a clearly visible phase error. In addition, the shape of the numerical solutions is deformed, in particular for . In contrast, the phase error of the conservative method is negligible. Nevertheless, the shape of the numerical solutions has changed a bit over these very long-time simulations, and small oscillations have appeared in the tails of the solitary wave. ( 1 ) Error of (non-conservative) Error of (conservative) Error of (non-conservative) Error of (conservative) Figure 8: Error growth in time for a numerically generated stegoton solution of the variablecoefficient -system (6.1), using (6.4) with Fourier collocation in space and fifth-order Runge-Kutta in time. The conservative method uses relaxation to enforce conservation of (6.3).
Although existing theory such as that reviewed in Section 2 cannot be applied to this system, it is perhaps not unreasonable to expect the behavior we have observed. Although stegotons are not translation invariant, they appear to be solutions of the form where˜ is periodic in its second argument; i.e.
where Ω = 1 is the period of the material coefficients ( ), ( ). This "periodic translation invariance" may fulfill the same role as simple translation invariance. Note also that the system considered here behaves, with a certain extreme choice of coefficients, like the discrete Toda lattice [46]; for the latter system, it is known that symplectic numerical integrators exhibit linear error growth [33, pp. 413-414].

The shallow water equations
We now turn to the shallow water equations in two space dimensions. Like the last example, this is a first-order hyperbolic system whose solutions generically develop shock  discontinuities. As with the last example, the constant-coefficient equations have no traveling wave solutions; however, for appropriate initial data and varying bathymetry, numerical experiments suggest that solitary wave solutions exist [57]. Nevertheless, this example differs in important ways from all the preceding ones and from all the results available in the literature on error growth for dispersive nonlinear waves. It is a twodimensional system, and its solitary wave solutions have a two-dimensional structure (localized in and periodic in ). Aside from the theoretical complication this involves, it also imposes a practical challenge. Since the exact form of the solitary wave solutions is not known, we will need to compute an approximate solitary wave initial condition, as we did in Section 6. However, in this case we cannot exhaustively resolve the solution and we expect that resulting initial condition is close to, but still different from, a solitary wave. The shallow water equations in two space dimensions with variable bathymetry ( , ) are ℎ + (ℎ ) + (ℎ ) = 0, where ℎ( , , ) is the water height and ( , , ), ( , , ) are the -and -components of velocity, respectively. The conserved energy is A two-parameter family of well-balanced, energy-conservative semidiscretizations of the shallow water equations with variable bathymetry was presented in [59,62]. The semidiscretization ℎ ℎ ℎ = − 1, ℎ ℎ ℎ − 1, ℎ ℎ ℎ , (7.3a) is a member of this family and was presented earlier in [31,81], see also [27] for related methods.
In the following, we use the gravitational constant = 9.8 and the smooth bottom topography in the domain [0, 20] × [−0.5, 0.5] with periodic boundary conditions. Results in [57] and additional computations we have conducted indicate that a wide range of initial conditions lead to solitary waves that propagate along the coordinate direction, transverse to the variation in bathymetry. These solitary waves have a non-trivial structure in both spatial dimensions.
To study the error growth in time, we construct a solitary wave solution as follows. Based on [57], we start with an initial condition given by where * = 0.75 and = 5 × 10 −2 . After propagating the initial condition up to a final time of = 340, we obtain a train of solitary waves, from which we isolate the largest wave. For more details see [67]. Next, we use the numerically constructed wave as initial condition for an energy-conservative Fourier pseudospectral discretization (7.3) with 2 10 nodes in and 2 6 nodes in . We solve the resulting ODE using the fifth-order Runge-Kutta method of [78] with adaptive time stepping and a tolerance of 10 −6 . The error growth is shown in Figure 10. In contrast to the examples considered previously, the error growth rates of both methods are sublinear. In agreement with previous results, the error of the conservative method still grows significantly slower, resulting in an error at the final time (after 15 periods) that is approximately an order of magnitude smaller. Snapshots of the numerical solutions at the final time are shown in Figures 11 and  12. The conservative method advects the profiles of the perturbations well, resulting in a numerical solution that is visually nearly indistinguishable from the initial profile. In contrast, the non-conservative method results in a slightly visible phase error. In addition, signs of nonlinear instabilities in the form of high-frequency oscillations start to manifest in the y-momentum and in reduced form also in the x-momentum and the total water height.
Possible reasons for error growth rates that are not linear/quadratic stem from differences relative to the other equations studied so far: this system is two-dimensional, there are no known exact solitary wave solutions, and the computed initial solitary wave is still affected by small numerical errors due to the cost of fully resolving the wave in two dimensions. Additionally, the spatial resolution may not be sufficient, so spatial discretization errors may be significant in these simulations. Although these numerical results presented for the shallow water equations do not match the linear/quadratic error growth rates of the other examples discussed previously, they still show the usefulness of conservative methods. In the current implementation, the CPU time of the conservative and non-conservative method are comparable while the accuracy differs by an order of magnitude. To obtain similarly accurate solutions using the non-conservative method would require an increased space/time resolution and hence more computational resources.

Summary and conclusions
We have studied the error growth in time of numerical solitary wave solutions of several systems of PDEs, focusing on the difference in behavior between conservative and non-conservative methods. In all cases, the general pattern is the same: conservative methods exhibit linear error growth, while non-conservative methods exhibit quadratic error growth. As one might expect based on this, the magnitude of the error itself is also much smaller for conservative methods.
We have shown that this phenomenon extends to a wide range of dispersive nonlinear wave equations, and that this behavior seems to arise when any nonlinear invariant is conserved, even if it is not quadratic (see Section 4.4). Furthermore, we have shown that this behavior extends to other classes of equations, including the -system example in Section 6. This is remarkable in that the equations in question are non-dispersive and the solutions are not traveling waves in the usual sense. Nevertheless, a similar advantage is obtained by using a conservative discretization.
This suggests that conservative methods may be much more efficient for solving wave problems that possess one or more conserved functionals. Efficiency depends also on the cost of the method, which we have not focused on here. In most previous works related to this phenomenon, energy conservation was achieved through the use of fullyimplicit Runge-Kutta time integration. With the relaxation Runge-Kutta approach, we can instead employ essentially explicit relaxation Runge-Kutta methods for non-stiff spatial discretizations, and diagonally implicit or linearly implicit (Rosenbrock) methods for stiff problems. These incur only a very small additional cost per time step (in order to  solve a scalar algebraic equation). For the methods we have compared here, at least, the conservative approach using relaxation is much more efficient (in terms of computational cost for a given level of error) than the corresponding non-conservative method. We expect that, for the nonlinear dispersive wave equations studied, a rigorous theoretical explanation of our results could be obtained by perturbation analysis similar to what has been done for other similar PDEs. It is much less clear how to obtain a theoretical justification for the results regarding variable-coefficient first-order hyperbolic PDEs. A first step in this direction might be a study focused on finite-dimensional non-autonomous Hamiltonian systems.
It is natural to ask whether the behavior studied here can be observed for more general solutions (not consisting of a single traveling wave). Some theoretical results are available for other classes of solutions; e.g. for multi-soliton solutions of KdV [4] and quasi-periodic solutions of KdV [21]. An extension of the current study to more general solutions is  Figure 11. the subject of ongoing work, and initial tests suggest a complex variety of behaviors. Nevertheless, these tests indicate that conservative methods are consistently much more accurate than their non-conservative counterparts.