Solving Parabolic and Hyperbolic Equations with Variable Coefficients Using Space-Time Localized Radial Basis Function Collocation Method

In this paper, we investigate the numerical approximation solution of parabolic and hyperbolic equations with variable coefficients and different boundary conditions using the space-time localized collocation method based on the radial basis function. The method is based on transforming the original d-dimensional problem in space into ðd + 1Þ-dimensional one in the space-time domain by combining the d-dimensional vector space variable and 1-dimensional time variable in one ðd + 1Þ-dimensional variable vector. The advantages of such formulation are (i) time discretization as implicit, explicit, θ-method, method-of-line approach, and others are not applied; (ii) the time stability analysis is not discussed; and (iii) recomputation of the resulting matrix at each time level as done for other methods for solving partial differential equations (PDEs) with variable coefficients is avoided and the matrix is computed once. Two different formulations of the d-dimensional problem as a ðd + 1Þ-dimensional space-time one are discussed based on the type of PDEs considered. The localized radial basis function meshless method is applied to seek for the numerical solution. Different examples in two and three-dimensional space are solved to show the accuracy of such method. Different types of boundary conditions, Neumann and Dirichlet, are also considered for parabolic and hyperbolic equations to show the sensibility of the method in respect to boundary conditions. A comparison to the fourth-order Runge-Kutta method is also investigated.


Introduction
Second-order parabolic and hyperbolic partial differential equations with variable coefficients are one of the most important problems in Mathematical and engineering fields. They attract the attention of a large number of scientists. The general heat transfer problem is considered one of the more interesting examples. Such initial boundary value problems (parabolic or hyperbolic) can be solved by coupling time-stepping algorithms with different numerical methods such as finite element, finite volume, boundary elements methods, meshless methods, fundamental solutions, and spectral and wavelet methods. In [1], Parzlivand and Shahrezaee presented a numerical technique based on a combination of collocation method and radial basis functions to solve parabolic partial differential equations with time-dependent coefficients. Dehghan and Shokri [2] have employed the thin-plate splines RBFs in the collocation RBF method to solve the two-space dimensional linear hyperbolic equation with variable coefficients, subject to appropriate initial and Dirichlet boundary conditions. A new version of the method of approximate particular solutions using radial basis functions has been proposed in [3] by Jiang et al. for solving elliptic partial differential equations with variable coefficients. In this paper [3], a comparison to Kansa's method and the method of fundamental solutions has also been investigated. Some other investigation of meshless methods to PDEs with variable coefficients is given in [4][5][6]. We can also mention the works of Gu et al. [7,8], where they used local versions of meshless methods leading to sparse matrices.
In most published works, these methods are based on first discretizing the time variable by applying any time time-stepping algorithms as implicit, explicit, Runge-Kutta, or others; and seeking the approximate solution at each instant t in a space domain problem.
The space-time methods have seen some investigation these last decades. Among them, we can mention the spacetime finite element method developed by Tayfun et al. in [9] for the computation of fluid-structure interaction problems. Klaij et al. have also developed the space-time discontinuous Galerkin finite element method for solving compressible Navier-Stokes equations in [10] and advection-diffusion problems in [11]. The technique has also been applied to shallow water flows by Ambati et al. in [12].
Although few works were published concerning the space-time meshless method for solving PDEs with independent variable coefficients [13][14][15][16][17][18], up to date and to our best knowledge, there is no investigation on the application of space-time meshless method for solving PDEs with either variable or time-dependent variable coefficients. In this paper, we investigate the application of space-time localized meshless collocation radial basis functions developed in [13] to a general second-order parabolic and hyperbolic problems, with variable coefficients.
The method is based on firstly transforming the considered d-dimensional evolutionary problem in space as a ðd + 1Þ-dimensional space-time one. The formulation starts by combining the d-dimensional vector space variable and 1-dimensional time variable in one ðd + 1Þ-dimensional variable vector and defining the space-time domain and its boundary. Then, the boundary conditions are sited on the global space-time boundary domain for PDEs with the first derivative with respect to time and on just one part of the boundary for PDEs with the second derivative with respect to time. The method does not need first to discretize all derivatives with respect to time and solve the problem in the space domain at each time level, as it is usually the case with many numerical methods. The developed formulation leads always to a square algebraic resulting system and benefit from the following advantages: (i) The time stability analysis is not discussed as it is the case for other time-stepping schemes as implicit, explicit, θ-method, etc.
(ii) Reducing the computational time as there is no need to recompute the matrix for the resulting algebraic system at each time level, unlike the case for others time integration methods used to solve PDEs with time-dependent coefficients (iii) Solving hyperbolic partial differential equations with variable coefficients as inverse problem, since just the hyperbolic equation with constant coefficients is less discussed in the literature and needs more sophisticated time integration technique The paper is organized as follows. In Section 2, we introduce the formulation of the parabolic and hyperbolic problem as space-time problem and the space-time localized RBF method implementation. Section 3 is devoted to the discussion of results obtained by solving different parabolic and hyperbolic examples in two-and three-dimensions in regular and irregular domains. A comparison of the given technique to the fourth-order Runge-Kutta method is also given in Section 4. We conclude in Section 5.

