Numerical Simulation of Partial Differential Equations via Local Meshless Method

Imtiaz Ahmad 1 , Muhammad Riaz 2, Muhammad Ayaz 2, Muhammad Arif 2 , Saeed Islam 2 and Poom Kumam 3,4,5,* 1 Department of Mathematics, University of Swabi, Khyber Pakhtunkhwa 23200, Pakistan; imtiazkakakhil@gmail.com 2 Department of Mathematics, Abdul Wali Khan University Mardan, Khyber Pakhtunkhwa 23200, Pakistan; riazmath83@gmail.com (M.R.); mayazmath@awkum.edu.pk (M.A.); marifmaths@awkum.edu.pk (M.A.); saeedislam@awkum.edu.pk (S.I.) 3 KMUTTFixed Point Research Laboratory, Room SCL 802 Fixed Point Laboratory, Science Laboratory Building, Department of Mathematics, Faculty of Science, King Mongkut’s University of Technology Thonburi (KMUTT), 126 Pracha-Uthit Road, Bang Mod, Thrung Khru, Bangkok 10140, Thailand 4 KMUTT-Fixed Point Theory and Applications Research Group, Theoretical and Computational Science Center (TaCS), 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 5 Department of Medical Research, China Medical University Hospital, China Medical University, Taichung 40402, Taiwan * Correspondence: poom.kum@kmutt.ac.th; Tel.: +66-2-4708-994


