Two-Dimensional Meshless Solution of the Non-Linear Convection-Diffusion-Reaction Equation by the Local Hermitian Interpolation Method

A meshless numerical scheme is developed for solving a generic version of the non-linear convection-diffusion-reaction equation in two-dimensional domains. The Local Hermitian Interpolation (LHI) method is employed for the spatial discretization and several strategies are implemented for the solution of the resulting non-linear equation system, among them the Picard iteration, the Newton Raphson method and a truncated version of the Homotopy Analysis Method (HAM). The LHI method is a local collocation strategy in 1 MSc. in Energetic Systems, carlos.bustamante@upb.edu.co, Universidad Pontificia Bolivariana, Medellín, Colombia. 2 PhD. in Computational Mechanics, henry.power@nottingham.ac.uk, University of Nottingham, Nottingham, UK. 3 PhD. in Computational Mechanics, whady.florez@upb.edu.co, Universidad Pontifica Bolivariana, Medellín, Colombia 4 MSc. in Technology Management, Universidad Pontifica Bolivariana, Medellín, Colombia. Universidad EAFIT 21| Two-dimensional meshless solution of the non-linear convection diffusion reaction equation by the LHI method which Radial Basis Functions (RBFs) are employed to build the interpolation function. Unlike the original Kansa’s Method, the LHI is applied locally and the boundary and governing equation differential operators are used to obtain the interpolation function, giving a symmetric and non-singular collocation matrix. Analytical and Numerical Jacobian matrices are tested for the Newton-Raphson method and the derivatives of the governing equation with respect to the homotopy parameter are obtained analytically. The numerical scheme is verified by comparing the obtained results to the one-dimensional Burgers’ and two-dimensional Richards’ analytical solutions. The same results are obtained for all the non-linear solvers tested, but better convergence rates are attained with the Newton Raphson method in a double iteration scheme.