Space-Time Localized RBF Method Formulation
To recall the space-time localized RBF method defined in [13], let Ω ⊂ ℝ d be a bounded domain with a sufficiently regular boundary ∂Ω and consider the following time-dependent boundary value problem or Bu where P ðx, tÞ is a differential operator of second order with variable coefficients of the form: and B is a linear boundary operator depending on the kind of problem treated. The given functions f , g, u 0 , and u 1 are assumed to be sufficiently regular. The formulation of the evolutionary problems given by (1-2-3) and (4-5-6-7) as a space-time one starts by combining the space variable x and the time variable t in one ðd + 1Þ-dimensional vector. The constructed variable vector belongs to the space-time domain Ω T = Ω × 0, T½ represented by Figure 1 for some simple cases. The boundary of the new formulated domain Ω T is given by As discussed in [13], the formulation of the d-dimensional problem as a ðd + 1Þ-dimensional space-time one depends on the type of equations considered. For the problem given by Eqs. Bu For the case of the hyperbolic equation, we do not need any boundary condition on the part on the part of the space-time boundary characterized by Ω × ft = Tg. The problem is solved as an ill-posed one, and the algebraic system is square (see [1]). Then, the system has the new form: Bu In general, to construct the approximate solution of this news system, we select a number of distinct centers points fx j = ðx j , t j Þg Ni j=1 and fx j = ðx j , t j Þg N j=Ni+1 in the domain Ω T and on its boundary ∂Ω T , respectively. The finitedimensional space V Φ that contains the RBF interpolation functions is defined by where m is the order of the used conditionally positive defined function Φ. The approximate function is given by We can mention that the space-time domain Ω T = Ω × 0, T½ satisfies the interior cone condition since the domain Ω and the interval 0, T½ are also satisfying the interior cone condition. So, the lemma of error interpolation cited in [19] is satisfied for the case of the space-time domain (see also [13]).
The general form of the system of Eqs. (9)-(12) or (13)- (16) can be written under the following form:
Following Li et al. [20] and by Yao et al. [21], a localized influence domain Ω s T for each pointx s ∈ Ω T is selected. It contains a number n s of neighboring points fx ½s k Þg n s k=1 ∈ Ω s T to the selected centerx s . (see Figure 3). Then, the proximate function in the subdomain Ω s T is given by the expansion where fα k g n s k=1 are the unknown coefficients, k·k is the Euclidean norm, and Φ is the chosen RBF.
Applying Eq. (20) to fx ½s k g n s k=1 ⊂ Ω s T set points, we get n s × n s linear equations given bŷ where Then, the problem of seeking the expansion coefficients fα k g n s k=1 is transformed into a determination of the values of the solutionû at each center points fx ½s k g n s k=1 ⊂ Ω s T by using the equation Forx s ∈ Ω s T , we apply the differential operator Lx to Eq. (20) to obtain the following equation whereû ½s = ½uðx s 1 Þ, uðx s 2 Þ,⋯,uðx s n s Þ and Y ½s = LxΦ ½s ðΦ ½s Þ −1 .
The vectorû = ½uðx 1 Þ, uðx 2 Þ,⋯,uðx N Þ is incorporated in the system of (23) by adding zeros at the proper locations based on the mapping ofû ½s toû, and considering the Y 1×N as the global expansion of Y ½s 1×n s . The global system of Eq. (23) is then written under the form In the same way, selecting a centerx s on the boundary ∂Ω T , we have where Θ s = BxΦ ½s ðΦ ½s Þ −1 . In a global form, the system is then written as where Θ is the expansion of Θ ½s by adding zeros.

