A Study of Traveling Wave Structures and Numerical Investigation of Two-Dimensional Riemann Problems with Their Stability and Accuracy

The Riemann wave system has a fundamental role in describing waves in various nonlinear natural phenomena, for instance, tsunamis in the oceans. This paper focuses on executing the generalized exponential rational function approach and some numerical methods to obtain a distinct range of traveling wave structures and numerical results of the two-dimensional Riemann problems. The stability of obtained traveling wave solutions is analyzed by satisfying the constraint conditions of the Hamiltonian system. Numerical simulations are investigated via the finite difference method to verify the accuracy of the obtained results. To extract the approximation solutions to the underlying problem, some ODE solvers in FORTRAN software are applied, and outcomes are shown graphically. The stability and accuracy of the numerical schemes using Fourier’s stability method and error analysis, respectively, to increase the reassurance are investigated. A comparison between the analytical and numerical results is obtained and graphically provided. The proposed methods are effective and practical to be applied for solving more partial differential equations (PDEs).


Introduction
Diverse phenomena in nature and technology are described using nonlinear evolution equations (NLEEs). In physics, for instance, the heat transfer and the traveling wave phenomena are successfully modeled by PDEs. In chemistry, the dispersion of a chemically reactive substance is controlled by PDEs. Also, PDEs are invoked to characterize population growth problems. Moreover, most physical incidents of shallow-water waves, quantum mechanics, plasma physics, electricity, fluid dynamics, and more others are investigated using PDEs. In order to better comprehend the qualitative characteristics of such equations, one can study their analytical solutions. A distinct range of traveling wave structures allows us to explain the mechanisms of immense complicated phenomena. As a result, some researchers have revealed numerous effective methods. Some of the powerful approaches are the Sine-Gordon expansion approach [1], the modified simple equation technique [2], the tanh-sech process [3,4], the trial equation method [5], the exp(−f (ζ ))-expansion principal [6,7] and the generalized exponential rational function approach [8]. More methods can be easily seen in references [9][10][11][12][13][14][15]. The Riemann wave system [16,17] is given by where α, β and γ are constants. System (1) describes a (2 + 1)-dimensional interaction for the propagation of Riemann waves in the x-axis and y-axis. The Riemann wave equations have been investigated by some scientists using various methods. For example, Shao [18] examined the stability of the solution of the quasi-linear hyperbolic systems under the influence of some BV perturbations. Shao [19] presented several solutions with satisfying some conditions. Chen et al. [20] used the generalized expansion approach of the Riccati equation to extract some soliton-like solutions for Eq. (1) when m = n = 4b, where b is a constant. The modified mapping principle was applied in [21] to establish a class of periodic wave solutions in terms of the Jacobi elliptic functions. Moreover, Abdelrahman et al. [22] utilized the singular manifold technique to construct two general solutions for system (1). Every obtained solution contains two arbitrary functions. Then, some periodic wave solutions were derived from the obtained general solutions. In [16], the tanh method was implemented to generate some traveling wave solutions for Eq. (1). The generalized Kudryashov approach was applied in [17] to systematically and graphically show the traveling wave solutions of Eq. (1). While, in this paper, we present analytical and numerical solutions to Eq. (1) using the generalized exponential rational function approach [8] and finite difference methods, respectively.
Since our knowledge of constructing exact solutions for the Riemann wave equations is basically based on few techniques, we utilize the generalized exponential rational function approach [8] to construct the exact solutions of the considered equations. This technique depends on Jacobi elliptic functions. As a result, various solitary traveling wave solutions can be simply generated in terms of trigonometric and hyperbolic functions. Regarding the numerical solution, some numerical simulations are presented using the finite difference method to emphasize that the obtained results are accurate. The numerical solutions of the considered system are obtained by approaching the equations on a mesh using the finite-difference notations. The domain is divided into a limited set of grids to achieve meshes for both independent variables x and y. It is well known that the wave of the solutions has areas with rapid spatial changes, for instance, steep fronts structures. In order to resolve these types of areas, fine numbers of grids, for x and y, are required. The step size, x = x i − x i−1 , y = y i − y i−1 should be extremely small to catch these regions. Computationally, it is intensive and expensive. Hence, I strive to achieve an alternative discretizing for the meshes to have a non-uniform mesh that manually sets more grids in where the solution varies rapidly and fewer grids outside these regions. And then, I used the stiff ODE solver DASPK [23] to solve the obtained ODEs of the semi-discretization of the system. This ODEs solver is an implicit iterative method, based on the Krylov subspace method. They are applied to solve the system of linearized equations. They also allow converting the Jacobian matrix to be LU factorized to make the calculations faster. I additionally studied, here, the stability and the error analysis for the numerical schemes.