Introduction
The Radial Basis Functions (RBFs) have been widely used in global and continuous interpolation of scattered data sets.Moreover, RBF collocation by using Multiquadric (MQ), Thin Plate Spline (TPS) and Inverse Multiquadric (IMQ) functions was considered as one of the best numerical techniques for multidimensional interpolation, in terms of accuracy and ease of implementation, among several schemes tested by Franke [1].It is fair to mention, that the TPS functions are optimal interpolants since they minimizes the functional ∫ ∂x i ∂x j d⃗ x for i = 1, • • • , n and j = 1, • • • , n.Recently, the RBFs have been employed as the base of meshless collocation approaches for solving partial differential equations (PDEs).The use of RBF interpolation technique has become the foundation of the RBF meshless collocation methods for the solution of PDEs, since the pioneer work on the Unsymmetric method by Kansa [2].Kansa use the MQ function to obtain an accurate meshless solution to the advection-diffusion and Poisson equations without employing any special treatment for the advection term (upwinding), due to the high order of the resultant scheme and the intrinsic relationship between governing equations and the interpolation.This strategy involved all the scattered nodes that cover the domain and therefore it produced a global fully populated matrix.With the aim of improving the Kansa's Method, Fasshauser [3] use Hermite interpolation to construct an RBF interpolating function which gives a non-singular symmetric collocation matrix.He concludes that the Hermitian (Symmetric) method performs slightly better than the Kansa (Unsymmetric) method.Jumarhon et al. [4] obtain a similar improvement using the Symmetric method and more recently Power and Barraco [5] attain better results by employing the Symmetric method for ing.cienc., vol.9, no.17, pp.21-51, enero-junio.2013.
a variety of problems including convection diffusion equation.Later, La Rocca et al. [6] implement the Symmetric method in the context of the Method of Fundamental Solutions (MFS) to solve transient convectiondiffusion problems showing better accuracy of the scheme in the cases of strong variable velocity field with respect to the traditional formulation.Also, Chinchapatnam et al. [7] solve the convection-diffusion equation using several RBFs by Symmetric and Unsymmetric method.Their results are in agreement with those presented by Fasshauser [3].
Full-domain RBF methods exhibit high-order convergence rates of the L 2 -norm of the difference between the interpolation function and a test (analytical) function in terms of the upper bound of the distances between points (see [8] for the MQ case), and are used by many authors to solve variety of problems (Laplace, Poisson, Helmholtz and Parabolic equations) showing better accuracy compared to traditional methods [9], [10], [11], [12].Nevertheless, the fully-populated matrix systems they produce lead to poor numerical conditioning as the size of the data-set increases.This problem is described by Shaback [13] as the uncertainty relation; better conditioning is associated with worse accuracy, and worse conditioning is associated with improved accuracy.As the system size is increased, this problem becomes more pronounced.Many techniques have been developed to reduce the effect of the uncertainty relation, such as RBF-specific preconditioners [14] and adaptive selection of data centres [15].However, at present the only reliable method of controlling numerical ill-conditioning and computational cost as problem size increases is through domain decomposition (see, for example, [16], [17], [18], [19]).
One of the first attempts in this direction is made by Lee et al. [20] who propose the local MQ approximation in which only the nodes inside the influence subdomain of one central node are used in the Unsymmetric method for solving the Poisson equation.Based on the above development Sarler and Vertnik [17] create an explicit scheme to solve transient diffusion equation.Later on, other schemes are implemented by modifying traditional methods using Symmetric and Unsymmetric method.Wright and Fornberg [19] generalise the Finite Difference Method (FDM) using RBF Hermite interpolation attaining high order local approximations.Moroney and Turner [21], [22] apply the RBF interpolation in the Finite Volume Method (FVM) to reconstruct gradients in a two-and three-dimensional non-linear diffusion problems by the Unsymmetric method.More recently, Orsini et al. [23] solve the diffusion convection equation by using the FVM and the Hermitian (Symmetric) and the Unsymmetric interpolation schemes.In the pure meshless RBF collocation strategies there are also recent developments.Divo and Kassab [18] solve the non-isothermal flow problem by implementing a localized Radial Basis method based on the formulation proposed in [17] and using a sequential algorithm.Afterwards, Stevens et al. [24] implement the Local Hermitian Interpolation Method (LHI) to solve transient and non-linear diffusion problems.Unlike the method proposed in [18], Stevens et al. [24] employ the Symmetric method considering the boundary operator at the local level.Stevens et al. [25] solve accurately the two and three dimensional convection diffusion equation by LHI method including the PDE operator and PDE centres in the approximation of the solution field for each local domain.The LHI method is the scheme used in the present work focusing in the solution of a generic version of the non-linear convection diffusion reaction equation.
Newton like methods have been widely use in the solution of nonlinear system of equations [26], [27].However, in the case of non-linear systems which arise from the discretization of non-linear PDEs other strategies are often used regarding a straightforward implementation as the Picard iteration scheme.For instance, Cui and Yue [28] develop a FDM scheme for systems of non-linear parabolic and hyperbolic differential equations by employing a Picard iteration strategy and Mohan et al. [29] employ Picard iteration in the solution of a non-linear diffusive transport equation.Despite of the ease of implementation, the Picard iteration strategy is limited by its slow convergence in terms of the L 2 -norm of the difference between the function value at present and previous iteration and its inability of dealing with strong nonlinearities unless variable transformations are available (for example see [30]).Nevertheless, such transformations require non-trivial operations in terms of computational cost and cannot be applied for general forms of the non-linear terms, i.e. each type of non-linear term needs a specific transformation to weaken its influence on the problem.
The Homotopy or continuation methods for the solution of nonlinear equation systems have been developed with the aim of reducing the influence of the initial guess in the validity of the approximation [31].For solving non-linear PDEs, Liao [32] proposes the Homotopy Analysis Method (HAM) with the aim of generalising the Boundary Element Method (BEM), obtaining accurate solutions for several linear and non-linear problems.Also, Liao [33] solves transient non-linear heat transfer in two-dimensional domains by using a fully implicit FDM for time discretization and the BEM for the spatial discretization in the sense of the HAM formulation.The numerical results obtained by Liao for a third order formulation are enough accurate in comparison to those given by an iterative BEM scheme showing that HAM is an efficient technique compared to Picard iterative schemes.Later on, the HAM in conjunction to numerical methods for the PDE spatial discretization has been applied to obtained numerical solution for a variety of non-linear problems (see [34], [35], [36]).
In this paper a meshless numerical scheme is proposed to solve a twodimensional generic convection-diffusion-reaction equation by employing the LHI method in conjunction to several strategies to solve the resulting set of non-linear equations.Among them the Picard iteration scheme, Newton Raphson method and a truncated version of the Homotopy Analysis Method.The main goal of this work is to evaluate which of the mentioned non-linear solvers performs better, in terms of computational efficiency, when solving the non-linear convection-diffusion-reaction equation discretized by the LHI method.In the following section the LHI solution of the generic convection-diffusion-reaction equation is detailed and the resulting non-linear equation system is presented.Then the application of the non-linear solvers to the resulting LHI discretization equation is explained starting with the Picard Iteration and following with the Newton Raphson method and the HAM.Afterwards the numerical results for the one-dimensional Burgers' equation and twodimensional Richards' equation are obtained and the convergence rates of the different solvers are compared.