Introduction
In present time, the use of PDEs are indispensable to many fields of modern physics, engineering and technology.Kortewege-de Vries-Burgers' (KdVB) equation is of high importance for describing different mechanisms derived by Su and Gardner [1] and is given as follows with initial and boundary conditions U(x, 0) = h 1 (x) U(a, t) = h 2 (t), U(b, t) = h 3 (t), where α, β and γ are arbitrary constants.There are various physical applications of the KdVB equation such as dissipative effects in nonlinear plasma waves [2], shallow water wave description in viscous fluids [3], and propagation of waves in elastic tubes filled with a viscous fluid [4].The KdVB equation has been solved by different numerical methods including meshless methods [5,6], the Chebyshev spectral collocation method [7], sinc-collocation method [8], and extended B-spline collocation method [9].
The coupled Drinfeld's-Sokolov-Wilson equations is an important nonlinear PDE model which can be written as follows with initial and boundary conditions V(a, t) = f 5 (t), V(b, t) = f 6 (t). ( Different numerical methods such as the homotopy analysis method [10], meshless method of lines [11] and radial basis functions (RBFs) meshfree algorithm based on scattered data interpolation [12] can be found in the literature to obtain the solution of this model.
The three asset American-put option under the Black-Scholes model [26] is considered as a free boundary-value problem and can be written as [27] ∂U ∂t The variables x, y, z are space variables showing prices of the underlying assets, r, D i , σ i , ρ ij , T and U represent risk free interest rate, dividend paid by the i-th asset, volatility of the i-th underlying asset, correlation between the i-th and j-th assets, expiry time and value of the option respectively.
Using the penalty term approach [28,29], Equation ( 10) is converted into a fixed domain model in the following form where The boundary conditions are given as The functions G 1 , G 2 and G 3 in the boundary are the solution of appropriate two asset American-put option problem.Various numerical methods such as [27,[29][30][31][32] are utilized for solving American options.
Recently, the meshless methods have gotten broad attention for the solution of different kinds of PDEs arising in several different disciplines of engineering.The meshless character is one for the most important reasons for rising the demand of such types of methods.These methods eliminate the construction of difficult mesh and are easily extendible for interpolating multi-dimensional data as compared to other conventional methods.The numerical results of RBFs based algorithms have demonstrated that it is truly meshfree, accurate and east to implement.Some interesting models can be found in [30,[33][34][35][36][37][38][39][40].
In the global meshless method (GMM), in case of shape parameter dependent RBFs, the selection of the optimum value of the shape parameter value c and dense ill-conditioned matrix are counted as a major deficiencies.To ward off the deficiencies in the GMM, the researchers [41,42] suggested local meshless methods.In the last few years, local meshless methods have been employed for the solution to complex PDE models (see [32,42,43]).
In the current work, the local meshless differential quadrature collocation method (LMM) is proposed for the numerical simulation of one dimensional Kortewege-de Vries-Burgers', one dimensional coupled Drinfeld's-Sokolov-Wilson, two dimensional linear diffusion and two dimensional regularized long wave equations.The LMM is an accurate and efficient numerical technique which requires two steps to approximate time-dependent PDEs.In the first step, space derivatives are approximated using RBFs which reduced the given PDE into a system of ODEs while in the second step, the system-obtained ODEs are solved by an appropriate ODE solver.
Contents of the rest of the paper are organized as follows; The suggested numerical method is highlighted in Section 2, time integration scheme is presented in Section 3, stability of the method is discussed in Section 4, implementation of the method to different test problems is presented in Section 5 and finally some conclusions are drawn in Section 6.

Implementation of Numerical Method
In this section a description of the LMM [33,43] is given for numerical solution of the PDE models given in Section 1.To pursue this, the derivatives of U(x, t) at the center x i are approximated by the function values at a set of nodes in the neighbourhood of {x i 1 , x i 2 , x i 3 , ..., Substituting RBF ψ( x − x l ) into Equation ( 12) to find the corresponding coefficient λ Equation ( 13) can be written as follows where Matrix notation of Equation ( 14) is where From Equation (15), we obtain λ (m) From Equation (12) and Equation ( 16), we get where

Local Meshless Method for KdVB Equation
Using the LMM to Equation ( 1) in space, we get The semi-discretized model Equation ( 17) with corresponding boundary conditions Equation ( 3) is given as follows where the symbol * represents element-wise multiplication of two vectors.
The corresponding initial condition Equation ( 2) is given as Up to now, the spatial derivatives are approximated by RBFs and the given PDE is reduced to a system of ODEs.Next, any ODE solver is needed to solve this system of ODEs.

Time Discretization
In this paper, three integration techniques, explicit Euler method (EEM), implicit method (IM) and the Crank-Nicolson method (CNM), are used to solve PDE models given in Section 1 and its explanation is given as:

A θ-Weighted Technique for 2D Diffusion Equation
Applying θ-weighted technique to approximate Equation (6) in time, we have After applying the LMM (discussed in Section 2), Equation ( 21) leads to where L is the weight matrix of the differential operator L and I is an identity matrix.Equation ( 22) reduces to EEM, CNM and IM for θ = 0, θ = 0.5 and θ = 1 respectively.

Stability Analysis
In this section, we used the matrix method for the stability of the local and global meshless methods as reported in [32,43].22), where Û(n) is the exact solution and U (n) is the approximate solution obtained by the LMM.The error E (n+1) can be expressed as where is the amplification matrix.ρ(W) denotes the spectral radius of W so the sufficient condition for numerical stability is ρ(W) ≤ 1, which leads to the following formula where λ is eigenvalue of the matrix L.
The inequality Equation ( 24) reduces to CNM for θ = 0.5, that is, In the case of complex eigenvalue λ = c + id, where c, d are any real numbers, the inequalities Equations ( 25)-( 27) take the following form For FI M For CNM 1 + 0.5dt(c + id) Similarly, for FEDF 1 provided Re(λ) ≤ 0 in all the inequalities in Equations ( 25)- (31).
The inequalities in Equations ( 25)-( 27) guarantee unconditional stability in case of implicit and Crank-Nicolson schemes.Conditional stability in case of explicit scheme is also shown in Figure 1.The stability of meshless methods also depends on the value of the shape parameter c.The GMM is sensitive to the value of c in comparison to the LMM and the system gets ill-conditioned quickly by varying c.We investigated this dependence numerically in Figure 2. It can be seen from the figure that the LMM satisfies the stability condition for a wide range of this parameter.Figure 3 shows sparsity pattern of the LMM corresponding to local sub-domain three in the 1D case and stencil size five in the 2D case.In comparison, the coefficient matrix of the global meshless method is dense.From the comparison of the coefficient matrices of both the local and the global meshless methods, it is clear that local meshless method produces sparse banded matrices which can be efficiently solved by sparse matrix solvers.In case of the global meshless method, the coefficient matrix is fully dense and often ill-conditioned.In such a scenario, a traditional solver, like Gauss elimination method, produces inaccurate results.Special procedure like truncated singular value decomposition (TSVD )is a more appropriate choice in such cases.

Numerical Analysis
The performance of the LMM is demonstrated against some problems in this section.The proposed scheme works for any radial basis functions but in the present work; multiquadric (MQ), Gaussian (GA), and inverse quadratic (IQ) RBFs are taken into account for spatial discretization.The proposed LMM is capable of approximating the solution on both the uniform and scattered nodal points.The size of a local sub-domain in one, two and three dimensional cases were taken to be three, five and seven, respectively.Computational efficiency and less sensitivity to the value of the shape parameter c were the major gains of the LMM, as GMM was found very sensitive to the value of shape parameter c.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 GHs 4 GB RAM.
Accuracy of the LMM were measured via L abs , L ∞ and L R error norms which are given as follows where Û, U denote exact and approximate solutions respectively.
Problem 1.The 1D KdVB Equation (1) with α = 1, having the following exact solution [5] U where In Table 1, numerical results of the LMM for Problem 1 using different values of β and time t over the space interval [−20, 20] are shown.We took α = 1, γ = 0.1, dt = 0.001, N = 81 and MQ RBF with c = 100.Better results of the LMM as compared to the result given in [5] were clear from the table.In Figures 4 and 5 2. It is clear from the table that the LMM was a more accurate than the method in reported in [8].Problem 2. The analytical solution [44] of the 1D coupled Drinfeld's-Sokolov-Wilson system (4) is given as The numerical simulations for Problem 2 are reported in Tables 3 and 4 of the coupled Drinfeld's-Sokolov-Wilson system Equation ( 4) on the spatial domain [−4, 4] from t = 0.1 to t = 0.5 with N = 81 and dt = 0.0001.Tables 3 and 4 show the L ∞ error norm of the LMM along with the comparison of the meshless method of lines based on RBFs [11] for α = 1, k = 0.001 and α = 0.0001, k = 0.01 respectively.One can guarantee from this comparison that the results of the LMM are more effective.Tables 5 and 6 show numerical results of the LMM and the GMM of the coupled Drinfeld's-Sokolov-Wilson system Equation (4) for N = 81, dt = 0.0001, t = 0.5, and different values of c.The accuracy of the global RBF methods is heavily dependent on c and this is evident from Table 6 as well.From the comparison shown in Tables 5 and 6, one can observe that the value of c was not an issue in the LMM as compared to the GMM.
The proposed LMM is employed on the test Problem Equation ( 3) and the results are shown in Figure 6 by having a time step length of dt = 0.0001, the shape parameter value of MQ RBF c = 100 and N = 20 × 20, at different times, up to t = 1.Furthermore, the results obtained by the LMM (Crank-Nicolson) have been compared in order to verify the accuracy of the LMM with the results reported in [13,14], which confirm the validity of the presented method.
Table 7 shows numerical results in terms of L ∞ error norm and CPU time for Problem 4 with v = 2.06, N = 50, dt = 0.5, [−80, 100] and up to time t = 20.From this table, it is evident that the proposed LMM is accurate and efficient.Eliminating the need of meshing and approximating the solution using uniform or scattered points in rectangular as well non-rectangular domains is one of the salient feature of the LMM.Numerical results of the RLW equation on non-rectangular domains are shown in Figures 14 and 15   Figure 16 (left) is evidence to the fact that, in case of LMM, the sensitivity to the shape parameter value is no longer a major problem.In Figure 16 (right) the condition number κ is also shown.Table 8 shows the numerical solution and CPU time (in s) of three asset American-put option for Problem 5 using LMM.The time step was dt = 0.001 and we had a fixed value of shape parameter c = 10, whereas Table 9 shows the L abs error at P(1, 1, 1) and P(1.1, 1.1, 1.1).The numerical solution of the American-put option are shown in Figure 18.

Conclusions
In the present work, the local meshless method based on radial basis functions is applied to one, two and three dimensional PDEs.To check the accuracy and efficacy of the proposed scheme on both rectangular and non-rectangular domains, different problems have been considered.Results of the local meshless method are compared with exact/approximate solutions available in the existent literature.The local meshless method has been found to be a flexible interpolation method as it removes sensitivity of the shape parameter and produces the coefficient matrix well-conditioned.

Figure 1 .
Figure 1.ρ(W) versus dt for the explicit scheme (left) and Implicit and Crank-Nicolson scheme (right) with N = 21 for Problem 3.

Figure 3 .
Figure 3. Matrix structure of the LMM of local sub-domain three in the 1D case (left), local sub-domain five in the 2D case (right).
solution profiles are shown for different values of β.

Figure 5 .
Figure 5. Numerical solution of the KdVB equation for β = 0.5 (left) β = 1 (right), for Problem 1.The numerical solution of the KdVB equation for different time t in Problem 1 with time step length dt = 0.02, α = 2, β = 0.005, γ = 0.1, N = 101, spatial domain [−100, 100] are shown in Table2.It is clear from the table that the LMM was a more accurate than the method in reported in[8].

Figure 6 .
Figure 6.Comparison of the 2D linear diffusion with the method reported in [13,14] for Problem 3. The numerical results of the 2D linear diffusion Equation (6) on non rectangular domains are shown in Figures 7-11 for Problem 3 using GA RBF.The numerical results shown in Figures 7-11 are performed with dt = 0.0001, t = 1 and α = 1.These figures show the efficiency of the LMM in non rectangular domains in term of absolute error L abs norm for Problem 3.

Figure 7 .
Figure 7. Computational domain, numerical solution and absolute error for Problem 3.

Figure 8 .
Figure 8. Computational domain, numerical solution and absolute error for Problem 3.

Figure 9 .Figure 10 .
Figure 9. Computational domain, numerical solution and absolute error for Problem

Figure 11 .
Figure 11.Computational domain and absolute error for Problem 3.

Problem 4 .
Finally, consider a two-dimensional nonlinear RLW Equation (8) with the following exact solution U

Figure 16 .Figure 17 .Problem 5 .
Figure 16.c versus κ (left), N versus κ (right) of the LMM for Problem 4.Comparison of CPU time of the LMM and the GMM are shown in Figure17by taking dt = 0.01 and t = 10.It is clear from this figure that LMM was more efficient then GMM.

Figure 18 .
Figure 18.Numerical solution of the LMM on the line x = y = z for Problem 5.

Table 8 .
CPU time and option value for Problem 5.