High-order compact difference methods for solving two-dimensional nonlinear wave equations

1.   Introduction

Nonlinear wave equations are used in various fields of science and engineering, such as quantum field theory, nonlinear optics, fluid mechanics, and so on ^[1,2,3]. As the mathematical models of many complex real-world phenomena involve coefficients that vary in time and space ^[4,5], this paper considers the following two-dimensional nonlinear wave equation with variable coefficients:

in which $\Omega = \left[{{a_1}, {b_1}} \right] \times \left[{{a_2}, {b_2}} \right]$ and ${a_i}, {b_i} \in {\mathbb R} \left({i = 1, 2} \right)$. The initial conditions are given as

and the boundary conditions are

or the periodic boundary conditions are

where $L_1 = b_1-a_1$, $L_2 = b_2-a_2$, ${v}\left({x, y} \right)$ is the acoustic velocity and $f(u, x, y, t)$ is a nonlinear expressions in terms of $u$. Assume that $\varphi(x, y)$, $\phi(x, y)$, $g_0(y, t)$, $g_1(y, t)$, $l_0(x, t)$ and $l_1(x, t)$ are smooth known functions.

When ${v}\left({x, y} \right) = 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\left({{\tau ^2} + {h^4}} \right)$ 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\left({{\tau ^2} + {h^4}} \right)$ for solving nonlinear wave equations containing a viscous term. Deng and Liang ^[6] presented two compact ADI schemes with truncation errors of $O\left({{\tau ^4} + {h^4}} \right)$ 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 energy-preserving 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.

