MESHLESS METHOD FOR THE STATIONARY INCOMPRESSIBLE NAVIER-STOKES EQUATIONS

Mathematical analysis is achieved on a meshless method for the stationary incompressible Stokes and Navier-Stokes equations. In particular, the Moving Least Square Reproducing Kernel(MLSRK) method is employed. The existence of discrete solution and its error estimate are obtained. As a numerical example for convergence analysis, we compute the numerical solutions for these equations to compare with exact solutions. Also we solve the driven cavity flow numerically as a test problem.

1. Introduction. The objective of this paper is to develop the numerical theory for the Galerkin formulation using the MLSRK(moving least square reproducing kernel) method especially for the stationary incompressible Stokes and Navier-Stokes equations.
Several methods for meshless approximations were proposed for various applications. We note that Smoothed Particle Hydrodynamics(SPH) by Gingold and Monaghan(1977) [1], Reproducing Kernel Particle Method(RKPM) by  [6,5], Diffuse Element Method(DEM) by Nayroles et al.(1992) [11], Element Free Galerkin Method(EFG) by Belytschko et al.(1994) [9] and Partition of Unity Finite Element Method(PUFEM) by Babuška and Melenk(1995) [10] were proposed. In particular, we are interested in the applications of Moving Least Square Reproducing Kernel Galerkin Method(MLSRK) proposed by Liu et al.(1996) [8] to the incompressible Navier-Stokes equations. One distinct advantage of this method over the standard finite element method is that it requires simple distribution of nodes, not the complex mesh generation depend on the geometry of the flow domain. Another advantage is that the desired regularity of the approximate solution can be readily achieved by introducing suitable window function with sufficient regularity.
Though there has been keen interest in developing meshless approximations for the Galerkin formulation of the partial differential equations in engineering, mathematical analysis on the existence and convergence criterion of discrete solution has not been made yet as far as we have known. In this paper, We have obtained the solvability and the convergence for successive approximation of solution of the Stokes and the Navier-Stokes equations, which results in the H 1 −error estimate of the velocity.
As a numerical example, we calculate the numerical solutions for the Stokes and the Navier-Stokes equations in two dimension and the errors of each components of the numerical solution are tabulated. Also the driven cavity flow is calculated numerically as a test case with non-zero boundary condition and the several plots for the numerical solution are shown in this paper.

Reproducing Formula.
Let Ω ⊂ R n be a bounded domain with smooth boundary ∂Ω and u(x) be a smooth function defined in Ω. Define the set of all basis polynomials of order less than or equal to m 1 · · · x αn n | |α| = α 1 + · · · + α n ≤ m}.
Here, the number of all components in P m (x) is (n+m)! n!m! . We choose a smooth nonnegative window function Φ(x) which has a compact support, say, supp Φ ⊂ B 1 (0). To describe the MLSRK approximation of u(x) with m−th order consistency, let us introduce a localized error residual functional x −x ρ and ρ > 0 is a dilation parameter. Minimizing the quadratic functional J(a(x)), we find that the minimizer a(x) satisfies here the matrix M (x) is so called the moment matrix defined by Since the polynomial basis P α (x)'s are linearly independent, M (x) is always invertible and det M (x) > 0. Now the local approximation of u(x) nearx is obtained as the following : For fixedx ∈ Ω, the manipulation is the standard weighted least square procedure.
Sincex is an arbitrary point in Ω for the weight least square procedure, we may choosex = x and we obtain a global approximation of u(x). More precisely the global approximation operator is defined by This formulation is so called reproducing kernel formulation by . For simplicity, we define the correction function as and the kernel function as then the global approximation (2) is written in a convolution form For this global approximation, any polynomial of order less than or equal to m satisfies and we call this property as m−th order consistency.
To show m−th order consistency for the above MLSRK approximation, let u(y) be a polynomial of order less than or equal to m. Then it is represented by the form where β = (β 1 , · · · , β n ) is a multi-index and x β = x β1 1 · · · x βn n for x ∈ R n . Here we note that the first coefficient c β0 (x) with β 0 = (0, · · · , 0) is the polynomial u(x), i.e., From the definition of moment matrix (1) and the reproducing kernel approximation where e β = [0, · · · , 1, · · · , 0] T is the β−th standard basis of dimension (n+m)! n! m! . Therefore, from the identity (5), we have the m−th order consistency of (4). Now we define the shape function to develop the theory of MLSRK for the Stokes and the Navier-Stokes equations. For a given set of nodes Λ = {x i |i = 1, · · · , NP }, 498 CHOE, KIM, KIM AND KIM employing discretized moment matrix we define discretized kernel function This set of functions will be used in this paper as MLSRK shape functions, which is called simply shape functions if there is no confusion. Also we will denote briefly In the finite element method, the mesh generation follows quasi-uniform or regular condition. Similarly, we have the following condition. Node set is not concentrated in some region of domain, and the support of each shape function overlap sufficiently many times with other shape function's support. The overlapping condition ensures the invertibility of the moment matrix defined in (6). In detail, we want each point x ∈ Ω is contained in the supports of at most L shape functions where L is some fixed positive integer independent of the number of nodes, i.e., for the support of shape function, the maximum number of intersection is bounded independent of the number of nodes. In this manner we define the followings.
where N (i, γ) be the number of nodes contained in B ρi (x i ). We also let ρ = γh be characteristic radius.
Remark 1. Let us call γ in Definition 2.1 as dilation ratio for regular node set. The parameter ρ i will be used as dilation parameter of shape function for all a α ∈ R NP , α = 1, ..., n.
Note that the above regular and admissibility conditions imply certain uniform condition for the node distance and support radius of shape functions. In making shape function, dilation parameter ρ i and node distance h i are depend on each x i , but with assuming regular node distribution, we may consider ρ i and h i as independent parameter with respect to each x i .

