Meshless Galerkin Method Based on RBFs and Reproducing Kernel for Quasi-Linear Parabolic Equations with Dirichlet Boundary Conditions

. The main aim of this paper is to present a hybrid scheme of both meshless Galerkin and reproducing kernel Hilbert space methods. The Galerkin meshless method is a powerful tool for solving a large class of multi-dimension problems. Reproducing kernel Hilbert space method is an extremely eﬃcient approach to obtain an analytical solution for ordinary or partial diﬀerential equations appeared in vast areas of science and engineering. The error analysis and convergence show that the proposed mixed method is very eﬃcient. Since the solution space spanned by radial basis functions do not directly satisfy essential boundary conditions, an auxiliary parameterized technique is employed. Theoretical studies indicate that this new method is very stable, though a parameterized problem is employed instead of the main problem.

1 Introduction various types of partial differential equations (PDEs). In these methods, system of algebraic equations is established for the whole PDEs by employing only a set of scattered nodes in the domain and on the boundary. Consequently, the method avoids difficulties arising from the mesh-based methods such as finite element method (FEM) and finite difference method (FDM). The main attractive features of RBfs is that they can be easily applied to higher dimensional cases due to the fact that they depend only on the distance between the nodal points.
The Galerkin method is a class of powerful tools to prove the existence of solutions for various problems and also to approximate solutions for linear and nonlinear PDEs. This method uses the weak-form of the underlying PDE which results in reducing the smoothness requirement of the approximation functions. Authors in [9,10] applied RBFs to Galerkin method for solving elliptic problems with Dirichlet boundary conditions by making an artificial Neumann boundary condition via some auxiliary parameters. They also proved the convergence and obtained an error estimate of their method in the weak form. Salehi et.al applied the so called moving least square method (MLS) to solving PDEs using polynomials and estimating smooth functions with an optimal accuracy [21]. An RBF meshless method was employed to solve two-dimensional magnetohydrodynamic (MHD) equations in [15] using variably scaled radial kernel. Also the RBF method based on collocation (spectral method) has been applied to some applications (see for example [5,16]). Dehghan et.al [14] solved the time fractional diffusion-wave equation by the method of lines using this technique. They first discretized the main problem by employing the Crank Nicolson method and then applied the meshless Galerkin method by the use of the auxiliary parameters technique presented in [9,10].
On the other hand, for about two decades, the reproducing kernel Hilbert space method (RKHS) which is an analytical approach, has been in progress for solving linear and nonlinear diverse problems [1,2,3,4,13]. Many hard and nonlinear problems such as system of boundary value problems [13], Bagley-Torvik and Painlevé equations [2], integro-differential equations of Fredholm operator type in the sense of the Atangana-Baleanu fractional operator [3] and ABC-Fractional Volterra integro-differential equations [4] have been successfully solved by this method.
We present a new method combining both Galerkin RBFs (GRBFs) and RKHS methods. For this purpose, we describe the method for the following special class of quasi-linear initial-boundary value problem.
x) + f (u(t, x)), in (t, x) ∈ (0, T ) × Ω, (1.1) Here Ω ⊂ R N , N = 2, 3, 4, . . ., is a convex and bounded polygon, I = (0, T ) is the time interval and g satisfies the required regularity conditions. Also f is a sufficiently smooth function in one variable without any growth conditions such as u 5 , e u or u 3 |u| as well as Fitz-Hugh-Nagumo or Allen-Cahn type nonlinearities like u 3 − αu with some positive α ∈ R. The existence, uniqueness and stability of the solution for the above equation have been studied in [7]. A general RKHS is defined as follows.
Definition 1. Let E be an arbitrary (non-empty) abstract set like a real Banach space and F (E) denote the set of all complex-valued functions on E. A RKHS on the set E is a Hilbert space H ⊆ F (E) with a function K : E×E → R that is called the reproducing kernel with reproducing property K p = K(., p) for all p ∈ E in which f (p) = (f, K p ) holds for all f ∈ H and p ∈ E. The Hilbert space H in which its corresponding reproducing kernel function is K, denote by H K .
We intend to improve the Galerkin method by combining it with an analytical method. To do so, we present a new hybrid method of both meshless based on RBFs and RKHS. Some suitable assumptions required in the error analysis and convergence are stated. There is an adapted auxiliary parameterized technique for elliptic problems [9,10]. We adapt this scheme for problem (1.1) and show that the solution of the auxiliary problem uniformly converges to the solution of problem (1.1) in the Sobolev norm. Further, the process of GRBFs is introduced for spatial semi-discretization. Then, a convergence theory is carried out for the approximate solution of the parameterized problem. The system of nonlinear ordinary differential equations of time variable will be dealt with by the RKHS method, analytically. Finally, we provide a convergence theory for numerical solution produced by hybrid GRBFs and the RKHS. The outline of this paper is as follows. The meshless Galerkin method and reproducing kernel Hilbert space will be presented in Section 2. The numerical analysis will be given in Section 3. Some numerical results will be exhibited in Section 4. Finally, we summarize the achievements in Section 5.