|26 2 The Local Hermitian Interpolation Method
The RBFs are characterised for its only dependence on the Euclidian distance between test and trial points and its radial symmetric, i.e. its value does not change when the coordinates system is rotated.Other important features are the monotonically increasing and decreasing of the function with the distance from the origin (test point), adjustable softness and excellent convergence properties.The RBFs more often used in interpolation are presented in the Table 1.Where r = ∥x − ξ j ∥ is the Euclidian distance between a domain node x and a trial point ξ j .

Name
Function The MQ function converges exponentially and, like TPS, is conditionally positive define of order m > 0, since it is necessary to add a polynomial term of order m − 1 to guarantee a non-singular resulting matrix.Besides, MQ has a shape parameter c which allows to change the function slope near the collocation point.For small values of c the function turns in a cone shape, that will flattens as the shape parameter increases.The TPS behaves well in the interpolation of multivariable functions and does not have adjustable parameter, although it converges linearly.The Gaussian and IMQ functions (β < 0) are positive define and hence requires the addition of a polynomial term.In the present work the MQ function is used considering its exponentially convergence.
A generic non-linear diffusion-convection-reaction problem is given by the governing PDE (1) and the boundary condition (2).The functions D,⃗ u and k are, respectively, the diffusive, convective and reactive coefficients and S is the source term.The coefficients a and b are functions of the position and according to its values a specific boundary condition is applied at each boundary point.The outward normal vector is denoted as ⃗ n.
To implement the LHI method it is necessary to cover the solution domain with two sets of scattered points as shown in Figure 1.The two sets of points are the nodes where the solution and the PDE will be enforced by interpolation; hence they are called, respectively, solution nodes and PDE centres.Each solution node has an associated stencil, which is the subdomain where a specific interpolation function is valid.The stencil (big circles in Figure 1) contains three kinds of points which are the solution nodes, the PDE centres and the boundary points when the stencil intersects the domain boundary.For more information about the effect of the stencil nodal distribution on the performance of the method, see [24], [25].In the LHI method, the dependent variable is locally approximated by a function in terms of RBFs.The Hermitian characteristic of the |28 method is a consequence of using the governing PDE and the boundary operators to define the local approximation.For solving linear problems the governing PDE operator can be used directly, while in the non-linear case some linearized version of the operator has to be applied.In the present work, an auxiliary linear operator L x is built by using known values of the dependent variable ϕ * to compute the equation coefficients as shown in the following expression: (3) Then, the dependent variable ϕ can be approximated by the expression where k refers to the stencil number, Ψ is the RBF and N 1 , N 2 and N 3 are, respectively, the number of solution centres, the boundary points and the PDE centres contained in the stencil k.N 4 refers to the number of polynomial terms which depends on the dimensions of the problem and the RBFs used.Like in the Symmetric method [3], [25], the collocation procedure is performed by evaluating the variable approximation (4) at the solution nodes, the approximated boundary condition B x (ϕ (k) (⃗ x)) at the boundary points and the governing PDE L x (ϕ (k) (⃗ x)) at the PDE centres.Considering the generic operator Ξ[] the collocation procedure can be expressed as: ing.cienc., vol.9, no.17, pp.21-51, enero-junio.2013.
with ϕ centre as the unknown variable value, g(⃗ x i ) as the known function of the boundary condition, S(ϕ * , ⃗ x i ) equivalent to the known source term function and (k) is obtained where the symmetric matrix A (k) and the column vector d (k) are defined as: The matrix A (k) is non-singular as long as Ψ is chosen appropriately, i.e.MQ, TPS, IMQ functions [6], and provided that no two data-centres sharing linearly dependent operators are placed at the same location [25].According to the equations ( 5) and ( 6), the variable ϕ at a point ⃗ x located inside the stencil k is given by the equation (7), with So far a local interpolation equation is obtained for the stencil k, considering the coefficient α (k) in terms of the variable values at the solution centres included in the stencil.Then the original operator of the problem N x is applied to the dependent variable expressed as (7) and the resulting equation is evaluated at the central solution point ⃗ x k , to obtain the non-linear expression (9) which relates the unknown variable values contained inside the stencil k.