Projection Error Estimate.
For the convergence analysis, we need interpolation error estimate between the solution space and the projection generated by the set of shape functions. We define discrete projection and find projection error estimate.
NP } be the admissible set of MLSRK shape functions generated by the window function Φ for the regular node set Λ = {x i |i = 1, · · · , NP }. Let u(x) ∈ C 0 (Ω) be a function and ρ i > 0 is a dilation parameter defined in the Definition 2.1. We define the discrete projection as (7) and Here, m denotes the order of generating polynomial basis P m , ρ is the characteristic dilation parameter and h stands for the distance between the nodes.
The following theorem implies that the projection error converges to zero as we enlarge the node set holding regular condition.
here ρ stands for the characteristic dilation parameter, i.e., some positive number satisfying min (Ω) be given. By definition of projection, we have and taking derivatives yields, for |β| ≤ m, Using the above expansion and the reproducing property of shape functions, we have the identity Therefore we have Considering scaling, we find Hence we find the following inequalities using the Minkowski inequality Therefore we have, for |β| = k ≤ m, For the discrete projection, we find the following Sobolev type embedding theorem which will be used for the analysis of the Navier-Stokes equations.
Theorem 2.2. We assume supp φ i ∩ Ω is convex. Suppose v ∈ C 1 0 (Ω) and the window function Φ ∈ C m 0 (R n ), m ≥ 1. If p > n holds, then we have the following inequalities sup Proof. Let ρ be given and Λ = {x i |i = 1, · · · , NP } be the regular node set. By the definition of the discrete projection R m ρ,h , we have where x ∈ Ω and x α , α = 1, · · · , n is component of x. By Taylor series expansion of the function v(x) for given y, we can write

CHOE, KIM, KIM AND KIM
Hence we obtain From the linear consistency of the shape functions, the first summation of the above equation vanishes, i.e., Also our construction of the shape function φ i (x) implies and each x is contained in the finite number of supports of φ i (x)'s from the overlapping assumption of the shape functions. We set Λ(x) = {x i ∈ Λ|x ∈ support of φ i (x)} and the number of the element of Λ(x) is bounded by for some fixed number L. Also we note that support of Therefore we have the following estimates, from the Hölder inequality and the Minkowski inequality, Note that the last inequality in (12), the finite intersection property of shape functions is crucial. Also remember our assumption p > n.
The first part of the above inequalities results from the usual Sobolev embedding and the second part have been proved at (12). Therefore we have proved our theorem.

