A novel high-order symmetric and energy-preserving continuous-stage Runge-Kutta-Nyström Fourier pseudo-spectral scheme for solving the two-dimensional nonlinear wave equation

1.   Introduction

The nonlinear partial differential equations have significant roles in a variety of fields in engineering and science (see, e.g., ^[4,39]), including quantum field theory, nonlinear optics, propagation of dislocations in crystals, nucleation, and solid state physics. In this paper, we consider the following nonlinear wave equation in two dimensions:

subject to the initial conditions

and the periodic boundary conditions

where $\kappa^2$ is a dimensionless positive parameter, $\varphi_0(x, y)$ and $\varphi_1(x, y)$ are the given $(L_1, L_2)$- periodic functions, and $L_1$ and $L_2$ are the basic positive periods. In the literature, many works have been made to explore the analytical solution for the nonlinear wave equations (see, e.g., ^[1,39]). However, it is difficult to obtain the general exact solutions for all the nonlinear wave equations. Therefore, the development of efficient and high-precision numerical methods for solving the two-dimensional nonlinear wave equations has became much more important. A great number of excellent numerical strategies have been proposed to study the nonlinear wave equations, including the finite difference methods (see, e.g., ^[12,17,18]), the finite element methods (see, e.g., ^[2,3,33]), the spectral methods ^[25], the domain decomposition methods ^[19], and the radial basis functions methods ^[11].

If the nonlinear function $f(u)$ is the negative derivative of a nonnegative function $F(u)$, i.e., $f(u) = -\frac{{\rm d}F(u)}{{\rm d}u}$, and the solution of (1.1) satisfies $\big(u, \frac{\partial u}{\partial t}\big)\in H^1(\Omega)\times L^2(\Omega)$, then the nonlinear wave Eqs (1.1)–(1.3) could conserve the energy

The energy conservation (1.4) is a significant property of the nonlinear wave equations, and plays prominent roles in investigating soliton theory. Under this case, the nonlinear wave equations like (1.1) are called nonlinear Hamiltonian wave equations. We know that the energy conservation along the exact flow is one most characteristic properties of the nonlinear Hamiltonian wave Eq (1.1). The energy-conserving numerical schemes usually yield correct physical phenomenons and numerical stability (see, e.g., ^[24,29,32]). Therefore, it will be meaningful to design suitable numerical schemes which could exactly preserve the discrete energy and symmetry of the nonlinear Hamiltonian wave Eq (1.1).

The development of energy-preserving numerical schemes for nonlinear Hamiltonian wave equations has garnered significant attention across various fields of mechanics. For example, Li et al. ^[21] proposed several finite difference schemes that preserve specific algebraic invariants of the nonlinear Klein-Gordon equations. Moreover, based on the concept of the discrete line integral method (see, e.g., ^[5,6]), L. Brugnano et al. developed the energy-preserving Hamiltonian boundary value methods (HBVMs) to solve the nonlinear Hamiltonian PDEs (see, e.g., ^[7,8]). The energy-preserving average vector field (AVF) method was initially developed for solving Hamiltonian ordinary differential equations (ODEs). Recently, the AVF method, when combined with appropriate spatial semi-discretization techniques, has been utilized to numerically investigate nonlinear Hamiltonian wave equations, thereby attracting significant attention from researchers. For instance, AVF finite element methods have been introduced to solve one-dimensional Hamiltonian wave equations (see ^[10,33]). In ^[31,32], combining the AVF method with the spatial fourth-order finite difference semidiscretisation, the authors developed energy-preserving methods for one- and two- dimensional Hamiltonian wave equations with Neumann boundary conditions. However, the previous schemes have only second-order accuracy in time. To enhance temporal accuracy, Hou et al. ^[17,18] integrated the fourth-order AVF temporal integrator with spatial compact finite difference (CFD) discretization to construct and analyze high-order energy-preserving schemes for solving one- and two-dimensional nonlinear wave equations with variable coefficients. This represents a significant advancement in the field of energy-preserving methods for nonlinear Hamiltonian wave equations. Building on these contributions, we aim to develop and analyze a high-order energy-preserving and symmetric scheme for two-dimensional nonlinear Hamiltonian wave equations by combining the continuous-stage Runge-Kutta-Nyström time integrator with Fourier pseudo-spectral spatial discretization.