Boundary condition Boundary condition
Initial condition (s) Boundary condition (s)

Modelling and Simulation in Engineering
By collocating at all centers points fx j g N j=1 using Eqs. (24) and (26), we get the following sparse linear system of equations Y where Y = : ð28Þ The approximate solution at the interpolation points fVðx j Þg N j=1 can be obtained by solving the above sparse linear system of equations.

Numerical Simulations
The investigated test examples in this section are two-and three-dimensional parabolic and hyperbolic equations with variable coefficients given as follows: Knowing the shape of the space-time domain Ω T , the distributed nodes can be done uniformly or randomly without making any distinction between space and time variables. Through these numerical simulations, n s denotes the number of neighboring points in an influence space-time domain Ω T . To balance between the approximation quality of PDEs solution and sparsity of the algebraic matrix, a reasonable value of n s is chosen to be at least 2 dim ðΩ T Þ + 1 when dealing with Dirichlet boundary condition. More points can be added in the case of the Neumann condition depending the on problem treated (see [13] for more details). Seeking optimal values of shape parameter is still an open subject in the area of numerical approximation of PDEs. Many techniques have been proposed in the literature without any theoretical analysis. The N = ðN d + N b Þ × N t is the total number of nodes used, where N b and N d are the number of boundary and interior nodes of space domain Ω, respectively, and N t is the number of nodes on the time axe which can be either uniform or random.
In all proposed simulations, Neumann and Dirichlet boundary conditions are considered for parabolic and hyperbolic equations to show the sensibility of the method in respect to different boundary conditions on the space domain. Although any RBF can be used to test the technique, the multiquadric radial basis function (MQ) is the function used for all given simulations. The MQ-RBF is defined by , where ε is the shape parameter and r is the distance between two nodes.
The validity and the accuracy of the presented technique are demonstrated by considering the following discrete error-norms; the maximum absolute error (MAE), the root mean squared error (RMSE), and the L 1 er relative errors are defined by where u j andû j are the exact and the approximate solutions at ðx j , t j Þ, respectively. 3.1. One-Dimensional Test and Comparison with Runge-Kutta Method. To show the efficiency of the space-time method in word of stability and CPU time, we compare it to the 4th order Runge-Kutta method. The investigated example is the one-dimensional convection-diffusion PDE with variable coefficients defined as follows: The initial and the boundary conditions are defined according to the analytical solution given by uðx, tÞ = exp ð2 x − t 2 Þ, and the space domain is Ω = ½0, 1.
As the coefficients are depending on t, the matrix construction is called 4 times at each time step using the 4th order Runge-Kutta method, which enlarges the CPU time. According to Table 1, the CPU computation time of the Runge-Kutta method is much larger than that of the space-time technique. We can also mention that time stability is not guaranteed in the Runge-Kutta method or any other time-stepping algorithms. Besides the implementation simplicity of the demonstrated technique, results show that our method is faster and more accurate than the used Runge-Kutta method. Table 1 shows that the space-time method is stable and accurate compared to the RK4 method. The space-time method converges for different values of the final time T, keeping the same value of N t . From the same Table 1, it can be remarked that the RK4 method diverges for the simulation with T = 1 and dt = 0:001, which means that a smaller value of dt is needed to get the convergence. Regarding to the Table 1 and incase the results of the space-time method is discussed according to dt, we can remark that dt is ranged from 0:05 to 0:025 and results still more accurate. All results are obtained using the same number of nodes on the space domain N x = 20 and the final time T ranges from 0:1 to 1. The shape parameter of the MQ-RBF is chosen to be ε = 0:1.