|30
By applying the equation ( 9) to each one of the N stencils related to the N solution centres scattered throughout the domain, a set of N nonlinear equations is obtained.It is important to notice that the solution attained after solving the non-linear system is not the final solution since the local interpolation is performed by using the values of the dependent variable at the last iteration.Hence a double iteration algorithm is required to achieve the final solution, i.e. an outer loop where the coefficients in the interpolation matrices, which are functions of ϕ and are contained in the PDE operator L[.], are updated, and an inner loop where the non-linear system of equations ( 9) is solved.The solution of the non-linear equation system, the double iteration algorithm and the possibility of solving the problem without the mentioned outer iteration are discussed in the next section.For a deeper analysis of the Hermitian collocation method as a spatial discretization strategy, see [4], [24], [25] where the order of convergence of the L 2 -norm of the error as a function of the nodal distance is studied as well as the computational complexity of the algorithm implemented, and [5], [24] where problems with nonhomogeneous Neumann boundary conditions are solved and the problem of instability in the solution is addressed.

Solution of the non-linear convection-diffusion-reaction equation
Four strategies for solving the non-linear equation system which arises from the LHI discretization of the generic non-linear convection diffusion reaction equation ( 1) are tested.The strategies are the Picard Iteration, Newton Raphson method with double iteration, Newton Raphson method with single iteration and the truncated Homotopy Analysis method with double iteration.The Picard iteration is the simplest strategy used here to solve the non-linear system.In this case, the linear problem represented by equation ( 3) is solved iteratively until convergence.It is important to mention that the equation ( 3) is equivalent to the expression (1) after substituting the coefficient and source term values with the ones obtained in the previous iteration.The solution of the linear problem is performed ing.cienc., vol.9, no.17, pp.21-51, enero-junio.2013.
as it is shown by Steven et al. [24].The stop criterion is the magnitude of the difference between the solution at the last and the current iteration

Newton-Raphson method
Newton-like methods has been widely used for solving non-linear equation systems.In this case, the Newton-Raphson method is employed with the aim of increasing the convergence rate and decreasing the computing times respect to the Picard iteration proposed above.The Newton Raphson method is a good option since it shows quadratic convergence rate when the initial guess is close to the problem solution [26] as long as the solution function is smooth.In this sense, the use of MQ functions for the approximation is appropriate since they are infinitely differentiable [1].For equation systems, the solution of the problem is given by the following equation: where [J] is the Jacobian matrix, [ϕ] is a column vector containing the solution at the iteration indicated by the super index and [F (ϕ n−1 )] is a column vector with the non-linear functions evaluated at the last iteration.In this case, the non-linear functions are the resulting equations of the LHI spatial discretization and are given by the expression (9).Both numerical and analytical Jacobian matrices are computed with the aim of verifying the performance of the scheme using the numerical one.The numerical Jacobian is obtained by means of a finite difference strategy according to the equation (11), where ∆ϕ is a small perturbation applied to the ϕ j variable.The perturbation ∆ϕ is found experimentally for each problem, i.e. several values are used and the greatest one with which the algorithm converges to an accurate solution is selected.

|32
The analytical Jacobian Ja ij is obtained by differentiating with respect to the variable ϕ j the non-linear function f i , which is the expression (9) applied to the stencil i. Regarding the interpolation matrix [A] (k) are not in terms of ϕ j but ϕ * j , the analytical Jacobian Ja ij is given by the expression where the Einstein summation convention is considered except for the terms with the index in parenthesis.The derivatives of ϕ (i) , its gradient and its Laplacian with respect to the variable ϕ j can be expressed as equations ( 13),( 14) and (15) according to (7).
Two algorithms have been implemented to solve the non-linear problem by the Newton-Raphson method.The first of them is the single iteration scheme in which the numerical Jacobian is computed with the interpolation matrix [A] as a function of the solution ϕ, i.e. no linearization of the differential operators employed in the interpolation is required.The single iteration algorithm can be described as follows.
2. Evaluate the non-linear equation ( 9) including the interpolation coefficients [α] in terms of ϕ 0 3. Find the numerical Jacobian according to the expression ( 11) by evaluating the non-linear set of function ( 9) in ϕ j +∆ϕ with j being the index of the solution nodes in the stencil.The interpolation coefficients are calculated for each component of the Jacobian and, in consequence, each local equation system has to be solved as many times as the number of solution nodes included in the stencil.
4. Compute the new value of the solution ϕ n+1 by using the equation (10).

