TY - JOUR
AU1 - Frasca-Caccia,, Gianluca
AU2 - Hydon, Peter, Ellsworth
AB - Abstract Conservation laws are among the most fundamental geometric properties of a partial differential equation (PDE), but few known finite difference methods preserve more than one conservation law. All conservation laws belong to the kernel of the Euler operator, an observation that was first used recently to construct approximations symbolically that preserve two conservation laws of a given PDE. However, the complexity of the symbolic computations has limited the effectiveness of this approach. The current paper introduces some key simplifications that make the symbolic–numeric approach feasible. To illustrate the simplified approach we derive bespoke finite difference schemes that preserve two discrete conservation laws for the Korteweg–de Vries equation and for a nonlinear heat equation. Numerical tests show that these schemes are robust and highly accurate compared with others in the literature. 1. Introduction The main goal of geometric integration is to reproduce, in a numerical approximation, key geometric properties of a given continuous differential problem (see Budd & Piggott, 2003; Hairer et al., 2006). For instance, Hamiltonian ordinary differential equations (ODEs) occur in applications from nanoscale molecular dynamics to the macro-scale of celestial mechanics (see Hairer et al., 2006; Brugnano & Iavernaro, 2016). They have two fundamental features: symplecticity of the flow in phase space and constancy of the Hamiltonian function on solutions. Consequently, geometric integration of Hamiltonian ODEs has followed two main approaches, preserving symplecticity and energy. Symplectic methods are obtained by requiring that the discrete map associated with a given numerical method is symplectic (Feng, 1985; Sanz Serna & Calvo, 1994; Sanz Serna, 1998; Leimkuhler & Reich, 2004; Hairer et al., 2006). Energy conservation has been achieved by using discrete line integral methods (see Brugnano et al., 2010, 2012, 2015d; Brugnano & Iavernaro, 2016), time finite element methods (Betsch & Steinmann, 2000; Tang & Chen, 2007; Tang & Sun, 2012) and discrete gradient methods (Gonzalez, 1996; McLachlan et al., 1999; Dahlby & Owren, 2011) such as the average vector field method (see Quispel & McLaren, 2008; Celledoni et al., 2009; Hairer, 2010). These structure-preserving approaches have been extended to Hamiltonian partial differential equations (PDEs) (Bridges, 1997; Leimkuhler & Reich, 2004; Bridges & Reich, 2006). A particularly powerful approach uses a multisymplectic reformulation of the equations (Bridges, 1997; Bridges & Reich, 2001, 2006; Islas et al., 2001; Chen et al., 2002; Sun & Qin, 2003, 2004; Islas & Schober, 2004; Leimkuhler & Reich, 2004; Ascher & McLachlan, 2004, 2005). Alternatively, the method of lines is used to create a semidiscretization, and the resulting Hamiltonian ODEs (in time) are integrated by a symplectic method (Qin & Zhang, 1990; Lu & Schmid, 1997; Bridges & Reich, 2001; Oliver et al., 2004; Ascher & McLachlan, 2005; Cano, 2006; Guan et al., 2009; Bambusi et al., 2013) or energy-conserving method (Furihata, 1999; Koide & Furihata, 2009; Dahlby & Owren, 2011; Frasca-Caccia, 2015; Guo & Xu, 2015; Brugnano et al., 2015a,b,c; Barletti et al., 2016, 2017, 2018; Brugnano & Iavernaro, 2016). For Hamiltonian PDEs McLachlan & Quispel (2014) made the useful observation that ‘if the semidiscretization has a semidiscrete energy conservation law, then a discrete gradient method applied to this semidiscretization will have a fully discrete energy conservation law’. The benefits of preserving global invariants have been examined for several Hamiltonian PDEs in De Frutos & Sanz-Serna (1997), Durán & Sanz-Serna (2000), Durán & López-Marcos (2003) and, in a more general context, Durán & Sanz-Serna (1998). The current paper introduces a simple bespoke approach to constructing finite difference schemes that preserve multiple conservation laws of a given PDE. Conservation laws are among the most fundamental features of the PDE, as their origin is topological. Our approach is a simplification of the symbolic–numeric strategy introduced in Grant (2011) and Grant & Hydon (2013) and developed in Grant (2015). There are three advantages to this approach. First, it does not require the PDE to have any special structure, so it is suitable for discretizing PDEs independently of whether or not they possess other geometric structures. Secondly, the discretizations obtained by using this strategy exactly preserve local discrete conservation laws. Conserving local features of the continuous PDE gives, in general, a stricter constraint than preserving the corresponding global features. Given suitable boundary conditions, the preservation of local conservation laws also ensures the conservation of the corresponding global invariants. Finally, our approach can be used to seek methods that preserve any number of conservation laws. However, imposing the preservation of more than two conservation laws can considerably increase the complexity of the scheme. For this reason, in this paper, we deal only with methods that preserve two conservation laws, as a reasonable compromize between reliability and complexity of the schemes. In Section 2 we review Grant’s symbolic–numeric approach and introduce the simplifications that we will use to construct new conservative finite difference schemes. A different strategy, the multiplier method, has been proposed in Wan et al. (2016) to construct conservative finite difference methods for ODEs and PDEs. We briefly discuss the two different approaches. In Section 3, the simplified symbolic–numeric approach is applied to the Korteweg–de Vries (KdV) equation. Several new schemes are constructed and numerical tests are presented to show their effectiveness by comparison with some known methods that preserve only one conservation law. In Section 4, we consider a nonlinear heat equation as an example of a non-Hamiltonian PDE having two conservation laws. A family of two-parameter methods preserving both conservation laws is introduced. (These are easily extended to the more general porous medium equation.) At the end of the section, we present numerical tests that show the conservative properties of the new schemes and comparisons with a standard second-order finite difference method. Some concluding remarks are given in Section 5. 2. How to preserve multiple conservation laws We begin this section with some basic results on conservation laws of PDEs. After reviewing the general symbolic–numeric strategy for preserving multiple conservation laws of a given scalar PDE, we introduce some simplifications that enable accurate schemes to be derived efficiently. We restrict attention to scalar PDEs with two independent variables; the approach generalizes to more variables, but more simplifications may be needed to make the symbolic computations tractable. Consider a PDE for |$u(x,t)$|⁠, \begin{equation} {\mathcal{A}}(x,t,[u])=0, \end{equation} (2.1) where |$[u]$| denotes |$u$| and finitely many of its derivatives. More generally, square brackets around a differentiable expression denote the expression and finitely many of its derivatives. To simplify the exposition we will assume that |${\mathcal{A}}$| is at most quadratic in |$[u]$|⁠; the generalization to PDEs that are polynomial in |$[u]$| is obvious. For an application of our approach to a PDE that is cubic in |$[u]$| see Frasca-Caccia (2018). A conservation law of (2.1) is a divergence expression, \begin{equation*}\mbox{Div}\,\textbf{F} \equiv D_x\{F(x,t,[u])\}+D_t\{G(x,t,[u])\},\end{equation*} which is zero on all solutions of (2.1), that is, \begin{equation} \mbox{Div}\,\textbf{F} =0\quad \mbox{when}\quad [{\mathcal{A}}=0]. \end{equation} (2.2) Here |$D_x$| and |$D_t$| are the total derivatives w.r.t. |$x$| and |$t$|⁠, respectively: \begin{eqnarray*} D_x&\equiv &\frac{\partial}{\partial x}+u_x\frac{\partial}{\partial u}+u_{xx}\frac{\partial}{\partial u_x}+u_{xt}\frac{\partial}{\partial u_t}+\cdots,\\ D_t&\equiv &\frac{\partial}{\partial t}+u_t\frac{\partial}{\partial u}+u_{tt}\frac{\partial}{\partial u_t}+u_{tx}\frac{\partial}{\partial u_x}+\cdots. \end{eqnarray*} The components |$F$| and |$G$| are commonly referred to as the flux and density, respectively. A conservation law (2.2) is trivial of the first kind if |$F$| and |$G$| are zero on solutions of (2.1). It is trivial of the second kind if the divergence in (2.2) is identically zero without any reference to the PDE. A conservation law is trivial if and only if it is a linear superposition of the two types of trivial conservation laws. Two conservation laws are equivalent if they differ by a trivial conservation law. If the conservation law (2.2) amounts to \begin{equation} \mbox{Div}\,\mathbf{F}={\mathcal{Q}}{\mathcal{A}}, \end{equation} (2.3) it is said to be in characteristic form and the multiplier |${\mathcal{Q}}$| is called a characteristic of the conservation law. Remark 2.1 If the PDE is in Kovalevskaya form, integrating any of its conservation laws by parts yields an equivalent conservation law in characteristic form (see Olver, 1993). A characteristic, |${\mathcal{Q}}$|⁠, is trivial if it vanishes on solutions of (2.1); two characteristics are equivalent if they differ by a trivial characteristic. If (2.1) is in Kovalevskaya form there is a one-to-one correspondence between equivalence classes of characteristics and equivalence classes of conservation laws (see Alonso, 1979; Olver, 1993). Therefore, characteristics can be used to test the equivalence of conservation laws. A crucial result, for our purposes, is the characterization of the kernel of the Euler operator, \begin{equation*} \mathcal E=\sum_{i,\,j}(-D_x)^i(-D_t)^{\,j}\frac{\partial}{\partial u_{x^{i}t^{\,j}}}\,,\qquad\textrm{where}\quad u_{x^it^{\,j}}=D_x^{\,i}D_t^{\,\,j}(u), \end{equation*} as the space of total divergences. Consequently, if |${\mathcal{Q}}$| is a function such that \begin{equation*}\mathcal E({\mathcal{Q}}{\mathcal{A}})\equiv 0,\end{equation*} then there exists |$\mathbf{F}$| s.t. |${\mathcal{Q}}{\mathcal{A}}=\mbox{Div}\,\mathbf{F}$|⁠, and therefore, |${\mathcal{Q}}$| is the characteristic of the corresponding conservation law. Conservation laws, being defined as divergences (2.2), are local features of (2.1). Their integrals over the spatial domain yield quantities that are globally conserved on solutions (provided that (2.1) is coupled with suitable boundary conditions). However, although the local preservation of (2.2) implies the preservation of the globally conserved quantities, the converse is not true. For this reason we seek finite difference schemes that preserve discrete analogues of continuous local conservation laws. For simplicity we consider only uniform discretizations of the PDE (2.1). Relative to a generic lattice point |$\textbf{n}=(m,n)$| the grid points are \begin{equation*}x_i=x(m+i)=x(m)+i\varDelta x,\qquad t_j=t(n+j)=t(n)+j\varDelta t,\end{equation*} and the approximated values of the dependent variable |$u\in \mathbb{R}$| at these points are \begin{equation*} u_{i,\,j}\approx u(x_i,t_j),\quad i,\,j\in\mathbb{Z}. \end{equation*} The forward shift operators |$S_m$| and |$S_n$| are defined on the lattice by \begin{equation*} S_m:(m,n)\mapsto (m+1,n),\qquad S_n:(m,n)\mapsto (m,n+1); \end{equation*} their action extends naturally to |$x_i,t_j$| and |$u_{i,\,j}$| as follows: \begin{equation*} S_m:(x_i,t_j,u_{i,\,j})\mapsto (x_{i+1},t_j,u_{i+1,\,j}),\qquad S_n:(x_i,t_j,u_{i,\,j})\mapsto (x_i,t_{j+1},u_{i,\,j+1}). \end{equation*} Combining |$S_m$| with the identity operator \begin{equation*} I:(m,n,x_i,t_j,u_{i,\,j})\to(m,n,x_i,t_j,u_{i,\,j}) \end{equation*} yields the forward difference, |$D_m$|⁠, and the forward average, |$\mu _m$|⁠, defined for all functions |$f$| by \begin{equation*} D_m(\,f)=\tfrac{1}{\varDelta x}\,(S_m-I)(\,f),\quad\ \mu_m(\,f)=\tfrac{1}{2}(S_m+I)(\,f). \end{equation*} To obtain backward versions of the above operators compose each with |$S_m^{-1}$|⁠. Similarly, \begin{equation*} D_n(\,f)=\tfrac{1}{\varDelta t}\,(S_n-I)(\,f),\quad \mu_n(\,f)=\tfrac{1}{2}(S_n+I)(\,f). \end{equation*} All of these operators commute with one another. Discretizing (2.1) by means of a suitable finite difference approximation for the derivatives of the dependent variable, one obtains a partial difference equation (P|$\varDelta $|E), \begin{equation} \widetilde{\mathcal{A}}(m,n,[u])=0. \end{equation} (2.4) Here |$[u]$| denotes |$u_{0,0}$| and a finite number of its shifts; more generally, square brackets around a difference expression denote the expression and finitely many of its shifts. We seek schemes with the following finite difference analogue of each preserved conservation law: \begin{equation} \mbox{Div}\,\widetilde{\textbf{F}}\equiv D_m\left(\,\widetilde{F}(m,n,[u])\right)+ D_n\left(\,\widetilde{G}(m,n,[u])\right)=0\quad\mbox{when}\quad [\,\widetilde{\mathcal{A}}=0\,], \end{equation} (2.5) where tildes represent discretizations of the corresponding continuous terms. The functions |$\widetilde{F}$| and |$\widetilde{G}$| are (respectively) the flux and the density of the conservation law (2.5). Just as in the continuous case a conservation law of (2.4) is trivial of the first kind if |$\widetilde{F}$| and |$\widetilde{G}$| vanish on solutions of (2.4), and trivial of the second kind if (2.5) is identically satisfied without any reference to (2.4) and its shifts (see Hydon, 2014). A difference conservation law is trivial if and only if it is a linear combination of trivial conservation laws of these two kinds. Two conservation laws are equivalent if they differ by a trivial conservation law. A conservation law of (2.4) is in characteristic form if \begin{equation*}\mbox{Div}\,\widetilde{\mathbf{F}}=\widetilde{\mathcal{Q}}(m,n,[u])\widetilde{\mathcal{A}}.\end{equation*} Here |$\widetilde{\mathcal{Q}}$| is the characteristic, which is trivial if it is zero on all solutions of (2.4); two characteristics are equivalent if their difference is a trivial characteristic. Remark 2.2 (Grant & Hydon, 2013; Hydon, 2014) P|$\varDelta $|Es that can be solved for a highest shift in one direction (such as explicit P|$\varDelta $|Es) admit a one-to-one correspondence between equivalence classes of characteristics and equivalence classes of conservation laws. Therefore, characteristics can be used to test equivalence for conservation laws of such P|$\varDelta $|Es. The key result that underpins the symbolic–numeric approach is due to Kupershmidt (1985): similarly, to the continuous case, the set of all divergence expressions (2.5) over |$\mathbb{Z}^2$| is precisely the kernel of the difference Euler operator, \begin{equation} \textsf{E}\equiv\sum_{i,\,j}S_m^{-i}S_n^{-j}\frac{\partial}{\partial u_{i,\,j}}\,. \end{equation} (2.6) (See Hydon & Mansfield, 2004, for the generalization of this result.) Thus, if a function |$\widetilde{\mathcal{Q}}$| satisfies \begin{equation*}\textsf{E}(\widetilde{\mathcal{Q}}\widetilde{\mathcal{A}})\equiv 0,\end{equation*} there exists |$\widetilde{\mathbf{F}}$| such that |$\widetilde{\mathcal{Q}} \widetilde{\mathcal{A}}=\mbox{Div}\,\widetilde{\mathbf{F}}$|⁠; therefore, |$\widetilde{\mathcal{Q}}$| is the characteristic of this difference conservation law. For consistency, restrict attention to discretizations |$\widetilde{{\mathcal{Q}}}$| of |${\mathcal{Q}}$|⁠; then the difference conservation law |$\widetilde{{\mathcal{Q}}}\widetilde{{\mathcal{A}}}$| is automatically a discretization of the continuous conservation law |${\mathcal{Q}}{\mathcal{A}}$|⁠. Grant’s basic symbolic–numeric approach is straightforward. Choose a stencil of points and consider the most general discretizations on the stencil, |$\widetilde{\mathcal{A}} $| of the PDE and |$\widetilde{\mathcal{Q}}$| of the characteristic of the desired conservation law. If the stencil is large enough, there will be some free parameters in the discretizations. To preserve the conservation law impose the condition |$\textsf{E}(\widetilde{{\mathcal{Q}}}\widetilde{{\mathcal{A}}})=0$|⁠. This condition amounts to a system of algebraic equations that express constraints on the parameters. The procedure can be iterated for multiple characteristics, |${\mathcal{Q}}_l$|⁠, provided that the corresponding system of algebraic equations admits a solution. Finally, consistency conditions are applied to ensure that |$\widetilde{{\mathcal{A}}}$| converges to |${\mathcal{A}}$|⁠, and each |$\widetilde{{\mathcal{Q}}}_l$| converges to |${\mathcal{Q}}_l$| as the stepsizes |$\varDelta t$| and |$\varDelta x$| tend to zero; these give further constraints on the free parameters. In this way bespoke finite difference schemes for a given PDE may be derived by symbolic computation. In more detail the basic method is as follows. Having chosen a stencil, the most general discretizations of the PDE (2.1) and the characteristics are based on Taylor series expansions of the grid function about the point |$(x(m),t(n))\equiv (x_0,t_0)$|⁠: \begin{align}\nonumber u_{i,\,j}&\approx u(x_i,t_j)=u(x_0+i\varDelta x,t_0+j\varDelta t)\\ &=u+i\varDelta x\,u_x+j\varDelta t\,u_t+\frac{(i\varDelta x)^2}{2!}u_{xx}+\frac{(\,j\varDelta t)^2}{2!}u_{tt}+i\varDelta x\,\,j\varDelta t\,u_{xt}+\cdots\Big\vert_{(x_0,t_0)}. \end{align} (2.7) For a rectangular stencil of points defined by |$i=A,\ldots , B$| and |$j=C,\ldots , D$|⁠, linear terms in |${\mathcal{A}}$| and |${\mathcal{Q}}$| are approximated by linear combinations, with undetermined coefficients, of terms of the form (2.7): \begin{equation} \frac{\partial^{r+s}}{\partial x^r\partial t^s}u\approx\frac{1}{\varDelta x^r}\frac{1}{\varDelta t^s}\sum_{i=A}^B\sum_{j=C}^D\alpha_{i,\,j}u_{i,\,j}, \end{equation} (2.8) where the coefficients |$\alpha _{i,\,j}$| depend on |$r$| and |$s$|⁠. If quadratic terms appear in |${\mathcal{A}}$| or in |${\mathcal{Q}}$| (as happens in our examples) we need to look at products of Taylor expansions: \begin{align}\nonumber u_{i,\,j}u_{k,l}=&\,u^2+(i+k)\varDelta x\,uu_x+(\,j+l\,)\varDelta t\,uu_t+ik\varDelta x^2u_x^2+(\,jk+il\,)\varDelta x\varDelta t\,u_xu_t+jl\varDelta t^2u_t^2\\ &+\frac{\varDelta x^2(i^2+k^2)}{2!}uu_{xx}+(ij+kl\,)\varDelta x\,\varDelta t\,uu_{xt}+\frac{\varDelta t^2(\,j^2+l^2)}{2!}uu_{tt}+\cdots\Big\vert_{(x_0,t_0)}. \end{align} (2.9) Just as for the linear terms quadratic quantities in |${\mathcal{A}}$| and |${\mathcal{Q}}$| are replaced by linear combinations, with undetermined coefficients, of terms of the form (2.9): \begin{equation} \frac{\partial^{m+n}u}{\partial x^m\partial t^n}\frac{\partial^{r+s}u}{\partial x^r\partial t^s}\approx\frac{1}{\varDelta x^{m+r}\varDelta t^{n+s}}\sum_{j=C}^D\sum_{i=A}^B\left(\sum_{k=i}^B\beta_{i,\,\,j,k,\,j}u_{i,\,j}u_{k,\,j}+\sum_{l=j+1}^D\sum_{k=A}^B\beta_{i,\,j,k,l}u_{i,\,j}u_{k,l}\right), \end{equation} (2.10) with the coefficients |$\beta _{i,\,j,k,l}$| depending on |$r$| and |$s$|⁠. For the right-hand sides of (2.8) and (2.10) to approximate the corresponding left-hand sides, we also need to impose a number of consistency conditions on the coefficients |$\alpha _{i,\,j}$| and |$\beta _{i,\,j,k,l}$|⁠. Remark 2.3 Terms involving higher powers of |$u_{i,\,j}$| may be added, provided that they vanish as |$\varDelta x$| and |$\varDelta t$| tend to zero. For simplicity, such terms are not included here. Having set |$\widetilde{{\mathcal{Q}}}$| and |$\widetilde{{\mathcal{A}}}$| to be the discretizations of |${\mathcal{Q}}$| and |${\mathcal{A,}}$| respectively, one must solve \begin{equation} \textsf{E}(\widetilde{{\mathcal{Q}}}\widetilde{{\mathcal{A}}})=0, \end{equation} (2.11) where |$\textsf{E}$| is the difference Euler operator (2.6). In general, this is not easy. Typically, even for a PDE (2.1) that is only quadratic in |$[u]$|⁠, (2.11) amounts to a large system of nonlinear algebraic equations (see Grant, 2011, 2015). In principle, such systems can be solved by finding a Groebner basis (Cox et al., 1992; Mansfield, 1992; Buchberger & Kauers, 2010, 2011). However, the calculation of the Groebner basis may take a huge amount of memory and a very long computation time. As the use of this approach is limited mainly by the cost of the symbolic computation, it is helpful to impose some additional assumptions. For instance, Grant (2015) describes some symmetry-based ansätze that can simplify the Groebner basis calculation. There is one other approach to constructing finite difference schemes that can preserve multiple conservation laws for systems of PDEs. This is the multiplier method, introduced in Wan et al. (2016), which is as follows. Given a scalar PDE, \begin{equation} {\mathcal{A}} =0, \end{equation} (2.12) that has one conservation law in the form \begin{equation} D_x F+ D_t G={\mathcal{Q}}{\mathcal{A}}, \end{equation} (2.13) let |$\widetilde{{\mathcal{Q}}}$|⁠, |$\widetilde{F}$| and |$\widetilde{G}$| be finite difference approximations of |${\mathcal{Q}}$|⁠, |$F$| and |$G$|⁠, respectively. Then, provided that |$\widetilde{{\mathcal{Q}}}^{-1}$| exists on the whole domain of definition of (2.12), \begin{equation} \widetilde{{\mathcal{A}}}\equiv\widetilde{{\mathcal{Q}}}^{-1}\left\{D_m\left(\,\widetilde{F}\,\right)+D_n\left(\,\widetilde{G}\,\right)\right\}=0 \end{equation} (2.14) is a finite difference approximation of (2.12) that preserves the conservation law (2.13). This method has been applied to find conservative schemes for the inviscid Burgers’ equation and momentum-preserving schemes for the KdV equation. It has also been applied to systems of |$r$| PDEs, such as the two-dimensional shallow-water equations, to preserve |$s$| conservation laws, with |$s\leqslant r$|⁠. The multiplier method has the advantage of being simple to use. However, there are two main disadvantages. First, it cannot find schemes preserving |$s$| conservation laws for systems of |$r$| PDEs with |$s>r$|⁠; in particular, it cannot preserve multiple conservation laws for a scalar PDE. Secondly, it requires the characteristic to be nonzero throughout the domain. For characteristics that involve the dependent variable, one cannot identify points where the characteristic is zero a priori. By contrast, Grant’s approach is in principle able to find all conservative finite difference methods on the chosen stencil, provided that one is able to solve (2.11). The procedure can be iterated to select discretizations that preserve further conservation laws, provided that the corresponding condition (2.11) admits a solution for each characteristic. Moreover, there is no need to choose a particular and arbitrary discretization of densities and fluxes, as these can be reconstructed from the characteristics (see Hydon, 2001). To simplify Grant’s approach we adopt a strategy that reduces the number of variables and the computational cost of solving the system of nonlinear equations. This is achieved by first looking for second-order accurate approximations only, building in consistency from the outset. If the stencil is as compact as possible, this immediately determines the discretizations of the highest-order derivatives. The problem can be further simplified by restricting the approximations of some terms in |$\widetilde{{\mathcal{A}}}$| and |$\widetilde{{\mathcal{Q}}}$| to use only points in a sub-stencil that is as compact as possible. In particular, by approximating nonlinear terms using as few points as possible, the number of variables may be considerably reduced to the point of being able to solve (2.11) with a fast symbolic computation that does not need a Groebner basis. Remark 2.4 Conservation laws of a given PDE are evaluated on hypersurfaces. In the discrete case, the smallest ‘surface’ on which a conservation law can be evaluated locally is the convex hull of the stencil (Hydon, 2014). Therefore, for all conservative schemes that are presented in this paper, the consistency conditions are imposed so as to approximate the conservation laws (and their characteristics) to second-order accuracy at the centre of the rectangular stencil. Remark 2.5 Obtaining second-order accurate approximations of the conservation laws at the centre |$(x,t)$| of the stencil is equivalent to finding second-order accurate approximations of the corresponding densities and fluxes at the points |$(x,t-\varDelta t/2)$| and |$(x- \varDelta x/2,t),$| respectively (see (2.5)). Figure 1 shows an example of a rectangular stencil. The circle denotes the centre, i.e., the point where we require second-order approximations of the PDE and of the characteristics (and hence of the conservation laws). The crosses denote the points where we require second-order approximations of the corresponding densities and fluxes. These are not necessarily lattice points. Fig. 1. View largeDownload slide Example of a rectangular stencil. Conservation laws are preserved to second-order accuracy at the central point |$(x, t)$| (circle), densities and fluxes, respectively, at |$(x, t-\varDelta t/2) $| and |$(x-\varDelta x/2, t)$| (crosses). Fig. 1. View largeDownload slide Example of a rectangular stencil. Conservation laws are preserved to second-order accuracy at the central point |$(x, t)$| (circle), densities and fluxes, respectively, at |$(x, t-\varDelta t/2) $| and |$(x-\varDelta x/2, t)$| (crosses). In the next two sections, we use the above simplifications to derive approximations that preserve two conservation laws of some well-known nonlinear wave equations, and show that these conservative methods can generate robust and highly accurate schemes. 3. KdV equation In this section we exploit the strategy introduced in Section 2 to develop conservative schemes for the KdV equation, \begin{equation} \mathcal{A} \equiv u_t+uu_x+u_{xxx}=0, \quad (x,t)\in\varOmega\equiv[a,b]\times[0,\infty). \end{equation} (3.1) These schemes are tested for two benchmark problems; they compare favourably with two well-known schemes that each preserve only one conservation law. Equation (3.1) has an infinite number of conservation laws. The first three, in increasing order, are \begin{align} D_t(G_1)+D_x(F_1) \equiv D_t(u)+D_x\left(\frac{1}2 u^2 + u_{xx}\right)=0, \end{align} (3.2) \begin{align} D_t(G_2)+D_x(F_2) \equiv D_t\left(\frac{1}2 u^2\right)+D_x\left(\frac{1}3u^3+uu_{xx}-\frac{1}2u_x^2\right)=0, \end{align} (3.3) \begin{align} & D_t(G_3)+D_x(F_3) \equiv D_t\left(\frac{u^3}3+uu_{xx}\right)+D_x\left(\frac{u^4}4+u_xu_t-uu_{xt}+u^2u_{xx}+u_{xx}^2\right)=0, \end{align} (3.4) which can be written in characteristic form (2.3) with characteristics \begin{equation} {\mathcal{Q}}_1=1,\qquad{\mathcal{Q}}_2=u,\qquad{\mathcal{Q}}_3=u^2+2u_{xx}, \end{equation} (3.5) respectively. For a water wave problem, conservation laws (3.2)–(3.4) describe the local conservation of mass, momentum and energy, respectively (see Drazin & Johnson, 1989). As these conservation laws have a physical meaning it seems particularly desirable to preserve them. When (3.1) is coupled with suitable (e.g., periodic or zero) boundary conditions, integrating (3.2)–(3.4) over the spatial domain gives the global conservation of, respectively, \begin{equation} \int G_1\,\textrm{d} x=\int u\,\textrm{d} x,\qquad \int G_2\,\textrm{d} x=\int \frac{1}{2}u^2\,\textrm{d} x, \qquad \int G_3\,\textrm{d} x=\int \frac{1}{3}u^3+uu_{xx}\, \textrm{d} x. \end{equation} (3.6) It is well known that (3.1) possesses the Hamiltonian structure \begin{equation*} u_t=D_x\frac{\delta}{\delta u}{\mathcal{H}}_1, \end{equation*} where |${\delta }/{\delta u}$| is the variational derivative and \begin{equation} {\mathcal{H}}_1=\int-\frac{1}2\left(\frac{1}3u^3-u_x^2\right)\textrm{d} x \end{equation} (3.7) is the Hamiltonian functional. Equivalently, one can use the alternative Hamiltonian functional \begin{equation} {\mathcal{H}}_1=\int-\frac{1}2\left(\frac{1}3u^3+uu_{xx}\right)\textrm{d} x=-\frac{1}{2}\int G_3\,\textrm{d} x. \end{equation} (3.8) With this choice of functional the conservation law (3.4) implies the preservation of |${\mathcal{H}}_1$|⁠. The KdV equation (3.1) can also be written in another Hamiltonian form (see Olver, 1993), \begin{equation*}u_t=\left(D_x^3+\frac{2}{3}uD_x+\frac{1}{3}u_x\right)\frac{\delta}{\delta u}{\mathcal{H}}_2,\end{equation*} with the Hamiltonian \begin{equation} {\mathcal{H}}_2=\int-\frac{1}{2}u^2\,\textrm{d}x=-\int G_2\,\textrm{d}x. \end{equation} (3.9) The conservation law (3.3) implies that |${\mathcal{H}}_2$| is preserved. These are special instances of the following well-known general result. Given a scalar Hamiltonian evolution equation for |$u(x,t)$|⁠, \begin{equation} u_t=\mathcal{D}\frac{\delta}{\delta u}{\mathcal{H}},\qquad{\mathcal{H}}=\int H([u]) \,\textrm{d}x, \end{equation} (3.10) where |$\mathcal{D}$| is a skew-adjoint differential operator (with respect to the |$L^2$| inner product), the Hamiltonian |${\mathcal{H}}$| is constant (provided that some technical conditions are satisfied). To discretize the KdV equation (3.1) one first needs to set the stencil. Having done this we will use our simplified version of Grant’s approach to construct two types of scheme: our energy-conserving schemes preserve discrete versions of the conservation laws (3.2) and (3.4), while our momentum-conserving schemes preserve discrete versions of (3.2) and (3.3). 3.1 Conservative methods for the KdV equation Eight-point schemes The most compact rectangular stencil consists of eight points, as shown in Fig. 2. Here and henceforth grid points are labelled with respect to the lattice point denoted with a square. From Remark 2.5 we seek second-order approximations of characteristics, densities and fluxes at |$(-1/2,1/2)$|⁠, |$(-1/2,0)$| and |$(-1,1/2)$|⁠, respectively. Fig. 2. View largeDownload slide The most compact rectangular stencil for (3.1). Conservation laws are preserved to second-order accuracy at the central point |$(-1/2, 1/2)$| (circle), densities and fluxes respectively at |$(-1/2,0) $| and |$(-1, 1/2)$| (crosses). Fig. 2. View largeDownload slide The most compact rectangular stencil for (3.1). Conservation laws are preserved to second-order accuracy at the central point |$(-1/2, 1/2)$| (circle), densities and fluxes respectively at |$(-1/2,0) $| and |$(-1, 1/2)$| (crosses). Energy-conserving schemes To simplify the symbolic computations for energy-conserving schemes on the eight-point stencil we use the approximations \begin{equation*} {u_{xx}}\approx\mu_n\left(D_m^2u_{-2,0}\right),\qquad{u_{xx}}\approx\mu_m\left(\mu_n\left(D_m^2u_{-2,0}\right)\right), \end{equation*} in |$\widetilde{F_1}$| and |$\widetilde{Q_3}$|⁠, respectively. The remaining terms in the approximations of |$G_1$| at |$(-1/2,0)$|⁠, |$F_1$| at |$(-1,1/2)$| and |${\mathcal{Q}}_3$| at |$(-1/2,1/2)$| are obtained from (2.8) and (2.10) by requiring that the coefficients |$\alpha _{i,\,j}$| and |$\beta _{i,\,j,k,l}$| satisfy all consistency conditions for second-order accuracy. This yields families of approximations that depend on just a few undetermined coefficients. The discretizations of the conservation laws (3.2) and (3.4) at |$(-1/2,1/2)$| are taken to be of the following form (setting |$\widetilde{{\mathcal{Q}}_1}=1$|⁠): \begin{equation*} \widetilde{{\mathcal{A}}}=D_m\left(\widetilde{F_1}\right)+D_n\left(\widetilde{G_1}\right)=0,\qquad\qquad \widetilde{{\mathcal{Q}}_3}\widetilde{{\mathcal{A}}}=0. \end{equation*} As |$\widetilde{{\mathcal{A}}}$| is defined to be a discrete conservation law, the condition |$\textsf{E}(\widetilde{{\mathcal{A}}})\equiv 0$| holds for any choice of the remaining coefficients in |$\widetilde{G_1}, \widetilde{F_1}$| and |$\widetilde{Q_3}$|⁠. We then find these undetermined coefficients by solving \begin{equation*} \textsf{E}(\widetilde{{\mathcal{Q}}}_3\widetilde{{\mathcal{A}}})\equiv 0. \end{equation*} This constraint determines all remaining coefficients. So only one scheme of this form preserves (3.2) and (3.4) to second-order accuracy at the centre of the stencil: \begin{equation} \mbox{EC}_8\equiv\widetilde{{\mathcal{A}}}=D_m\left(\widetilde{F_1}\right)+D_n\left(\widetilde{G_1}\right)=0, \end{equation} (3.11) where \begin{equation} \widetilde{F_1}=\frac{1}{3}\,\mu_n\left(\left(\mu_m^2u_{-2,0}\right)^2\right)+\frac{1}{6}\left(\mu_m^2u_{-2,0}\right)\left(\mu_m^2u_{-2,1}\right)+\mu_nD_m^2u_{-2,0},\qquad \widetilde{G_1}=\mu_mu_{-1,0}. \end{equation} (3.12) The scheme |$\mbox{EC}_8$| preserves the following discrete version of the conservation law (3.4): \begin{align} \widetilde{{\mathcal{Q}}_3}\widetilde{{\mathcal{A}}}=& D_m\widetilde{F_3}+ D_n\widetilde{G_3}=0, \end{align} (3.13) where \begin{align}\nonumber \widetilde{{\mathcal{Q}}_3}=&\,2\mu_m\widetilde{F_1},\\ \nonumber \widetilde{F_3}=&\,\widetilde{F_1}^2+\left(\mu_nD_m\mu_mu_{-2,0}\right)\left(D_n\mu_m^2u_{-2,0}\right)-\left(\mu_n\mu_m^2u_{-2,0}\right)\left(D_nD_mu_{-2,0}\right)\\ \nonumber &\,+\frac{\varDelta x^2}{6}(\mu_n\mu_m^2u_{-2,0})\left\{(2\mu_mD_nu_{-1,0})(\mu_nD_m\mu_mu_{-2,0})-D_n(\mu_mu_{-1,0}D_m\mu_mu_{-2,0})\right\}\!,\\ \nonumber \widetilde{G_3}=&\,\frac{1}{3}{(\mu_mu_{-1,0})\mu_m\left(\left(\mu_m^2u_{-2,0}\right)^2\right)}+(\mu_mu_{-1,0})(\mu_mD_m^2u_{-2,0}). \end{align} The last term in the flux |$\widetilde{F_3}$| vanishes as the spatial stepsize tends to zero and does not correspond to an expression in the continuous flux. The scheme |$\mbox{EC}_8$| is equivalent to one introduced in Grant (2015). When (3.1) is coupled with zero or periodic boundary conditions, the scheme |$\mbox{EC}_8$| preserves at each time step the following discretization of the Hamiltonian (3.8): \begin{equation} \widetilde{{\mathcal{H}}_1}(\,j)=-\frac{\varDelta x}2\sum_i\left(\frac{1}{3}{(\mu_mu_{i-1,\,j})\mu_m\left(\left(\mu_m^2u_{i-2,\,j}\right)^2\right)}+(\mu_mu_{i-1,\,j})\big(\mu_mD_m^2u_{i-2,\,j}\big)\right). \end{equation} (3.14) In general, |$\mbox{EC}_8$| fails to preserve the conservation law (3.3), even to first order: given any approximation, \begin{equation*}\widetilde{{\mathcal{Q}}_2}=\sum_{i=-2}^1\sum_{j=0}^1\eta_{i,\,j}u_{i,\,j},\qquad \widetilde{{\mathcal{Q}}_2}= {\mathcal{Q}}_2 +{\mathcal{O}}(\varDelta x,\varDelta t),\end{equation*} the condition |$E(\widetilde{{\mathcal{Q}}_2}\widetilde{{\mathcal{A}}})=0$| cannot be satisfied when |$\widetilde{{\mathcal{A}}}$| is given by (3.11). Momentum-conserving schemes Grant (2015) introduced several momentum-conserving schemes, including a one-parameter family obtained by using the most compact second-order approximation of |$u_t$| in (3.1). This amounts to \begin{equation*} \mbox{MC}_8(\alpha)\equiv\widetilde{{\mathcal{A}}}=D_m\left(\widetilde{F_1}\right)+D_n\left(\widetilde{G_1}\right)=0, \end{equation*} where \begin{align*}\nonumber \widetilde{F_1}=&\,\frac{1}{6}(\mu_nu_{-1,0})\mu_n(u_{-2,0}+u_{-1,0}+u_{0,0})+\mu_nD_m^2u_{-2,0}\\\nonumber &\,+\alpha\varDelta x^2\left\lbrace(\mu_nu_{0,0})\big(D_m^2\mu_nu_{-2,0}\big)+D_m((\mu_nu_{-2,0})(D_m\mu_nu_{-2,0}))\right\rbrace\!,\\ \widetilde{G_1}=&\,\mu_mu_{-1,0}. \end{align*} For any value of |$\alpha $| these methods preserve the discrete momentum conservation law \begin{equation*}\widetilde{{\mathcal{Q}}_2}\widetilde{{\mathcal{A}}}= D_m\left(\widetilde{F_2}\right)+ D_n\left(\widetilde{G_2}\right)=0,\end{equation*} with \begin{align*} \widetilde{{\mathcal{Q}}_2}=&\,\mu_m\mu_nu_{-1,0}\,,\\ \widetilde{F_2}=&\, \frac{1}{3}(\mu_m\mu_nu_{-2,0})(\mu_m\mu_nu_{-1,0})\mu_n(u_{-1,0}+2\alpha\varDelta x^2D_m^2u_{-2,0})\\ &\,+\big(\mu_n\mu_m^2u_{-2,0}\big)\big(\mu_nD_m^2u_{-2,0}\big)-\frac{1}{2}\mu_m\left((\mu_nD_mu_{-2,0})^2\right)\!,\\ \widetilde{G_2}=&\,\frac{1}{2}(\mu_mu_{-1,0})^2. \end{align*} For zero or periodic boundary conditions these schemes preserve the following discretization of the Hamiltonian (3.9) at each time step: \begin{equation} \widetilde{{\mathcal{H}}_2}(\,j)=-\frac{\varDelta x}2\sum_i\left(\mu_mu_{i-1,\,j}\right)^2. \end{equation} (3.15) The local truncation error of the scheme |$\mbox{MC}_8(\alpha )$| is |${\mathcal{O}}(\varDelta x^2)+{\mathcal{O}}(\varDelta t^2)$|⁠. Restricting attention to the case |$\varDelta t \ll \varDelta x$| the truncation error can be reduced considerably by choosing |$\alpha $| optimally. No choice of |$\alpha $| eliminates the second-order terms identically, so the optimal value will depend on the particular problem. Ten-point schemes To find new schemes that preserve two conservation laws one must use a wider stencil. This is beyond what can be tackled in full generality, but the symbolic computations are made tractable (indeed, fast) by the simplifications that we have introduced. Adding one further pair of nodes in the spatial direction gives the 10-point stencil in Fig. 3. Hence, according to Remark 2.5, our goal is to find second-order approximations of characteristics, densities and fluxes at |$(0,1/2)$|⁠, |$(0,0)$| and |$(-1/2,1/2)$|⁠, respectively. Fig. 3. View largeDownload slide The 10-point rectangular stencil for (3.1). Conservation laws are preserved with higher order at the central point |$(0, 1/2)$| (circle), densities and fluxes, respectively, at |$(0,0) $| and |$(-1/2, 1/2)$| (crosses). Fig. 3. View largeDownload slide The 10-point rectangular stencil for (3.1). Conservation laws are preserved with higher order at the central point |$(0, 1/2)$| (circle), densities and fluxes, respectively, at |$(0,0) $| and |$(-1/2, 1/2)$| (crosses). Energy-conserving schemes Just as for the eight-point schemes let \begin{equation*} \widetilde{\mathcal{A}}=D_m\left(\widetilde{F_1}\right)+D_n\left(\widetilde{G_1}\right)\!, \end{equation*} so that (3.2) is preserved for any choice of the undetermined coefficients. The preservation of (3.4), obtained by requiring that \begin{equation} \textsf{E}\left(\widetilde{{\mathcal{Q}}_3}\widetilde{{\mathcal{A}}}\,\right)\equiv 0, \end{equation} (3.16) is simplified by setting the sub-stencils for |$\widetilde{G_1}$| and the quadratic term in |$\widetilde{{\mathcal{Q}}_3}$| to be as compact as possible, given that the approximations must be second order: \begin{align} \widetilde{G_1}&=u_{0,0},\\ \nonumber \widetilde{{\mathcal{Q}}_3}&=\xi\left(u_{0,0}^2+u_{0,1}^2\right)+(1-2\xi)u_{0,0}u_{0,1}+\frac{1}{\varDelta x^2}\sum_{i=-2}^2\sum_{j=0}^1\gamma_{i,\,j}u_{i,\,j},\quad \xi\in\mathbb{R}. \end{align} (3.17) The undetermined coefficients in |$\widetilde{F_1}$| and |$\widetilde{{\mathcal{Q}}_3}$| are obtained by solving (3.16). This yields a one-parameter family of schemes, \begin{equation*} \widetilde{{\mathcal{A}}}(\lambda)=D_m\left(\widetilde{F_1}\right)+D_n\left(\widetilde{G_1}\right)=0, \end{equation*} with \begin{equation*} \widetilde{F_1}=\mu_m\varphi_{-1,0},\qquad \widetilde{G_1}=u_{0,0}, \end{equation*} where \begin{equation*}\varphi_{-1,0}=\frac{u_{-1,1}^2+u_{-1,0}^2+u_{-1,0}u_{-1,1}}{6}+D_m^2\mu_nu_{-2,0}+\lambda D_nD_m\mu_mu_{-2,0},\end{equation*} and |$\lambda ={\mathcal{O}}(\varDelta x^2,\varDelta t^2)$|⁠. These schemes preserve \begin{align*} \widetilde{{\mathcal{Q}}_3}\widetilde{{\mathcal{A}}}=& D_m\widetilde{F_3}+ D_n\widetilde{G_3}=0, \end{align*} where \begin{align*} \widetilde{{\mathcal{Q}}_3}=&\,2\varphi_{0,0}\,,\\[2pt]\nonumber \widetilde{F_3}=&\,\varphi_{0,0}\varphi_{-1,0}+(D_m\mu_nu_{-1,0})(D_n\mu_mu_{-1,0})-(\mu_m\mu_nu_{-1,0})(D_nD_mu_{-1,0})\\[2pt] &\,+\lambda(D_nu_{0,0})(D_nu_{-1,0}),\\[2pt]\nonumber \widetilde{G_3}=&\,\frac{1}{3}u_{0,0}^3+u_{0,0}D_m^2u_{-1,0}\,. \end{align*} For zero or periodic boundary conditions these schemes preserve at each time step \begin{equation} \widetilde{{\mathcal{H}}_1}(\,j)=-\frac{\varDelta x}2\sum_i\left(\frac{1}3u_{i,\,j}^3+u_{i,\,j}\,D_m^2u_{i-1,\,j}\right), \end{equation} (3.18) but none of them preserves the conservation law (3.3). Assuming for simplicity that |$\varDelta t\ll \varDelta x$|⁠, the leading term in the local truncation error amounts to |${\mathcal{O}}(\varDelta x^2)$|⁠. This suggests that by setting |$\lambda =\alpha \varDelta x^2$| one may be able to remove at least part of this error by choosing |$\alpha \in \mathbb{R}$| optimally. However, Taylor expansion shows that no choice of |$\alpha $| will give a higher order method. Indeed, the optimal value depends on the initial conditions. In the results section, we write the one-parameter family of schemes as \begin{equation*}\mbox{EC}_{10}(\alpha)\equiv\widetilde{{\mathcal{A}}}(\alpha\varDelta x^2).\end{equation*} The scheme |$\mbox{EC}_{10}(0)$| was originally found by the discrete variational derivative method (see Furihata, 1999). Dahlby & Owren (2011) proved that the Furihata scheme can also be derived by the average vector field method, which approximates (3.10) by |$D_nu_{0,0}=\widetilde{\mathcal{D}}(\widetilde{\delta }\widehat{\mathcal{H}})$|⁠, where the operator |$\widetilde{\mathcal{D}}$| is skew-adjoint with respect to the |$\ell ^2$| inner product. In this case, \begin{equation*} \widetilde{\mathcal{D}}=D_m\mu_mS_m^{-1},\qquad \widetilde{\delta}\widehat{\mathcal{H}}=-\frac{1}6\left\{u_{0,1}^2+u_{0,1}u_{0,0}+u_{0,0}^2\right\}-D_m^2\mu_nu_{-1,0}\,\!. \end{equation*} None of the other |$\mbox{EC}_{10}(\alpha )$| schemes can be derived in this way. Momentum-conserving schemes To simplify the derivation of momentum-conserving schemes on the 10-point stencil in Fig. 3, use the following approximations in which |$\widetilde{G_1}$|⁠, |$\widetilde{{\mathcal{Q}}_2}$| and the quadratic term in |$\widetilde{F_1}$| are compact:1 \begin{align} \widetilde{G_1}&\,=u_{0,0}, \end{align} (3.19) \begin{align} \widetilde{{\mathcal{Q}}_2}&\,=\mu_nu_{0,0},\\ \nonumber \widetilde{F_1}&\,=\sum_{j=0}^1\sum_{i=-1}^0\sum_{k=i}^0\beta_{i,\,j,k,\,j}u_{i,\,j}u_{k,\,j}+\sum_{i=-1}^0\sum_{k=-1}^0\beta_{i,0,k,1}u_{i,0}u_{k,1}+\frac{1}{\varDelta x^2}\sum_{i=-2}^1\sum_{j=0}^1\alpha_{i,\,j}u_{i,\,j}. \end{align} (3.20) Proceeding as before one obtains a two-parameter family of momentum-conserving methods: \begin{equation*} \widetilde{\mathcal{A}}(\lambda,\nu)=D_m\left(\widetilde{F_1}\right)+D_n\left(\widetilde{G_1}\right)=0, \end{equation*} with |$\widetilde{G_1}$| as defined in (3.19) and \begin{equation*} \widetilde{F_1}=\,\frac{(\mu_nu_{-1,0})^2+(\mu_nu_{0,0})^2+(\mu_nu_{-1,0})(\mu_nu_{0,0})}6+\mu_nD_m^2\mu_mu_{-2,0} +D_nD_m(\lambda u_{-1,0}+\nu D_m^2u_{-2,0}); \end{equation*} here |$\lambda ={\mathcal{O}}(\varDelta x^2,\varDelta t^2)$| and |$\nu ={\mathcal{O}}(\varDelta x^2,\varDelta t^2)$|⁠. These schemes preserve \begin{equation*} \widetilde{{\mathcal{Q}}_2}\widetilde{{\mathcal{A}}}= D_m\widetilde{F_2}+ D_n\widetilde{G_2}=0, \end{equation*} with |$\widetilde{{\mathcal{Q}}_2}$| as given in (3.20) and \begin{align*} \widetilde{F_2}=&\,\frac{1}{3}(\mu_nu_{0,0})(\mu_nu_{-1,0})(\mu_n\mu_mu_{-1,0}) +\frac{1}2\left\lbrace(\mu_nu_{0,0})\big(D_m^2\mu_nu_{-2,0}\big)+(\mu_nu_{-1,0})\big(D_m^2\mu_nu_{-1,0}\big)\right\rbrace\\ &\!-\frac{1}2(D_m\mu_nu_{-1,0})^2+\lambda\left\lbrace(\mu_m\mu_nu_{-1,0})(D_nD_mu_{-1,0})-\frac{1}{2}D_n(\mu_mu_{-1,0}D_mu_{-1,0})\right\rbrace\\\nonumber &\!+\nu\Big\lbrace\!(\mu_nu_{0,0})\big(D_nD_m^3u_{-2,0}\big)\!-\!(D_m\mu_nu_{-1,0})\big(D_nD_m^2u_{-1,0}\big)\!+\!\frac{1}2D_n\!\left[\!D_mu_{-1,0}D_m^2u_{-1,0}-u_{0,0}D_m^3u_{-2,0}\!\right]\!\!\Big\rbrace\!,\\ \widetilde{G_2}=&\,\frac{1}2u_{0,0}^2+\frac{1}{2}u_{0,0}D_m^2\big(\lambda u_{-1,0}+\nu D_m^2u_{-2,0}\big). \end{align*} For suitable boundary conditions the |$\widetilde{{\mathcal{A}}}(\lambda ,\nu )$| scheme preserves at each time step \begin{equation} \widetilde{{\mathcal{H}}_2}(\,j)=-\frac{\varDelta x}2\sum_i\left({u_{i,\,j}^2}+u_{i,\,j}D_m^2\left(\lambda u_{i-1,\,j}+\nu D_m^2u_{i-2,\,j}\right)\right). \end{equation} (3.21) To simplify the main sources of local truncation error we will consider only |$\varDelta t\ll \varDelta x$|⁠; then the leading term is |${\mathcal{O}}(\varDelta x^2)$|⁠, so we define the two-parameter family of schemes \begin{equation*}\mbox{MC}_{10}(\alpha,\beta)\equiv\widetilde{\mathcal{A}}(\alpha\varDelta x^2,\beta\varDelta x^2).\end{equation*} Again, it is not possible to obtain higher-order methods for any choice of the free parameters; the optimal values depend on the particular problem. It turns out that the |$\mbox{MC}_{10}(0,0 )$| scheme can also be derived by the average vector field method, approximating the right-hand side of (3.10) by \begin{equation*} \widetilde{\mathcal{D}}\widetilde{\delta}\widehat{\mathcal{H}}=D_m^3S_m^{-2}\mu_m(\widetilde{\delta}\widehat{\mathcal{H}})+\frac{2}3S_m^{-1}\mu_m\left\{(\mu_m\mu_nu_{0,0})D_m(\widetilde{\delta}\widehat{\mathcal{H}})\right\}+\frac{1}3D_m(\mu_m\mu_nu_{-1,0})\,\widetilde{\delta}\widehat{\mathcal{H}}, \end{equation*} where \begin{equation*}\widetilde{\delta}\widehat{\mathcal{H}}=-\mu_n \, u_{0,0}. \end{equation*} In this case, the skew-adjoint operator |$\widetilde{\mathcal{D}}$| is not constant. 3.2 Numerical tests In this subsection, two benchmark solutions are used to illustrate the effectiveness of the schemes developed in Section 3.1 by comparison with two well-known schemes that each preserve a single conservation law. These are the multisymplectic scheme proposed in Ascher & McLachlan (2004, 2005), which we rewrite as \begin{equation} D_m\left(\frac{1}{2}\mu_m(\mu_m\mu_nu_{-2,0})^2+D_m^2\mu_nu_{-2,0}\right)+D_n\left(\mu_m^3 u_{-2,0}\right)=0 \end{equation} (3.22) and the narrow box scheme, defined in the same references, which amounts to \begin{equation} D_m\left(\frac{1}{2}(\mu_n u_{-1,0})^2+D_m^2\mu_nu_{-2,0}\right)+D_n\left(\mu_m u_{-1,0}\right)=0. \end{equation} (3.23) Both the multisymplectic scheme (3.22) and the narrow box scheme (3.23) are defined on the eight-point stencil in Fig. 2, and preserve a discrete version of the mass conservation law (3.2). Each scheme considered in this section is solved by using the Newton method, simplified by using a ‘frozen’ Jacobian. This procedure is computationally attractive because the inversion of the Jacobian is performed just once for a single instance of the iterative method. The iterations are run until the error reaches full machine accuracy (up to rounding errors) in double precision. For each of our numerical experiments the computational cost is approximately the same for all of the schemes. In the following we consider (3.1) subject to periodic boundary conditions. We evaluate the error in the solution at the final time |$t=T$| as \begin{equation} \left.\frac{\|u-u_{\textrm{exact}}\|}{\|u_{\textrm{exact}}\|}\right\vert{}_{t=T}. \end{equation} (3.24) For a grid with |$M$| points in space and |$N$| points in time, the errors in the invariants (3.6) are \begin{equation} \textrm{Err}_\ell=\varDelta x\max_{j=1,\ldots,N}\left|\sum_{i=1}^M \left(\widetilde{G_\ell}(x_i,t_j)-\widetilde{G_\ell}(x_i,t_1)\right)\right|,\qquad \ell=1,2,3. \end{equation} (3.25) Where some of the discrete densities |$\widetilde{G_1}$|⁠, |$\widetilde{G_2}$| and |$\widetilde{G_3}$| are undefined because the considered scheme does not preserve the corresponding conservation laws, we instead evaluate the corresponding errors as \begin{align}\nonumber \textrm{Err}_1=&\,\varDelta x\max_{j=1,\ldots,N}\left|\sum_{i=1}^M (v_{i,\,j}-v_{i,1})\right|,\\ \textrm{Err}_2=&\,\varDelta x\max_{j=1,\ldots,N}\left|\frac{1}2\sum_{i=1}^M \big(v_{i,\,j}^2-v_{i,1}^2\big)\right|,\\ \nonumber \textrm{Err}_3=&\,\varDelta x\max_{j=1,\ldots,N}\left|\sum_{i=1}^M \left( \frac{v_{i,\,j}^3}3+v_{i,\,j}D_m^2(v_{i-1,\,j})-\frac{v_{i,1}^3}3-v_{i,1}D_m^2(v_{i-1,1})\right)\right|. \end{align} (3.26) Here |$v_{i,\,j}=u_{i,\,j}$| for the 10-point stencil in Fig. 3 and |$v_{i,\,j}=\mu _mu_{i-1,\,j}$| for the 8-point stencil in Fig. 2, where |$u_{i,\,j}\simeq u(a+i\varDelta x, j\varDelta t)$|⁠; subscripts denote shifts from the point |$(x,t)=(a,0)$|⁠. Note that |$\textrm{Err}_2$| shows how well each scheme preserves the corresponding discretization of the Hamiltonian |${\mathcal{H}}_2$| because \begin{equation*}\max_{j} |\widetilde{{\mathcal{H}}_2}(\,j)-\widetilde{{\mathcal{H}}_2}(1)|=\textrm{Err}_2.\end{equation*} Similarly, |$\textrm{Err}_3$| shows how well |${\mathcal{H}}_1$| is preserved because \begin{equation*}\max_{j} |\widetilde{{\mathcal{H}}_1}(\,j)-\widetilde{{\mathcal{H}}_1}(1)|=\textrm{Err}_3/2.\end{equation*} As a first numerical test we consider equation (3.1) for |$t\in [0,2]$|⁠, with periodic boundary conditions over the interval |$[-20,20]$| and the initial condition \begin{equation} u(x,0)=3c\,\textrm{sech}^2\left(\frac{\sqrt c}{2}(x+d)\right). \end{equation} (3.27) The exact solution of (3.1) with initial condition (3.27) on an infinite domain is a single soliton, \begin{equation} u_{\textrm{exact}}(x,t)=3c\,\textrm{sech}^2\left(\frac{\sqrt c}{2}(x-ct+d)\right). \end{equation} (3.28) Each scheme is solved for the parameters |$c=d=5$| with stepsizes |$\varDelta x=0.1$| and |$\varDelta t=0.01$|⁠. For this problem the solution errors for |$\mbox{MC}_8(\alpha _1)$|⁠, |$\mbox{MC}_{10}(\alpha _2,\beta _2)$| and |$\mbox{EC}_{10}(\alpha _3)$| are minimized when |$\alpha _1=-0.069$|⁠, |$(\alpha _2,\beta _2)=(0.39,0.04)$| and |$\alpha _3=0.12$|⁠. When the solution error cannot be evaluated because the exact solution is unknown, another criterion is needed to optimize the free parameters. For instance, to minimize the error in the nonpreserved conservation law, the approximate parameter values are |$\alpha _1=-0.073$|⁠, |$(\alpha _2,\beta _2)=(0.21,0.03)$| and |$\alpha _3=0.17$|⁠. Table 1 shows that the MC and EC schemes described in Section 3.1 preserve two conservation laws to machine accuracy. The most accurate of these schemes is |$\mbox{EC}_{10}(0.12)$|⁠. Minimizing the error in the nonpreserved conservation law does not optimize the numerical solution, but nevertheless yields a solution error that is comparable to (for MC|$_{10}$|⁠) or smaller than (for MC|$_{8}$| and EC|$_{10}$|⁠) the errors in the Furihata, multisymplectic and narrow box schemes. Table 1 Errors in conservation laws and solution at the final time |$T=2$|⁠, when solving (3.1) with initial condition (3.27) (for |$c=d=5$|⁠) and periodic boundary conditions over |$[-20,20]$|⁠, using various schemes with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.01$|⁠. The nonzero values of the free parameter minimize the error indicated by an asterisk Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\mbox{EC}_8$| 9.24e–14 0.0019 6.71e–12 0.0964 |$\mbox{MC}_8(0)$| 1.24e–13 6.82e–13 0.0308 0.0584 |$\mbox{MC}_8(-0.069)$| 1.24e–13 1.05e–12 0.0028 0.0052|$^*$| |$\mbox{MC}_8(-0.073)$| 8.17e–14 6.25e–13 0.0014|$^*$| 0.0063 Furihata; |$\mbox{EC}_{10}(0)$| 8.17e–14 9.93e–04 2.16e–12 0.0217 |$\mbox{EC}_{10}(0.12)$| 7.46e–14 3.40e–04 2.61e–12 0.0020|$^*$| |$\mbox{EC}_{10}(0.17)$| 7.46e–14 7.29e–05|$^*$| 2.27e–12 0.0095 |$\mbox{MC}_{10}(0,0)$| 6.75e–14 3.69e–13 0.0033 0.0335 |$\mbox{MC}_{10}(0.39,0.04)$| 7.46e–14 3.41e–13 0.0127 0.0033|$^*$| |$\mbox{MC}_{10}(0.21,0.03)$| 6.39e–14 4.26e–13 6.07e–04|$^*$| 0.0224 Multisymplectic 1.24e–13 7.03e–04 0.0436 0.0385 Narrow box 1.24e–13 0.0033 0.0325 0.0235 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\mbox{EC}_8$| 9.24e–14 0.0019 6.71e–12 0.0964 |$\mbox{MC}_8(0)$| 1.24e–13 6.82e–13 0.0308 0.0584 |$\mbox{MC}_8(-0.069)$| 1.24e–13 1.05e–12 0.0028 0.0052|$^*$| |$\mbox{MC}_8(-0.073)$| 8.17e–14 6.25e–13 0.0014|$^*$| 0.0063 Furihata; |$\mbox{EC}_{10}(0)$| 8.17e–14 9.93e–04 2.16e–12 0.0217 |$\mbox{EC}_{10}(0.12)$| 7.46e–14 3.40e–04 2.61e–12 0.0020|$^*$| |$\mbox{EC}_{10}(0.17)$| 7.46e–14 7.29e–05|$^*$| 2.27e–12 0.0095 |$\mbox{MC}_{10}(0,0)$| 6.75e–14 3.69e–13 0.0033 0.0335 |$\mbox{MC}_{10}(0.39,0.04)$| 7.46e–14 3.41e–13 0.0127 0.0033|$^*$| |$\mbox{MC}_{10}(0.21,0.03)$| 6.39e–14 4.26e–13 6.07e–04|$^*$| 0.0224 Multisymplectic 1.24e–13 7.03e–04 0.0436 0.0385 Narrow box 1.24e–13 0.0033 0.0325 0.0235 View Large Table 1 Errors in conservation laws and solution at the final time |$T=2$|⁠, when solving (3.1) with initial condition (3.27) (for |$c=d=5$|⁠) and periodic boundary conditions over |$[-20,20]$|⁠, using various schemes with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.01$|⁠. The nonzero values of the free parameter minimize the error indicated by an asterisk Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\mbox{EC}_8$| 9.24e–14 0.0019 6.71e–12 0.0964 |$\mbox{MC}_8(0)$| 1.24e–13 6.82e–13 0.0308 0.0584 |$\mbox{MC}_8(-0.069)$| 1.24e–13 1.05e–12 0.0028 0.0052|$^*$| |$\mbox{MC}_8(-0.073)$| 8.17e–14 6.25e–13 0.0014|$^*$| 0.0063 Furihata; |$\mbox{EC}_{10}(0)$| 8.17e–14 9.93e–04 2.16e–12 0.0217 |$\mbox{EC}_{10}(0.12)$| 7.46e–14 3.40e–04 2.61e–12 0.0020|$^*$| |$\mbox{EC}_{10}(0.17)$| 7.46e–14 7.29e–05|$^*$| 2.27e–12 0.0095 |$\mbox{MC}_{10}(0,0)$| 6.75e–14 3.69e–13 0.0033 0.0335 |$\mbox{MC}_{10}(0.39,0.04)$| 7.46e–14 3.41e–13 0.0127 0.0033|$^*$| |$\mbox{MC}_{10}(0.21,0.03)$| 6.39e–14 4.26e–13 6.07e–04|$^*$| 0.0224 Multisymplectic 1.24e–13 7.03e–04 0.0436 0.0385 Narrow box 1.24e–13 0.0033 0.0325 0.0235 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\mbox{EC}_8$| 9.24e–14 0.0019 6.71e–12 0.0964 |$\mbox{MC}_8(0)$| 1.24e–13 6.82e–13 0.0308 0.0584 |$\mbox{MC}_8(-0.069)$| 1.24e–13 1.05e–12 0.0028 0.0052|$^*$| |$\mbox{MC}_8(-0.073)$| 8.17e–14 6.25e–13 0.0014|$^*$| 0.0063 Furihata; |$\mbox{EC}_{10}(0)$| 8.17e–14 9.93e–04 2.16e–12 0.0217 |$\mbox{EC}_{10}(0.12)$| 7.46e–14 3.40e–04 2.61e–12 0.0020|$^*$| |$\mbox{EC}_{10}(0.17)$| 7.46e–14 7.29e–05|$^*$| 2.27e–12 0.0095 |$\mbox{MC}_{10}(0,0)$| 6.75e–14 3.69e–13 0.0033 0.0335 |$\mbox{MC}_{10}(0.39,0.04)$| 7.46e–14 3.41e–13 0.0127 0.0033|$^*$| |$\mbox{MC}_{10}(0.21,0.03)$| 6.39e–14 4.26e–13 6.07e–04|$^*$| 0.0224 Multisymplectic 1.24e–13 7.03e–04 0.0436 0.0385 Narrow box 1.24e–13 0.0033 0.0325 0.0235 View Large The upper part of Fig. 4 shows the initial condition and the numerical solution |$\mbox{EC}_{10}(0.12)$| at the final time |$T=2$|⁠. The lower plot shows only the top of the soliton, comparing the exact solution (3.28) with the numerical solutions given by |$\mbox{EC}_{10}(0.12)$|⁠, the multisymplectic, narrow box and Furihata (⁠|$\mbox{EC}_{10}(0)$|⁠) schemes. The |$\mbox{EC}_{10}(0.12)$| solution is the closest to the exact solution, reflecting the results in Table 1. Fig. 4. View largeDownload slide One-soliton solution for the KdV equation (3.1) with initial condition (3.27) and |$c=d=5$|⁠; the stepsizes are |$\varDelta x=0.1$| and |$\varDelta t=0.01$|⁠. Top: initial condition (dashed curve) and numerical solution |$\mbox{EC}_{10}(0.12)$| at |$T=2$| (solid curve). Bottom: top of the soliton at |$T=2$|⁠; exact solution and numerical solutions from |$\mbox{EC}_{10}(0.12)$|⁠, |$\mbox{EC}_{10}(0)$|⁠, multisymplectic and narrow box schemes. Fig. 4. View largeDownload slide One-soliton solution for the KdV equation (3.1) with initial condition (3.27) and |$c=d=5$|⁠; the stepsizes are |$\varDelta x=0.1$| and |$\varDelta t=0.01$|⁠. Top: initial condition (dashed curve) and numerical solution |$\mbox{EC}_{10}(0.12)$| at |$T=2$| (solid curve). Bottom: top of the soliton at |$T=2$|⁠; exact solution and numerical solutions from |$\mbox{EC}_{10}(0.12)$|⁠, |$\mbox{EC}_{10}(0)$|⁠, multisymplectic and narrow box schemes. The second benchmark test is the interaction between two solitons. The exact solution on the infinite line is \begin{equation} u_{\textrm{exact}}(x,t)=\frac{12(c_1-c_2)\left(c_1\cosh^2\xi_2+c_2\sinh^2\xi_1\right)}{\left((\sqrt{c_1}-\sqrt{c_2})\cosh(\xi_1+\xi_2)+(\sqrt{c_1}+\sqrt{c_2})\cosh(\xi_1-\xi_2)\right)^2}, \end{equation} (3.29) where \begin{equation} \xi_1=\frac{\sqrt{c_1}}2(x+d_1-c_1t),\qquad \xi_2=\frac{\sqrt{c_2}}2(x+d_2-c_2t). \end{equation} (3.30) Again, we use stepsizes |$\varDelta x=0.1$| and |$\varDelta t=0.01$| over the spatial domain |$[-20,20]$| with periodic boundary conditions, on the temporal interval |$[0,2]$|⁠. The initial condition is obtained by evaluating (3.29)–(3.30) at |$t=0$|⁠, using the parameters \begin{equation} c_1=10, \quad c_2=5, \quad d_1=12, \quad d_2=10. \end{equation} (3.31) For this problem the values |$\alpha _1=-0.099$|⁠, |$(\alpha _2,\beta _2)=(-0.011,-0.031)$| and |$\alpha _3=0.23$| minimize the solution error for |$\mbox{MC}_8(\alpha _1)$|⁠, |$\mbox{MC}_{10}(\alpha _2,\beta _2)$| and |$\mbox{EC}_{10}(\alpha _3)$|⁠, respectively. The values |$\alpha _1=-0.084$| and |$\alpha _3=0.26$|⁠, which minimize the error in the nonpreserved conservation law, both produce fairly accurate solutions. The error in the energy conservation law for |$\mbox{MC}_{10}(\alpha _2,\beta _2)$| is minimized (but remains large) for each tested |$\beta _2$| by a large negative value of |$\alpha _2$|⁠; this produces a large solution error. So minimizing the remaining conservation law is a poor criterion for selecting between the schemes |$\mbox{MC}_{10}(\alpha _2,\beta _2)$|⁠. Table 2 shows the solution error (3.24) and the error in the three conservation laws according to (3.25) or, for nonpreserved conservation laws, (3.26). The table includes the error in the phase shift for the fastest soliton at the final time, \begin{equation*}\mbox{Err}_\phi=(x_{\textrm{max}}-\tilde{x}_{\textrm{max}}){\vert_{t=2}}, \end{equation*} Table 2 Errors in conservation laws, solution and phase shift at the final time |$T=2$| for the two-soliton problem with parameters (3.31) over |$[-20,20]$| with periodic boundary conditions using different methods with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.01$|⁠. The nonzero values of the free parameter minimize the error indicated by an asterisk Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\mbox{Err}_\phi $| |$\mbox{EC}_8$| 2.27e–13 0.0201 5.64e–11 0.4561 0.36 |$\mbox{MC}_8(0)$| 3.69e–13 5.17e–12 30.3753 0.3338 0.26 |$\mbox{MC}_8(-0.099)$| 2.56e–13 3.87e–12 4.9224 0.0301|$^*$| −0.04 |$\mbox{MC}_8(-0.084)$| 2.70e–13 5.00e–12 0.3785|$^*$| 0.0625 0.06 Furihata; |$\mbox{EC}_{10}(0)$| 1.71e–13 1.3595 2.18e–11 0.1706 0.16 |$\mbox{EC}_{10}(0.23)$| 1.99e–13 0.1725 2.00e–11 0.0213|$^*$| −0.04 |$\mbox{EC}_{10}(0.26)$| 1.85e–13 0.0187|$^*$| 2.55e–11 0.0301 −0.04 |$\mbox{MC}_{10}(0,0)$| 1.42e–13 1.36e–12 27.7429 0.2391 0.16 |$\mbox{MC}_{10}(-0.011,-0.031)$| 1.85e–13 1.71e–12 28.6356 0.0253|$^*$| −0.04 Multisymplectic 1.85e–13 0.4373 28.1328 0.2557 0.26 Narrow box 1.71e–13 0.8633 17.9549 0.0255 0.06 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\mbox{Err}_\phi $| |$\mbox{EC}_8$| 2.27e–13 0.0201 5.64e–11 0.4561 0.36 |$\mbox{MC}_8(0)$| 3.69e–13 5.17e–12 30.3753 0.3338 0.26 |$\mbox{MC}_8(-0.099)$| 2.56e–13 3.87e–12 4.9224 0.0301|$^*$| −0.04 |$\mbox{MC}_8(-0.084)$| 2.70e–13 5.00e–12 0.3785|$^*$| 0.0625 0.06 Furihata; |$\mbox{EC}_{10}(0)$| 1.71e–13 1.3595 2.18e–11 0.1706 0.16 |$\mbox{EC}_{10}(0.23)$| 1.99e–13 0.1725 2.00e–11 0.0213|$^*$| −0.04 |$\mbox{EC}_{10}(0.26)$| 1.85e–13 0.0187|$^*$| 2.55e–11 0.0301 −0.04 |$\mbox{MC}_{10}(0,0)$| 1.42e–13 1.36e–12 27.7429 0.2391 0.16 |$\mbox{MC}_{10}(-0.011,-0.031)$| 1.85e–13 1.71e–12 28.6356 0.0253|$^*$| −0.04 Multisymplectic 1.85e–13 0.4373 28.1328 0.2557 0.26 Narrow box 1.71e–13 0.8633 17.9549 0.0255 0.06 View Large Table 2 Errors in conservation laws, solution and phase shift at the final time |$T=2$| for the two-soliton problem with parameters (3.31) over |$[-20,20]$| with periodic boundary conditions using different methods with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.01$|⁠. The nonzero values of the free parameter minimize the error indicated by an asterisk Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\mbox{Err}_\phi $| |$\mbox{EC}_8$| 2.27e–13 0.0201 5.64e–11 0.4561 0.36 |$\mbox{MC}_8(0)$| 3.69e–13 5.17e–12 30.3753 0.3338 0.26 |$\mbox{MC}_8(-0.099)$| 2.56e–13 3.87e–12 4.9224 0.0301|$^*$| −0.04 |$\mbox{MC}_8(-0.084)$| 2.70e–13 5.00e–12 0.3785|$^*$| 0.0625 0.06 Furihata; |$\mbox{EC}_{10}(0)$| 1.71e–13 1.3595 2.18e–11 0.1706 0.16 |$\mbox{EC}_{10}(0.23)$| 1.99e–13 0.1725 2.00e–11 0.0213|$^*$| −0.04 |$\mbox{EC}_{10}(0.26)$| 1.85e–13 0.0187|$^*$| 2.55e–11 0.0301 −0.04 |$\mbox{MC}_{10}(0,0)$| 1.42e–13 1.36e–12 27.7429 0.2391 0.16 |$\mbox{MC}_{10}(-0.011,-0.031)$| 1.85e–13 1.71e–12 28.6356 0.0253|$^*$| −0.04 Multisymplectic 1.85e–13 0.4373 28.1328 0.2557 0.26 Narrow box 1.71e–13 0.8633 17.9549 0.0255 0.06 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\mbox{Err}_\phi $| |$\mbox{EC}_8$| 2.27e–13 0.0201 5.64e–11 0.4561 0.36 |$\mbox{MC}_8(0)$| 3.69e–13 5.17e–12 30.3753 0.3338 0.26 |$\mbox{MC}_8(-0.099)$| 2.56e–13 3.87e–12 4.9224 0.0301|$^*$| −0.04 |$\mbox{MC}_8(-0.084)$| 2.70e–13 5.00e–12 0.3785|$^*$| 0.0625 0.06 Furihata; |$\mbox{EC}_{10}(0)$| 1.71e–13 1.3595 2.18e–11 0.1706 0.16 |$\mbox{EC}_{10}(0.23)$| 1.99e–13 0.1725 2.00e–11 0.0213|$^*$| −0.04 |$\mbox{EC}_{10}(0.26)$| 1.85e–13 0.0187|$^*$| 2.55e–11 0.0301 −0.04 |$\mbox{MC}_{10}(0,0)$| 1.42e–13 1.36e–12 27.7429 0.2391 0.16 |$\mbox{MC}_{10}(-0.011,-0.031)$| 1.85e–13 1.71e–12 28.6356 0.0253|$^*$| −0.04 Multisymplectic 1.85e–13 0.4373 28.1328 0.2557 0.26 Narrow box 1.71e–13 0.8633 17.9549 0.0255 0.06 View Large where |$x_{\textrm{max}}$| and |$\tilde{x}_{\textrm{max}}$| denote the location of the peak of the fastest soliton in the exact and numerical solution, respectively. Table 2 shows that |$\mbox{MC}_{8}(-0.099)$|⁠, |$\mbox{EC}_{10}(0.23)$|⁠, |$\mbox{MC}_{10}(-0.011,-0.031)$| and the narrow box scheme give the most accurate solutions. The schemes obtained by choosing the free parameter in |$\mbox{MC}_8$| and |$\mbox{EC}_{10}$| to minimize the error in the nonpreserved conservation law are more accurate than the Furihata and multisymplectic schemes. The upper part of Fig. 5 shows the initial condition (dashed line) and the numerical solution |$\mbox{EC}_{10}(0.23)$| at time |$T=2$| (solid line). The lower plot shows the exact solution (3.29) and the numerical solutions from |$\mbox{EC}_{10}(0.23)$|⁠, the Furihata (⁠|$\mbox{EC}_{10}(0)$|⁠), multisymplectic and narrow box schemes. The narrow box and |$\mbox{EC}_{10}(0.23)$| schemes are the most accurate and give similar results. Fig. 5. View largeDownload slide Two-soliton solutions of KdV with the parameters (3.31) and stepsizes |$\varDelta x=0.1$| and |$\varDelta t=0.01$|⁠. Top: initial condition (dashed curve) and numerical solution given by |$\mbox{EC}_{10}(0.23)$| at |$T=2$| (solid curve). Bottom: the top of the fastest soliton at |$T=2$|⁠; exact solution and numerical solutions from the |$\mbox{EC}_{10}(0.23)$|⁠, |$\mbox{EC}_{10}(0)$|⁠, multisymplectic and narrow box schemes. Fig. 5. View largeDownload slide Two-soliton solutions of KdV with the parameters (3.31) and stepsizes |$\varDelta x=0.1$| and |$\varDelta t=0.01$|⁠. Top: initial condition (dashed curve) and numerical solution given by |$\mbox{EC}_{10}(0.23)$| at |$T=2$| (solid curve). Bottom: the top of the fastest soliton at |$T=2$|⁠; exact solution and numerical solutions from the |$\mbox{EC}_{10}(0.23)$|⁠, |$\mbox{EC}_{10}(0)$|⁠, multisymplectic and narrow box schemes. As a last numerical test we solve the two-soliton problem on a coarser time grid, setting |$\varDelta x=0.1$| and |$\varDelta t=0.02$|⁠. For these stepsizes, the values |$\alpha _1=-0.22$|⁠, |$(\alpha _2,\beta _2)=(0.815,-0.001)$| and |$\alpha _3=0.66$| minimize the solution error for |$\mbox{MC}_8(\alpha _1)$|⁠, |$\mbox{MC}_{10}(\alpha _2,\beta _2)$| and |$\mbox{EC}_{10}(\alpha _3)$|⁠. The values |$\alpha _1=-0.16$| and |$\alpha _3=0.53$| yield the minimal error in the nonpreserved conservation law. Most results in Table 3 are qualitatively similar to their counterparts in Table 2, though with larger solution and phase errors. However, the narrow box scheme is far less accurate for the larger time step. The solution error in the most accurate scheme, |$\textrm{EC}_{10}(0.66)$|⁠, is around |$3.5$| times that of the most accurate |$\textrm{EC}_{10}$| scheme for the smaller time step |$\varDelta t=0.01$|⁠. This growth in solution error is slightly greater than those of the Furihata (⁠|$\textrm{EC}_{10}(0)$|⁠) and multisymplectic schemes (whose phase errors also grow more slowly). Even so, |$\textrm{EC}_{10}(0.66)$| is by far the most accurate of the schemes (see Fig. 6). Table 3 Errors in conservation laws, solution and phase shift at the final time |$T=2$|⁠, for the two-soliton KdV problem with parameters (3.31) over |$[-20,20]$| with periodic boundary conditions, using different methods with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.02$|⁠. The nonzero values of the free parameter minimize the error indicated by an asterisk Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\textrm{Err}_\phi $| |$\mbox{EC}_8$| 2.13e–13 1.2732 6.37e–11 0.6704 0.66 |$\mbox{MC}_8(0)$| 1.99e–13 1.99e–12 55.4259 0.6281 0.56 |$\mbox{MC}_8(-0.22)$| 3.98e–13 5.97e–12 15.8171 0.1065|$^*$| 0.06 |$\mbox{MC}_8(-0.16)$| 2.27e–13 4.21e–12 5.997|$^*$| 0.2176 0.16 Furihata; |$\mbox{EC}_{10}(0)$| 1.99e–13 2.7234 2.82e–11 0.4501 0.36 |$\mbox{EC}_{10}(0.66)$| 1.42e–13 0.5161 2.00e–11 0.0749|$^*$| 0.06 |$\mbox{EC}_{10}(0.53)$| 1.14e–13 0.1326|$^*$| 2.73e–11 0.1245 0.06 |$\mbox{MC}_{10}(0,0)$| 1.56e–13 1.82e–12 54.2290 0.5621 0.46 |$\mbox{MC}_{10}(0.815,-0.001)$| 1.56e–13 1.25e–12 59.5767 0.0947|$^*$| 0.06 Multisymplectic 2.27e–13 0.4306 53.2081 0.5678 0.46 Narrow box 2.13e–13 0.8481 10.3316 0.3860 0.36 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\textrm{Err}_\phi $| |$\mbox{EC}_8$| 2.13e–13 1.2732 6.37e–11 0.6704 0.66 |$\mbox{MC}_8(0)$| 1.99e–13 1.99e–12 55.4259 0.6281 0.56 |$\mbox{MC}_8(-0.22)$| 3.98e–13 5.97e–12 15.8171 0.1065|$^*$| 0.06 |$\mbox{MC}_8(-0.16)$| 2.27e–13 4.21e–12 5.997|$^*$| 0.2176 0.16 Furihata; |$\mbox{EC}_{10}(0)$| 1.99e–13 2.7234 2.82e–11 0.4501 0.36 |$\mbox{EC}_{10}(0.66)$| 1.42e–13 0.5161 2.00e–11 0.0749|$^*$| 0.06 |$\mbox{EC}_{10}(0.53)$| 1.14e–13 0.1326|$^*$| 2.73e–11 0.1245 0.06 |$\mbox{MC}_{10}(0,0)$| 1.56e–13 1.82e–12 54.2290 0.5621 0.46 |$\mbox{MC}_{10}(0.815,-0.001)$| 1.56e–13 1.25e–12 59.5767 0.0947|$^*$| 0.06 Multisymplectic 2.27e–13 0.4306 53.2081 0.5678 0.46 Narrow box 2.13e–13 0.8481 10.3316 0.3860 0.36 View Large Table 3 Errors in conservation laws, solution and phase shift at the final time |$T=2$|⁠, for the two-soliton KdV problem with parameters (3.31) over |$[-20,20]$| with periodic boundary conditions, using different methods with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.02$|⁠. The nonzero values of the free parameter minimize the error indicated by an asterisk Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\textrm{Err}_\phi $| |$\mbox{EC}_8$| 2.13e–13 1.2732 6.37e–11 0.6704 0.66 |$\mbox{MC}_8(0)$| 1.99e–13 1.99e–12 55.4259 0.6281 0.56 |$\mbox{MC}_8(-0.22)$| 3.98e–13 5.97e–12 15.8171 0.1065|$^*$| 0.06 |$\mbox{MC}_8(-0.16)$| 2.27e–13 4.21e–12 5.997|$^*$| 0.2176 0.16 Furihata; |$\mbox{EC}_{10}(0)$| 1.99e–13 2.7234 2.82e–11 0.4501 0.36 |$\mbox{EC}_{10}(0.66)$| 1.42e–13 0.5161 2.00e–11 0.0749|$^*$| 0.06 |$\mbox{EC}_{10}(0.53)$| 1.14e–13 0.1326|$^*$| 2.73e–11 0.1245 0.06 |$\mbox{MC}_{10}(0,0)$| 1.56e–13 1.82e–12 54.2290 0.5621 0.46 |$\mbox{MC}_{10}(0.815,-0.001)$| 1.56e–13 1.25e–12 59.5767 0.0947|$^*$| 0.06 Multisymplectic 2.27e–13 0.4306 53.2081 0.5678 0.46 Narrow box 2.13e–13 0.8481 10.3316 0.3860 0.36 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| |$\textrm{Err}_3$| Solution error |$\textrm{Err}_\phi $| |$\mbox{EC}_8$| 2.13e–13 1.2732 6.37e–11 0.6704 0.66 |$\mbox{MC}_8(0)$| 1.99e–13 1.99e–12 55.4259 0.6281 0.56 |$\mbox{MC}_8(-0.22)$| 3.98e–13 5.97e–12 15.8171 0.1065|$^*$| 0.06 |$\mbox{MC}_8(-0.16)$| 2.27e–13 4.21e–12 5.997|$^*$| 0.2176 0.16 Furihata; |$\mbox{EC}_{10}(0)$| 1.99e–13 2.7234 2.82e–11 0.4501 0.36 |$\mbox{EC}_{10}(0.66)$| 1.42e–13 0.5161 2.00e–11 0.0749|$^*$| 0.06 |$\mbox{EC}_{10}(0.53)$| 1.14e–13 0.1326|$^*$| 2.73e–11 0.1245 0.06 |$\mbox{MC}_{10}(0,0)$| 1.56e–13 1.82e–12 54.2290 0.5621 0.46 |$\mbox{MC}_{10}(0.815,-0.001)$| 1.56e–13 1.25e–12 59.5767 0.0947|$^*$| 0.06 Multisymplectic 2.27e–13 0.4306 53.2081 0.5678 0.46 Narrow box 2.13e–13 0.8481 10.3316 0.3860 0.36 View Large Fig. 6. View largeDownload slide Two-soliton solutions for (3.1) with the parameters (3.31) and stepsizes |$\varDelta x=0.1$| and |$\varDelta t=0.02$|⁠. Top of the fastest soliton at time |$T=2$|⁠: exact solution and numerical solutions from the |$\mbox{EC}_{10}(0.66)$|⁠, |$\mbox{EC}_{10}(0)$|⁠, multisymplectic and narrow box schemes. Fig. 6. View largeDownload slide Two-soliton solutions for (3.1) with the parameters (3.31) and stepsizes |$\varDelta x=0.1$| and |$\varDelta t=0.02$|⁠. Top of the fastest soliton at time |$T=2$|⁠: exact solution and numerical solutions from the |$\mbox{EC}_{10}(0.66)$|⁠, |$\mbox{EC}_{10}(0)$|⁠, multisymplectic and narrow box schemes. 4. A nonlinear heat equation In this section we consider the nonlinear heat equation \begin{equation} {\mathcal{A}}\equiv u_t-u_x^2-uu_{xx}=0,\qquad (x,t)\in\varOmega~\equiv~[a,b]\times[0,\infty), \end{equation} (4.1) coupled with suitable initial and boundary conditions: \begin{equation} u(x,0)=\psi(x),\qquad u(a,t)=\varphi_1(t),\qquad u(b,t)=\varphi_2(t). \end{equation} (4.2) Equation (4.1) has only two independent (equivalence classes of) conservation laws: \begin{eqnarray} {\mathcal{A}} &=& D_t(G_1)+D_x(F_1)\equiv D_t(u) + D_x(-uu_x), \end{eqnarray} (4.3) \begin{eqnarray} x{\mathcal{A}} &=& D_t(G_2)+D_x(F_2)\equiv D_t(xu) + D_x\left({u^2}/{2}-xuu_x\right), \end{eqnarray} (4.4) with characteristics \begin{equation} {\mathcal{Q}}_1=1,\quad{\mathcal{Q}}_2=x, \end{equation} (4.5) respectively (see Ibragimov, 1994). To construct finite difference schemes that preserve a discrete version of (4.3) and (4.4) we use the following general results. Theorem 4.1 Any PDE of the form \begin{equation} D_t\left(g[u]\right)+D_x^{\,k}\left(\,f[u]\right)=0,\qquad k\in\mathbb{N}, \end{equation} (4.6) where |$g[u]$| and |$f[u]$| are smooth functions of |$u$| and its derivatives, has |$k$| conservation laws whose characteristics are |$Q_i=x^{i-1},\, i=1,\ldots ,k.$| Any scheme of the form \begin{equation} D_n\left(\widetilde{g[u]}\right)+D_m^{\,k}\left(\,\widetilde{f[u]}\right)=0, \end{equation} (4.7) where |$\widetilde{g[u]}$| and |$\widetilde{f[u]}$| are finite difference approximations to |$g[u]$| and |$\,f[u],$| respectively, has |$k$| conservation laws whose characteristics are |$\widetilde{Q_i}=x_0^{i-1},\, i=1,\ldots ,k$|⁠. Proof. On solutions of (4.6), using integration by parts, \begin{equation*}x^{i-1}\left( D_t(g[u])+D_x^{\,k}(\,f[u])\right)=D_t\left(x^{i-1}g[u]\right)+(-1)^kD_x^{\,k}(x^{i-1})\,f[u]+D_x(h[u]),\end{equation*} for some function |$h$| of |$u$| and its derivatives. As |$i-1<k$| this simplifies to \begin{equation*}x^{i-1}\left( D_t(g[u])+D_x^{\,k}(\,f[u])\right)=D_t\left(x^{i-1}g[u]\right)+D_x(h[u]),\end{equation*} which is a divergence and, therefore, a conservation law with characteristic |$Q_i=x^{i-1}.$| The proof for the discrete case is similar. Summation by parts gives \begin{equation*}x_0^{i-1}\left( D_n\left(\widetilde{g[u]}\right)+D_m^{\,k}\left(\,\widetilde{f[u]}\right)\right)=D_n\left(x_0^{i-1}\widetilde{g[u]}\right)+\left(-S_m^{-1}D_m\right)^k(x_0^{i-1})\,\widetilde{f[u]}+D_m\left(\widetilde{h[u]}\right)\end{equation*} on solutions of (4.7), for some function |$\widetilde{h[u]}$| of |$u$| and its shifts. Again, |$i-1<k$| so \begin{equation*}x_0^{i-1}\left( D_n\left(\widetilde{g[u]}\right)+D_m^{\,k}\left(\,\widetilde{f[u]}\right)\right)=D_n\left(x_0^{i-1}\widetilde{g[u]}\right)+D_m\left(\widetilde{h[u]}\right)\!,\end{equation*} which is a discrete conservation law with characteristic |$\widetilde{Q_i}=x_0^{i-1}.$| 4.1 Conservative methods for the nonlinear heat equation The nonlinear heat equation (4.1) is of the form (4.6), with |$k=2$|⁠. Therefore, according to Theorem 4.1, both conservation laws can be preserved by finding suitable finite difference approximations of |$g[u]=u$| and |$f[u]=-{u^2}/{2}$|⁠. This can be achieved to second order on the most compact rectangular stencil for (4.1), which consists of six points (see Fig. 7). Fig. 7. View largeDownload slide The six-point rectangular stencil for (4.1). Conservation laws are preserved to second order at the central point |$(0, 1/2)$| (circle), densities and fluxes, respectively, at |$(0,0) $| and |$(-1/2, 1/2)$| (crosses). Fig. 7. View largeDownload slide The six-point rectangular stencil for (4.1). Conservation laws are preserved to second order at the central point |$(0, 1/2)$| (circle), densities and fluxes, respectively, at |$(0,0) $| and |$(-1/2, 1/2)$| (crosses). This results in a two-parameter family of schemes, \begin{equation} \mbox{CS}(\alpha,\beta)\equiv\widetilde{{\mathcal{A}}}\equiv D_m\widetilde{F_1}+D_n\widetilde{G_1}=0, \end{equation} (4.8) with \begin{align} \widetilde{F_1}&=D_m\left(\,\widetilde{f[u]}\right),\quad\textrm{where}\quad\widetilde{f[u]}=-\frac{1}2u_{-1,0}u_{-1,1}+\varDelta t^2\beta(D_nu_{-1,0})^2, \end{align} (4.9) \begin{align}\widetilde{G_1}=\widetilde{g[u]}=u_{0, 0}+\alpha \varDelta x^2D_m^2u_{-1, 0}\,. \end{align} (4.10) These schemes preserve a discrete version of (4.3) with |$\widetilde{Q_1}=1$| and of (4.4), namely, \begin{equation} \widetilde{Q_2}\widetilde{\mathcal{A}}=D_m\widetilde{F_2}+D_n\widetilde{G_2}, \end{equation} (4.11) where \begin{equation*} \widetilde{Q_2}=x_0\,,\qquad \widetilde{F_2}=\mu_m(x_{-1})\widetilde{F_1}-\mu_m\left(\,\widetilde{f[u]}\right)\!,\qquad \widetilde{G_2}=x_0\widetilde{G_1}. \end{equation*} Except for |$\mbox{CS}(0,0)$|⁠, the approximated densities and fluxes include derivative terms that do not appear in the corresponding continuous quantities. However, these vanish as the stepsizes approach zero. The following schemes are particularly straightforward. The scheme |$\mbox{CS}(0,-1/4)$| has perhaps the most obvious discretization of |$\widetilde{F_1}$| and |$\widetilde{G_1}$|⁠, namely, \begin{equation*} \widetilde{F_1}=D_m\left(-\tfrac{1}{2}\mu_nu_{-1,0}^2\right)\!,\qquad \widetilde{G_1}=u_{0, 0}. \end{equation*} The components of the second conservation law are \begin{equation*} \widetilde{F_2}=\mu_m(x_{-1})\widetilde{F_1}+\tfrac{1}2\mu_m\left(\mu_nu_{-1,0}^2\right)\!,\qquad \widetilde{G_2}=x_0u_{0,0}. \end{equation*} The scheme |$\mbox{CS}(0,-1/8)$| amounts to \begin{equation*} D_n(u_{0,0})+D_m^2\left(-\tfrac{1}{2}\left(\mu_nu_{-1,0}\right)^2\right)=0, \end{equation*} which is obtained by applying the implicit midpoint rule to the following simple semidiscretization of (4.1), with |$U_i(t)$| approximating |$u(x_i,t)$|⁠: \begin{equation*}U_0^{\prime}+D_m^2\left(-\tfrac{1}2(U_{-1})^2\right)=0.\end{equation*} The densities and fluxes of the discrete conservation laws preserved by |$\mbox{CS}(0,-1/8)$| are \begin{equation*} \widetilde{F_1}=D_m\left(-\tfrac{1}{2} \left(\mu_nu_{-1,0}\right)^2\right)\!,\qquad \widetilde{G_1}=u_{0, 0}; \end{equation*} \begin{equation*} \widetilde{F_2}=\mu_m(x_{-1})\widetilde{F_1}+\tfrac{1}2\mu_m\left(\left(\mu_nu_{-1,0}\right)^2\right)\!,\qquad \widetilde{G_2}=x_0u_{0,0}. \end{equation*} The scheme with the fewest terms is |$\mbox{CS}(0,0)$|⁠, whose densities and fluxes are \begin{equation*} \widetilde{F_1}=D_m\left(-\tfrac{1}{2}{u_{-1,0}u_{-1,1}}\right)\!,\qquad \widetilde{G_1}=u_{0, 0}; \end{equation*} \begin{equation*} \widetilde{F_2}=\mu_m(x_{-1})\widetilde{F_1}+\tfrac{1}2\mu_m\left({u_{-1,0}u_{-1,1}}\right)\!,\qquad \widetilde{G_2}=x_0u_{0,0}. \end{equation*} This scheme can be solved explicitly at the |$k$|th integration step, provided that the matrix \begin{equation*} \frac{1}{\varDelta t}\,I-\frac{1}{2\varDelta x^2}\,M\circ V_{k-1} \end{equation*} is invertible; here |$I$| is the identity matrix, |$\circ $| is the Hadamard product,2 \begin{equation*}M=\left(\!\!\begin{array}{ccccc} -2& 1 & & &\\ 1 & -2 & 1 & &\\ & 1 & -2 & 1 & \\ & & \ddots & \ddots & \ddots \end{array}\!\!\right)\!,\qquad V_{k-1}=\textbf{e}\,u_{k-1},\end{equation*} where |$\textbf{e}$| is the column vector whose entries are all |$1$| and |$u_{k-1}$| is the row vector whose entries are the approximation (from the previous step) of |$u$| at the spatial grid points. Remark 4.2 The nonlinear heat equation (4.1) is a special case of the porous medium equation, \begin{equation} D_t(u)+D_x^2\left(-\,\frac{u^s}{s}\right)=0,\qquad (x,t)\in\varOmega~\equiv~[a,b]\times[0,\infty), \end{equation} (4.12) where |$s\in \mathbb{N}\backslash \{1\}$|⁠. For each |$s$| this equation has only two conservation laws, with characteristics |$1$| and |$x$|⁠. The approach that we have used for the nonlinear heat equation can be used to obtain conservative schemes for (4.12) with |$s>2$|⁠. On the six-point stencil in Fig. 7 this yields an |$s$|-parameter family of second-order methods. 4.2 Numerical tests In this section three benchmark numerical tests for the problems (4.1)–(4.2) are used to show the effectiveness of the methods developed in Section 4.1. We compare the results from several |$\mbox{CS}(\alpha ,\beta )$| schemes, which preserve both conservation laws, with those from the following second-order scheme that, in general, does not preserve either conservation law: \begin{align} D_nu_{0,0}-(D_m\mu_m\mu_nu_{-1,0})^2-\mu_n(u_{0,0})D_m^2(\mu_nu_{-1,0})=0. \end{align} (4.13) We call the scheme ML/IM, as it is obtained by applying the implicit midpoint method to the following standard second-order method-of-lines semidiscretization of (4.1): \begin{align*} U_0^{\prime}-(D_m\mu_mU_{-1})^2-U_0D_m^2U_{-1}=0. \end{align*} Again, implicit methods are solved by a simplified Newton method with frozen Jacobian, run until the error reaches machine accuracy. The (relative) solution error at the final time |$t=T$| is evaluated as \begin{equation} \frac{\|u-u_{\textrm{exact}}\|}{\|u_{\textrm{exact}}\|}\bigg\vert_{t=T}. \end{equation} (4.14) The errors in the discrete conservation laws (4.8) and (4.11) are evaluated respectively as3 \begin{align} \mbox{Err}_1=\,\varDelta x\max_{j=1,\ldots, N-1} &\left|\sum_{i=2}^{M-1}D_n\!\left(u_{i,\,j}+\alpha\varDelta x^2D_m^2u_{i-1,\,j}\right)-D_m\left\{\left[\frac{u_{r,\,j}u_{r,\,j+1}+2\beta\varDelta t^2(D_nu_{r,\,j})^2}{2\varDelta x}\right]_{r=1}^{M-1}\right\}\right|\,, \end{align} (4.15) \begin{align}\nonumber \mbox{Err}_2=\,\varDelta x\max_{j=1,\ldots, N-1} &\left|\sum_{i=2}^{M-1}x_i\,D_n\!\left(\!u_{i,\,j}+\alpha\varDelta x^2D_m^2u_{i-1,\,j}\right)\!-\!\left[\!\frac{\mu_m(x_r)D_m(u_{r,\,j}u_{r,\,j+1})-\mu_m(u_{r,\,j}u_{r,\,j+1})}{2\varDelta x}\!\right]_{r=1}^{M-1}\right.\\[5pt] &\left.+\,\beta\varDelta t^2\left[\frac{\mu_m(x_r)D_m\!\big((D_nu_{r,\,j})^2\big)-\mu_m\big((D_nu_{r,\,j})^2\big)}{\varDelta x}\right]_{r=1}^{M-1}\,\, \right|\,, \end{align} (4.16) where |$u_{i,\,j}\simeq u(a+i\varDelta x, j\varDelta t)$|⁠, so that subscripts denote shifts with respect to |$(x,t)=(a,0)$|⁠. To evaluate the error in the conservation laws resulting from ML/IM we use (4.15) and (4.16), setting |$\alpha =\beta =0$|⁠. The first benchmark problem is (4.1) with the initial and boundary conditions \begin{align} u(x,0)=\left(1-\frac{x^2}{6}\right)_{\!\!+}, \quad -6\leqslant x\leqslant 6,\qquad u(-6,t)=u(6,t)=0,\quad t\in[0,4], \end{align} (4.17) where |$f_+=\max (\,f,0)$|⁠. These conditions yield the Barenblatt solution of the porous medium equation (4.12) with |$s=2$|⁠, which is \begin{equation*} u_{\textrm{exact}}(x,t)=(t+1)^{-1/3}\left(1-\frac{x^2}{6(t+1)^{2/3}}\right)_{\!\!+}. \end{equation*} For all |$t>0$| this solution has compact support with the interface moving outward at a finite speed. The Barenblatt solution is a (weak) energy solution, but is not a classical solution as it is not differentiable at the interface points. Such solutions cause difficulties in numerical simulation. Standard finite element methods can create oscillations close to the interface, but negative values have no meaning physically (see Zhang & Wu, 2009). Here we show that, by contrast, various conservative finite difference schemes |$\textrm{CS}(\alpha ,\beta )$| are effective for nonsmooth solutions. For simplicity, we will consider only the one-parameter family obtained by setting |$\alpha =0$|⁠. Table 4 shows the errors in the conservation laws for various |$\textrm{CS}(0,\beta )$| given the stepsizes |$\varDelta x=0.25$| and |$\varDelta t=0.333$|⁠. These schemes locally preserve both conservation laws to machine accuracy. The solution error at the final time |$T=4$|⁠, evaluated according to (4.14), is minimized by setting |$\beta =0.21$|⁠. Nevertheless, the explicitly solved scheme |$\mbox{CS}(0,0)$|⁠, though slightly less accurate, is a better option because of its low computation time. Table 4 Errors in solution and conservation laws at the final time |$T=4$|⁠, when solving (4.1) with the conditions (4.17), using CS|$(0,\beta )$| schemes with |$\varDelta x=0.25$|⁠, |$\varDelta t=0.333$| and the scheme ML/IM with |$\varDelta x=0.25$|⁠, |$\varDelta t=0.0267$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 8.11e–16 5.66e–16 0.0038 0.032 |$\mbox{CS}(0,-1/8)$| 4.55e–16 3.92e–16 0.0035 0.026 |$\mbox{CS}(0,0)$| 5.95e–15 5.41e–15 0.0032 0.002 |$\mbox{CS}(0,0.21)$| 8.61e–16 1.03e–15 0.0028 0.030 ML/IM 0.0671 5.30e–15 0.0307 0.221 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 8.11e–16 5.66e–16 0.0038 0.032 |$\mbox{CS}(0,-1/8)$| 4.55e–16 3.92e–16 0.0035 0.026 |$\mbox{CS}(0,0)$| 5.95e–15 5.41e–15 0.0032 0.002 |$\mbox{CS}(0,0.21)$| 8.61e–16 1.03e–15 0.0028 0.030 ML/IM 0.0671 5.30e–15 0.0307 0.221 View Large Table 4 Errors in solution and conservation laws at the final time |$T=4$|⁠, when solving (4.1) with the conditions (4.17), using CS|$(0,\beta )$| schemes with |$\varDelta x=0.25$|⁠, |$\varDelta t=0.333$| and the scheme ML/IM with |$\varDelta x=0.25$|⁠, |$\varDelta t=0.0267$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 8.11e–16 5.66e–16 0.0038 0.032 |$\mbox{CS}(0,-1/8)$| 4.55e–16 3.92e–16 0.0035 0.026 |$\mbox{CS}(0,0)$| 5.95e–15 5.41e–15 0.0032 0.002 |$\mbox{CS}(0,0.21)$| 8.61e–16 1.03e–15 0.0028 0.030 ML/IM 0.0671 5.30e–15 0.0307 0.221 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 8.11e–16 5.66e–16 0.0038 0.032 |$\mbox{CS}(0,-1/8)$| 4.55e–16 3.92e–16 0.0035 0.026 |$\mbox{CS}(0,0)$| 5.95e–15 5.41e–15 0.0032 0.002 |$\mbox{CS}(0,0.21)$| 8.61e–16 1.03e–15 0.0028 0.030 ML/IM 0.0671 5.30e–15 0.0307 0.221 View Large ML/IM does not converge on such a coarse grid. Only by reducing the time step so that |$\varDelta t<\varDelta x^2$| can this scheme be made to converge. Reducing the time step to |$\varDelta t=0.0267$| the solution error is still larger than those of the |$\mbox{CS}(0,\beta )$| methods, which converge even when |$\varDelta t>\varDelta x$|⁠. Note that ML/IM preserves the conservation law (4.4); this is a consequence of the reflectional symmetry of the scheme and the boundary conditions. Tables 5 and 6 show the outcomes of solving the same problem with various |$\mbox{CS}(0,\beta )$| on the finer grids |$\varDelta x=0.1$|⁠, |$\varDelta t=0.133$| and |$\varDelta x=0.025$|⁠, |$\varDelta t=0.03$|⁠. On these grids the values |$\beta =0.07$| and |$\beta =0.05,$| respectively, minimize the solution error. The explicit scheme |$\mbox{CS}(0,0)$| is by far the most efficient and has a low solution error. Of the implicit schemes the optimized scheme is the fastest in each case. The errors in the conservation laws are tiny, but grow as the grid is refined due to the accumulation of rounding errors. Again, ML/IM requires smaller time steps for convergence; even then, the solution error is still far greater than those of the conservative methods. Table 5 Errors in solution and conservation laws at the final time |$T=4$|⁠, when solving (4.1) with the conditions in (4.17), using CS|$(0,\beta )$| methods with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.133$| and the scheme ML/IM with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.005$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 1.66e–15 2.57e–15 0.0013 0.15 |$\mbox{CS}(0,-1/8)$| 1.58e–15 1.83e–15 0.0012 0.10 |$\mbox{CS}(0,0)$| 4.88e–14 4.67e–14 0.0011 0.014 |$\mbox{CS}(0,0.07)$| 3.70e–15 3.36e–15 9.77e–04 0.09 ML/IM 0.0232 2.37e–14 0.0126 6.48 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 1.66e–15 2.57e–15 0.0013 0.15 |$\mbox{CS}(0,-1/8)$| 1.58e–15 1.83e–15 0.0012 0.10 |$\mbox{CS}(0,0)$| 4.88e–14 4.67e–14 0.0011 0.014 |$\mbox{CS}(0,0.07)$| 3.70e–15 3.36e–15 9.77e–04 0.09 ML/IM 0.0232 2.37e–14 0.0126 6.48 View Large Table 5 Errors in solution and conservation laws at the final time |$T=4$|⁠, when solving (4.1) with the conditions in (4.17), using CS|$(0,\beta )$| methods with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.133$| and the scheme ML/IM with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.005$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 1.66e–15 2.57e–15 0.0013 0.15 |$\mbox{CS}(0,-1/8)$| 1.58e–15 1.83e–15 0.0012 0.10 |$\mbox{CS}(0,0)$| 4.88e–14 4.67e–14 0.0011 0.014 |$\mbox{CS}(0,0.07)$| 3.70e–15 3.36e–15 9.77e–04 0.09 ML/IM 0.0232 2.37e–14 0.0126 6.48 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 1.66e–15 2.57e–15 0.0013 0.15 |$\mbox{CS}(0,-1/8)$| 1.58e–15 1.83e–15 0.0012 0.10 |$\mbox{CS}(0,0)$| 4.88e–14 4.67e–14 0.0011 0.014 |$\mbox{CS}(0,0.07)$| 3.70e–15 3.36e–15 9.77e–04 0.09 ML/IM 0.0232 2.37e–14 0.0126 6.48 View Large Table 6 Errors in solution and conservation laws at the final time |$T=4$|⁠, when solving (4.1) with the conditions in (4.17), using CS|$(0,\beta )$| methods with |$\varDelta x=0.025$|⁠, |$\varDelta t=0.03$| and the scheme ML/IM with |$\varDelta x=0.025$|⁠, |$\varDelta t=0.000267$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 3.16e–14 1.78e–14 6.51e–05 14.35 |$\mbox{CS}(0,-1/8)$| 2.18e–14 1.50e–14 5.87e–05 8.99 |$\mbox{CS}(0,0)$| 3.22e–14 2.69e–14 5.48e–05 0.38 |$\mbox{CS}(0,0.05)$| 3.20e–14 4.49e–14 5.42e–05 6.71 ML/IM 0.0095 2.50e–13 0.0035 1612.61 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 3.16e–14 1.78e–14 6.51e–05 14.35 |$\mbox{CS}(0,-1/8)$| 2.18e–14 1.50e–14 5.87e–05 8.99 |$\mbox{CS}(0,0)$| 3.22e–14 2.69e–14 5.48e–05 0.38 |$\mbox{CS}(0,0.05)$| 3.20e–14 4.49e–14 5.42e–05 6.71 ML/IM 0.0095 2.50e–13 0.0035 1612.61 View Large Table 6 Errors in solution and conservation laws at the final time |$T=4$|⁠, when solving (4.1) with the conditions in (4.17), using CS|$(0,\beta )$| methods with |$\varDelta x=0.025$|⁠, |$\varDelta t=0.03$| and the scheme ML/IM with |$\varDelta x=0.025$|⁠, |$\varDelta t=0.000267$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 3.16e–14 1.78e–14 6.51e–05 14.35 |$\mbox{CS}(0,-1/8)$| 2.18e–14 1.50e–14 5.87e–05 8.99 |$\mbox{CS}(0,0)$| 3.22e–14 2.69e–14 5.48e–05 0.38 |$\mbox{CS}(0,0.05)$| 3.20e–14 4.49e–14 5.42e–05 6.71 ML/IM 0.0095 2.50e–13 0.0035 1612.61 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 3.16e–14 1.78e–14 6.51e–05 14.35 |$\mbox{CS}(0,-1/8)$| 2.18e–14 1.50e–14 5.87e–05 8.99 |$\mbox{CS}(0,0)$| 3.22e–14 2.69e–14 5.48e–05 0.38 |$\mbox{CS}(0,0.05)$| 3.20e–14 4.49e–14 5.42e–05 6.71 ML/IM 0.0095 2.50e–13 0.0035 1612.61 View Large The upper part of Fig. 8 shows the initial condition (dashed line) and the numerical solution given by |$\mbox{CS}(0,0)$| for equation (4.1), with conditions in (4.17), setting |$\varDelta x=0.025$| and |$\varDelta t=0.03$|⁠. The method does not produce any spurious oscillations close to the interface. Magnifying the left interface, as shown at bottom of Fig. 8, one can see that the solution of |$\mbox{CS}(0,0)$| is closer to the exact solution at time |$T=4$| than the solution given by ML/IM, even though the time step used to advance ML/IM is much smaller. Furthermore, the interface of the numerical solution has moved at the correct speed and overlaps the interface of the exact solution. The solutions given by |$\mbox{CS}(0,\beta )$| for the optimal value of |$\beta $| overlap the |$\mbox{CS}(0,0)$| solution, so we omit the corresponding figures. Fig. 8. View largeDownload slide Results for (4.1) with conditions in (4.17). Top: initial condition (dashed curve) and numerical solution from |$\mbox{CS}(0,0)$| (solid curve), setting |$\varDelta x=0.025$| and |$\varDelta t=0.03$| over |$[-6,6]$| at time |$T=4$|⁠. Bottom: left interface: exact solution (solid curve), numerical solutions |$\mbox{CS}(0,0)$| with |$\varDelta x=0.025$| and |$\varDelta t=0.03$| (crosses) and ML/IM with |$\varDelta x=0.025$| and |$\varDelta t=0.000267$| (diamonds). Fig. 8. View largeDownload slide Results for (4.1) with conditions in (4.17). Top: initial condition (dashed curve) and numerical solution from |$\mbox{CS}(0,0)$| (solid curve), setting |$\varDelta x=0.025$| and |$\varDelta t=0.03$| over |$[-6,6]$| at time |$T=4$|⁠. Bottom: left interface: exact solution (solid curve), numerical solutions |$\mbox{CS}(0,0)$| with |$\varDelta x=0.025$| and |$\varDelta t=0.03$| (crosses) and ML/IM with |$\varDelta x=0.025$| and |$\varDelta t=0.000267$| (diamonds). The second benchmark problem is (4.1) with |$(x,t)\in \varOmega =[0,15]\times [0,10]$| and the following initial and boundary conditions: \begin{align}\begin{aligned} &u(x,0)=0,\quad x\in[0,15],\\ u(0,t)=t,\quad t\in&\,[0,10],\qquad u(15,t)=0,\quad t\in[0,10]. \end{aligned}\end{align} (4.18) The exact solution of this problem (see Ibragimov, 1994) is again not smooth: \begin{equation*} u_{\textrm{exact}}(x,t)=\begin{cases} t-x,&\quad 0\leqslant x \leqslant t,\\ 0,&\quad x> t; \end{cases} \end{equation*} this is a wave travelling with unit speed into an undisturbed medium. Table 7 shows the errors for various |$\mbox{CS}(0,\beta )$| with |$\varDelta x=0.375$| and |$\varDelta t=0.333$|⁠. The value |$\beta =-0.14$| gives the minimum solution error. As ML/IM does not converge on this grid, we use the finer time step |$\varDelta t=0.00667$| for this method; the problem is not reflectionally symmetric, so neither conservation law is preserved. The results are similar to those for the first benchmark problem. Again, the sub-optimal scheme |$\mbox{CS}(0,0)$|⁠, solved explicitly, is convenient because of its low computation time. Table 7 Errors in solution and conservation laws at the final time |$T=10$|⁠, when solving (4.1) with conditions in (4.18), using CS|$(0,\beta )$| schemes with |$\varDelta x=0.375$|⁠, |$\varDelta t=0.333$| and ML/IM with |$\varDelta x=0.375$|⁠, |$\varDelta t=0.00667$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 2.13e–14 3.73e–14 0.0013 0.06 |$\mbox{CS}(0,-1/8)$| 2.40e–14 2.13e–14 9.26e–04 0.04 |$\mbox{CS}(0,0)$| 3.60e–14 5.86e–14 0.0035 0.002 |$\mbox{CS}(0,-0.14)$| 2.40e–14 2.66e–14 9.13e–04 0.04 ML/IM 0.0845 0.7938 0.0114 1.76 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 2.13e–14 3.73e–14 0.0013 0.06 |$\mbox{CS}(0,-1/8)$| 2.40e–14 2.13e–14 9.26e–04 0.04 |$\mbox{CS}(0,0)$| 3.60e–14 5.86e–14 0.0035 0.002 |$\mbox{CS}(0,-0.14)$| 2.40e–14 2.66e–14 9.13e–04 0.04 ML/IM 0.0845 0.7938 0.0114 1.76 View Large Table 7 Errors in solution and conservation laws at the final time |$T=10$|⁠, when solving (4.1) with conditions in (4.18), using CS|$(0,\beta )$| schemes with |$\varDelta x=0.375$|⁠, |$\varDelta t=0.333$| and ML/IM with |$\varDelta x=0.375$|⁠, |$\varDelta t=0.00667$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 2.13e–14 3.73e–14 0.0013 0.06 |$\mbox{CS}(0,-1/8)$| 2.40e–14 2.13e–14 9.26e–04 0.04 |$\mbox{CS}(0,0)$| 3.60e–14 5.86e–14 0.0035 0.002 |$\mbox{CS}(0,-0.14)$| 2.40e–14 2.66e–14 9.13e–04 0.04 ML/IM 0.0845 0.7938 0.0114 1.76 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 2.13e–14 3.73e–14 0.0013 0.06 |$\mbox{CS}(0,-1/8)$| 2.40e–14 2.13e–14 9.26e–04 0.04 |$\mbox{CS}(0,0)$| 3.60e–14 5.86e–14 0.0035 0.002 |$\mbox{CS}(0,-0.14)$| 2.40e–14 2.66e–14 9.13e–04 0.04 ML/IM 0.0845 0.7938 0.0114 1.76 View Large Similar results are obtained in Table 8 by solving the same problem using |$\mbox{CS}(0,\beta )$| schemes on the finer grid |$\varDelta x=0.05$|⁠, |$\varDelta t=0.025$|⁠, for which the value |$\beta =0.34$| minimizes the solution error. As ML/IM does not converge on these grids we reduce the time step to |$\varDelta t=8\text{e-}05$|⁠. Again, the conservative methods are more accurate. Table 8 Errors in solution and conservation laws at the final time |$T=10$|⁠, when solving (4.1) with conditions in (4.18), using CS|$(0,\beta )$| methods with |$\varDelta x=0.05$|⁠, |$\varDelta t=0.025$| and ML/IM with |$\varDelta x=0.05$|⁠, |$\varDelta t=8$|e-05 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 9.55e–13 1.79e–12 1.16e–04 4.92 |$\mbox{CS}(0,-1/8)$| 7.92e–13 1.34e–12 9.99e–05 3.67 |$\mbox{CS}(0,0)$| 9.96e–13 2.06e–12 8.16e–05 0.19 |$\mbox{CS}(0,0.34)$| 2.60e–12 5.37e–12 2.94e–05 6.90 ML/IM 0.0114 0.1131 0.0017 1400.70 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 9.55e–13 1.79e–12 1.16e–04 4.92 |$\mbox{CS}(0,-1/8)$| 7.92e–13 1.34e–12 9.99e–05 3.67 |$\mbox{CS}(0,0)$| 9.96e–13 2.06e–12 8.16e–05 0.19 |$\mbox{CS}(0,0.34)$| 2.60e–12 5.37e–12 2.94e–05 6.90 ML/IM 0.0114 0.1131 0.0017 1400.70 View Large Table 8 Errors in solution and conservation laws at the final time |$T=10$|⁠, when solving (4.1) with conditions in (4.18), using CS|$(0,\beta )$| methods with |$\varDelta x=0.05$|⁠, |$\varDelta t=0.025$| and ML/IM with |$\varDelta x=0.05$|⁠, |$\varDelta t=8$|e-05 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 9.55e–13 1.79e–12 1.16e–04 4.92 |$\mbox{CS}(0,-1/8)$| 7.92e–13 1.34e–12 9.99e–05 3.67 |$\mbox{CS}(0,0)$| 9.96e–13 2.06e–12 8.16e–05 0.19 |$\mbox{CS}(0,0.34)$| 2.60e–12 5.37e–12 2.94e–05 6.90 ML/IM 0.0114 0.1131 0.0017 1400.70 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 9.55e–13 1.79e–12 1.16e–04 4.92 |$\mbox{CS}(0,-1/8)$| 7.92e–13 1.34e–12 9.99e–05 3.67 |$\mbox{CS}(0,0)$| 9.96e–13 2.06e–12 8.16e–05 0.19 |$\mbox{CS}(0,0.34)$| 2.60e–12 5.37e–12 2.94e–05 6.90 ML/IM 0.0114 0.1131 0.0017 1400.70 View Large Figure 9 compares numerical solutions of problem (4.1) with (4.18) on the finer grid, which has a mesh point at the interface |$x=10$|⁠. Although ML/IM is qualitatively correct it produces a slight lag near to the interface. The scheme |$\mbox{CS}(0,0)$| is very accurate except at the interface, and does not produce spurious oscillations in the solution. The most accurate scheme, |$\mbox{CS}(0,0.34)$|⁠, models the moving interface extremely well, but produces a very small oscillation in the error at nearby points. Fig. 9. View largeDownload slide Results at |$T=10$| for equation (4.1) subject to (4.18), setting |$\varDelta x=0.05$|⁠. Exact solution (solid curve), |$\textrm{CS}(0,\beta )$| with |$\varDelta t=0.025$| (crosses) and ML/IM with |$\varDelta t=8$|e-05 (diamonds). Top: |$\beta =0$|⁠. Bottom: |$\beta =0.34$|⁠. Fig. 9. View largeDownload slide Results at |$T=10$| for equation (4.1) subject to (4.18), setting |$\varDelta x=0.05$|⁠. Exact solution (solid curve), |$\textrm{CS}(0,\beta )$| with |$\varDelta t=0.025$| (crosses) and ML/IM with |$\varDelta t=8$|e-05 (diamonds). Top: |$\beta =0$|⁠. Bottom: |$\beta =0.34$|⁠. The final benchmark problem is (4.1) with the initial and boundary conditions \begin{align} \begin{aligned} &u(x,0)=0, \quad 0\leqslant x\leqslant 5,\\ &u(0,t)=(\tilde{t}-t)^{-1}\left[1-\left(1-t/\tilde t\,\right)^{2/3}\right],\qquad u(5,t)=0,\quad t\in[0,T]\subset[0,\tilde t\,), \end{aligned}\end{align} (4.19) where |$\tilde{t}$| is a positive constant (the time of existence of the solution). The exact solution of this problem is (see Galaktionov & Posashkov, 1988; Ibragimov, 1994) \begin{equation*} u_{\textrm{exact}}(x,t)=\begin{cases} (\tilde{t}-t)^{-1}\left[\left(1-x/\sqrt{6}\right)^2-\left(1-t/\tilde t\,\right)^{2/3}\right],&\quad 0\leqslant x \leqslant \tilde x(t),\\ 0,&\quad \tilde x(t)<x\leqslant 5, \end{cases} \end{equation*} where \begin{equation*}\tilde{x}(t)=\sqrt{6}(1-(1-t/\tilde{t}\,)^{1/3}),\qquad 0\leqslant t\leqslant T<\tilde{t}.\end{equation*} This is another nonsmooth solution of equation (4.1); it exists on the finite time interval |$[0,\tilde t\,)$|⁠, blowing up when |$t$| approaches |$\tilde t$|⁠. Tables 9 and 10 summarize the numerical solutions of the nonlinear heat equation over the time interval |$[0,10]$|⁠, with the initial and boundary conditions in (4.19) for |$\tilde{t}=11$|⁠. For the coarser grid with |$\varDelta x=0.25$| the optimal value of |$\beta $| is |$0.06$|⁠; for the finer grid, |$\varDelta x=0.1$|⁠, |$\beta =0.05$| is optimal. On both grids we have chosen |$\varDelta t>\varDelta x$| for the |$\mbox{CS}(0,\beta )$| schemes to show that this does not produce instability. This contrasts markedly with ML/IM, which requires |$\varDelta t<\varDelta x^2$| for convergence. Once again |$\mbox{CS}(0,0)$| shows itself to be a highly efficient, reasonably accurate scheme. Table 9 Errors in solution and conservation laws at the final time |$T=10$|⁠, when solving (4.1) with conditions in (4.19), using CS|$(0,\beta )$| methods with |$\varDelta x=0.25$|⁠, |$\varDelta t=0.333$| and ML/IM with |$\varDelta x=0.25$|⁠, |$\varDelta t=0.04$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 3.05e–16 1.39e–17 0.0162 0.013 |$\mbox{CS}(0,-1/8)$| 8.33e–17 1.39e–17 0.0099 0.013 |$\mbox{CS}(0,0)$| 2.78e–17 1.39e–17 0.0036 0.002 |$\mbox{CS}(0,0.06)$| 5.55e–17 2.08e–17 0.0017 0.011 ML/IM 0.0099 0.0117 0.0165 0.076 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 3.05e–16 1.39e–17 0.0162 0.013 |$\mbox{CS}(0,-1/8)$| 8.33e–17 1.39e–17 0.0099 0.013 |$\mbox{CS}(0,0)$| 2.78e–17 1.39e–17 0.0036 0.002 |$\mbox{CS}(0,0.06)$| 5.55e–17 2.08e–17 0.0017 0.011 ML/IM 0.0099 0.0117 0.0165 0.076 View Large Table 9 Errors in solution and conservation laws at the final time |$T=10$|⁠, when solving (4.1) with conditions in (4.19), using CS|$(0,\beta )$| methods with |$\varDelta x=0.25$|⁠, |$\varDelta t=0.333$| and ML/IM with |$\varDelta x=0.25$|⁠, |$\varDelta t=0.04$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 3.05e–16 1.39e–17 0.0162 0.013 |$\mbox{CS}(0,-1/8)$| 8.33e–17 1.39e–17 0.0099 0.013 |$\mbox{CS}(0,0)$| 2.78e–17 1.39e–17 0.0036 0.002 |$\mbox{CS}(0,0.06)$| 5.55e–17 2.08e–17 0.0017 0.011 ML/IM 0.0099 0.0117 0.0165 0.076 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 3.05e–16 1.39e–17 0.0162 0.013 |$\mbox{CS}(0,-1/8)$| 8.33e–17 1.39e–17 0.0099 0.013 |$\mbox{CS}(0,0)$| 2.78e–17 1.39e–17 0.0036 0.002 |$\mbox{CS}(0,0.06)$| 5.55e–17 2.08e–17 0.0017 0.011 ML/IM 0.0099 0.0117 0.0165 0.076 View Large Table 10 Errors in solution and conservation laws at the final time |$T=10$|⁠, when solving (4.1) with conditions in (4.19), using CS|$(0,\beta )$| methods with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.133$| and ML/IM with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.00667$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 1.78e–16 6.66e–17 0.0030 0.05 |$\mbox{CS}(0,-1/8)$| 2.22e–16 4.44e–17 0.0018 0.04 |$\mbox{CS}(0,0)$| 1.33e–16 1.33e–16 7.01e–04 0.01 |$\mbox{CS}(0,0.05)$| 5.33e–16 4.44e–17 5.23e–04 0.04 ML/IM 0.0036 0.0044 0.0088 0.91 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 1.78e–16 6.66e–17 0.0030 0.05 |$\mbox{CS}(0,-1/8)$| 2.22e–16 4.44e–17 0.0018 0.04 |$\mbox{CS}(0,0)$| 1.33e–16 1.33e–16 7.01e–04 0.01 |$\mbox{CS}(0,0.05)$| 5.33e–16 4.44e–17 5.23e–04 0.04 ML/IM 0.0036 0.0044 0.0088 0.91 View Large Table 10 Errors in solution and conservation laws at the final time |$T=10$|⁠, when solving (4.1) with conditions in (4.19), using CS|$(0,\beta )$| methods with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.133$| and ML/IM with |$\varDelta x=0.1$|⁠, |$\varDelta t=0.00667$| Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 1.78e–16 6.66e–17 0.0030 0.05 |$\mbox{CS}(0,-1/8)$| 2.22e–16 4.44e–17 0.0018 0.04 |$\mbox{CS}(0,0)$| 1.33e–16 1.33e–16 7.01e–04 0.01 |$\mbox{CS}(0,0.05)$| 5.33e–16 4.44e–17 5.23e–04 0.04 ML/IM 0.0036 0.0044 0.0088 0.91 Method |$\textrm{Err}_1$| |$\textrm{Err}_2$| Solution error Computational time |$\mbox{CS}(0,-1/4)$| 1.78e–16 6.66e–17 0.0030 0.05 |$\mbox{CS}(0,-1/8)$| 2.22e–16 4.44e–17 0.0018 0.04 |$\mbox{CS}(0,0)$| 1.33e–16 1.33e–16 7.01e–04 0.01 |$\mbox{CS}(0,0.05)$| 5.33e–16 4.44e–17 5.23e–04 0.04 ML/IM 0.0036 0.0044 0.0088 0.91 View Large Fig. 10. View largeDownload slide Results for equation (4.1) with conditions in (4.19). Top: numerical solution given by method |$\mbox{CS}(0,0)$|⁠, setting |$\varDelta x=0.1$| and |$\varDelta t=0.133$| over |$[0,5]$| at time |$T=10$|⁠. Bottom: exact solution (solid curve), numerical solutions close to the interface from |$\mbox{CS}(0,0)$| with |$\varDelta x=0.1$| and |$\varDelta t=0.133$| (crosses) and ML/IM with |$\varDelta x=0.1$| and |$\varDelta t=0.00667$| (diamonds). Fig. 10. View largeDownload slide Results for equation (4.1) with conditions in (4.19). Top: numerical solution given by method |$\mbox{CS}(0,0)$|⁠, setting |$\varDelta x=0.1$| and |$\varDelta t=0.133$| over |$[0,5]$| at time |$T=10$|⁠. Bottom: exact solution (solid curve), numerical solutions close to the interface from |$\mbox{CS}(0,0)$| with |$\varDelta x=0.1$| and |$\varDelta t=0.133$| (crosses) and ML/IM with |$\varDelta x=0.1$| and |$\varDelta t=0.00667$| (diamonds). Figure 10 compares the numerical solutions from |$\mbox{CS}(0,0)$| and ML/IM on the finer grid. Again, |$\mbox{CS}(0,0)$| is very close to the true solution (as are the other CS schemes in the tables), whereas ML/IM has a small lag close to the interface. 5. Conclusions and discussion Motivated by the basic principle of geometric integration that numerical schemes should preserve key structural features of the approximated problem to the extent that is possible, we have presented a strategy for developing finite difference methods that preserve two local conservation laws. This new strategy simplifies the approach introduced in Grant & Hydon (2013) and developed in Grant (2015). Depending on the stencil Grant’s method can have a very long symbolic computation time (typically several days on a fast PC for a 10-point stencil), which is a strong limitation. However, this difficulty can be overcome by restricting attention to schemes that are second order, with key terms that are as compact as possible. Such schemes can be determined by hand or by a short symbolic computation (of no more than a few minutes), even for larger stencils. We have developed new parameterized families of conservative numerical schemes for the solution of the KdV equation and a nonlinear heat equation. These schemes seem to be more robust and efficient than other well-known methods that do not preserve multiple conservation laws, perhaps due to topological and analytic advantages. Conservation laws have a topological origin as cohomology classes in the restricted variational bicomplex; schemes that preserve these retain discrete analogues of what may be essential topological features. Furthermore, parameters typically multiply terms that regularize the approximation of |${\mathcal{A}}$| in some way. By using benchmark problems and optimizing the parameters with respect to the solution error, we have found members of each family that are highly accurate. In practice the exact solution to a given problem is not usually known. Nevertheless, one can optimize the parameters numerically in order to achieve the best regularization for a given problem. Depending on the problem being approximated, it may be advantageous to choose the parameters in a way that best preserves other geometric structures, such as symplecticity, symmetries or further conservation laws. Acknowledgements We are grateful to our colleague John Pearson, University of Edinburgh, and to the anonymous referees, for their insightful remarks and constructive suggestions which have helped to improve this paper. Footnotes 1 A family of second-order schemes depending on eight free parameters can be found by removing the compactness assumption on the quadratic term in |$\widetilde{F_1}$|⁠. 2 The Hadamard product gives |$(A\circ B)_{i,\,j}=A_{i,\,j}B_{i,\,j},$| for each |$i,\,j.$| 3 If periodic or zero boundary conditions apply, (4.15) and (4.16) can be replaced with \begin{equation*}{\varDelta x}\max_{j=1,\ldots, N-1} \sum_{i=1}^{M}D_n\left(u_{i,1}+\alpha\varDelta x^2D_m^2u_{i-1,\,j}\right),\qquad{\varDelta x}\max_{j=1,\ldots, N-1} \sum_{i=1}^{M}{x_i}D_n\left(u_{i,\,j}+\alpha\varDelta x^2D_m^2u_{i-1,\,j}\right),\end{equation*} which measure the error in the conservation of the global invariants \begin{equation*}\int G_1\,\mathrm d x = \int u \,\mathrm d x \quad \mbox{and} \quad \int G_2\,\mathrm d x = \int xu \,\mathrm d x.\end{equation*} References Alonso , L. M. ( 1979 ) On the Noether map . Lett. Math. Phys. , 3 , 419 – 424 . Google Scholar Crossref Search ADS Ascher , U. M. & McLachlan , R. I. ( 2004 ) Multisymplectic box schemes and the Korteweg–de Vries equation . Appl. Numer. Math. , 48 , 255 – 269 . Google Scholar Crossref Search ADS Ascher , U. M. & McLachlan , R. I. ( 2005 ) On symplectic and multisymplectic scheme for the KdV equation . J. Sci. Comput. , 25 , 83 – 104 . Google Scholar Crossref Search ADS Bambusi , D. , Faou , E. & Grébert , B. ( 2013 ) Existence and stability of ground states for fully discrete approximations of the nonlinear Schrödinger equation . Numer. Math. , 123 , 461 – 492 . Google Scholar Crossref Search ADS Barletti , L. , Brugnano , L. , Frasca-Caccia , G. & Iavernaro , F. ( 2016 ) Recent advances in the numerical solution of Hamiltonian partial differential equations . AIP Conf. Proc. , 1776 , 020002 . Barletti , L. , Brugnano , L. , Frasca-Caccia , G. & Iavernaro , F. ( 2017 ) Solving the nonlinear Schrödinger equation using energy conserving Hamiltonian Boundary Value Methods . AIP Conf. Proc. , 1863 , 160002 . Barletti , L. , Brugnano , L. , Frasca-Caccia , G. & Iavernaro , F. ( 2018 ) Energy-conserving methods for the nonlinear Schrödinger equation . Appl. Math. Comput. , 318 , 3 – 18 . Betsch , P. & Steinmann , P. ( 2000 ) Inherently energy conserving time finite elements for classical mechanics . J. Comput. Phys. , 160 , 88 – 116 . Google Scholar Crossref Search ADS Bridges , T. J. ( 1997 ) Multisymplectic structures and wave propagation . Math. Proc. Camb. Philos. Soc. , 121 , 147 – 190 . Google Scholar Crossref Search ADS Bridges , T. J. & Reich , S. ( 2001 ) Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity. Phys. Lett . A. , 284 , 184 – 193 . Bridges , T. J. & Reich , S. ( 2006 ) Numerical methods for Hamiltonian PDEs . J. Phys. A. , 39 , 5287 – 5320 . Google Scholar Crossref Search ADS Brugnano , L. , Frasca-Caccia , G. & Iavernaro , F. ( 2015a ) Energy conservation issues in the numerical solution of the semilinear wave equation . Appl. Math. Comput. , 270 , 842 – 870 . Brugnano , L. , Frasca-Caccia , G. & Iavernaro , F. ( 2015b ) Energy conservation issues in the numerical solution of Hamiltonian PDEs . AIP Conf. Proc. , 1648 , 020002 . Brugnano , L. , Frasca-Caccia , G. & Iavernaro , F. ( 2015c ) Recent advances in the numerical solution of Hamiltonian PDEs . AIP Conf. Proc. , 1648 , 150008 . Brugnano , L. , Iavernaro , F. & Trigiante , D. ( 2015d ) Analysis of Hamiltonian boundary value methods (HBVMs): a class of energy-preserving Runge–Kutta methods for the numerical solution of polynomial Hamiltonian systems . Commun. Nonlinear Sci. Numer. Simul. , 20 , 650 – 667 . Google Scholar Crossref Search ADS Brugnano , L. & Iavernaro , F. ( 2016 ) Line Integral Methods for Conservative Problems . Monograph and Research Notes in Mathematics . Boca Raton, FL : CRC Press . Brugnano , L. , Iavernaro , F. & Trigiante , D. ( 2010 ) Hamiltonian boundary value methods (energy preserving discrete line integral methods) . JNAIAM. J. Numer. Anal. Ind. Appl. Math. , 5 , 17 – 37 . Brugnano , L. , Iavernaro , F. & Trigiante , D. ( 2012 ) A simple framework for the derivation and analysis of effective one-step methods for ODEs . Appl. Math. Comput. , 218 , 8475 – 8485 . Buchberger , B. & Kauers , M. ( 2010 ) Groebner basis . Scholarpedia , 5 , 7763 . Google Scholar Crossref Search ADS Buchberger , B. & Kauers , M. ( 2011 ) Buchbergers’s algorithm . Scholarpedia , 6 , 7764 . Google Scholar Crossref Search ADS Budd , C. J. & Piggott , M. D. ( 2003 ) Geometric integration and its applications . Handb. Numer. Anal. , 11 , 35 – 139 . Cano , B. ( 2006 ) Conserved quantities of some Hamiltonian wave equations after full discretization . Numer. Math. , 103 , 197 – 223 . Google Scholar Crossref Search ADS Celledoni , E. , McLachlan , R. I. , McLaren , D. I. , Owren , B. , Quispel , G. R. W. & Wright , W. M. ( 2009 ) Energy-preserving Runge–Kutta methods . M2AN Math. Model. Numer. Anal. , 43 , 645 – 649 . Google Scholar Crossref Search ADS Chen , J. B. , Qin , M. Z. & Tang , Y. F. ( 2002 ) Symplectic and multi-symplectic methods for the nonlinear Schrödinger equation . Comput. Math. Appl. , 43 , 1095 – 1106 . Google Scholar Crossref Search ADS Cox , D. , Little , J. & O’Shea , D. ( 1992 ) Ideals, Varieties and Algorithms. An Introduction to Computational Algebraic Geometry and Commutative Algebra . New York : Springer . Dahlby , M. & Owren , B. ( 2011 ) A general framework for deriving integral preserving numerical methods for PDEs . SIAM J. Sci. Comput. , 33 , 2318 – 2340 . Google Scholar Crossref Search ADS De Frutos , J. & Sanz-Serna , J. M. ( 1997 ) Accuracy and conservation properties in numerical integration: the case of the Korteweg–de Vries equation . Numer. Math. , 75 , 421 – 445 . Google Scholar Crossref Search ADS Drazin , P. G. & Johnson , R. S. ( 1989 ) Solitons: an Introduction . Cambridge : Cambridge University Press . Durán , A. & López-Marcos , M. A. ( 2003 ) Conservative numerical methods for solitary wave interactions . J. Phys. A Math. Gen. , 36 , 7761 – 7770 . Google Scholar Crossref Search ADS Durán , A. & Sanz-Serna , J. M. ( 1998 ) The numerical integration of relative equilibrium solutions. Geometric theory . Nonlinearity , 11 , 1547 – 1567 . Google Scholar Crossref Search ADS Durán , A. & Sanz-Serna , J. M. ( 2000 ) The numerical integration of relative equilibrium solutions . The nonlinear Schrödinger equation. IMA J. Numer. Anal. , 20 , 235 – 261 . Google Scholar Crossref Search ADS Feng , K. ( 1985 ) On difference schemes and symplectic geometry . Proceedings of the 1984 Beijing Symposium on Differential Geometry and Differential Equations . Beijing : Sci. Press Beijing , pp. 42 – 58 . Frasca-Caccia , G. ( 2015 ) A new efficient implementation for HBVMs and their application to the semilinear wave equation. Ph.D . Thesis. Università degli Studi di Firenze , Italy . Frasca-Caccia , G. ( 2018 ) Bespoke finite difference methods that preserve two local conservation laws of the modified KdV equation . AIP Conf. Proc , (in press) arXiv:1808.09370 . Furihata , D. ( 1999 ) Finite difference schemes for |$\partial u/\partial t={\left (\partial /\partial x\right )}^{\alpha }\delta G/\delta u$| that inherit energy conservation or dissipation property . J. Comput. Phys , 156 , 181 – 205 . Google Scholar Crossref Search ADS Galaktionov , V. A. & Posashkov , S. A. ( 1988 ) A method for investigating unbounded solutions of quasilinear parabolic equations . Zh. Vychisl. Mat. Mat. Fiz. , 28 , 842 – 854 . Gonzales , O. ( 1996 ) Time integration and discrete Hamiltonian systems . J. Nonlinear Sci. , 6 , 449 – 467 . Google Scholar Crossref Search ADS Grant , T. J. ( 2011 ) Characteristics of conservation laws for difference equations . Ph.D Thesis , University of Surrey , UK . Grant , T. J. ( 2015 ) Bespoke finite difference schemes that preserve multiple conservation laws. LMS . J. Comput. Math. , 18 , 372 – 403 . Grant , T. J. & Hydon , P. E. ( 2013 ) Characteristics of conservation laws for difference equations. Found. Comput . Math. , 13 , 667 – 692 . Guan , H. , Jiao , Y. , Liu , J. & Tang , Y. ( 2009 ) Explicit symplectic methods for the nonlinear Schrödinger equation . Commun. Comput. Phys. , 6 , 639 – 654 . Guo , L. & Xu , Y. ( 2015 ) Energy conserving local discontinuous Galerkin methods for the nonlinear Schrödinger equation with wave operator . J. Sci. Comput. , 65 , 622 – 647 . Google Scholar Crossref Search ADS Hairer , E. ( 2010 ) Energy-preserving variant of collocation methods . JNAIAM. J. Numer. Anal. Ind. Appl. Math. , 5 , 73 – 84 . Hairer , E. , Lubich , C. & Wanner , G. ( 2006 ) Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations , 2nd edn. Berlin : Springer . Hydon , P. E. ( 2001 ) Conservation laws of partial difference equations with two independent variables . J. Phys. A. , 34 , 10347 – 10355 . Google Scholar Crossref Search ADS Hydon , P. E. ( 2014 ) Difference Equations by Differential Equation Methods . Cambridge : Cambridge University Press . Hydon , P. E. & Mansfield , E. L. ( 2004 ) A variational complex for difference equations . Found. Comput. Math. , 4 , 187 – 217 . Google Scholar Crossref Search ADS Ibragimov , N. H. ( 1994 ) CRC Handbook of Lie Group Analysis of Differential Equations , vol. 1. Boca Raton, FL : CRC Press . Islas , A. L. , Karpeev , D. A. & Schober , C. M. ( 2001 ) Geometric integrators for the nonlinear Schrödinger equation . J. Comput. Phys. , 173 , 116 – 148 . Google Scholar Crossref Search ADS Islas , A. L. & Schober , C. M. ( 2004 ) On the preservation of phase space structure under multisymplectic discretization . J. Comput. Phys. , 197 , 585 – 609 . Google Scholar Crossref Search ADS Koide , S. & Furihata , D. ( 2009 ) Nonlinear and linear conservative finite difference schemes for regularized long wave equation . Japan J. Indust. Appl. Math. , 26 , 15 – 40 . Google Scholar Crossref Search ADS Kupershmidt , B. A. ( 1985 ) Discrete Lax Equations and Differential-Difference Calculus . Astérisque, 123. Paris : Société Mathématique de France . Leimkuhler , B. & Reich , S. ( 2004 ) Simulating Hamiltonian Dynamics , vol. 14. Cambridge : Cambridge University Press . Lu , X. & Schmid , R. ( 1997 ) A symplectic algorithm for wave equations . Math. Comput. Simulation , 43 , 29 – 38 . Google Scholar Crossref Search ADS Mansfield , E. L. ( 1992 ) Differentaial Groebner bases. Ph.D Thesis, University of Sydney , Australia . McLachlan , R. I. & Quispel , G. R. W. ( 2014 ) Discrete gradient methods have an energy conservation law . Discrete Contin. Dyn. Syst. , 34 , 1099 – 1104 . Google Scholar Crossref Search ADS McLachlan , R. I. , Quispel , G. R. W. & Robidoux , N. ( 1999 ) Geometric integration using discrete gradient . R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci. , 357 , 1021 – 1045 . Google Scholar Crossref Search ADS Oliver , M. , West , M. & Wulff , C. ( 2004 ) Approximate momentum conservation for spatial semidiscretization of semilinear wave equations . Numer. Math. , 97 , 493 – 535 . Google Scholar Crossref Search ADS Olver , P. J. ( 1993 ) Application of Lie Groups to Differential Equations , 2nd edn. Graduate Texts in Mathematics , vol. 107. New York : Springer . Qin , M. Z. & Zhang , M. Q. ( 1990 ) Multi-stage symplectic schemes of two kinds of Hamiltonian systems for wave equations . Comput. Math. Appl. , 19 , 51 – 62 . Quispel , G. R. W. & McLaren , D. I. ( 2008 ) A new class of energy-preserving numerical integration methods . J. Phys. A , 41 , 045206 . Sanz Serna , J. M. ( 1998 ) Runge–Kutta schemes for Hamiltonian systems . BIT , 28 , 877 – 883 . Google Scholar Crossref Search ADS Sanz Serna , J. M. & Calvo , M. P. ( 1994 ) Numerical Hamiltonian Problems . London : Chapman & Hall . Sun , J. Q. & Qin , M. Z. ( 2003 ) Multi-symplectic methods for the coupled 1D nonlinear Schrödinger system . Comput. Phys. Comm. , 155 , 221 – 235 . Google Scholar Crossref Search ADS Sun , Y. J. & Qin , M. Z. ( 2004 ) A multi-symplectic scheme for RLW equation . J. Comput. Math. , 22 , 611 – 621 . Tang , Q. & Chen , C. ( 2007 ) Continuous finite element methods for Hamiltonian systems . Appl. Math. Mech. , 28 , 1071 – 1080 . Google Scholar Crossref Search ADS Tang , W. & Sun , Y. ( 2012 ) Time finite elements methods: a unified framework for numerical discretizations of ODEs . Appl. Math. Comput. , 219 , 2158 – 2179 . Wan , A. , Bihlo , A. & Nave , J. C. ( 2016 ) The multiplier method to construct conservative finite difference schemes for ordinary and partial differential equations . SIAM J. Numer. Anal. , 54 , 86 – 119 . Google Scholar Crossref Search ADS Zhang , Q. & Wu , Z. ( 2009 ) Numerical simulation for porous medium equation by local discontinuous Galerkin finite element method . J. Sci. Comput. , 38 , 127 – 148 . Google Scholar Crossref Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Simple bespoke preservation of two conservation laws
JF - IMA Journal of Numerical Analysis
DO - 10.1093/imanum/dry087
DA - 2018-12-14
UR - https://www.deepdyve.com/lp/oxford-university-press/simple-bespoke-preservation-of-two-conservation-laws-0iVUp7p9rI
SP - 1
VL - Advance Article
IS - 
DP - DeepDyve
ER -