2.   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 $\left[{{a_1}, {b_1}} \right] \times \left[{{a_2}, {b_2}} \right]\times(0, T]$, in which $\left[{{a_1}, {b_1}} \right]$ is divided into $M_x$ uniform intervals, $a_1 = \mathop x\nolimits_0 < \mathop x\nolimits_1 < \mathop x\nolimits_2 < \cdot \cdot \cdot < \mathop x\nolimits_{M_x} = b_1$, $\left[{{a_2}, {b_2}} \right]$ is divided into $M_y$ uniform intervals, $a_2 = \mathop y\nolimits_0 < \mathop y\nolimits_1 < \mathop y\nolimits_2 < \cdot \cdot \cdot < \mathop y\nolimits_{M_y} = b_2$, and $(0, T]$ is divided into $N$ uniform intervals. The spatial steps are $h_x = \frac{b_1-a_1}{M_x}$ and $h_y = \frac{b_2-a_2}{M_y}$, and the time step is $\tau = \frac{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\tau$, $i = 0, 1, \cdots, M_x$, $j = 0, 1, \cdots, M_y$, $n = 0, 1, \cdots, N$.

We denote $\Omega_{\tau} = \{t_n|0\leq n\leq N\}$, $\Omega_h = \{(x_i, y_j)|0\leq i\leq M_x, 0\leq j\leq M_y\}$, $\Omega_{h\tau} = \Omega_h\times \Omega_{\tau}$. Define the grid function $v = \{v_{i, j}^n|0\leq i\leq M_x, 0\leq j\leq M_y, 0\leq n\leq N\}$ on $\Omega_{h\tau}$, and denote

Lemma 1. ^[6] Suppose that $q\left(t \right) \in {C^4}\left(I \right)$, $t \in I$. Then, we have $q\left({{t_{n + 1}}} \right) = 4q\left({{t_n}} \right) - 6q\left({{t_{n - 1}}} \right) + 4q\left({{t_{n - 2}}} \right) - q\left({{t_{n - 3}}} \right) + O\left({\Delta {t^4}} \right)$, where $\Delta t = {t_{n + k}} - {t_{n + k - 1}}, \left({k = - 2, - 1, 0, \; or\; 1} \right)$.

This section focuses on the construction of HOC difference schemes for the problem with Dirichlet boundary conditions (1.3).

2.1.   HOC scheme for the start-up time step

Let $\vartheta \left({x, y} \right) = {v^2}\left({x, y} \right)$. Equation (1.1) can be rewritten as

First, using the Taylor series expansion at the grid nodes $\left({{x_i}, {y_j}, {t_0}} \right)$ yields

According to Eqs (1.1) and (1.2), ${\frac{{{\partial ^2}u}}{{\partial {t^2}}}}\left({{x_i}, {y_j}, 0} \right)$, ${\frac{{{\partial ^3}u}}{{\partial {t^3}}}}\left({{x_i}, {y_j}, 0} \right)$, and ${\frac{{{\partial ^4}u}}{{\partial {t^4}}}}\left({{x_i}, {y_j}, 0} \right)$ can be easily calculated. For simplicity, we denote $f\left({u_{i, j}^n, {x_i}, {y_j}, {t_n}} \right): = f\left({u_{i, j}^n} \right)$. Considering $\varphi \left({{x_i}, {y_j}} \right)$, $\psi \left({{x_i}, {y_j}} \right)$ and (2.2), the following high-order approximation for $u\left({{x_i}, {y_j}, \tau } \right)$ can be obtained:

in which ${f_t}\left({\varphi _{i, j}} \right)$ and ${f_{tt}}\left({\varphi _{i, j}} \right)$ can be obtained by the chain rule as

2.2.   HOC difference scheme for the other time steps

For the discretization of ${\frac{{{\partial ^2}u}}{{\partial {t^2}}}}$, the central difference with truncated remainders yields

where $\left({R_t} \right)_{i, j}^n = O\left({{\tau ^4}} \right)$. According to Eq (2.1), we have

Substituting (2.7) into (2.6) and discretizing the spatial derivatives with the central difference operator, we obtain

in which $\left({R_1} \right)_{i, j}^n = O\left({{\tau ^4} + {\tau ^2}h_x^2 + {\tau ^2}h_y^2} \right)$.

Substituting ${\frac{{{\partial ^2}u}}{{\partial {t^2}}}}$ into the above formula with Eq (2.1) gives

Then, $\frac{{{\partial ^2}f\left({u_{i, j}^n} \right)}}{{\partial {t^2}}}$ 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-Ⅰ for convenience. Letting ${\lambda _x} = \frac{\tau }{{{h_x}}}$ and ${\lambda _y} = \frac{\tau }{{{h_y}}}$, HOC-Ⅰ scheme can be rewritten as

Remark 1: For nonlinear wave equations with periodic boundaries (1.4), there is ${U^n_{i, j}} = {U^n_{i+M_x, j}}$ for $i \le 0$, ${U^n_{i, j}} = {U^n_{i, j+M_y}}$ for $j \le 0$, ${U^n_{i, j}} = {U^n_{i-M_x, j}}$ for $i > M_x$ and ${U^n_{i, j}} = {U^n_{i, j-M_y}}$ for $j > M_y$. At the same time, the grid point range $\left\{ {\left({{x_i}, {y_j}, {t_n}} \right)|1 \le i \le M_x, 1 \le j \le M_y, 1 \le n \le N - 1} \right\}$ for high-order scheme (2.11) calculations.

Obviously, we note that the calculation of HOC-Ⅰ requires the values of $\left({\frac{{{\partial ^2}U}}{{\partial {x^2}}}} \right)_{i, j}^n$ and $\left({\frac{{{\partial ^2}U}}{{\partial {y^2}}}} \right)_{i, j}^n$ 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.

Remark 2: For nonlinear problems with periodic boundaries (1.4), $\left({\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right)_{i, j}^n$ and $\left({\frac{{{\partial ^2}u}}{{\partial {y^2}}}} \right)_{i, j}^n$ 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-Ⅰ difference scheme with a truncation error of $O\left({{\tau ^4} + {\tau ^2}h_x^2 + {\tau ^2}h_y^2 + h_x^4 + h_y^4} \right)$ for solving Eqs (1.1)–(1.4). This scheme achieves fourth-order accuracy in both the temporal and spatial dimensions. We can apply HOC-Ⅰ using Algorithm 1.

Remark 3: In the Step 2 of Algorithm 1, $\left({\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right)_{i, j}^n$ and $\left({\frac{{{\partial ^2}u}}{{\partial {y^2}}}} \right)_{i, j}^n$ can be calculated directly through (2.12) and (2.13) for the problems with periodic boundaries.

2.3.   Linearized HOC difference scheme

As the HOC-Ⅰ 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\left({u_{i, j}^{n + 1}} \right)$ of HOC-Ⅰ to obtain a linear HOC scheme. According to Lemma 1, we have

Substituting (2.18) into (2.10) gives the following linear HOC difference scheme:

which we call HOC-Ⅱ for convenience. We can apply HOC-Ⅱ using Algorithm 2.

Remark 4: In the Step 3 of Algorithm 2, $\left({\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right)_{i, j}^n$ and $\left({\frac{{{\partial ^2}u}}{{\partial {y^2}}}} \right)_{i, j}^n$ can be calculated directly through (2.12) and (2.13) for the problems with periodic boundaries.

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.

3.   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]:

in which $\Omega = \left[{{a_1}, {b_1}} \right] \times \left[{{a_2}, {b_2}} \right]$ and $\beta, \alpha > 0$ are constants. The initial conditions are

and the boundary conditions are

We use the same derivation method as for HOC-Ⅰ to obtain the following difference scheme:

The terms $u_{i, j}^1$ and $w_{i, j}^1$ can be derived as in Section 2.

In addition, we generalize HOC-Ⅱ for solving the above equations. The corresponding calculation procedures are similar to the above, and are omitted for brevity.

4.   Numerical experiments

This section first presents three experiments to confirm the effectiveness and accuracy of the two proposed HOC schemes. We use ADI-Ⅰ and ADI-Ⅱ to represent the ADI schemes in ^[6], and ADI-Ⅲ and ADI-Ⅳ 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 single-ring 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.60$GHz $2.59$ GHz and $16$ GB of RAM.

Definition 1. The $L_\infty$ 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_{i, j}^N$ represents the numerical solution at point $(x_i, y_j, t_N)$. In contrast, we introduce ${L_\infty } = \max \{ {L_\infty }\left(u \right), {L_\infty }\left(w \right)\}$ and ${L_2 } = \max \{ {L_2 }\left(u \right), {L_2 }\left(w \right)\}$ in Tables 7 and 8.

4.1.   The sine–Gordon equation with constant coefficients ^[6,19,21]

First, we consider the following nonlinear sine–Gordon equation:

The analytical solution is $u\left({x, y, t} \right) = 4\arctan \left({{e^{x + y - t}}} \right)$.

Tables 1–3 present the results obtained by computing Problem 4.1 with the various difference schemes. From 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-Ⅰ are obviously larger than those of HOC-Ⅰ. This is because the start-up level calculations for ADI-Ⅰ 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-Ⅰ scheme becomes significantly better than ADI-Ⅰ in terms of computational efficiency. From Table 2, we can obtain the same conclusion. In addition, we find that HOC-Ⅱ requires only half the CPU time of HOC-Ⅰ to compute Problem 4.1 with the same stride length. This illustrates the importance of constructing the linearized difference scheme HOC-Ⅱ. 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_{\infty}$ error and $L_{2}$ error convergence order graphs in Figure 1. For both the $L_{\infty}$ 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-Ⅰ and ADI-Ⅱ 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.

4.2.   The sine–Gordon equation with periodic boundary condition ^[22]

Next, we consider the sine–Gordon equation with periodic boundary condition:

in which ${v^2}\left({x, y} \right) = 1 + 0.1x + 0.1y$, $f\left({x, y, t} \right) = 2{\pi ^2}t\left({1 + 0.1x + 0.1y} \right)\sin \left({\pi x} \right)\sin \left({\pi y} \right) + \sin \left[{t\sin \left({\pi x} \right)\sin \left({\pi y} \right)} \right]$, and the analytical solution is $u\left({x, y, t} \right) = t\sin \left({\pi x} \right)\sin \left({\pi y} \right)$.

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 ${\tau = h/2}$ and ${\tau = 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-Ⅰ and HOC-Ⅱ, and both have fourth-order convergence. In addition, the HOC-Ⅱ scheme requires about half the CPU time of HOC-Ⅰ. This demonstrates the significance of the HOC-Ⅱ scheme described in this paper.

Figure 2 shows the surfaces and the contour plots of the computational results using the HOC-Ⅰ scheme when $\tau = 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.

4.3.   Coupled sine–Gordon equations ^[6,9,10,35]

Consider the coupled sine–Gordon equations (3.1)–(3.6) in a square domain $\left[{- 2, 2} \right] \times \left[{- 2, 2} \right]$ and $t \in \left({0, T} \right]$. The set of analytical solutions is

in which

To compare with the calculation results in Ref. ^[6], we take the same parameters $\alpha = \sqrt {3}$, $\beta = 0.5$, $\eta = 1.5$. The initial conditions and boundary conditions are given by the analytical solutions.

We solve this problem using the difference schemes HOC-Ⅰ and HOC-Ⅱ at $T = 1$. For comparison, we also derive solutions using ADI-Ⅰ and ADI-Ⅱ ^[6]. Table 7 presents the numerical results calculated by HOC-Ⅰ and ADI-Ⅰ ^[6]. The results using HOC-Ⅰ are smaller than the values in the literature for both the $L_{\infty}$ norm error and the $L_{2}$ norm error. Additionally, our schemes achieve fourth-order convergence, which is consistent with the theoretical accuracy. ADI-Ⅰ 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-Ⅰ 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-Ⅱ ^[6], the calculation results of HOC-Ⅱ are better. We can also see that HOC-Ⅱ is more computationally efficient than HOC-Ⅰ, suggesting the superiority of the linearized HOC-Ⅱ scheme.

Figure 3 shows the surfaces of the numerical solutions $u_{i, j}^n$, $w_{i, j}^n$ 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.

4.4.   Evolution of single ring soliton ^[6,25,26]

Consider the sine–Gordon equation in the square field $\Omega = \left[{- 15, 15} \right] \times \left[{- 15, 15} \right]$:

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-Ⅱ.

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-Ⅱ in solving such problems.

4.5.   Evolution of double solitons in a layered medium ^[22]

Finally, in a layered medium covering $[-8, 8] \times [-8, 8]$, we research the following Klein–Gordon equation:

in which

We consider periodic boundary conditions and initial conditions of

The HOC-Ⅱ 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.

5.   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-Ⅰ, and the other is a linearized fourth-order compact scheme, denoted as HOC-Ⅱ. In the linear HOC-Ⅱ 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-Ⅱ is significantly better than HOC-Ⅰ and several previously developed difference schemes in terms of computational efficiency. Additionally, HOC-Ⅱ 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 three-dimensional 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.

Acknowledgments

We would like to thank the editors and the referees whose constructive comments and suggestions are helpful to improve the quality of this paper. This work is partially supported by National Natural Science Foundation of China (12161067), Natural Science Foundation of Ningxia (2022AAC02023, 2020AAC03059), National Youth Top-notch Talent Support Program of Ningxia, and the First Class Discipline Construction Project in Ningxia Universities: Mathematics.

Conflicts of nterest

The authors declare no conflict of interest.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Appendix. Linear stability analysis

For simplicity, the stability of the difference scheme for the linear problem is analyzed by using Fourier method. Here, we consider x- and y-directions taking the same step size, with the corresponding difference scheme as follows.

Lemma 2. ^[41] The sufficient and necessary condition for the roots of the quadratic equation ${\mu ^2} - b\mu - c = 0$ with real coefficients to be less than or equal to 1 is $\left| c \right| \le 1, \left| b \right| \le 1 - c$.

Theorem 1. Assume that there is no error in the boundary and the source term, then the stability condition of difference scheme (1) is

where ${v_{\max }} = \mathop {\max }\limits_{ 0 \le i \le {M_x} \atop 0 \le j \le {M_y} } \;\;\;\;\left| {{v_{i, j}}} \right|$ and $\lambda = \frac{\tau }{h}$.

Proof. Let $U_{i, j}^n = {\zeta ^n}{e^{I({\sigma _1}{x_i} + {\sigma _2}{y_j})}}$, $\left({{U_{xx}}} \right)_{i, j}^n = {\mu ^n}{e^{I({\sigma _1}{x_i} + {\sigma _2}{y_j})}}$ and $\left({{U_{yy}}} \right)_{i, j}^n = {\kappa ^n}{e^{I({\sigma _1}{x_i} + {\sigma _2}{y_j})}}$, where $\zeta$, $\mu$ and $\kappa$ are amplitudes, $\sigma_1$ and $\sigma_2$ are wavenumber and $I = \sqrt { - 1}$. According to (2.12) and (2.13) we have

Let $W_{i, j}^{n + 1} = U_{i, j}^n$, $\alpha = \mathop {\max }\limits_{ 0 \le i \le {M_x}\atop 0 \le j \le {M_y}} {\vartheta _{i, j}} = \mathop {\max }\limits_{ 0 \le i \le {M_x}\atop 0 \le j \le {M_y}} v_{i, j}^2$ and assume that there is no error in the source term, (1) can be written in the following matrix form

Then, we can obtain

where $G = \left[{\begin{array}{*{20}{c}} A & { - 1}\\ 1 & 0 \end{array}} \right]$ is the error propagation matrix, and $A = 2BC{\alpha ^2}{\lambda ^4} + 12B\alpha {\lambda ^2} + 2$, $B = \frac{{\left[{\cos \left({{\sigma _1}h} \right) - 1} \right]}}{{\left[{\cos \left({{\sigma _1}h} \right) + 5} \right]}} + \frac{{\left[{\cos \left({{\sigma _2}h} \right) - 1} \right]}}{{\left[{\cos \left({{\sigma _2}h} \right) + 5} \right]}}$, $C = \left[{\cos \left({{\sigma _1}h} \right) + \cos \left({{\sigma _2}h} \right) - 2} \right]$. The characteristic equation corresponding to the error propagation matrix is

According to Lemma 2, the stability condition of the difference scheme is $\left| {BC{\alpha ^2}{\lambda ^4} + 6B\alpha {\lambda ^2} + 1} \right| \le 1$, i.e., $- 2 \le BC{\alpha ^2}{\lambda ^4} + 6B\alpha {\lambda ^2} \le 0$.

First, consider the inequality

It can be easily obtained

Next consider the inequality

Likewise, inequality (7) holds true when $\cos \left({{\sigma _1}h} \right) = \cos \left({{\sigma _2}h} \right) = 1$. Otherwise, let quadratic function $F\left(\omega \right) = {\omega ^2}B C + 6 B \omega + 2$. Since $B C > 0$, the two roots of quadratic function are

Further analysis shows that the smallest root of the quadratic equation is 1/2. Hence, the solution satisfying the inequality (7) is

Combining (6) and (8), the stability condition of the difference scheme (1) is $\left[{0, \frac{{\sqrt 2 }}{2}} \right]$. This completes the proof.