Numerical Simulation of PDEs by Local Meshless Differential Quadrature Collocation Method

Imtiaz Ahmad 1 , Muhammad Ahsan 1,2, Iltaf Hussain 2, Poom Kumam 3,4,* and Wiyada Kumam 5,* 1 Department of Mathematics, University of Swabi, Swabi 23430, Pakistan; imtiazkakakhil@gmail.com (I.A.); ahsankog@uoswabi.edu.pk (M.A.) 2 Department of Basic Sciences, University of Engineering and Technology, Peshawar 25000, Pakistan; iltafhussain@uetpeshawar.edu.pk 3 KMUTTFixed Point Research Laboratory, Room SCL 802 Fixed Point Laboratory, Department of Mathematics, Science Laboratory Building, Faculty of Science, King Mongkut’s University of Technology Thonburi (KMUTT), 126 Pracha-Uthit Road, Bang Mod, Thrung Khru, Bangkok 10140, Thailand 4 Department of Medical Research, China Medical University Hospital, China Medical University, Taichung 40402, Taiwan 5 Program in Applied Statistics, Department of Mathematics and Computer Science, Faculty of Science and Technology, Rajamangala University of Technology Thanyaburi (RMUTT), Thanyaburi, Pathumthani 12110, Thailand * Correspondences: poom.kum@kmutt.ac.th (P.K.); wiyada.kum@rmutt.ac.th (W.K.)


Introduction
The Klein-Gordon (KG) equation can describe various vital phenomena in chemical and physical sciences.The KG equation can be written as: along with initial and boundary conditions: where α, β, and γ are the parameters.
The nonlinear KG equation has applications in numerous fields such as nonlinear optics, quantum mechanics, solid state physics, and in mathematical physics [1,2].Various numerical techniques are employed for solving the KG equation such as the symplectic finite difference method [3], the spectral method [4], the meshless RBF method [5], the finite-difference collocation method [6], the differential quadrature method [7], the Haar wavelet method [8], the decomposition method [2], the Laplace transformation and Legendre wavelets method [9], the lattice Boltzmann method [10], and the multiquadric quasi-interpolation method [11].
The coupled Burgers' equations are related to many physical problems including traffic flow, acoustic transmission, flow of a shock wave traveling in a viscous fluid, airfoil flow theory, supersonic flow, and in turbulence phenomena (see [12][13][14][15] and the references therein for details).
Recently, meshless methods have seen broad attention for solving different kinds of PDE model areas in almost all disciplines of engineering.The meshless character is one of the most important reasons for the rising demand of such types of methods.Meshless methods reduce the complexity caused due to dimensionality to a large extent, which is faced in the carrying out of conventional methods like the finite-element and finite-difference procedures.Meshing in the case of complicated geometries is another cause for the growing demand of meshless methods.The numerical results of RBF-based algorithms have demonstrated that they are truly meshless, accurate, and easy to implement.Some interesting models can be found in [26][27][28][29][30].
It is noted that the global meshless method (GMM), which is based on the global interpolation paradigm, has faced the problems of dense ill-conditioned matrices and finding the optimum value of the shape parameter.To avoid the limitations of the GMM, a local meshless method, which is based on local interpolation in the sub-domains, is used as a substitute to get a stable and accurate solution for the PDE models (see [31][32][33][34][35][36][37]).
In the current work, the local meshless differential quadrature collocation method (LMM) based on radial basis functions (RBFs) is proposed for the numerical simulation of 1D nonlinear KG, 2D coupled Burgers' equations, and the 2D RLW equation.The LMM is an accurate and efficient numerical technique, which requires two steps to approximate a time-dependent PDEs.Firstly, the spatial derivatives are approximated by using the RBFs, which could convert the given PDE into a system of ODEs.Then, the obtained ODEs will be solved by suitable ODE solvers.
The rest of the paper is organized as follows: the suggested numerical method is highlighted in Section 2; the implementation of the method with different test problems is presented in Section 3; and finally, some concluding remarks are given in Section 4.