The rest of the paper is organized as follows: In Section 2, the two-dimensional nonlinear wave Eq (1.1) will first be reformulated as an abstract infinite-dimensional separable Hamiltonian ODE system in an appropriate function space. Then, the application of a continuous-stage Runge-Kutta-Nyström time integrator to the yielded ODE's system to derive the time-stepping scheme is presented. The energy preservation and symmetry of the proposed time-stepping scheme will be investigated. Furthermore, by approximating the spatial differential operator using the two-dimensional Fourier pseudo-spectral method, we derive a fully discrete scheme. A rigorous analysis of the energy conservation properties of this scheme is then conducted. The error analysis demonstrates that the proposed scheme achieves sixth-order accuracy in the relatively low regularity function space $C^2\big([t_0, T], \mathcal{B}\big)$. Numerical experiments are presented in Section 4. Lastly, a concise conclusion is provided in Section 5.

2.   Temporal semi-discretisation: energy-preserving and symmetric time-stepping scheme

In this section, we will first represent the two-dimensional nonlinear wave Eqs (1.1)–(1.3) as an abstract nonlinear ODE on an appropriate infinite-dimensional Hilbert space. Then, we will develop and analyze a novel energy-conserving time-stepping scheme for the abstract ODE.

2.1.   Abstract Hamiltonian ODE's system

According to the analysis in references ^[30,37], by defining the mapping

we can express the nonlinear wave Eqs (1.1)–(1.3) as the following abstract ODE (e.g., ^[29,30,38]):

where $\mathcal{A}$ is the linear differential operator

and $\mathcal{B}$ is the infinite-dimensional Hilbert space

For any $\phi(x, y), \psi(x, y)\in \mathcal{B}$, we introduce the inner product and the norms

Then, by taking the inner product of the abstract ODE's system (2.1) with $u^{\prime}(t)$, we are able to find that system (2.1) can preserve the separable energy

where the kinetic energy part $\mathcal{H}_1[u^{\prime}(t)]$ and the potential energy part $\mathcal{H}_2[u(t)]$ are

respectively. Obviously, the energy $\mathcal{H}[u(t), u^{\prime}(t)]$ of the abstract ODE (2.1) is the same as the energy $E(t)$ of the two-dimensional nonlinear wave Eqs (1.1)–(1.3). Moreover, the abstract ODE's system (2.1) is actually a Hamiltonian system

where $v(t) = u'(t)$ and $\mathcal{S}$ is a skew-adjoint operator

In light of the definition of the variational derivatives, we are able to check that

The main purpose of this work is to design a suitable time-stepping scheme for the two-dimensional nonlinear wave Eqs (1.1)–(1.3), which could exactly preserve the separable energy $\mathcal{H}[u(t), u^{\prime}(t)]$ or $E(t)$. To achieve this purpose, the temporal discretization strategy will be first considered for the abstract ODE (2.1) in the infinite-dimensional function space.

2.2.   Formulation of the energy-preserving time-stepping scheme

For any positive integer $N$, we define the temporal mesh grid as

with time step size $\Delta t = (T-t_0)/N$, and introduce the following approximations:

Then, applying the energy-preserving integrators, which are proposed for the second-order Hamiltonian ordinary differential systems (see ^[22,26]), to the abstract ODEs (2.1), we can establish the time-stepping scheme for the two-dimensional nonlinear wave Eqs (1.1)–(1.3).

Definition 2.1. For any one temporal single step $t_n$ to $t_{n+1}$, a continuous-stage Runge-Kutta-Nyström (RKN) time-stepping scheme for the abstract ODE (2.1) is defined as

where the weight function $P_{3, 2}(\tau, \sigma)$ is a cubic binary polynomial of the form

Remark 2.1. In ^[22,26], the authors introduced a framework for energy-preserving continuous-stage RKN methods for solving second-order Hamiltonian ODEs. Drawing upon the methodologies proposed in ^[22,26], we develop a novel energy-preserving time-stepping scheme utilizing the weight function $P_{3, 2}(\tau, \sigma)$, and extend this scheme to the two-dimensional nonlinear wave Eqs (1.1)–(1.3). Furthermore, it is important to emphasize that the selection of the weight function $P_{3, 2}(\tau, \sigma)$ is not unique. Different choices of weight functions can result in numerical methods exhibiting varying accuracy.

Now, we focus on verifying the energy conservation of the continuous-stage RKN time-stepping scheme defined in Definition 2.1 for the two-dimensional nonlinear wave Eqs (1.1)–(1.3).