Test Function Space and Boundary Transformation.
To develop the theory using the Galerkin formulation based on MLSRK method, it is necessary to clarify the test function space. In this subsection we define the test function space denoted by V h 0 and describe how to construct this space. Also we prove the projection error estimate for this space.
Let Λ = {x i |i = 1, · · · , NP } be a regular node set and A = {φ i |i = 1, · · · , NP } be a set of admissible shape functions which is generated from the smooth window function Φ. First we divide node set Λ to three parts and re-indexing as the following for convenience. Λ Γ = {x 1 , · · · , x NΓ } is the set of boundary node points. Λ Γ * = {x NΓ+1 , · · · , x N Γ * } is the set of boundary influence node points, i.e., shape functions associated with these node points have non-vanishing values at some boundary node points. Λ 0 = {x N Γ * +1 , · · · , x NP } is the set of interior node points such that shape function correspond to this node point has compact support in Ω. We define the discrete function space V h by Since it is complicated to describe the function space V h 0 with shape function set A, we consider the following linear transformation.
We call this transformation as boundary transformation. Note that . This is closely related to the shape of window function Φ and the dilation ratio γ. Basically, the shape function is similar to the window function. Hence we need to design the window function Φ and to choose the dilation ratio γ so that the matrix Here φ i 's satisfy m-th order consistency, and dilation ratio is γ. We define Φ as proper window function of m-th order and γ m as proper dilation ratio if where n is some positive constant. Choose the proper window function of second order as where f is defined by Also choose the proper dilation ratio as As a summary, using boundary transformation, we defined V h 0 (Ω) from V h (Ω). Now we define the test function space for MLSRK method as The next thing to study is projection for C 0 (Ω), the space of continuous function with compact support in Ω. For this, we define new projection.
NP } is the transformed shape functions using the boundary transformation(13). Define a discrete projection of function in C 0 (Ω) as Here m, ρ and h represent order of generating polynomial basis, dilation parameter for φ i and node distance.
As a Corollary of Theorem 2.1, we find the following projection error estimate for R m ρ,h .

Then the following interpolation estimate holds
Proof. Since we have and from Theorem 2.1, it is enough to show Note that G −1 is bounded. From the assumption of proper window function, smallness of off-diagonal terms for G is guaranteed. Hence using the consistency of shape functions, we observe where G * is a matrix with G * = O (1). Now from the Sobolev inequality, we have For the second part of (16), In above inequalities, the inequality , is justified from the assumptions and considering Gram-Schmidt process for making d i j . Hence by following similar argument in Theorem 2.1, we have for all 0 ≤ k ≤ m, and this completes the proof.
There is a Sobolev type embedding theorem for R m ρ,h also, and we state it as a corollary of Theorem 2.2.
If p > n holds, then we have the following inequalities sup The above corollary can be proved directly from Theorem 2.2 by considering the Or one can prove it using Corollary 2.1 and the Poincaré inequality.