Methodology
The generalized exponential rational function approach is comprehensively summarized in this section, as presented in [8]. Assume that is a given nonlinear PDE in the unknown functions V = V (x, y, t), U = U(x, y, t) and = (x, y, t). P is a polynomial in V , U and and their partial derivatives. Suppose that ζ = μx + δy − wt, where μ, δ and w are unknown parameters that will be computed later. Thus, Then, Eq. (2) is converted into where = d dζ . According to the proposed method, the formal solution of Eq. (4) is written as where p 1 p 2 , p 5 , p 6 , q 1 , q 2 , q 3 , and q 4 are complex (or real) constants so that Eq. (2) is expressed as The balance principle is used to determine the value of N appearing in Eq. (6). Moreover, the coefficients A 0 , A k and B k (k = 1, 2, . . . , N) is evaluated such that Eq. (6) satisfies Eq. (4). Inserting Eq. (6) into Eq. (4) gives a polynomial from which one can obtain an algebraic system solved using Mathematica or Maple. The values of the above-mentioned constants are included in the solutions of this system. Substituting these values into Eq. (6) yields the exact solutions of Eq. (2).

Traveling Wave Solution
In this section, we study the exact solutions of system (1) using the generalized exponential rational function approach. Substitute Eq. (3) into system (1) to have −wv ζ + μδ 2 αψ ζ ζ ζ + βδψv ζ + γ μuv ζ = 0, We now integrate each equation in system (7) with respect to ζ to obtain Balancing the highest order of v ζ ζ with the nonlinear term v 2 in system (8) evaluates the value of N given by N = 2. As a result, the exact solutions are shown as follows: where ϒ(ζ ) is presented in Eq. (5). We insert systems (9) into (8) to obtain the traveling wave solutions of (1) which are written as follows.
Family 1: Substituting Eq. (10) into the first equation of system (9) and then we insert the result into the first equation of system (8) give an algebraic system whose solutions are given as follows: Substituting Eq. (11) into system (9) and Eq. (10) gives the exact solutions of system (1) which are illustrated as follows: Case 2: Inserting Eq. (13) into system (9) and Eq. (10) yields the exact solutions of system (1) which are expressed as Case 3: Plugging Eq. (15) into system (9) and Eq. (10) gives the traveling wave solutions of system (1) which are given by Case 4: Putting Eq. (11) into system (9) and Eq. (10) leads to the exact solutions of system (1) which are shown as Case 1: Substituting Eq. (20) into system (9) and Eq. (19) gives the exact solutions of system (1) which are shows as follows: Case 2: Inserting Eq. (22) into system (9) and Eq. (19) yields the traveling wave solutions of system (1) which are Case 3: Putting Eq. (24) into system (9) and Eq. (19) gives the exact solutions of system (1) which are Case 4: Substituting Eq. (26) into system (9) and Eq. (19) leads to the exact solutions of system (1) which are expressed as

Stability of the Analytical Solution
We introduce the Hamiltonian system in this section. Hamiltonian system is applied on analytical solutions to test their stability on a specific interval. The Hamiltonian system [24,25] is expressed by where indicates the momentum function. Furthermore, w presents the wave speed and v(ζ ) is the considered analytical solution. The sufficient condition for the stability is When we apply Eqs. (28) and (29) where the parameters A 2 = 1.2, β = −0.5, α = 2.70, and γ = −2.20. As a result, the analytical solutions are unconditionally stable.

Finite Difference Semi-Discretization Scheme on a Fixed Mesh
This section is devoted to study the numerical solutions of system (1)  Here, x and y denote the width of the sub-intervals in x and y directions, respectively. System (1) is now converted into some equations of ODEs by implementing finite differences on spatial derivatives. We keep the temporal derivative continuous. Completing this gives where δ 2 y | k+1 n,m = k+1 n,m+1 − 2 k+1 n,m + k+1 n,m−1 .
We compute the boundary conditions as follows: These boundary conditions allow us to employ the fictitious points in computing the spatial derivatives at the boundaries of the domain. It is worth noting that the initial conditions are established by evaluating the exact solution at t = 0, as follows: where In order to extract the numerical solutions of the considered equation, we implement the finite difference approach which depends on a standard ODE solver in FORTRAN software, DASPK solver [25]. The standard backward differentiation operators are invoked to approximate the time derivatives in this solver. In addition, the Jacobian matrix of the linearised system is approximated by applying LU factorization. To have less bandwidth for this matrix, we use a unique system numbering for the unknowns V 1,1 , V 2,1 , V 3,1 , . . . , V N+1,M+1 .

Stability of the Numerical Solution
This section investigates the stability of the numerical solution using a Fourier's stability technique. From the second and third equations of system (8), we have Substituting these equations into the first equation of system (1) yields where α 0 = μ δ α, β 0 = μ δ β and γ 0 = δ μ γ . Since the Fourier stability approach is applied on linear equations, we linearise Eq. (34) as follows: where L 0 and L 1 are constants quantity defined by The boundary conditions are ignored. Consider the point (x n , y m , t k ), where x n = n x, y m = m y and t k = k t. Let V k n,m = λ k e i ξ π n x e i r π m y and then V k+1 n,m = λV k n,m , n = 1, 2, . . . , N, and m = 1, 2, . . . , M. (36) plugging Eq. (36) into scheme (35) yields Dividing both sides of the result by V k n,m gives Assume that a = sin(rπ y) t y 2α Then, Eq. (38) becomes Hence, According to the Fourier stability, the stability of the considered scheme occurs if the absolute value of λ does not exceed one. This constrained condition is perfectly satisfied in our analysis. It is clear from Eq. (39) that the absolute value of λ is less than one. Consequently, the numerical scheme is unconditionally stable.