Theorem 2.1. The continuous-stage Runge-Kutta-Nyström time-stepping scheme (2.8) can exactly preserve the energy $\mathcal{H}[u(t), v(t)]$ of the two-dimensional nonlinear wave Eqs (1.1)–(1.3) or the infinite-dimensional abstract ODE's system (2.1), that is,

Proof. Noticing the form of the separable energy (2.3) and (2.4), we have

It follows from inserting the expression of $v^{n+1}$ into (2.10) and after a careful calculation that

Moreover, it is evident that

Substituting the expressions of $U_{\tau}^n$ into (2.12) leads to

Noticing the form of the weight function $P_{3, 2}(\tau, \sigma)$, we have

Therefore, Eq (2.13) can be simplified as

Comparing (2.11) with (2.14), we obtain

The conclusion of the theorem is confirmed. □

The symmetric time integration method usually exhibits superior long time computational behavior along the numerical flows (see Chapter V in ^[16]). We know that the two-dimensional nonlinear wave Eqs (1.1)–(1.3) are temporal reversible (see, e.g., ^[24,28,29]). Therefore, it will be significant to investigate the symmetry of the energy-preserving continuous stage RKN time-stepping scheme.

Theorem 2.2. The energy-preserving continuous-stage Runge-Kutta-Nyström time-stepping scheme (2.1) for solving the two-dimensional nonlinear wave Eqs (1.1)–(1.3) is temporal reversible.

Proof. According to the concept of the time reversible integration method (see Chapter V in ^[16]), and applying the following transformations

to the time-stepping scheme (2.1), we obtain the adjoint scheme

After a careful calculation, the last two equations of the adjoint scheme (2.15) can be rewritten as

Inserting Eq (2.16) into the first equation of (2.15), we obtain

The integral transformation $\tau = 1-\xi$ yields that

Therefore, we see that the adjoint scheme (2.15) is the same as scheme (2.8). That means the energy-preserving continuous stage RKN time-stepping scheme is symmetric or temporal reversible. □

Utilizing the variation-of-constants formula on the infinite-dimensional abstract ODE system (2.1), the exact solution of the abstract system (2.1) can be expressed as

Furthermore, it is easy to obtain from Eq (2.19) that

Inserting the exact solution $u(t)$ of the infinite-dimensional abstract ODE system (2.1) into the time-stepping scheme (2.8), we have

where the residual $R^n(\tau)$ is a function of $\tau\in[0, 1]$. Applying the Taylor expansion with integral remainder

to the nonlinear integrands appearing in (2.19) and the first equation of (2.21) results in

and

Comparing the formulae (2.23) with (2.24), and noticing $g(u) = -\mathcal{A}u+f(u)$, we can approximate the local residuals $R^n(\tau), \; \; 0\leq \tau\leq 1,$ in the following theorem.

Theorem 2.3. Suppose that the exact solution $u$ of the abstract ODE's system (2.1) satisfies $u\in C^2\big([t_0, T], \mathcal{B}\big)$ and the nonlinear function $f'\in L^{\infty}\big([t_0, T], \mathcal{B}\big)$. Then, the remainder $R^n(\tau)$ satisfies the estimations

where $C$ is a constant and independent of $\Delta t$.

Proof. Subtracting (2.24) from (2.23) and noticing $g(u) = -\mathcal{A}u+f(u)$, we obtain

where

and

It follows from the definition of the bilinear polynomial weight function $P_{3, 2}(\tau, \sigma)$ that

Therefore, it is easy to check that

Utilizing the Taylor expansion of $f(\cdot)$, i.e.,

in (2.28) and recalling Eq (2.29), we have

Inserting the results (2.30) and (2.31) into (2.26), and taking the $L^2$ norms on both sides of (2.29), it is easy to verify the estimated result of the theorem.

□

3.   Structure-preserving fully discrete scheme

In this section, by combing the Fourier pseudo-spectral spatial approximation with the continuous-stage RKN time-stepping scheme (2.8), we will construct the energy-preserving fully discrete scheme for the two-dimensional nonlinear wave Eqs (1.1) and (1.2).

Choose $M_1$ and $M_2$ to be any even integers, and define $\Delta x: = \frac{L_1}{M_1}$ and $\Delta y: = \frac{L_2}{M_2}$ as the spatial steps. Then, the temporal-spatial grid points are denoted as $\Omega^{N}_{M} = \Omega_{M}\times\Omega^{N}$, where the temporal grid $\Omega^{N}$ is given by (2.7), and the spatial grid $\Omega_{M}$ is defined as