Check the tolerance criterion ||
End, if it is satisfied, or go to the step 2 with ϕ 0 = ϕ n+1 , otherwise.
The second algorithm implemented requires a double iteration or two loop strategy since the non-linear problem is solved with constant interpolation coefficients [α].Thus the problem is split in a local linear problem that is solved in an external loop and a global non-linear problem which is solved by the Newton-Raphson method in an internal loop.The algorithm is given by the following sequence.

External loop: Compute the interpolation coefficients for each
stencil by evaluating the interpolation matrix [A] at the current solution value ϕ n (or ϕ 0 at the first iteration) 3. Internal loop: Evaluate the non-linear equation ( 9) in terms of ϕ n to obtain [F (ϕ n )] and calculate the numerical or the analytical Jacobian according to the expressions (11) or (12), respectively.
4. Compute the new value of the solution ϕ n+1 by using the equation (10).

Check the tolerance criterion for the linear or external loop ||
is the final solution, otherwise go to step 2 with ϕ 0 = ϕ n+1 .

Truncated Homotopy Analysis Method
The Homotopy Analysis Method (HAM) proposed by Liao [32] as a generalisation of the BEM is employed here for the solution of the nonlinear convection diffusion reaction problem.The HAM solves the nonlinear problem as a combination of linear auxiliary solutions.The nonlinear solution is expressed as the homotopy Θ which is constructed by performing a continuous deformation of the initial linear solution ϕ 0 .The zeroth order deformation equation is given by expression (16), where L x and N x refer, respectively, to the linear and the non-linear differential operators as they are defined in the linearized equation ( 3) and the governing PDE (1).The homotopy Θ is a function of the domain and the parameter p ∈ [0, 1] and, according to the deformation equation, The homotopy function Θ is considered smooth enough to assure the continuity of its derivatives respect to the parameter p. Then it is possible to define Θ as a expansion in Taylor series according to the expression (17), where 0 < λ < ρ < 1 with λ being a necessary parameter when the convergence ratio ρ is less than the unit.For brevity the m−order derivative of Θ respect to p for p = 0 is written as Θ m 0 (⃗ x).
The iterative formulation given by the following equation is obtained from ( 17) when the series is truncated at m = M , being M the order of the formulation.
Hence, the iterative solution of the non-linear problem is achieved by computing the mth-order derivatives of the homotopy Θ m 0 (m = 1, . . ., M ).Thus differentiating the zeroth order deformation equation ( 16) with respect to the parameter p, an expression for the m-th order homotopy derivative is obtained.The generalised mth-order deformation equation can be expressed as with the function f m given by equation ( 20) when m = 1 and ( 21) when m > 1.
The equation ( 19) is solved for Θ m 0 as a boundary value problem with the following boundary condition obtained from the original condition (2) as in the formulation proposed by Liao [33].
In this work, the homotopy series is truncated at M = 2. Hence, just the first derivative ∂Nx(Θ) ∂p has to be calculated.The equation ( 23) results from differentiating the generic expression (1) with respect to p and evaluating it at p = 0.After solving the first order deformation equation for Θ 1 0 , the equation ( 23) is evaluated by reconstructing the Θ 0 0 and Θ 1 0 gradients and Laplacians at the stencil centre using the LHI discretization.
To implement the HAM, a double iteration algorithm is used with the aim of updating both the local interpolation operators and the auxiliary linear operator L x .The HAM solution is achieved by performing the following algorithm.

External loop: Compute the interpolation coefficients for each
stencil by evaluating the interpolation matrix [A] at the current solution value ϕ n (or ϕ 0 at the first iteration).
3. Internal loop: Sequentially, calculate the right hand side of the mth-order deformation equation ( 19) and solve the resulting linear boundary value problem for each m, with m = 1, . . ., M .

4.
Compute the new value of the solution ϕ n+1 by using the equation (18).

Check the tolerance criterion for the linear or external loop ||
is the final solution, otherwise go to step 2 with ϕ 0 = ϕ n+1 .

