Hyperbolic B-Spline Function-Based Differential Quadrature Method for the Approximation of 3D Wave Equations

: We propose a differential quadrature method (DQM) based on cubic hyperbolic B-spline basis functions for computing 3D wave equations. This method converts the problem into a system of ODEs. We use an optimum ﬁve-stage and order four SSP Runge-Kutta (SSPRK-(5,4)) scheme to solve the obtained system of ODEs. The matrix stability analysis is also investigated. The accuracy and efﬁciency of the proposed method are demonstrated via three numerical examples. It has been found that the proposed method gives more accurate results than the existing methods. The main purpose of this work is to present an accurate, economically easy-to-implement, and stable technique for solving hyperbolic partial differential equations.


Introduction
Some of the existing numerical approaches for solving wave equations as well as fractional wave equations involve the temporal extrapolation method, finite difference method (FDM), finite volume method (FVM), finite element method (FEM), and boundary element method (BEM) [1][2][3][4][5]. It has been noticed that the most popularly used numerical approaches for solving 3D wave equations are based on the FDM [6,7]. In short, the FDM is utilized to handle the time derivative, and the space derivatives are discretized by other numerical techniques. In particular, the radial basis and B-spline basis functions-based collocation methods are extensively applied for solving 3D wave equations. Ranocha et al. [8] set up fully discrete conservative techniques for various disseminative wave equations. Recently, Wang et al. [9] presented radial basis function-based single-step mesh free technique for 2D variable coefficients wave equation. Bakushinsky and Leonov [10] presented a fast Fourier transform-based algorithm to solve the 3D wave inverse problem in a cylindrical system.
In recent decades, wave equations have been approximated by numerous researchers. Dehghan [11] approximated the solution of 1D hyperbolic PDEs with nonlocal boundary specifications, while in [12], the author used ADI, fully implicit, fully explicit FD methods, and the Barakat and Clark type explicit formulae to approximate the 2D Schrodinger equation. Mohanty and Gopal [13] presented an off-step discretization-based technique for the approximation of 3D wave equations. Titarev and Toro [14] implemented fourthorder accurate ADER schemes for 3D hyperbolic systems. Zhang et al. [15] proposed an improved element-free Galerkin (EFG) method, while EFG method and Meshless local Petrov-Galerkin (MLPG) method have been proposed by Shivanian [16]. Recently, Shukla et al. [17] proposed an Expo-MCBDQM to approximate the aforementioned equations.
Bellman et al. [18] introduced DQM. DQM based on various basis functions has been presented to solve several PDEs such as sinc DQM [19], Fourier expansion based DQM [20], harmonic DQM [21], quintic B-spline DQM [22] and many more. The authors of [23][24][25][26] proposed a cubic B-spline (CB-spline) based DQM for the Burgers', coupled with and with boundary conditions where Ω = {(x, y, z) : x, y, z ∈ [0, 1]} is the problem domain and ∂Ω is its boundary. The functions f , ξ 1 , ξ 2 and ζ are known whereas u = u(x, y, z, t) is to be determined. The terms α 1 , α 2 , δ are real constants. The rest of the paper is organized as follows. In Section 2, the procedure of the cubic hyperbolic B-spline DQM is described. Section 3 examines the stability analysis of the proposed method. Section 4 demonstrates the computational results. Finally, Section 5 presents the conclusion of our study.
In DQM, we approximate u xx , u yy and u zz as given below: where u ijk = u x i , y j , z k , t , and a (p) kr are the weighting coefficients corresponding to The cubic hyperbolic B-spline functions are given as [39]: where The cubic hyperbolic B-spline functions {H f 0 , H f 1 , . . . , H f M , H f M+1 } form a basis over Ω. Table 1 shows the values of cubic hyperbolic B-spline functions with derivatives at the knots, where Preserving the matrix system remains diagonally dominant, we modify the hyperbolic B-spline functions as:Ĥ which can be written in the form of tridiagonal system as follows: where T is the unknown vector and We solve the system (10) to find the weighting coefficient vector a (1) ij . Similarly, by fixing x and z in ∂u ijk ∂y and using the modified cubic hyperbolic B-spline functionsĤ f r (y), r = 1, 2, . . . , M x in Equation (5), and fixing x and y in ∂u ijk ∂z and using the modified cubic hyperbolic B-spline functionsĤ f r (z), r = 1, 2, . . . , M x in Equation (6), we can compute b ij , for i, j = 1, 2, . . . , M z . (13) Axioms 2022, 11, 597 5 of 12 Now, using u t = w, u tt = w t and substituting the approximated u xx , u yy and u zz by the cubic hyperbolic B-spline DQM in problem (1) and (2), we have and where and where, Finally, we apply SSPRK-(5,4) scheme [41] to solve the above systems.