The grid function space $\mathbb{V}_{M}$ defined on $\Omega_{M}$ is given by

For any $u = (u_{j, k})\in\mathbb{V}_{M},$ we can reformulate it as the vector form

Therefore, the vector space of the grid functions $\mathcal{V}_{M}$, which is identical to $\mathbb{V}_{M}$, can be presented as

In addition, the corresponding discrete inner product and norm are defined as

3.1.   Spatial semi-discretisation: Fourier pseudo-spectral method

Define the interpolation space as

where $g_j(x)$ and $g_k(y)$ are trigonometric polynomials

with $\mu_1 = \frac{2\pi}{L_1}, \; \; \mu_2 = \frac{2\pi}{L_2}$, and

Therefore, for any periodic function $u(x, y)\in L_p^2(\Omega)$, the interpolation operator $I_M : L_p^2(\Omega)\to \mathcal{S}_M^p$ is

where the Fourier coefficients $\tilde{u}_{k_1, k_2}$ are

Moreover, it is simple to check that

It follows from applying the differential operator $\mathcal{A}$ to the interpolation $I_Mu(x, y)$ that

where $I_{M_i}, i = 1, 2$ are the unity matrices, and $D_2^x = F_{M_1}^{H}\Lambda_1F_{M_1}$ and $D_2^y = F_{M_2}^{H}\Lambda_2F_{M_2}$ are the spectral differential matrices. Here, we should point out that $F_{M}$ is the discrete Fourier transform matrix with elements $(F_{M})_{j, k} = \frac{1}{\sqrt{M}}e^{\frac{-2\pi{\rm \textbf{i}}(j)(k)}{M}}, \; j, k = 0, 1, \dots, M-1,$ $F_{M}^H$ is the conjugate transformation matrix of $F_{M}$, and $\Lambda_1, \Lambda_2$ are the frequency matrices with entries

Thus, the spectral differential matrix $A$ for approximating the 2D differential operator $\mathcal{A}$ can be expressed as (see, e.g., ^[13,14,40])

Actually, the spectral differential matrix $A$ is a symmetric and semi-positive matrix, and $A\mathbf{u}$ can be fast computed by the two-dimensional FFT function ${\rm ifft2}\big(\Lambda.*{\rm fft2}(\rm\mathbf{u})\big)$ built in MATLAB.

Using the spectral differential matrix $A$ to approximate the differential operator $\mathcal{A}$, the two-dimensional nonlinear wave Eqs (1.1) and (1.2) or the abstract ODE system (2.1) can be converted into the semi-discrete system

where $\mathbf{u}(t)\in\mathcal{V}_M$ and $g\big(\mathbf{u}(t)\big) = -A\mathbf{u}(t)+f\big(\mathbf{u}(t)\big)$ with

Taking the discrete inner product on both sides of the semi-discrete system (3.7) with $\mathbf{u}^{\prime}(t)$, we have

which means

Therefore, we can conclude that the semi-discrete system (3.7) is energy-conserving.

Theorem 3.1. Suppose that $\mathbf{u}(t)\in\mathcal{V}_M$ is the solution of the semi-discrete system (3.7). Then, the semi-discrete system (3.7) can conserve the discrete energy

where the discrete kinetic energy $H_1[\mathbf{u}^{\prime}(t)]$ and discrete potential energy $H_2[\mathbf{u}(t)]$ are given by

Proof. The proof process is straightforward along the above analysis. Therefore, we omit the details. □

Remark 3.1. Actually, the energies (3.8) and (3.9) the discrete versions of the energies (2.3) and (2.4) of the the two-dimensional nonlinear wave Eqs (1.1) and (1.2) or the abstract system (2.1). Therefore, we will explore the fully discrete scheme which can exactly preserve the discrete energies (3.8) and (3.9) in this work.

3.2.   Energy-preserving fully discrete scheme

The main strategy of the construction of the energy-preserving fully discrete scheme is to approximate the differential operator $\mathcal{A}$ in the continuous-stage RKN time-stepping scheme (2.8) by the spectral differential matrix $A$. Therefore, the following theorem will show that the fully discrete scheme could preserve the discrete energy (3.8) and (3.9) exactly.