Implementation of the Numerical Method
The LMM [26,33] is extended to the PDE models discussed in Section 1.The derivatives of U(x, t) at the center x i are approximated by the function values at a set of nodes in the neighborhood of {x i1 , x i2 , x i3 , . . ., N n , where i = 1, 2, . . ., N n .For the one-dimensional case, n = 1 and x = x, and for the two-dimensional case, n = 2 and x = (x, y).Now, for the 1D case, we have: To find the corresponding coefficient λ k , radial basis function ψ( x − x l ) can be substituted in Equation (7) as follows: where ψ( x ik − x l ) = 1 + (c x ik − x l ) 2 and ψ( x ik − x l ) = (1 + (c x ik − x l ) 2 ) −1 in the case of multiquadric (MQ) and inverse quadratic (IQ) radial basis functions, respectively.The matrix notation of Equation ( 8) is: where: From Equation (9), we obtain: From Equations ( 7) and (10), we get: where: For the 2D case, the derivatives of U(x, y, t) with respect to x are approximated in a similar way as stated above and can be written as: and the corresponding coefficients γ (m) k , k = 1, 2, . . ., n i can be found as follows: Similarly, the derivatives of U(x, y, t) with respect to y and its corresponding coefficients η k , k = 1, 2, . . ., n i can be found as follows: The proposed LMM is combined with a technique based on the local supported domain, called an upwind technique, in the case of convection-dominated PDE models.This technique has the ability to avoid spurious oscillatory solutions.Two types of local supported domains are used, i.e., central and upwind, as shown in Figures 1 and 2.

Implementation of LMM for the KG Equation
Now, using the above meshless method for Equation (1) in space, we get the second order ODE, which is reduced to the first order ODE by substituting U t (x, t) = V(x, t) as follows: with initial and boundary conditions: Now, applying the LMM to Equation ( 15): Equations ( 17) and ( 18) in matrix form: with the following corresponding initial condition:

Implementation of LMM for the 2D Model Equations
Applying the LMM to 2D model equations in space along with the prescribed boundary and initial conditions at each nodal point, we get the following form of an initial value problem: where A represents the sparse coefficient matrix of order N n × N n (n = 1 in the 1D case and n = 2 in the 2D case).The matrix A is obtained once U and its spatial derivatives are discretized by the LMM.The vector h denotes the boundary conditions of the problem, and the vector fis a vector of the corresponding initial condition of the problem.Orders of the vectors h and f are N n × 1, where n = 1; 2 for one-and two-dimensional PDEs, respectively.

Numerical Analysis
To test the accuracy and applicability of the proposed local meshless method, several problems have been considered.The proposed scheme works for any radial basis functions, but in this study, multiquadric and inverse quadratic RBFs were taken into account for spatial discretization.
The accuracy of the LMM was measured via the L abs , L 2 , and L ∞ error norms, which are given as follows: where Û and U represent the exact and approximate solutions, respectively.The following formulas were used to compute the numerical rate of convergence: where U ∆x i and U ∆t i represent the numerical solutions with spatial step size ∆x i and time step size ∆t i , respectively.The proposed LMM is truly meshless and capable of approximating the solution on both uniform and scattered nodal points.The size of the local sub-domain in the one-dimensional case was taken as three, whereas in the two-dimensional case, it was taken as five.In all numerical simulations, the multiquadric RBF with the value of the shape parameter c = 10 was used for the Klein-Gordon equation, and the inverse quadric RBF with shape parameter c = 20 was used for the coupled Burgers', as well as for regularized long wave equations.The forward Euler difference formula (FEDF) was used as the time integrator throughout the numerical simulation.The central processing unit (CPU) time was calculated in seconds in all cases.All the computations were performed using MATLAB (R2012a) on a Dell PC Laptop (Windows 7, 64 bit) with an Intel (R) Core(TM)i5-240M CPU 2.50 GHz 2.5 GHs4 GB RAM.
with initial conditions: Using the substitution U t = V, we get: The numerical results presented in Table 1 were obtained by the local meshless method at coarse grids by using nodal points N = 11, time step size ∆t = 0.0001, spatial domain [0, 1], and up to final time t = 1.The numerical results of the global meshless method of lines [38] versus the suggested local meshless method are shown in Table 1.It is clear from Table 1 that the LMM gives better accuracy compared to the numerical procedures reported in [38].In Table 2, the spatial convergence rate in terms of L ∞ and L 2 error norms for the MQ radial basis function and condition number κ are given for N = 11, 21, 31, 41, ∆t = 0.0001, t = 1.It can be seen from Table 2 that the condition number κ, as well as the convergence rate both increased with the increase in the number of collocation points N. The table also shows the second order of convergence rate of the LMM.Table 3 shows the time convergence rate in terms of the L ∞ and L 2 error norms for the time step sizes ∆t = 0.1, 0.05, 0.01, 0.005, N = 11, and t = 1.The results in Table 3 show that the FEDF has the first order convergence rate.The numerical results obtained by the local and global meshless method for a long range of shape parameter value c, with N = 11 and t = 1, are shown in Figure 4 for Test Problem 1. Less sensitivity to the selection of the shape parameter c in the case of the LMM, in comparison to the GMM, can be observed from Figure 4. Figure 5 illustrates the convergence of the proposed method, which shows that the error decreased (L ∞ and L 2 error norms) with the decrease of both the time step size ∆t, as well as the distance between nodes ∆x.Tables 2 and 3 show the convergence rate corresponding to Test Problem 1.A faster spatial convergence rate as compared to the time convergence rate can easily be seen from the tables.Test Problem 2. Consider the second test problem, by taking h(x, t) = sin(x) sin(t), α = −1, β = −2, and γ = 0 with the exact solution given in [38,39].
with initial and boundary conditions: Using the substitution U t = V, we get: The results of Test Problem 2, in terms of the L ∞ error norm at various times t, are shown in Table 4 for spatial domain [0, π/2], nodal points N = 11, and time step size ∆t = 0.0001.The results obtained by the LMM are more accurate than the results reported in [38,39].

Table 4.
Comparison of the numerical results of the global meshless method of lines [38] and the spline collocation method [39] with the LMM in terms of the L ∞ error norm for Test Problem 2.

Method
In Table 5, the condition number κ and spatial convergence rate are given for N = 11, 21, 31, 41, ∆t = 0.0001, and t = 1 for Test Problem 2. Table 5 shows that the increase in N caused the increase in both the condition number κ and convergence rate.One can observe from the table that the LMM had almost second order convergence.One of the drawbacks of meshless methods using shape parameter-dependent RBFs is the sensitivity of the shape parameter value.The comparison of both the local and global version of the meshless method for Test Problem 2 is shown in Figure 7.The results were obtained for MQ RBF using N = 41 and t = 1.It is clear from Figure 7 that the LMM gave stable results for a long range of shape parameter value c as compared to the GMM. Figure 8 illustrates the convergence rate of the proposed meshless method in which one can see that the error decreased with the decrease of both the time step size ∆t, as well as the distance between nodes ∆x.Tables 5 and 6 show the convergence rate corresponding to Test Problem 1.A faster spatial convergence rate as compared to the time convergence rate is verified in this case, as well.Test Problem 3. In the third test problem, consider Equation (1) with h(x, t) = −x cos(t) + x 2 cos 2 (t) and α = −1 β = 0, γ = 1 with the exact solution given in [38].
with initial conditions: Using the substitution U t = V, we get: Numerical results for different values of t, [−1, 1], N = 11, ∆t = 0.0001 are presented in Table 7, which shows the better performance of the LMM in comparison to the results reported in [5,10,11].Lattice Boltzmann method [10], N = 100 TPSRBFs method [5], N = 100 MQ quasi-interpolation scheme [11], N = 10 2.0694×10 −5 3.7065×10 −5 6.9684×10 −5 8.1943×10 −5 . . .In Figure 10, we compare the sensitivity of shape parameter value c for both meshless methods, local and global.It is clear from the figure that the LMM gave stable results in the range c ∈ (0, 100), but on the other hand, the GMM gave stable results only in the range c ∈ (0, 0.18).Figure 11 shows the condition number κ versus shape parameter c and number of nodal points N, which is self explanatory.Test Problem 4. The exact solutions for the 2D coupled Burgers' Equation (3) are written as [12,14]: In Table 8, the numerical results were obtained by the LMM in terms of the L ∞ error norm for Test Problem 4 using the IQ radial basis function with c = 20.We have used different N and Re and ∆t = 0.0001, t = 0.5.In this table, we have compared the numerical results of the proposed Local meshless method with the local method of approximate particular solutions (LMAPS) [15] and the local RBF-based DQmethod (LDQ) [15].A full agreement between the results of the LMM and the methods reported in [15] has been observed.
In Table 9, the performance of the LMM for Test Problem 4 is shown against the other methods [12-14] at selected points by letting Re = 100 and ∆t = 0.0001, t = 2, N = 20 × 20.These results show the accuracy and stable performance of the LMM.
Numerical solutions of the LMM corresponding to Test Problem 4 are shown in Figures 12-15 at different Reynolds numbers.It is clear from the figures that up to Re = 300, the proposed methods can handle the coupled Burgers' equation, but when we increased Reynolds number, the result became unstable.In contrast, the LMM combined with the upwind technique showed stable results at Re = 1000, and this phenomena can be seen in Figure 15.To show the computational efficiency of the LMM over the GMM, we have done a CPU time (in seconds) comparison, which is shown in Figure 16.Test Problem 5. Finally, consider the 2D nonlinear RLW Equation ( 5) with the following exact solution: U(x, y, t) = q 2 sech 2 √ q 2p (x + y − vt − x 0 − y 0 ) , (36) where q = 6(v−2) 2 and p = √ 6v.
The numerical results and CPU time for Test Problem 5 for v = 2.06, N = 50, ∆t = 0.5, [−80, 100], and up to time t = 20 are shown in Table 10.This table makes evident that the LMM is an accurate and efficient.In Figure 17, the L abs error is shown at t = 10.Eliminating the requirement of meshing and approximating the solution using uniform or scattered points in rectangular, as well non-rectangular domains is one of the salient features of the LMM.Numerical results of the RLW equation on non-rectangular domains are shown in Figures 18 and 19

Figure 1 .
Figure 1.Schematics of the local supported domain in 2D geometry for n i = 5 [34].

Figure 2 .
Figure 2. Schematic of the local supported domain in 2D geometry for n i = 5 in the upwind technique [34].

Figure 3 .
Figure 3. Numerical and exact solutions for t = 1 with N = 11 for Test Problem 1.

Figure 4 .
Figure 4. c versus the L ∞ error norm of the local meshless method (LMM) (left) and the global meshless method (GMM) (right) for Test Problem 1.

Figure 5 .
Figure 5. Convergence in space (left) and convergence in time (right) for Test Problem 1.

Figure 7 .
Figure 7. c versus the L ∞ error norm of the LMM (left) and GMM (right) for Test Problem 2.

Figure 8 .
Figure 8. Convergence in space (left) at t = 1 and convergence in time (right) at t = 0.1 for the test problem 2.

Figure 10 .
Figure 10.c versus the L ∞ error norm of the LMM (left) and GMM (right) for Test Problem 3.

Figure 11 .
Figure 11.c versus κ (left) and the number of nodal points N versus κ (right) of the LMM for Test Problem 3.

1 VFigure 14 .
Figure 14.Numerical solution of the 2D coupled Burgers' combined with the upwind technique for Re = 500, t = 1 for Test Problem 4.

1 VFigure 15 .
Figure 15.Numerical solution of the 2D coupled Burgers' combined with the upwind technique for Re = 1000, t = 1 for Test Problem 4.

Figure 16 .
Figure 16.Comparison of the CPU time of the LMM and the GMM with ∆t = 0.0001, t = 0.5 for Test Problem 4.

Figure 17 .
Figure 17.L abs error norm of the LMM for Test Problem 5.

Table 1 .
Numerical comparison of the L ∞ error norm for Test Problem 1.

Table 2 .
Spatial convergence rate of the L ∞ and L 2 error norms for Test Problem 1.

Table 3 .
Time convergence rate of the L ∞ and L 2 error norms for Test Problem 1.Comparisons of the exact and numerical solutions using the FEDF are plotted in Figure3.

Table 5 .
Spatial convergence rate of the L ∞ and L 2 error norms for Test Problem 2. Time convergence rate for Test Problem 2 is shown in Table 6 with time step sizes ∆t = 0.1, 0.05, 0.01, 0.005, N = 11, and t = 0.1.First order convergence is evident from the table.

Table 6 .
Time convergence rate of the L ∞ and L 2 error norms for Test Problem 2.

Table 7 .
Comparison of the L ∞ and L 2 error norms for Test Problem 3.

Table 8 .
Comparison of the numerical results of the 2D coupled Burgers' equation for Test Problem 4.

Table 9 .
Comparison of the 2D coupled Burgers' equation by using the FEDF for Test Problem 4.

Table 10 .
L ∞ error norm and CPU time for Test Problem 5.