Stability Analysis
For stability, we rewrite Equations (16) and (17) by choosing α 1 , α 2 ≥ 0 and α 1 > δ f as follows: where, , N and I are null and identity matrices, respectively. We have ij and c (2) ij , respectively, and given as follows: Axioms 2022, 11, 597 6 of 12 The order of null matrices N y and N z are (M y − 2)(M z − 2) and (M z − 2), respectively which is same as the order of identity matrices I x and I z .
Now, we suppose that λ A be an eigenvalue of A associated with the eigenvector (X 1 , X 2 ) tr , where the order of each component vector is which implies that This illustrates that λ A (λ A + α 1 − δ f ) is the eigenvalue of F as follows: The eigenvalues of F 1 with h x = h y = h z = 0.2, 0.1 and 0.05 are represented in Figure 1, where real and negative eigenvalues are observed. Equation (25) implies that all eigenvalues of B are real and negative.
From Equation (26), we conclude that x < , we conclude that th A λ will be negative. Therefore, the proposed method is stable for the 3D w discretized system.

Computational Results
Now, we consider three examples of the 3D wave equation (1) to check and efficiency of the proposed method. The ˆ( , , , ) p x y z t is described appropriately. We choose 0.1, 0.05 h = a Table 2 shows the comparison between the proposed method and the ex terms of RMS error norms. It can be noted that the present solutions are m than the solutions presented in [16] by EFG and MLPG methods, and Ex [17]. Figure 2 illustrates the absolute error norms for fixed z = 0.5 with h while Figure 3 shows the behavior of the solutions. From Figures 2 and 3, o that the absolute error norms are very small, and analytical and numerica

Now, let λ
y is negative and real, that is, From Equation (26), we conclude that Since α 1 > δ f and so x < 0, we conclude that the real part of λ A will be negative. Therefore, the proposed method is stable for the 3D wave equations discretized system.