Theorem 3.2. By following the continuous stage RKN Fourier pseudo-spectral scheme

with $g(u) = -Au+f(u)$ and $\Delta t$ the time step size, the discrete energies (3.8) and (3.9) are conserved, i.e.,

Proof. We calculate the separable energy

Inserting $\mathbf{v}^{n+1}$ into $H_1[\mathbf{v}^{n+1}]$ and keeping that $A$ is a symmetric matrix in mind gives

On the other hand, we have

It follows from inserting

into (3.14) that

Combining the results of (3.12), (3.13), and (3.15), we have the desired result. □

Remark 3.2. Similar as the proof process of Theorem 2.2, it can be verified that the continuous stage RKN Fourier pseudo-spectral scheme (3.10) is time reversible. Moreover, we have noticed that the authors in ^[40] considered the integrating factor technique and the 4th-order (2-stage) Gauss-Legendre Runge-Kutta scheme to propose a symplectic time integration method for three-dimensional nonlinear water waves. This method can sufficiently use the oscillation generated by the spatial discretisation. Moreover, we have observed that the authors in ^[40] employed the integrating factor technique and the fourth-order (two-stage) Gauss-Legendre Runge-Kutta scheme to develop a symplectic time integration method for three-dimensional nonlinear water waves. This approach efficiently utilizes the oscillations resulting from spatial discretization, and typically yields accurate results at a reasonable computational cost. Perhaps, the combination of our proposed energy-preserving time integrator with the integrating factor technique could lead to a more efficient energy-preserving scheme. This will be considered in our future research.

Remark 3.3. In practice, the integrals in the fully discrete scheme usually cannot be easily calculated. Therefore, the $s$-point Gauss-Legendre formula $(b_i, c_i)_{i = 1}^s$ will be used to evaluate the integrals

Since the $s$-point GL quadrature formula is symmetric, the formula (3.16) is also symmetric.

To date, we have developed an energy-preserving fully discrete scheme for solving the two-dimensional nonlinear wave Eqs (1.1) and (1.2). This was achieved by initially semidiscretizing the temporal derivatives using a continuous-stage RKN method, followed by applying the Fourier spectral differential matrix $A$ to approximate the spatial differential operator $\mathcal{A}$. It has been observed that the Fourier pseudo-spectral method for approximating spatial derivatives can achieve spectral precision of order $\mathcal{O}(M^{-r})$ provided that the spatial regularity conditions are adequately satisfied.

Assume that $u(t)$ and $v(t)$ represent the exact solution and its derivative of the abstract ODE's system (2.1), while $\mathbf{u}(t)$ and $\mathbf{v}(t)$ denote the exact solution and its derivative of the semi-discrete system (3.7). Additionally, $\mathbf{u}^n$ and $\mathbf{v}^n$ signify the numerical solutions obtained from the continuous stage RKN Fourier pseudo-spectral scheme (3.10). It follows from inserting the exact solution $\mathbf{u}(t)$ of the semi-discrete system (3.7) into the continuous stage RKN Fourier pseudo-spectral scheme (3.10) that

where $\mathbf{R}^n(\tau)\in\mathcal{V}_{M}$ is the temporal local truncation error. Similar to the analysis of Theorem 2.3, we obtain the estimation for the residual $\mathbf{R}^n(\tau)$ in the following theorem.

Theorem 3.3. Suppose that the semi-discrete system (3.7) is well-posed and satisfies $\mathbf{u}(t)\in C^2\big([t_0, T]\big)$ and the nonlinear function $f'\in L^{\infty}\big([t_0, T]\big)$. Then, the local truncation error $\mathbf{R}^n(\tau)$ could be estimated as

where $\tilde{C}$ is a constant and independent of $\Delta t$.

Proof. The details of the proof are similar to the process of Theorem 2.3, so we omit the details. □

Letting

and subtracting (3.10) from (3.17), we obtain

We suppose the two-dimensional nonlinear wave Eqs (1.1) and (1.2) are well-posed. Subsequently, we present the error estimation for the fully discrete scheme (3.10) as detailed in the following theorem.

Theorem 3.4. If the exact solution $u(x, y, t)$ of the two-dimensional nonlinear wave Eqs (1.1) and (1.2) satisfies $u(x, y, t)\in C^2\big([t_0, T], \mathcal{B}\big)$, and the nonlinear function $f(\cdot)$ satisfies $f'\in L^{\infty}\big([t_0, T], \mathcal{B}\big)$, then under the limitation of $0 < \Delta t\leq h_0$ with a sufficiently small $h_0$ such that $h_0\|A\| < 1$, we obtain the error bounds