Hybrid meshless Galerkin method and reproducing kernel Hilbert space
In this section, first, we implement the proposed method and then present the convergence analysis showing that the approximate solution converges to the weak solution of (1.1). All over in this work we use the following spaces for any integer number k, m, with its norm for u belonging to it, as follows dt.
When k = 0, this space is denoted with L 2 [0, T ]; H m (Ω) which its norms is that the elements of this space have finite value in its norm.

Implementation of new method
In this subsection, the meshless Galerkin method based on RBFs is introduced for solving problem (1.1). As is well known, in the weak form of the main problem (1.1), an integral curve of directional derivative of the unknown function along the boundary ∂Ω appears in the formulation which is the challenging part of the problem. One way to avoid this difficulty is apply an auxiliary parametrized technique to convert a Dirichlet boundary condition to a Neumann one. To do so, problem (1.1) is converted to the following parametrized quasi-linear parabolic equation: where n is the outward unit vector of ∂Ω and β is an auxiliary parameter that it should be created enough big in our numerical computation. As well known, the corresponding variational problem of (2.1) finds the solution u β ∈ H 2 [0, T ]; H 1 (Ω) , such that for all v ∈ H 1 (Ω), where u(t) is the abbreviation of u(t, x) and We discretize the weak form in (2.2) by choosing a finite dimensional subspace x − x i )s radial basis function (RBF) create by Φ(x) = φ( x ) that is positive definite or m-order conditionally positive definite on R n with respect to the set of polynomials Π n m having total degree m − 1 or less when all nonzero a ∈ R n satisfying N i=1 a i p(x i ) = 0, for all p ∈ Π n m , we have that [5,8,11,15,16,18,21,23]. We then approximate u β by the discrete u β,M given by the expansion Inserting the expansion into (2.2) and setting v = Φ(x − x i ), provides the following system of ODEs, whose solution determines the coefficients u β,i , where A, B and C are symmetric positive definite matrices with ) and U M β,i (t) are vectors as follows, By solving the system of ODEs in (2.3) the unknown coefficients u β,i can be determined. Since in the current work, the compactly supported RBFs are employed, the above coefficient matrices are sparse. The Matlab 2017b is used to evaluate the integrals. We now obtain an one variable unknown vector function of the system ODE in (2.3) based on reproducing kernel Hilbert space.
In this space, one can precisely express the solution of a nonlinear differential equation by series of elementary functions [1,2,3,4,13]. We rewrite the ODE in (2.3) as follows where L β is a continuous linear operator from T ] : such that u(0) = 0 similar to that in [1]. Also H 1 and H 2 are Hilbert spaces with inner products , respectively (see [6,13] where i = 1, 2, 3, . . . , e j for j = 1, 2, . . . , M are the standard basis on R M and L * β is the adjoint operator of L β . Also Ψ i,j (t), i = 1, 2, . . . , j = 1, 2, . . . , M , are a complete basis functions for H 1 that they are independent [1,13]. The normal orthogonal basis functions Ψ i,j (t) for H 1 can be constructed by the Gram-Schmidt orthogonalization process on Ψ i,j (t) to example one can see Algorithm 1 of [1] to make the normal orthogonal basis functions and corresponding their coefficients, where β i,j k,l s are called the Gram-Schmidt coefficients of Ψ i,j . Similar with Algorithm 2 of [1], the method can be implemented by making a sequence convergent to the exact solution of (2.4) by the following process =0. The convergence of the above sequence in H 1 has been proven in [1,13]. This is due to converging the above sequence, W M β,n (t), to the exact solution, W M β (t), of problem (2.4) in H 1 .

Numerical analysis
In currently subsection, we intend to survey the growth rate of convergence of the approximate solution to the weak solution (1.1). We know that the ordinary differential equation (2.2) has a unique solution, u β (t; v), depending continuously on v ∈ B = v ∈ H 1 (Ω) : v Ω ≤ 1 . Since f is a Lipschitz function, there exists a positive and continuous function By boundary regularity property u • vanishes and using trace theorem gives The weak solution (2.2) can be obtained for v ∈ B, but B is a weakly compact set of H 1 (Ω) by the Banach-Alaoglu theorem [20]. Let L = sup L(B) < ∞ then the inequality (3.1) is obtained for the weak solution u β (t) as follows Now we obtain a upper bound for the weak solution of problem (2.1). Substi- Integrating both sides of inequality (3.2) yields Using Grönwall's lemma results in the following inequality where C(t) is a positive function bounded by e Lt . Also the following inequality is a straightforward result of the last two relations Integrating both sides over [0, t] gives The inequalities (3.3), for t = 0, and (3.4), for t = T , can result the following inequality . (3.5) By boundary condition (2.1), one can demonstrate that with respect to (3.6) and multiplication both sides of the inequality (3.5) by β, then we have .
Using (3.5) and (3.6) and applying Trace theorem, one can conclude that u β (t) converges to the weak solution of problem (1.1).
Applying trace theorem and using the bound of e β (t) ∂Ω,k in (3.5) yields Integrating both sides of the above inequality, gives (3.11) Addition of inequalities (3.10) and (3.11) for C > 0, results in For the case f = 0, the problem is reduced to a linear case for which one can show the following inequality. we have established the following inequality for Galerkin method applied to problem (2.1), where R β M u(t) is called the elliptic or Ritz projection on V N . We define it as the orthogonal projection with the inner product (∇., ∇.) Ω + β(., .) ∂Ω [6], hence (3.12) Using Cea's lemma [6] and definition of Ritz projection, we obtain the following inequality for all Theorem 4.1 of [22], yields Since in the current work, 1 β ≤ O(h), we obtain the error bound or, equivalently, by trace theorem The following theorem states a priori error estimate for GRBFs.
Proof. We write the error in the quasi-linear parabolic problem (2.1) as a sum of two terms, The first term θ(t) will be the main object of the analysis. In order to find an upper bound for θ(t), we note that by our definition of u β,M (t) and attention to (3.12) for all φ ∈ V M . From (3.14) we have Because the variables are independent, we have used the easily established fact that the operator R β M commutates with time differentiation. Since θ belongs to V M , we may choose φ = θ and conclude Here the second term is nonnegative. Also by the locally Lipschitz of f , we have d dt e β (t) Ω ≤ ∂ρ(t) ∂t Ω + L e β (t) Ω and integrating the above inequality over the interval [0, t], yields From Grönwall's lemma [6], Theorem 2 and the inequality (3.13), we can conclude From the proof process of Theorem 2, we can result C(n • , t) ≤ e Lt where L is the Lipschitz constant of f depending on the above bound of u. The following theorem states the growth rate of convergence for reproducing kernel Hilbert space method.
We obtain an upper bound for W M β H1 . One can obtain the solution of ODE (2.4) as follows where E(t) = |P| t 0 e − D1+βD2 s |P t C −1 |ds is an increasing bounded linear operator, D 1 and D 2 are the positive diagonal matrices, and where P is a unitary matrix and |.| is absolute function on the vectors and matrices. By using Grönwall's lemma, we get By consistency of matrix norms, it can be shown that, for any integer number n, there are two positive constants C 1 and C 2 such that ∀U ∈ H 1 The following inequality is derived from (3.16), (3.17) and (3.18), LetΦ(x) be a column vector of Φ(x − x i ). We know that . By Theorem 2 we find, a 2 i and Cauchy-Schwartz inequality, we get .
Integrating both sides of the above inequality over [0, T ] gives Thus we have Now the proof is complete.

Numerical investigation
In this section, we present some numerical results and compare them with the analytical solution to demonstrate the performance of the proposed method.
We cos(π(x + y)) − cos(π(x − y) w i (x, y) where where u h,δt,β denotes the numerical solution when selecting the parameters h, δt and β. The numerical errors for the unit square are shown in Figure 2 in three parts each showing the effect of one parameter, h, δt, or β on the r.l.s.e using log-log plot. The red charts in each part of the Figure demonstrate the predictions of the error via the theoretical results, and the blue curves denote the error based on the observed numerical results. Figure 2a indicates the graph of the r.l.s.e against the ratio of the spatial parameter h over the initial value h 1 while Figure  2b and c, respectively, show the error values with respect to the time and the auxiliary parameter β with the initial values δt 1 and β 1 . Also, Figure 2a and b show that the numerical results are in a good agreement with the theoretical conclusion. In Figure 2a, the r.l.s.e is plotted against the spatial parameter h/h 1 with constant values of the other parameters β = 10 6 , δt = 0.01 and the same plot is shown in Figure 2b for the time parameter δt/δt 1 with the parameter values β = 10 6 , h = 1 40 . Also in Figure 2c δt = 0.01, h = 1 40 . As is observed, by decreasing the parameters h and δt the accuracy is increased which confirms the theoretical results. But in the implementation of our method, the condition numbers of the stiff and mass matrices are growing up as h is declining and this causes serious computational challenges for spatial discretization ( see [17,18]). Since the closeness of the time trial points leads to the correlated basis functions of the produced RKHS method, the practical use of this method increases the error of calculating the Gram-Smith coefficients. Figure 2c shows the relationship between the auxiliary parameter β and r.l.s.e. We observe that the convergence rate is faster than the theoretical prediction. The green curve is a linear regression curve made by the r.l.s.e values excluding the last two error values which are farther from the others. This curve shows that the computational results may be more accurate than that of the theoretical expectation. It turns out that the error analysis stated in this article can be improved with respect to the parameter β (see [10]). The numerical errors shown in parts b and c are dominated by the mixed errors produced by two other parameters. By Theorem 4, we can guess the behavior of the upper bound of r.l.s.e. This analysis will be mentioned in details in the conclusion section.
A similar discussion can be made for the other regions but we do not present it in here. Some more results computed for different parameters h, β and δt = 0.01 are given in Tables 1 and 2 with uniform and non-uniform grid points, respectively. The accuracy of numerical results is apparently increasing when the boundary is smooth.
We construct the non-uniform grid points as x i = X i + d i cos(α i ) and y i = Y i + d i sin(α i ) where the point (X i , Y i ) is the i−th point in the uniform grid points, d i and α i are stochastic variables whose probability distribution are discrete uniform distribution with their support set [0, 1 π 2 ] and [0, 2π], respectively.

Conclusions
A mixed method of meshless Galerkin radial basis functions and reproducing kernel Hilbert space methods was proposed in this work. Also, the stability of the solution for the underlying problem was signified in the L p , Sobolev and Hölder spaces. In addition, the uniqueness and boundedness of the solution were considered. An auxiliary problem with Neumann boundary condition, related to the main problem, was obtained using a parametrization technique, and this boundary condition was imposed on the weak form of the problem (2.1). Moreover, we showed that the weak solution of the auxiliary problem converges to the weak solution of the main problem as β → ∞. An approximate solution was obtained for the auxiliary problem using the Galerkin method for spatial variables. We theoretically demonstrated that the new method was more accurate than the other meshless methods. The reproducing kernel Hilbert space method was employed to compute a sequence of approximate solutions for the system of nonlinear ordinary differential equations applying the Galerkin method to the auxiliary problem. It seems that the new method has the advantages of both methods. The approximate solution has the main properties of the Galerkin method based on the radial basis functions and also it is independent from the dimension of the problem. A convergence growth rate of O δt 2 + h k + 1 √ β was shown for the method . Finally, it turns out that if 1 √ β ≤ O(δt 2 + h k ), then the convergence depends only on the spatial-time approximation method without any dependence on the auxiliary parameter which seems to be a general property for any spatial-time discretization method.