Error Analysis
Taylor series is used in this section to examine the order of the accuracy of scheme (31). We evaluate the truncation error to obtain the order. Suppose that where e k+1 n,m denotes the error, V k+1 n,m and V (x n , y m , t k+1 ) represent the approximation solution and the analytical solution at (x n , y m , t k+1 ), respectively. We now insert Eq. where T k n,m is the truncation error which is expressed as Hence, The leading terms mentioned above are known as the fundamental part of the local truncation error, and we have accepted the truth that u(x, y, t) is the solution of the underlying system. Therefore, The truncation error, which is generated in every step, is given by O( t, x 2 , y 2 ).

Convergence of the Numerical Schemes
Now consider that a sequence of computations is carried out using given initial data, with the refinement of three meshes, so that x → 0, y → 0 and t → 0. Then, the numerical scheme is said to be convergent if, for each fixed point ( We have established above that the implicit schema is unconditional stability. So, we will show that the implicit schemes are unconditional convergence. Suppose that the error e is given by Now, V k n,m satisfies the scheme (Eq. (31)) exactly, while V (x n , y m , t k ) omits the error indicated by the truncation error tT k n,m . Hence by using the fact that V y = U x , and V x = y and subtraction, we have where a 0 = α t 2 x 2 y , a 1 = β t 2 y , and a 2 = γ t 2 x and T k n,m is truncation error (see Eq. (42)). If we suppose that the maximum error for time step is given by Since the given initial data is used we can identify E 0 = 0. Hence the inequality is given by But T k n,m → 0 as x, y, t → 0, then Hence, the scheme (Eq. (31)) is convergent as x, y, t → 0.

Results and Discussion
In this section, we discuss the results shown in this work. We extract a distinct range of traveling wave structures of the two-dimensional Riemann problems via the generalized exponential rational function method. The obtained solutions are presented in terms of trigonometric and hyperbolic functions. We examine the stability of Eq. (18)  The numerical solutions are discussed by employing the finite difference method to convert the underlying problems into the system of ODEs with keeping time derivatives continuous. Then, I solve the resulted ODEs system using the DASPK solver. This method gives reliable and powerful results. This can be clearly seen in the graphical comparisons presented in the above-mentioned figures. For instance, Figs. 1 and 2 illustrate the behavior of the analytical and numerical solutions for t = 5 and t = 10, respectively. It can be easily observed from these figures that the solutions nearly have the same behavior. In Fig. 3, the exact and numerical solutions approximately have the same behavior when t = 20. Fig. 4 presents the time evolution of V(x, y, t) to the traveling wave structures (a) the analytical and (b) the numerical solutions with parameter values A 2 = 1.2, β = −0.5, α = 2.70, and γ = −2.20. Fig. 4b also illustrates the performance of the used approaches. Moreover, Fig. 5 shows acceptable performance for the used numerical technique when a massive number of meshes is used. For example, when we use N = 120, the error is high. Nevertheless, the numerical solutions approximately approach the exact solution (blue line) for N = 3000. The stability of the numerical results is investigated using the Fourier technique. Since |λ| ≤ 1 in Eq. (40), the numerical scheme is unconditionally stable. Furthermore, the accuracy of the numerical scheme is of O( t, x 2 , y 2 ).     Table 1 illustrates L 2 errors and the CPU times taken to reach t = 20. The error decreases for large N but the method takes more time to give a small error. We begin with N = 120, x = 0.50 and y = 0.01. The L 2 error stands at 5.37×10 −4 during 0.74 min. When we use N = 3000 with x = 0.02 and y = 0.01, the L 2 error dramatically decays and stands at 1.073 × 10 −6 during 20.63 min. Fig. 6 shows the decay in the L 2 error as N increases.  Table 1 10 Conclusion We have favorably implemented an accurate finite difference method on a uniform mesh and the generalized exponential rational function method for the two-dimensional Riemann problems. The main advantage of the results is to show the traveling wave structures and prove their accuracy using numerical methods. By comparing the exact solutions using the generalized exponential rational function method with those from the numerical scheme, it was said that the numerical results are almost identical to the analytical results. It is well known that the numerical scheme is stable and allows a meaningful reduction in memory requirements. Using a fine mesh in both spatial variables x and y permits us to resolve the wave-like structures. We can conclude that the used methods can be efficiently applied to more nonlinear evolution models to construct their exact and numerical solutions.
Funding Statement: The author received no specific funding for this study.

Conflicts of Interest:
The author declares that they have no conflicts of interest to report regarding the present study.