Here, we should point out that $A\lesssim B$ means there is a constant $C$ such that $A\leq CB$, and $M = M_1 = M_2$ is the spatial grid scale. Moreover, the constant $C$ depends on $T$, but is independent of $M, \|A\|$, and $\Delta t$.

Proof. The concept of the temporal-spatial error splitting method suggests that

Therefore, to obtain the accuracy of the fully discrete scheme, it is essential to concentrate on the analysis of temporal accuracy. Taking norms on both sides of (3.20) leads to

Then, under the restriction of the time step size $\Delta t\leq h_0$ with sufficiently small $h_0$ satisfying $h_0\|A\| < 1$, the first inequality of (3.23) results in

Summing up the last two inequalities of (3.23), we have

Moreover, the third inequality of (3.23) results in

Combining (3.25) and (3.26), we obtain

The mathematical induction will be an efficient approach to prove the result of the theorem.

Step Ⅰ. Letting $n = 0$ in (3.24) and noticing $\mathbf{e}^0 = 0, \mathit{\boldsymbol{\eta}}^0 = 0$, we have

Furthermore, noticing the limitation of the time step size again, the inequality (3.25) leads to

Step Ⅱ. Now, we assume that the estimation (3.21) is valid for $1\leq n\leq m-1$. That is,

Then, by mathematical induction, we only need to verify that the estimation (3.21) is still valid for $n = m$. Setting $n = m-1$ in (3.27) and using the above assumptions leads to

Therefore, the proof of (3.21) is completed. □

Remark 3.4. The conclusion of Theorem 3.3 indicates that the continuous-stage RKN Fourier pseudo-spectral scheme has at least fourth-order accuracy in the temporal domain. Owing to the temporal reversibility of the scheme (3.10), the forthcoming numerical experiments demonstrate that the proposed energy-preserving continuous-stage RKN Fourier pseudo-spectral scheme (3.10) can achieve sixth-order convergence in time.

Remark 3.5. In general, the preservation of energy typically ensures the stability of the fully discrete scheme. The analysis process of Theorem 3.2 demonstrates that the proposed energy-preserving scheme (3.10) is unconditionally stable. However, according to the result presented in Theorem 3.4, it can be concluded that the continuous stage RKN Fourier pseudo-spectral scheme (3.10) exhibits convergence under the condition $0 < \Delta t \leq h_0$ with $h_0\|A\| < 1$. In fact, the constraint $h_0\|A\| < 1$ corresponds to the CFL condition, as the differential matrix is intrinsically linked to the spatial discretization scale.

4.   Numerical experiments

We observe that the weight function $P_{3, 2}(\tau, \sigma)$ of the the continuous stage RKN Fourier pseudo-spectral scheme defined in Definition 2.1 is a cubic binary polynomial. Hereafter, the continuous stage RKN Fourier pseudo-spectral scheme will be denoted as EP3-FP. In this section, we will calculate the two-dimensional Klein-Gordon equation and the two-dimensional sine-Gordon equation to verify the precision, the efficiency, and the energy preservation of the derived EP3-FP scheme. Additionally, the following energy-preserving time integrators are chosen for comparison:

• AVF: the energy-preserving second-order averaged vector field method (see, e.g., ^[15,16]);

• HEP3: the symmetric sixth-order energy-preserving integrator constructed by Hairer in ^[15];

• SRKN3: the continuous-stage symplectic RKN method of order six (see, e.g., ^[34,35,36]).

The fully discrete scheme is obtained after discretizing the spatial derivatives with the Fourier pseudo-spectral method. We compute the temporal convergence rate by the following formula:

where the global error ${\rm GE}(h)$ is the difference of the exact solution $U(T; h)$ with the numerical solution $u(T; h)$ at time $T$ with step $h$. Moreover, it is known that the exact solution of the two-dimensional sine-Gordon equation could not be explicitly represented. Therefore, we will use the posterior error (see, e.g., ^[9,20]) to calculate the convergence rate, i.e.,

Furthermore, it is important to emphasize that the energy-preserving continuous-stage RKN time-stepping method introduced in this work, as well as the numerical methods selected for comparison, are closely associated with nonlinear integrals. To approximate these nonlinear integrals, the four-point Gauss-Legendre quadrature formula will be employed in numerical simulation.