Applications of the MLSRK Method to Incompressible Flows.
In this section, we study the stationary incompressible flow with zero velocity on the boundary. Throughout this section, we suppose the following. Let Ω be a bounded domain in R n with smooth boundary. Λ V = {x V i | i = 1, · · · , NP } is a regular node set in Ω, and A V = {φ i | i = 1, · · · , NP } is an admissible shape function set with m-th order consistency. Shape function set A V is made by proper window function Φ(x) ∈ C m+1 0 (R n ) of m-th order and proper dilation ratio γ. After re-indexing, let Λ V = {x V i | i = 1, · · · , NP } be a set of interior node and boundary influence node of Λ V , and is a test function space as defined in the previous section. Λ V and A V will be used for velocity approximation. Note that even transformed shape function φ i ∈ A V does not vanish on the entire boundary ∂Ω, it has zero value only on the boundary node points. But the L ∞ norm of φ i ∈ A V is sufficiently small on the boundary of Ω so that φ i L∞(∂Ω) ≤ cρ 2m , as long as the proper window function is used for shape functions.
For the pressure, let Λ P = {x P j | j = 1, · · · , MP } be a regular node set in Ω, and A P = {ψ j | j = 1, · · · , MP } be an admissible shape function set with m-th order consistency from window function Ψ(x) ∈ C m 0 (R n ). Here we assume that the characteristic distance for pressure node is compatible to the characteristic distance for velocity node, and so is the characteristic support radius of shape function. Hence the discrete solution space for velocity and pressure will be as followings: In addition to the assumption on the node distribution, we assume the technical hypothesis. It is about the coercivity of the bilinear form which is induced from the viscous term of governing equations, i.e. we assume that there is a positive constant β 0 independent of ρ such that for all a ∈ R NP −NΓ , where NP is the number of nodes, N Γ is the number of boundary nodes and n is the dimension of space.
We denote the discrete projection for velocity as R m ρ,h and for pressure as S m ρ,h , in detail, we have For further analysis, we define the following inf-sup condition.
Definition 3.1. We say the pair of shape function sets ( for all P ∈ M h (Ω).
We find a sufficient condition for ( A V , A P ) to satisfy the inf-sup condition.
where α 0 > 0 is a constant and Φ i and Ψ i are the index sets defined as Then ( A V , A P ) satisfies the inf-sup condition.
Proof. Let {p 1 , · · · , p MP } be the pressure coefficients so that the discrete pressure is given. Let S 0 be the index set of pressure nodes. Choose p a1 such that and let S a2 = {j | supp ψ j ∩ supp ψ a2 = ∅}. Continuing this process, we have a subset of pressure coefficients such that From our assumption of finite intersection property, we have where H is the maximum number of elements of S aj for all j = 1, · · · , L. The maximum intersection number H is independent of the dilation parameter ρ when we refine or coarsen nodes. Thus, we obtain  Now, we choose the value of u α j adequately where u α j is to be the α-th component of the coefficient of the velocity shape function φ j for j = 1, · · · , NP. From our choice of a k 's, we note that all of the index sets Φ a k 's are mutually disjoint. For each k = 1, · · · , L, let us choose u α j 's such that, for all j ∈ Φ a k , If we normalize the chosen velocity coefficients by MP k=1 |p k | 2 1 2 and use (31), then Therefore, it completes the proof.

Remark 3. Theorem 3.1 implies that if every support of velocity shape function is located at the place where the pressure shape function is steep then the inf-sup condition is satisfied.
One can suggest to make node distributions to satisfy inf-sup condition as the following. First, make pressure node distribution x P i . For each pressure node point x P i , consider adjacent node points which are included in the support of shape functions at x P i . Now consider mid-points of line segments between x P i and adjacent node points. Adding all of those mid points to pressure node points, velocity node distribution is made.

Let us define
Then we obtain the system of the discrete Stokes equations Note that B T is the transpose of B andū = [ū 1 , · · · ,ū n ] T . From the hypothesis (21), we have for a positive constant β 0 , where a = [a 1 , · · · , a n ] with each a α ∈ R NP . Let b = (b 1 , · · · , b MP ) ∈ R MP be an arbitrary vector such that Then, the inf-sup condition, the hypothesis (21) and the admissibility of shape functions for pressure imply for some λ > 0 and β 1 > 0 independent of ρ, hence , from (34) and (36), there exists some β 2 > 0 independent of ρ such that for any vector b satisfying (35). Therefore then the pair (ū,p) satisfies our discrete Stokes equations (33). For the uniqueness, we consider two discrete solutions (ū 1 ,p 1 ) and (ū 2 ,p 2 ). Then the residual pair (v,q) = (ū 1 −ū 2  Since M is positive definite, we havev = 0. From the inequality (36), which impliesq = 0. Therefore the uniqueness holds.
Let (u, p) be the solution of the Stokes equations (32), and (U, P ) be a discrete solution of (33). Then we have error equations Adding and subtracting ∇u to ∇(U − R m ρ,h u), we have From the Hölder inequality, we have Using the divergence free condition (38) of U in the discrete sense and ∇·u = 0 in the continuous sense, we obtain Again using the Hölder inequality, Now using ∇·u = 0 in (32), we have From the error equation (37) and the inf-sup condition, we can estimate the pressure term in (41) as from the assumption of proper window function and the coercivity assumption (21). Likewise, we have Therefore we have the following pressure estimate From the above inequality (44) with (41), we obtain Now for the estimate of IV , suppose U − R m ρ,h u = i α i φ i , α = (α 1 , · · · , α NP ). Then we have, from the assumption of proper window function, the coercivity assumption (21). Likewise, we have Applying interpolation theorem (9) to the pressure term p − S m ρ,h p L 2 and applying (15) to ∇ u − R m ρ,h u L 2 , and combining all the above estimates (39), (40), (45), (46) and (47) we obtain the following H 1 estimate:

Now from the coercivity (8), (21) and interpolation theorem (15),
Also from the pressure estimate (44) and the interpolation inequality for the pressure we get Therefore we proved the following Theorem.
the solution of the discrete Stokes problem (33). If m > n 2 − 1, then the following error estimate holds Remark 4. When we apply the meshfree method which use the Galerkin formulation to the Dirichlet problem, terms involving boundary integral are not exactly zero. But as we have seen in the proof of the above theorem, those terms involving boundary integral can be controlled by choosing the window function and dilation parameter. Indeed, estimates for terms involving boundary integrals are somehow trivial, but we include for the completeness.

Navier-Stokes Problem.
In this subsection we study the stationary incompressible Navier-Stokes equations with zero boundary condition for the space dimension n = 2 and 3. Let (u, p) ∈ H 1 0 (Ω) × L 2 (Ω)/R be the solution of the Navier-Stokes equations −ν∆u + (u · ∇)u + ∇p = f The necessary regularity of f will be assumed. We will assume the same hypotheses for the velocity and the pressure nodes as in previous section, also we follow all notations and assumptions. The discrete solution (U, For the existence of discrete solution we prepare the following. Indeed, with the admissibility condition of shape functions and Sobolev type inequality (2.2), we have an L 2 inverse type inequality which provides a fixed point theorem.
for given ε > 0 and for some constant c(ε) depending on ε.
Proof. Note that F (x) is represented by For fixed ε, from the Sobolev type inequality (19) and finite overlapping property of shape function, we have 516 CHOE, KIM, KIM AND KIM Theorem 3.4. For the space dimension n = 2 or 3, let w ∈ H 1 0 (Ω) be divergence free in the discrete sense, i.e., Ω ∇·w ψ j dx = 0, j = 1, · · · , MP and F ∈ V h 0 (Ω) be a given function. Then we have for any ε ∈ (0, 1), where ρ > 0 is the dilation parameter of the shape function.
. Assume S m ρ,h is the discrete projection in terms of the shape functions ψ j 's. Since ∇·w = 0 in the discrete sense, and thus we can write Using Hölder inequality and interpolation inequality (9), we obtain Now from the previous lemma, we have and hence we have Theorem 3.5. There exists ρ 0 > 0 depending on f H −1 and ν such that if 0 < ρ < ρ 0 , then there is a discrete solution (U, P ) to the discrete Navier-Stokes equations (50) for space dimension n = 2, 3.
Proof. We only prove the three dimensional case. Indeed the two dimensional case is easier and we omit the proof. Let us introduce a solution map whose fixed point is a solution to the discrete Navier-Stokes equations. We define the finite dimensional function space Let W ∈ W be fixed. We define a map L : W → W such that U = L(W) satisfying First we claim the following. For given ν > 0 and γ > 0, there is ρ 0 > 0 so that if ρ < ρ 0 then for each w ∈ W(γ) the matrix M (w) is invertible. Here the matrix M (w), that is associated with the map L, is defined by Indeed, for any given Since ∇·w = 0 in the discrete sense, from the estimate (51), the following is obtained Hence, if ρ is chosen so small that then we have Since c α is chosen arbitrarily, M (w) is invertible. If we define the assembled matrices M(w) and B as in the case of the Stokes equations (see the proof of Theorem 3.2 for the definition of B), we obtain Since M(w) is invertible, the above matrix equation has a unique solution as long as B satisfies the inf-sup condition (23). Now we have to prove U ∈ W(γ) for sufficiently small γ > 0. From the energy estimates, if C ρ Thus the following is derived

CHOE, KIM, KIM AND KIM
The inequality (54) is derived from the estimates (52) and the following identities is an into map. To complete our proof, we have to show that the map L is continuous. Let W 1 and W 2 be arbitrary two functions in W(γ) and TakingŪ as a test function and using (52), the following estimate is obtained If 0 < ρ < ρ 0 , from the Sobolev embedding theorem, we obtain Consequently, we have shown the continuity of the map L by Now by the Leray-Schauder fixed point theorem, we prove our existence theorem for the incompressible Navier-Stokes equations.
Remark 5. Although we could show the existence of discrete solutions to the discrete Navier-Stokes equations without any restrictions on the size of the external force f as long as ρ is small, to prove uniqueness, we need an intrinsic smallness condition for the size of H −1 norm of the external force f for the three dimension. Now let us consider the uniqueness of solution for the discrete Navier-Stokes equations in the case n = 3. Suppose U 1 and U 2 be two solutions of the discrete Navier-Stokes equations. If we assume that 2 f H −1 ν 2 is sufficiently small, then, from the energy estimate, we have The above inequality implies that the solution is unique if Let (U, P ) is a solution of the discrete Navier-Stokes equations (50), and (u, p) is a solution of the Navier-Stokes equations (49). Then we have the following error for all V ∈ [V h 0 (Ω)] n and Q ∈ M h (Ω). Taking U − R m ρ,h u as a test function to the discrete error equation (55), we have the equality The first term I is estimated as follows Also we can obtain the following estimates for II, III and V , Terms V and V I can be estimated as in the case of the Stokes equations, we refer (46) and (47). For the pressure term IV , we apply the same procedure as in the Stokes problem, i.e., use the inf-sup condition which is the hypothesis for the shape functions φ's and ψ's. At first, we can obtain the followings We want to find the estimate of P −S m ρ,h p L 2 . As in the case of the Stokes problem, the inf-sup condition implies that Thus from the estimates (56) and (57), we obtain Note that the followings are true Combining the estimates of I, II, III and IV all together, we finally obtain the following L 2 -error estimates for ∇(U − u) such that Then, considering interpolation theorems Theorem 2.1, Corollary 2.1 and L 2 estimate for ∇(U − u), we obtain the L 2 estimates for u and p. Details of the proof is exactly same as the case of the Stokes equations. Hence we state the stability theorem of the MLSRK scheme for the Navier-Stokes equations.
then the following error estimate is obtained where C 0 depends only on ν.
It is interesting that the smallness of H −1 norm of f is necessary for the error estimates.

Numerical Examples.
Using the MLSRK method discussed in the previous sections, the numerical experiments are performed for the Stokes and the Navier-Stokes equations in two dimensional case. The results show good agreement with stability theorems that we proved.
The domain considered in our examples is the unit square Ω = [0, 1]×[0, 1] ⊂ R 2 . The zero boundary conditions for the velocity are imposed on the boundary of Ω. Node distributions for the velocity approximation and the pressure approximation are shown in Fig. 1. About the window function, we use the function defined by where S(t) is the cubic spline function defined below : The shape functions of the velocity and the pressure associated with the window function Φ(x, y) are drawn in Fig. 2. It turns out that such a distribution of velocity and pressure nodes satisfies the inf-sup condition of the definition 3.

Pressure Node
Support of ϕ i ψ j Support of Velocity Node Then its corresponding force f can be exactly calculated for each case of the Stokes and the Navier-Stokes equations. Under this situation, we made the numerical solutions U and P of the Stokes flow and the Navier-Stokes flow problems using MLSRK scheme. Fig. 3 shows the logarithmic scale plots of decay rates of relative errors for the Stokes problem and the navier-Stokes problem. As a reference, relative L 2 −, H 1 errors for U = (U, V ) and relative L 2 − error for P are tabulated in Table 1, Table  2 and Table 3  Though we discussed the problems only with the zero velocity condition on the boundary of Ω in this paper, the generalization for the problem of non-zero boundary condition is readily obtained. As an example of the problem of non-zero boundary condition, the driven cavity flow is calculated numerically. We assume  To impose boundary condition, we made boundary transformation for the shape function in the numerical code as we did for the previous example. In numerical  Table 3. Relative L 2 − errors of P − p experiment, we put the Reynolds number Re = 1, Re = 100 for the Stokes problem and the Navier-Stokes problem. The stream lines and the pressure are shown in Fig. 4, respectively. Here the number of velocity nodes is 6561 and that of pressure nodes is 1681. In Fig. 4, the symmetry of our numerical solution is well presented, which is intrinsic property of the Stokes flow in a symmetric domain. The skewness of the solution for the Navier-Stokes flow is also well illustrated.

5.
Conclusions. The meshless method for the Stokes and the Navier-Stokes problems are analyzed rigorously and we have performed several numerical experiments successfully. The numerical results show a good agreement with the theoretic analysis for the convergence of the discrete solution. The MLSRK method seems recommendable for the incompressible Navier-Stokes flow and the Stokes flow. In the MLSRK method for the incompressible viscous flow problems, only if we distribute the nodes in the computational domain, we can calculate the numerical solution as smooth as we want. Furthermore, this method has the advantage of mesh adaptation without doubt. Therefore, we may conclude that the MLSRK must be one of the promising methods recently discovered and much more intense researches are required for the more general applications to various directions.