High-order compact di ff erence methods for solving two-dimensional nonlinear wave equations

: Nonlinear wave equations are widely used in many areas of science and engineering. This paper proposes two high-order compact (HOC) di ff erence schemes with convergence orders of O (cid:16) τ 4 + h 4 x + h 4 y (cid:17) that can be used to solve nonlinear wave equations. The first scheme is a nonlinear compact di ff erence scheme with three temporal levels. After calculating the second-order spatial derivatives of the previous time level using the Pad´e scheme, numerical solutions of the next time level are obtained through repeated iterations. The second scheme is a three-level linearized compact difference scheme. Unlike the first scheme, iterations are not required and it obtains numerical solutions through an explicit calculation. The two proposed schemes are applied to solutions of the coupled sine-Gordon equations. Finally, some numerical experiments are presented to confirm the e ff ectiveness and accuracy of the proposed schemes.

When v (x, y) = 1 and the nonlinear term f (u, x, y, t) is a polynomial with respect to u, −sinh(u), or −sin(u), Eq (1.1) can be simplified to the corresponding Klein-Gordon equation, sinh-Gordon equation, or sine-Gordon equation [6]. These equations embody different nonlinear physical phenomena. Soliton waves, which is one of the most conspicuous solutions to the sine-Gordon and Klein-Gordon equations, occur in many physical processes [7]. In addition, Khusnutdinova and Pelinovsky [8] derived the coupled sine-Gordon equations reflecting certain physical and biological phenomena [9][10][11].
The widespread use of nonlinear wave equations has led to strong interest in their solutions. In the past decades, analytical solutions to problems with special initial boundary value conditions have been obtained using various analytical methods [12][13][14][15]. In addition, researchers have investigated the corresponding Cauchy problem and the properties of soliton solutions [16][17][18]. However, it is not possible to obtain analytical solutions to problems with arbitrary initial boundary value conditions. Thus, different numerical methods have been developed to solve nonlinear wave equations, such as finite difference methods [6,[19][20][21][22], differential quadrature methods [3], finite element methods [23,24] and predictor-corrector schemes [25,26], etc.
In recent years, various high-order compact schemes have been used to solve nonlinear wave equations. In [19], two high-order compact alternating direction implicit (ADI) schemes with truncation errors of O τ 2 + h 4 have been proposed for the generalized sine-Gordon equation. One of them is a nonlinear scheme and the other is a linear scheme. A three-level compact ADI scheme with second-order convergence in the time direction and fourth-order convergence in space has been used to solve nonlinear wave equations [21], and the time accuracy can be increased to the fourth order through an extrapolation formula. Two-and three-level compact ADI schemes have been presented for solving the two-dimensional telegraph equations with nonlinear forcing [20]. Deng and Zhang [27] derived a compact multi-step ADI scheme to solve the nonlinear viscous wave equation by combining the second-order backward differential formula and the fourth-order Padé approximation, and finally applied a three-level Richardson extrapolation algorithm to reach fourth-order accuracy in the time direction. Deng [28] proposed two ADI schemes with convergence orders of O τ 2 + h 4 for solving nonlinear wave equations containing a viscous term. Deng and Liang [6] presented two compact ADI schemes with truncation errors of O τ 4 + h 4 for nonlinear wave equations. One of them is nonlinear and involves three temporal levels, while the other is a five-level linear scheme. In [22], an energypreserving average vector field compact difference scheme for nonlinear wave equations is described. This achieves fourth-order convergence in different directions. For the coupled nonlinear sine-Gordon equations, various numerical methods have been proposed in recent years. A numerical simulation method is proposed in [29], while Deng [30] established a high-order compact ADI scheme by combining the fourth-order compact difference operator and the Crank-Nicolson method. The coupled nonlinear sine-Gordon equations are solved by two compact ADI schemes in [6].
The existing finite difference schemes for solving two-dimensional nonlinear wave equations are typically based on the ADI method. At present, some schemes only have second-order convergence in the time direction, while other schemes have fourth-order convergence in both the temporal and spatial dimensions, but their computational efficiency is relatively poor. Explicit difference schemes offer high computational efficiency in solving partial differential equations [31][32][33]. Therefore, it is a challenging task to construct numerical methods with high temporal and spatial accuracy, and higher computational efficiency. So, in order to improve computational efficiency, this paper presents two high-order compact difference schemes based on the explicit difference approach. The remainder of this paper is structured as follows. In Section 2, two high-order compact (HOC) difference schemes for solving the nonlinear wave equation are presented. In Section 3, two high-order difference schemes are employed to simulate the coupled sine-Gordon equations. In Section 4, we present five numerical examples to verify the effectiveness and accuracy of the proposed HOC schemes. The conclusions to this study are presented in Section 5.

HOC difference schemes
In this section, two HOC difference schemes for simulating the nonlinear system of (1.1)-(1.4) are derived. First, we give some definitions and a lemma that will be useful in constructing the HOC difference schemes. We consider the domain [a 1 is divided into M y uniform intervals, a 2 = y 0 < y 1 < y 2 < · · · < y M y = b 2 , and (0, T ] is divided into N uniform intervals. The spatial steps are h x = b 1 −a 1 M x and h y = b 2 −a 2 M y , and the time step is τ = T N . Denoting each mesh point as (x i , y j , t n ), x i = a 1 + ih x , y j = a 2 + jh y , t n = nτ, i = 0, 1, · · · , M x , j = 0, 1, · · · , M y , n = 0, 1, · · · , N.
This section focuses on the construction of HOC difference schemes for the problem with Dirichlet boundary conditions (1.3).

HOC scheme for the start-up time step
Let ϑ (x, y) = v 2 (x, y). Equation (1.1) can be rewritten as (2.1) First, using the Taylor series expansion at the grid nodes x i , y j , t 0 yields According to Eqs (1.1) and (1.2), ∂ 2 u ∂t 2 x i , y j , 0 , ∂ 3 u ∂t 3 x i , y j , 0 , and ∂ 4 u ∂t 4 x i , y j , 0 can be easily calculated. For simplicity, we denote f u n i, j , x i , y j , t n := f u n i, j . Considering φ x i , y j , ψ x i , y j and (2.2), the following high-order approximation for u x i , y j , τ can be obtained: in which f t φ i, j and f tt φ i, j can be obtained by the chain rule as

HOC difference scheme for the other time steps
For the discretization of ∂ 2 u ∂t 2 , the central difference with truncated remainders yields Substituting (2.7) into (2.6) and discretizing the spatial derivatives with the central difference operator, we obtain ∂t 2 into the above formula with Eq (2.1) gives Then, is discretized using the central difference scheme. Rearranging (2.9) and neglecting the truncation error, we have the HOC difference scheme which we call HOC-I for convenience. Letting λ x = τ h x and λ y = τ h y , HOC-I scheme can be rewritten as Remark 1: For nonlinear wave equations with periodic boundaries ( Obviously, we note that the calculation of HOC-I requires the values of ∂ 2 U ∂x 2 n i, j and ∂ 2 U ∂y 2 n i, j to be known. These can be calculated by the following fourth-order Padé schemes [34]: where the boundaries can be calculated exactly by the following formulas: Noting that the two coefficient matrices constituted by (2.12)-(2.17) are all tridiagonal matrices, so we can efficiently solve these linear systems using the Thomas algorithm. can be calculated directly through (2.12) and (2.13) without calculating boundaries (2.14)-(2.17). The coefficient matrices can be written as Similarly, we can solve these linear systems using Thomas algorithm for quasi-tridiagonal matrices.
In this way, we obtain the HOC-I difference scheme with a truncation error of . This scheme achieves fourth-order accuracy in both the temporal and spatial dimensions. We can apply HOC-I using Algorithm 1.
Algorithm 1: Step 3 : Step 4 : Let n ← n + 1, repeat the above S tep 2 and S tep 3 till the time reaches the f inal moment, and the calculation is terminated. can be calculated directly through (2.12) and (2.13) for the problems with periodic boundaries.

Linearized HOC difference scheme
As the HOC-I scheme is nonlinear, multiple iterations are required at each time level, which increases the computation time. To avoid iterations and improve computational efficiency, we linearize the nonlinear term f u n+1 i, j of HOC-I to obtain a linear HOC scheme. According to Lemma 1, we have (2.18) Substituting (2.18) into (2.10) gives the following linear HOC difference scheme: which we call HOC-II for convenience. We can apply HOC-II using Algorithm 2.
Remark 5: For the stability of the compact difference scheme, we give the stability range for the corresponding linear problem using discrete Fourier analysis in Appendix.

For the coupled sine-Gordon equations
This section discusses the generalization of the two HOC schemes derived in Section 2 for solving the following coupled sine-Gordon equations [8][9][10][11]:  4) and the boundary conditions are

(3.6)
We use the same derivation method as for HOC-I to obtain the following difference scheme: The terms u 1 i, j and w 1 i, j can be derived as in Section 2. In addition, we generalize HOC-II for solving the above equations. The corresponding calculation procedures are similar to the above, and are omitted for brevity.

Numerical experiments
This section first presents three experiments to confirm the effectiveness and accuracy of the two proposed HOC schemes. We use ADI-I and ADI-II to represent the ADI schemes in [6], and ADI-III and ADI-IV are the difference schemes in [19] and [21], respectively. Note that the CPU time and numerical errors of these ADI schemes are calculated using the PHOEBE Solver (www.phoebesolver.com) on the same computer. Then, we simulate the motion states of a singlering soliton and a soliton in a layered medium. The following examples are programmed in Fortran 90 and executed on a laptop with an Intel Core i7-10750H CPU@2.60GHz 2.59 GHz and 16 GB of RAM.
Definition 1. The L ∞ and L 2 norm errors and Rate are defined as follows: where u(x i , y j , t N ) represents the analytical solution at point (x i , y j , t N ), u N i, j represents the numerical solution at point (x i , y j , t N ). In contrast, we introduce L ∞ = max{L ∞ (u) , L ∞ (w)} and L 2 = max{L 2 (u) , L 2 (w)} in Tables 7 and 8.

The sine-Gordon equation with constant coefficients [6, 19, 21]
First, we consider the following nonlinear sine-Gordon equation:  Table 1, we can see that when the spatial steps are 1/2, 1/2 2 , 1/2 3 and 1/2 4 , the errors calculated by the two difference schemes are essentially the same. When the spatial steps are 1/2 5 and 1/2 6 , the calculation errors of ADI-I are obviously larger than those of HOC-I. This is because the start-up level calculations for ADI-I have only third-order accuracy in the time direction, which does not match the fourth-order accuracy of the main difference scheme. This affects the calculation results. Furthermore, with the reduction in the spatial step size, the HOC-I scheme becomes significantly better than ADI-I in terms of computational efficiency. From Table 2, we can obtain the same conclusion. In addition, we find that HOC-II requires only half the CPU time of HOC-I to compute Problem 4.1 with the same stride length. This illustrates the importance of constructing the linearized difference scheme HOC-II. Table 3 presents the calculation results using the difference schemes in [21] and [19]. From the table, we find that the error produced by these difference schemes is greater than that calculated by the proposed method when the spatial step size is the same. The two HOC schemes described in this paper have better computational efficiency than the previous difference schemes considered here. According to the data in these tables, we plot the corresponding L ∞ error and L 2 error convergence order graphs in Figure 1. For both the L ∞ norm and the L 2 norm, the schemes proposed in this paper have the same fourth-order convergence speed as the previously developed schemes. This also means that the proposed schemes achieve fourth-order convergence in both the temporal and spatial dimensions, which accords with their theoretical accuracy. However, the convergence accuracy of ADI-I and ADI-II decreases in the fine mesh, reflecting the effects of the convergence accuracy of the start-up level. Furthermore, we substitute the exact solution into the source term to degenerate Problem 4.1 into a linear problem to verify the stability of difference scheme (1), and numerical results are listed in Table  4. It can be seen from Table 4 that the stability conditions we obtained through numerical experiments are consistent with the theoretical analysis.

The sine-Gordon equation with periodic boundary condition [22]
Next, we consider the sine-Gordon equation with periodic boundary condition: in which v 2 (x, y) = 1 + 0.1x + 0.1y, f (x, y, t) = 2π 2 t (1 + 0.1x + 0.1y) sin (πx) sin (πy) + sin t sin (πx) sin (πy) , and the analytical solution is u (x, y, t) = t sin (πx) sin (πy). Table 5 presents the numerical results using the HOC schemes proposed in this paper and the difference scheme in [22] under different spatial step sizes, where the time steps are τ = h/2 and τ = h, respectively. From the table, we can clearly see that the numerical errors calculated by the two difference schemes have little difference. Moreover, the computational results show that the proposed difference scheme achieves fourth-order convergence in both temporal and spatial dimensions for problems with variable coefficients. This agrees with our theoretical derivation and analysis. Table 6 presents the calculation results with the two difference schemes constructed in this paper. From the table, we can see that the calculated errors are almost the same for HOC-I and HOC-II, and both have fourthorder convergence. In addition, the HOC-II scheme requires about half the CPU time of HOC-I. This demonstrates the significance of the HOC-II scheme described in this paper. Table 5. L 2 norm error and convergence rate of different schemes at T = 2 for Problem 4.2.
EP-AVF(4)-CFD Scheme [22] HOC  Figure 2 shows the surfaces and the contour plots of the computational results using the HOC-I scheme when τ = 0.01 and h = 0.02. The numerical results do not oscillate. In addition, the exact solutions and numerical solutions have a high level of agreement, which means that the difference schemes in this paper have good simulation ability. [6,9,10,35] Consider the coupled sine-Gordon equations (3.1)-(3.6) in a square domain [−2, 2] × [−2, 2] and t ∈ (0, T ]. The set of analytical solutions is

Coupled sine-Gordon equations
To compare with the calculation results in Ref. [6], we take the same parameters α = √ 3, β = 0.5, η = 1.5. The initial conditions and boundary conditions are given by the analytical solutions.
We solve this problem using the difference schemes HOC-I and HOC-II at T = 1. For comparison, we also derive solutions using ADI-I and ADI-II [6]. Table 7 presents the numerical results calculated by HOC-I and ADI-I [6]. The results using HOC-I are smaller than the values in the literature for both the L ∞ norm error and the L 2 norm error. Additionally, our schemes achieve fourth-order convergence, which is consistent with the theoretical accuracy. ADI-I does not achieve fourth-order accuracy when the spatial step sizes are 1/2 5 and 1/2 6 , because its initial time stage has only third-order convergence instead of fourth-order convergence. In addition, the computational efficiency of HOC-I is significantly better than that of the difference scheme in the literature. Table 8 presents the calculation results for the two linearized difference schemes. Similarly, compared with ADI-II [6], the calculation results of HOC-II are better. We can also see that HOC-II is more computationally efficient than HOC-I, suggesting the superiority of the linearized HOC-II scheme. Figure 3 shows the surfaces of the numerical solutions u n i, j , w n i, j and the corresponding errors. From these numerical solution surfaces, we can see that the calculation results of the proposed difference scheme do not produce numerical oscillations. Therefore, the difference schemes constructed in this paper have good simulation ability.
The exact solution to this problem is unknown. We take time and space steps of 0.01 and 0.1, respectively, and then solve Problem 4.4 using HOC-II. Similar to Refs. [6,25,26], the surfaces of sin(U n /2) at different times are shown in Figure 4, and the contours at different time are plotted in Figure 5. From Figures 4 and 5, it is apparent that the ring soliton begins to shrink at t = 0, before radiating and oscillating outward in a proper sequence until another ring soliton is formed at about t = 10.5. The next cycle of evolution then begins. There are no oscillations in these groups of figures, which proves the effectiveness and accuracy of HOC-II in solving such problems.
The HOC-II scheme with time and space steps of 0.01 and 0.1 is applied to solve Problem 4.5. Figures 6 and 7 show the propagation of waves in a two-layered medium at different times. There are very complex interactions, such as contraction, radiation, and collisions between the two solitons. Furthermore, as the propagation velocity of the wave in the left medium is greater than that in the right medium, this allows us to see that the wave propagates faster in the left medium than in the right. Additionally, the greatest rate of change in u occurs at the collision center of the two solitons.

Conclusions
This paper has described the derivation of two HOC difference schemes with fourth-order accuracy for both temporal and spatial dimensions. These schemes are suitable for solving systems of the form of Eqs (1.1)-(1.4). One scheme is a nonlinear fourth-order compact scheme with three temporal levels, denoted as HOC-I, and the other is a linearized fourth-order compact scheme, denoted as HOC-II. In the linear HOC-II scheme, the function value of the next time level can be directly and explicitly calculated on the premise that the spatial derivative of the previous time level has been calculated.  Therefore, we can obtain numerical solutions of acceptable accuracy at a reasonable time cost. The two fourth-order compact difference schemes proposed in this paper have been applied to simulations of the coupled sine-Gordon equations. The accuracy and efficiency of the proposed schemes have been confirmed through a number of numerical examples, and it has been shown that HOC-II is significantly better than HOC-I and several previously developed difference schemes in terms of computational efficiency. Additionally, HOC-II has been used to simulate the motion of ring solitons and solitons in a layered medium.
In recent years, several researchers have proposed numerical methods for solving the threedimensional wave equations and different forms of the nonlinear fourth-order wave [36][37][38][39][40]. The method in this paper can be extended to solve such nonlinear wave equations, and related conclusions and results will be reported in the future.