Problem 1. Consider the two-dimensional nonlinear periodic Klein-Gordon equation (see, e.g., ^[18])

with the right-hand function $f(x, y, t) = \cos\big(\pi(x+y-t)\big)+\cos^3\big(\pi(x+y-t)\big)$. The exact solution is given by

We choose the parameters as $a = 1, b = 1,$ and $c^2 = \frac{1}{2}$. In Table 1, we list the errors and the corresponding convergence rates of the proposed EP3-FP scheme by varying the spatial and temporal step sizes. In Figure 1, we set the spatial grid scales as $M = M_x = M_y = 64$, and the logarithms of the global errors $\log_{10}(GE)$ against different time steps and the CPU times are plotted in Figure 1(a) and Figure 1(b), respectively. The logarithms of the energy errors of the EP3-FP scheme are plotted in Figure 1(b), which show that the proposed scheme is energy-preserving. The numerical results in Table 1 and Figure 1(a) illustrate that the proposed EP3-FP scheme achieves sixth-order temporal accuracy. Figure 1(c) shows that the EP3 time-stepping scheme has better computational efficiency than the chosen numerical methods.

Problem 2. Consider the two-dimensional sine-Gordon equations (see, e.g., ^[26,27])

with the dimensionless parameter $\kappa = 1/20$, and the initial conditions

Suppose that the two-dimensional sine-Gordon Eq (4.4) is equipped with periodic boundary conditions. Some snapshots of the numerical solution by the EP3-FP scheme are shown in Figure 2. These results demonstrate that the proposed EP3-EP scheme can efficiently simulate the two-dimensional sine-Gordon Eq (4.4) in a relatively long time domain. Moreover, it can be clearly observed from Figure 2 that the ring soliton shrinks during the initial stage, and oscillations and radiations have emerged by $t = 34.2$. Furthermore, the graphs also illustrate that the pulse simulated by the 2D sine-Gordon equation exhibits periodic oscillation. These phenomena are indeed valid, as other numerical methods have been employed to simulate this problem and exhibit similar phenomena. Here, we do not display the graphs obtained by other numerical methods. The numerical data listed in Table 2 demonstrates the convergence rate of the proposed EP3-FP scheme by varying the spatial and temporal step sizes. In Figure 3, after discretizing the spatial derivatives by the Fourier pseudo-spectral method with the fixed spatial mesh grid scales $M = M_x = M_y = 64$, the problem is calculated by a different time-stepping scheme. These phenomena further validate the accuracy, efficiency, and long-term energy conservation of the EP3-FP scheme.

5.   Conclusions

In this paper, based on the blend of the energy-preserving continuous-stage RKN time integrator with the Fourier pseudo-spectral spatial discretization, we presented a novel energy-preserving and symmetric fully discrete scheme for solving the two-dimensional nonlinear wave equations. The discrete energy of the two-dimensional nonlinear wave Eqs (1.1) and (1.2) is well conserved by the proposed scheme. Meanwhile, another significant discovery is that the derived EP3-FP scheme can achieve sixth-order temporal accuracy under the low regularity assumption $u\in C^2\big([t_0, T], \mathcal{B}\big)$. Numerical experiments are presented to illustrate the accuracy, efficiency, and long-term energy conservation of the derived EP3-FP scheme.

In light of a similar process, the derived EP3-FP scheme could be generalized to investigate other Hamiltonian PDEs, including the fractional nonlinear Hamiltonian wave equation, the Klein-Gordon equation with weak nonlinearity, the Klein-Gordon equation in the nonrelativistic limit regime, and the Klein-Gordon-Zakharov system.

Author contributions

Dongjie Gao wrote the main manuscript. Peiguo Zhang prepared numerical experiments. Longqin Wang, Zhenlong Dai and Yonglei Fang made the corrections. All the authors have agreed and given their consent for the publication of this research paper.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The work is supported by the Natural Science Foundation of China under Grant 11801280 and 12071419, the Natural Science Foundation of the Jiangsu Higher Education Institutions under Grant 21KJD110002, the Natural Science Foundation of Shandong Province under Grant ZR2024MA056, the foundation of innovative science and technology for youth in universities of Shandong Province under Grant 2023KJ278.

Conflict of interest

The authors declare no conflict of interest.