Numerical results
The non-linear LHI scheme with the algorithms presented in the last section is verified by solving one and two-dimensional problems involving non-linear PDEs.The four non-linear solvers are the Newton Raphson method with single iteration and numerical Jacobian (NR1INJ), the Newton Raphson method with double iteration and analytical Jacobian (NR2IAJ), the Newton Raphson method with double iteration and numerical Jacobian (NR2INJ) and the Homotopy Analysis method with M = 2 (HAM2).All cases are solved in a PC INTEL CORE 2 DUO, 2.66 GHz, 2 GB RAM, using a FORTRAN 90 code with the IMSL libraries for linear algebra calculations.
The first problem is the one-dimensional steady Burgers' equation which can be obtained from the generic form (1).In the second case, the diffusive coefficient is an exponential function of the dependent variable obtaining a non-linear PDE which represent the two dimensional steady Richards' equation.In all problems, the shape parameter value is a function c = dh, where the average distance between the central point and the rest of the points in a stencil is given by h and d is a factor specified for each problem solved.
The analytical solution of the Burgers equation with Dirichlet boundary conditions ϕ(0) = u 0 and ϕ(L) = 0 is reported in [37] as the expression where the parameter R = Lu 0 µ is analogous to the Reynolds number and η is the root of the non-linear algebraic expression η−1 η+1 = exp(−ηR), which is obtained by an independent non-linear scalar equation solver.

One-dimensional Burgers' equation
The one-dimensional solution is obtained by means of the developed twodimensional code by setting Neumann boundary conditions ∂ϕ ∂x i n i = 0 on the extreme x 2 −constant lines (x 2 = 0.0 and x 2 = 0.2).Numerical solutions for different values of the parameter R are presented in Table 2 and Figure 2. According to the numerical solutions presented in the Figure 2 it is clear that as the R value is increased the ϕ-gradient near to x 1 = 1.0 is higher and, in consequence, obtaining adequate numerical solutions for high R depends on the spatial accuracy of the discretization scheme.In the Table 2, the global iterations and the ratio between the CPU time spending of the employed nonlinear solver and the Picard iteration scheme are shown.A homogeneous nodal distributions are employed and the shape parameter factor d is set experimentally.The boundary conditions are applied in the local approximation i.e. no global equation is evaluated at the boundary points.The RMS error is computed by the following equation: where ϕ num and ϕ are the numerical and the analytical solution, respectively, and N the total number of nodes.
The assessment of the non-linear solvers is based on the number of global iterations and the CPU time knowing the RMS error is the same for all the computations since they are obtained by means of the same LHI spatial discretization.Therefore the same solution is attained by all the solvers tested as can be seen in Figure 3    The NR1INJ strategy is the most expensive in time (2.50 to 3.16 times of the CPU time spending using Picard iteration) despite showing the highest global convergence rate and, in consequence, the fewest number of global iterations for R ≤ 100.For Re = 1000 the number of global iteration in the NR1INJ strategy is the highest due to the ill-conditioning of the local interpolation matrices which causes an inaccurate response to the small perturbations needed to compute the numerical Jacobian matrix.The double iteration strategies perform better than the single iteration one considering the CPU time spending and the stability of the method as the parameter R is increased.As shown in Figures (3) b. and (4) b. the convergence rate of the double iteration strategies are almost the same while in the NR1INJ case the rate is higher for R = 100 and smaller for R = 1000.As is expected, the Picard iteration strategy gives the lowest convergence rate.and NR2IAJ) present a high convergence rate similar in trend to the quadratic convergence rate expected from Newton like methods [26].Besides, the behaviour is the same for both the analytical and numerical Jacobian cases showing that the Jacobian obtained by finite difference scheme is a suitable strategy to solve non-linear problems without the need of computing an analytical expression for the Jacobian.The HAM2 strategy requires more iterations and CPU time than the double iteration Newton Raphson schemes.The convergence rate of the HAM2 can be improved by taking a higher value of M , nevertheless the CPU time will be higher.