Computational Results
Now, we consider three examples of the 3D wave Equation (1) to check the accuracy and efficiency of the proposed method. Example 1. We consider Equation (1) for δ = 0 and α i = 2, i = 1, 2 with the analytical solution u(x, y, z, t) = e −2t sinhx sinhy sinhz.
Thep(x, y, z, t) is described appropriately. We choose h = 0.1, 0.05 and ∆t = 0.01. Table 2 shows the comparison between the proposed method and the existing ones in terms of RMS error norms. It can be noted that the present solutions are more accurate than the solutions presented in [16] by EFG and MLPG methods, and Expo-MCBDQM [17]. Figure 2 illustrates the absolute error norms for fixed z = 0.5 with h = 0.05 at t = 1 while Figure 3 shows the behavior of the solutions. From Figures 2 and 3, one can notice that the absolute error norms are very small, and analytical and numerical solutions are very close each other which shows the accuracy of the proposed method.     The function ˆ( , , , ) p x y z t is described appropriately and we choose the parameters  Example 2. Next, we consider Equation (1) for α 1 = δ = κ, α 2 = 0 and f (u) = u 2 with the analytical solution u(x, y, z, t) = sin πx sin πy sin πz e −κt .
The functionp(x, y, z, t) is described appropriately and we choose the parameters h = 0.05, 0.1, ∆t = 0.01, and κ = 3 for the numerical approximation of this example. The proposed method is again compared with EFG [16] and MLPG [16] methods, and Expo-MCBDQM [17] in terms of RMS error norms and shown in Table 3. It is noticed that the proposed method shows better solutions than the solutions presented in [16,17]. Figure 4 illustrates the absolute error norms for z = 1, h = 0.05 at t = 1 while Figure 5 demonstrates the comparison of the analytical and numerical solutions, where a close agreement is noticed between analytical and numerical solutions. From Figure 4, one can notice that the absolute error norms are very small in 10 −19 , which shows that the proposed method provides very accurate results.
The functionp(x, y, z, t) is chosen appropriately. We choose the parameters h = 0.05 and 0.1, ∆t = 0.01. The present solutions are compared with the solutions obtained by EFG [16] and MLPG [16] methods, and Expo-MCBDQM [17] in terms of RMS error norms and shown in Table 4. Again, it is noticed that the proposed method provides better solutions than the existing methods. Figure 6 illustrates the absolute error norms for z = 0.5, h = 0.05 at t = 1. Figure 7 shows the comparison of the analytical and numerical solutions. MCBDQM [17] in terms of RMS error norms and shown in Table 3. It is noticed that the proposed method shows better solutions than the solutions presented in [16,17]. Figure 4 illustrates the absolute error norms for z = 1, h = 0.05 at t = 1 while Figure 5 demonstrates the comparison of the analytical and numerical solutions, where a close agreement is noticed between analytical and numerical solutions. From Figure 4, one can notice that the absolute error norms are very small in 10 −19 , which shows that the proposed method provides very accurate results.     The function ˆ( , , , ) p x y z t is chosen appropriately. We choose the parameters h = 0.05 and 0.1, t Δ = 0.01. The present solutions are compared with the solutions obtained by EFG [16] and MLPG [16] methods, and Expo-MCBDQM [17] in terms of RMS error norms and shown in Table 4. Again, it is noticed that the proposed method provides better solutions than the existing methods. Figure 6 illustrates the absolute error norms for Figure 7 shows the comparison of the analytical and numerical solutions.

Computational Complexity
The system of equations, where the coefficient matrix is tridiagonal, is solved by using the Thomas algorithm which takes 3 M subtractions, 3 M multiplications, and 2 M + 1 divisions. Therefore, the algorithm needs a total of 8 M + 1 simple arithmetic operations. Therefore, the Thomas algorithm requires ( ) O n operations. The computational cost of the SSPRK-(5,4) technique is same as the cost of traditional ODE solvers. Hence, the proposed technique is not too complex from the computational point of view.

Conclusions
This work proposed a differential quadrature method based on cubic hyperbolic Bspline functions together with SSPRK-(5,4) scheme to solve 3D wave equations. The numerical examples show that the proposed method provides more accurate solutions than those discussed in [16,17]. The matrix stability analysis is also investigated, and we found that the proposed method is stable. Additionally, the method is economically easy-to-implement for solving hyperbolic partial differential equations. Moreover, the computational complexity shows that the proposed technique is not too complex from the computational point of view.

Computational Complexity
The system of equations, where the coefficient matrix is tridiagonal, is solved by using the Thomas algorithm which takes 3M subtractions, 3M multiplications, and 2M + 1 divisions. Therefore, the algorithm needs a total of 8M + 1 simple arithmetic operations. Therefore, the Thomas algorithm requires O(n) operations. The computational cost of the SSPRK-(5,4) technique is same as the cost of traditional ODE solvers. Hence, the proposed technique is not too complex from the computational point of view.

Conclusions
This work proposed a differential quadrature method based on cubic hyperbolic B-spline functions together with SSPRK-(5,4) scheme to solve 3D wave equations. The numerical examples show that the proposed method provides more accurate solutions than those discussed in [16,17]. The matrix stability analysis is also investigated, and we found that the proposed method is stable. Additionally, the method is economically easy-to-implement for solving hyperbolic partial differential equations. Moreover, the computational complexity shows that the proposed technique is not too complex from the computational point of view.