Two-Dimensional Parabolic Equations.
For examples treated herein, we consider the considered as computation space domain Ω, a two-planar domain given by Ω = fðx, yÞ ∈ ℝ 2 : kðx, yÞk ≤ 1g. Two different distributions of nodes, uniform and circular, illustrated in Figure 4, are considered. For uniform distribution of nodes, we take N d = 326 and N b = 50 and N d = 320 and N b = 50 for circular one.
Example 1. As a first simulation, we consider the parabolic example with the analytical solution given by where aðx, y, tÞ = 1/1 + x 2 + y 2 and bðx, y, tÞ = ðx + 2Þðy + 2Þ.  Modelling and Simulation in Engineering The problem was solved using different values of both N t and shape parameter ε. Considering Dirichlet boundary condition on ∂Ω, Tables 2 and 3 illustrate results for N t = 10 and N t = 20, respectively. They also show the influence of shape parameter ε on the approximate solution. Using ε = 0:25 and ε = 0:43 for uniform and circular distributions of nodes, respectively, Figure 5 shows the absolute error of the approximate solution at t = 1 on the computation space-domain Ω for Dirichlet boundary condition on ∂Ω × ½0, T. Using uniform distribution of nodes, the obtained maximum absolute error and the root mean square error are MAE = 2:71 × 1 0 −5 and RMSE = 1:09 × 10 −5 , respectively. For the circular distributed nodes, we obtained MAE = 5:47 × 10 −5 and RMSE = 2:37 × 10 −5 . As the Ω T ⊂ ℝ 3 and the boundary condition are Dirichlet type, the used value of the number of neighboring points is n s = 7.
To demonstrate the influence of the boundary condition on the approximate solution and then on the technique, mixed Neumann-Dirichlet boundary condition is considered on ∂Ω. Dirichlet condition was used on the part of the space domain boundary defined by Γ 1 = fðx, yÞ ∈ ∂Ω, y > 0g and Neumann condition on the other part Γ 2 = ∂Ω − Γ 1 . For this where aðx, y, tÞ = ð1 + x 2 Þð2 − tÞ and bðx, y, tÞ = xyt + 1.
The computational domain and boundaries are the same as for Example 1. The simulation of this second example shows some difficulties compared to Example 1 because of the pseudosingularity that present at points with t = 2. Keeping the same neighboring points value n s = 7, Tables 6 and 7     Figure 7, we can remark that at t = 1, we obtain MAE = 2:51 × 10 −4 with ε = 0:53 for the uniform distributed nodes and MAE = 4:13 × 10 −4 with ε = 0:26 for the circular distributed nodes.
Following the same step of simulation done before, Example 2 was then solved by considering mixed Neumann-Dirichlet boundary condition: Neumann condition on Γ 2 and Dirichlet conditions on Γ 1 . On the boundary with the Neumann condition, the number of selected neighboring point is set to be n s = 9. Tables 8 and 9 give the found results found.
Example 3. The investigated problem is the widely used convection-diffusion flow problem of a Gaussian pulse given by the equation: In our simulation, we take σ 2 = 0:001 and ν = 0:01. The computational space-time domain Ω T is ½0, 1 × ½0, 1 × ½0, π /2. In Figure 8, we show the numerical solution at final time t = T face to the analytical solution, the total number of nodes is N = 30 3 . Table 10 displays the absolute maximum error obtained for different values of N; it can be remarked that the rate of convergence is near quadric. In all simulations for Example 3, the shape parameter is chosen to be ε = 5.

Two-Dimensional Hyperbolic Equations.
Although it has been shown that the technique gives good results solving the parabolic equation with different boundary conditions, the hyperbolic equation is still an important and interesting equation to be tested using the described methodology. In all simulations, the problem is treated as an ill-posed one as it has mentioned before and the algebraic matrix obtained is square.
Taking aðx, y, tÞ = xye −2t , the functions u 0 , h, and f are chosen according to the analytical solution given by uðx, tÞ = sin x sin y sin t.
Applying the "ill-posed-problem" technique [22], we obtain accurate results by setting the number of neighboring points to be n s = 9 in the entire domain. Table 11 shows results obtained of different errors, taking ε = 0:1 and using the space-time domain 0, 1½ 3 . Figure 9 shows that the MAE error decreases when increasing the number of nodes.
As in the example of the one-dimensional problem, the convergence rate of the technique using different uniformly distributed sources points and ε = 0:1 is given in Table 12 showing that the order of convergence is quadratic.