Two-dimensional Richards' equation
The two-dimensional Richards' equation for the steady state is obtained from the general form (1) by setting D = e αϕ , u 1 = 0, u 2 = −αe αϕ , k = 0 and S = 0.By doing some algebra, the Richards' equation is reduced to the expression which is solved by the non-linear LHI strategies implemented.As the parameter α increases, the non-linearity becomes stronger and the gradients near to the boundaries higher, as in shown in the numerical results obtained for different α values (Figure 5).The two-dimensional domain is the square and the analytical solution is given by equation ( 28) applying the Dirichlet boundary conditions ( 29) and (30) [38].
ing.cienc., vol.9, no.17, pp.21-51, enero-junio.2013.The results presented in Table 3 are obtained using a homogeneous nodal distribution.The RMS error is the same for all the solvers employed and the shape parameter factor d is found experimentally as in the Burgers' equation

|44
Ingeniería y Ciencia case.The boundary conditions are applied at the local level.The number of global iterations and the CPU time increase with the parameter α but, unlike in the solution of Burgers' equation, the RMS error are much higher and it is not possible to attained a solution by using some solvers for α = 3 and α = 4.When α = 3 the NR1INJ strategy is unable to solve the problem whatever the initial values are and the HAM2 spends approximately 4 times the CPU time required by the Picard iteration strategy.In case of NR1INJ, the ill-conditioning of the interpolation matrices causes poor accuracy in the numerical Jacobian computation.The equivalence among the attained solutions is verified in the Figure 6 a. which shows the absolute error at line x 2 = 0.5 for α = 3.The Figure 6 b. allows to have a qualitatively idea of the high global convergence rates which the double iteration strategies present with respect to the Picard iterative scheme.The Figure 6 c. shows the slow convergence rate of the HAM strategy in the solution of the non-linear equation system while the double iteration Newton Raphson schemes present a similar convergence rate such as in the solution of the Burgers' equation (Figure 3 c and 4 c).A small difference can be noticed between the NR2INJ and NR2IAJ being the last (analytical Jacobian) better in terms of number of iterations and CPU time.However the analytical Jacobian computation increases the input data (coefficients plus their derivatives with respect to the dependent variable) making difficult the code implementation for the solution of a generic-convection-diffusion reaction equation compared to the numerical Jacobian strategy.Finally, in the case of α = 4 the HAM2 solver fails to obtained a solution and the Newton Raphson results are attained spending less CPU time compared to the Picard iteration case.After solving the one-dimensional Burgers' and the two-dimensional Ri-chards' equations the better non-linear equation solver to be used in conjunction to the LHI discretization is the double iteration Newton Raphson schemes NR2INJ and NR2IAJ.Although the NR2IAJ performs better than the NR2INJ, the implementation of the numerical Jacobian strategy (NR2INJ) is straightforward compared to the NR2IAJ since no analytical effort is required to compute the Jacobian.Hence the numerical results obtained above prove that the finite difference is a suitable scheme to compute the Jacobian.

Conclusions
A meshless procedure for the solution of the generic two-dimensional convection-diffusion-reaction equation is implemented.Several non-linear solvers for the solution of the equation system arising from the LHI spatial discretization are tested.The one-dimensional Burgers' and the two dimensional Richards' equations are solved to assess the developed non-linear solvers.Based on the number of global iterations and the CPU time, the best strategies are the Newton Raphson method in a double iteration scheme with numerical (NR2INJ) and analytical (NR2IAJ) Jacobians.The single iteration Newton Raphson (NR1INJ) scheme is more sensitive to the ill-conditioning of the local interpolation matrices and in consequence fails when the nonlinearity becomes strong.The convergence rate of the truncated Homotopy Analysis Method (HAM) with M = 2 depends on the strength of the nonlinearity, making the method unsuitable for strongly non-linear problems.The NR2INJ requires more CPU time than the NR2IAJ but its implementation for the solution of non-linear PDE systems is simpler.Therefore from the non-linear solvers tested the NR2INJ scheme is the most appropriate for the solution of generic convection diffusion reaction problems.

Figure 1 :
Figure 1: Representation of the stencils.Solution centres are denoted as •, boundary nodes as • and PDE centres as ×

Figure 2 :
Figure 2: One-dimensional Burgers numerical solution: a. Re = 1; b.Re = 100; c.Re = 1000 a. for Re = 100 and Figure 4 a. for Re = 1000, in which the absolute errors on the line x 2 = 0.1 are presented.

Table 1 :
Radial Basis Functions

Table 2 :
One-dimensional Burgers' equation solution by several non-linear solvers

Table 3 :
Two-dimensional Richards' equation solution by several non-linear solvers