Three-Dimensional Parabolic and Hyperbolic Examples.
In this section, we show the robustness and the accuracy of the proposed method implemented in high dimension problems. The traditional three-dimensional problems are transformed into 4D problems.
Example 5. The problem investigated is the 3D heat transfer problem given by the equation: x, y, z, t ð Þ− a x, y, z, t ð ÞΔu x, y, z, t ð Þ= f x, y, z, t ð Þ: ð38Þ  With aðx, y, z, tÞ = xyze t , the initial and boundary conditions are chosen in such a way that the analytical solution is given by the function uðx, y, z, tÞ = e −t sin x sin y sin z. The space-time domain is chosen Ω = 0, 1½ 4 and the right-hand side is f ðx, y, z, tÞ = sin ðxÞ sin ðyÞ sin ðzÞe −t ð−1 + 3xyze t Þ.
Taking n s = 9, Table 13 shows the errors found in two simulations. The total number of nodes is N = N x × N t × N y × N z = 5 4 and N = 9 4 .
We can remark that good accuracy is obtained even with weak number of nodes. We note that similar accuracy is found using different values of the shape parameter ε in the interval ½0,1:5 Example 6. The considered hyperbolic 3D problem is given by We use aðx, y, z, tÞ = xyze −t and ε = 0:12. The computational space-time domain is Ω = 0, 1½ 4 . The boundary and initial conditions and the function f are chosen according to the solution uðx, y, z, tÞ = sin x sin y sin z sin t. In this example, the number of neighboring points is set to n s = 11 on the entire domain. The total number of nodes is N = N x × N t × N y × N z = 9 4 . The errors obtained for this example are MAE = 1:54 × 10 −3 , RMSE = 9:23 × 10 −5 , and L 1 er = 6:00 × 10 −4 . The same ill-posed technique is used with the squared algebraic linear system.

Conclusion
The local space-time RBF collocation method developed in [13] is extended to solve parabolic and hyperbolic equations with variable coefficients in the space-time domain. Following the same technique used in [13], the parabolic problem is solved using the governing equation as boundary condition on the boundary characterized by the final time T, and the hyperbolic equation is solved as an ill-posed problem with incomplete data boundary condition on Ω × ft = Tg. The first one-dimensional tested numerical simulation concerns the comparison of the 4th order Runge-Kutta method to the presented one. It has been shown that although the Runge-Kutta method did not need the inversion of the matrix and it is considered as an explicit technique, it needs more CPU time than the proposed one.
The numerical results obtained for other two-and threedimensional parabolic and hyperbolic problems show the performance of the technique and its accuracy for solving different examples. It has been demonstrated that the technique is simple, straightforward. The formulation of the hyperbolic    11 Modelling and Simulation in Engineering equation as an ill-posed problem and also its application to three-dimensional problems and quasisingular problems with variables coefficients are the most important ideas of the paper. We can also remark that the main advantages of the technique which are (1) The time stability analysis is not discussed as it is the case for other time-stepping algorithms as implicit, explicit, θ-method, etc.
(2) Reducing the computational time as there is no need to recompute the matrix for the resulting algebraic system at each time level, unlike the case for other    (3) Solving hyperbolic partial differential equations with variable coefficients as an inverse problem, since just the hyperbolic equation with constant coefficients is less discussed and needs more sophisticated time integration technique are more representative for partial differential equation with variable coefficients than with constant one. Further work will focus on the application of the developed technique to nonlinear equations. Stability analysis of the technique should also be investigated in further work.

Further Work
Further work will focus on applying the method to more real PDEs like the equation of the telegraph and the transient heat conduction problem in an anisotropic medium, all this PDEs will be investigated with variable coefficients.

Data Availability
The data used to support the study is available within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest