Generalized finite difference solution for the Motz Problem

<!-- metadata commented in wiki content
==Generalized finite difference solution for the Motz Problem==

'''C. Chávez-Negrete<math>^2</math>, F.&nbsp;J.&nbsp;Domínguez-Mota<math>^{1,2}</math>, J.S Lucas-Martínez<math>^1</math>

J.G. Tinoco-Ruiz<math>^1</math>, D.&nbsp;Santana-Quinteros<math>^1</math> 

<math>^1</math>Facultad de Ciencias Físico-Matemáticas, U.M.S.N.H., México 
<math>^2</math>Facultad de Ingeniería Civil,U.M.S.N.H., México'''
-->
==Abstract==

In this paper, it is presented a formulation of a generalized finite difference scheme to solve the Motz problem. It is based on a general difference scheme defined by an optimality condition, which has been developed to solve Poisson-like equations whose domains are approximated by a wide variety of grids over general regions. Numerical examples showing second-order accuracy of the calculated solutions are presented.

'''Keywords''': Generalized finite difference, Motz problem

==1. The Motz problem==

The Motz problem is a benchmark for the Laplace equation problem that is often used for testing various special numerical methods; it was proposed in the literature for the solution of elliptic boundary value problems with boundary singularities. The boundary value problem is stated as follows (Figure [[#img-1|1]])

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{\partial ^2 u}{\partial x^2}+\frac{\partial ^2 u}{\partial y^2}=0,\qquad -1\le x\le 1,\quad 0\le y \le 1</math>
|-
| style="text-align: center;" | <math> u|_{\overline{OD}}=0</math>
|-
| style="text-align: center;" | <math> u|_{\overline{AB}}=500</math>
|-
| style="text-align: center;" | <math> \frac{\partial u}{\partial n}|_{\overline{OA}\cup \overline{BC}\cup \overline{CD}}=0. </math>
|}
|}

A singularity similar to that at the origin appears analogously in several hydrological applications, often related to an abrupt change in the boundary condition, particularly in groundwater infiltration in the alternation of pervious and impervious boundary zones. 

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
|-
|style="padding:10px;"|[[Image:Draft_Dominguez Mota_636399092-motz.png|498px|Motz problem domain and boundary conditions.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figure 1'''. Motz problem domain and boundary conditions
|}


The closed-form solution of Motz problem is given as a series
<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u=\sum _{i=0}^{\infty } A_i r^{(i+\frac{1}{2})}\cos \left[\left(i+\frac{1}{2}\right)\theta \right], </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">r</math> and <math display="inline">\theta </math> are the polar coordinates; its convergence ratio is 2, and therefore valid over the whole domain  (Figure [[#img-1|1]]).

Several papers discuss the accurate numerical calculation of the coefficients <math display="inline">A_i</math>. Lu <span id='citeF-1'></span>[[#cite-1|[1]]] calculated them with 17 significant digits using a collocation Treftz method. Li ''et al.'' <span id='citeF-2'></span>[[#cite-2|[2]]] used a special  boundary approximation method to derive the leading coefficients of the asymptotic solution expansion. Similar approaches were followed by Li <span id='citeF-3'></span>[[#cite-3|[3]]] and Arad ''et al.'' <span id='citeF-4'></span>[[#cite-4|[4]]], who  employed boundary methods combined with least-squares and penalty method techniques to handle the boundary information, and by Georgiou  <span id='citeF-5'></span>[[#cite-5|[5]]], who employed Lagrange multipliers.

Some other authors focus on solving the Motz problem using different numerical techniques. Bernal and Kindelan <span id='citeF-6'></span>[[#cite-6|[6]]] thoroughly discussed the performance of the radial basis function method in the solution of the Motz problem, including the global radial basis function method due to Kansa <span id='citeF-7'></span>[[#cite-7|[7]]], as well as  Shu's local  differential quadrature method based on radial basis functions <span id='citeF-8'></span>[[#cite-8|[8]]]. They showed that the exponential convergence obtained by increasing resolution and increasing the shape parameter is lost due to the presence of the singularity. The accuracy of the solution can be increased by using collocation at some boundary nodes to restore the convergence of the method, although it is necessary to use functions similar to the terms in equation ([[#eq-1|1]]), which are capable of capturing the behaviour of the solution near the discontinuity.

To produce an accurate solution based on the flexibility of meshfree methods, Hong <span id='citeF-9'></span>[[#cite-9|[9]]] addressed a singular basis function enrichment technique in a meshless method, using the leading terms of the local series expansion at the boundary singularity,  as enrichment functions for a local approximation space constructed to calculate the solution. The enrichment technique is derived from the Generalized Finite Element Method (GFEM) and produces accurate results, although it requires smooth global support basis functions and basis enrichment to  produce the numerical solution.

The aforementioned papers give an insight into the different  approaches considered in the problem. Finite elements and radial basis functions can be certainly used to produce accurate results near the singularity,  although it is also possible to consider the little restrictive and more versatile alternative given by the Generalized Finite Differences Method (GFDM), which has recently  been used to produce satisfactory approximations to solve differential equations in a wide number of applications modeled over irregular regions; for instance, on structured grids, a Generalized Crank Nicolson scheme is presented in <span id='citeF-10'></span>[[#cite-10|[10]]].  Chávez ''et al'' showed a stable calculation of the solution of nonlinear Richards equation which avoids oscillations in high gradient regions in <span id='citeF-11'></span>[[#cite-11|[11]]], and Tinoco ''et al'' addressed the conditions for stability in general schemes in <span id='citeF-12'></span>[[#cite-12|[12]]]. 

The GFDM is not limited to a grid structure, and can also be applied on unstructured grids and even as a meshless method in clouds of points: Benito ''et al'' studied its use on parabolic and hyperbolic equations <span id='citeF-13'></span>[[#cite-13|[13]]], advection-diffusion equations <span id='citeF-14'></span>[[#cite-14|[14]]], seismic wave problems <span id='citeF-15'></span>[[#cite-15|[15]]] and nonlinear elliptic equations <span id='citeF-16'></span>[[#cite-16|[16]]] . Fan ''et al'' have applied general finite differences on Burgers equation <span id='citeF-17'></span>[[#cite-17|[17]]], obstacle problems <span id='citeF-18'></span>[[#cite-18|[18]]], shallow water equations <span id='citeF-19'></span>[[#cite-19|[19]]] and natural double convection problems <span id='citeF-20'></span>[[#cite-20|[20]]].

In this paper, bearing in mind the engineering applications where the use of GFDM has been successfully applied, this paper shows an accurate numerical approximation to the solution of Motz problem.

It is organized as follows: generalized finite differences are presented in Section [[#2 General finite differences|2]]. The description of the meshes used is given in Section [[#3 Grids and stencils|3]]. Numerical tests are described in Section [[#4 Numerical tests|4]]; discussion and conclusions in Section [[#5 Results and discussion|5]], followed by conclusions in Section [[#6 Conclusions|6]].

==2. General finite differences==

A very common and widely used discretization procedure to approximate the solutions of differential equations is the finite difference method. This is because it can be easily applied in block-rectangular regions.

The well-known classic schemes for the first and second-order derivatives can be easily deduced from basic analytical considerations or by straightforward application of the Taylor theorem. Despite the basic idea of the theorem is quite simple, its use leads to more complicated schemes on non-rectangular regions when the geometry associated with the scheme changes. Nevertheless, when symmetry is involved, very interesting approaches can be considered; like those due to  Liszka <span id='citeF-21'></span>[[#cite-21|[21]]] for regular hexagonal geometries and to Benito <span id='citeF-22'></span><span id='citeF-13'></span>[[#cite-22|[22,13]]] for slightly perturbed clouds of nodes.

In the following paragraphs, we discuss the use of a general difference scheme that can be calculated by considering a  set of nodes <math display="inline">p_0=(x_0,y_0),p_1=(x_1,y_1),p_2=(x_2,y_2),...,p_k=(x_k,y_k)</math> where the solution of a differential equation is evaluated or approximated. For these nodes it is required to find coefficients <math display="inline">\Gamma _0,\Gamma _1,...,\Gamma _k</math> in such a way that consistency condition is fulfilled. In other words, the local truncation error to approximate the second-order linear operator

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L(u,A(x,y),B(x,y),C(x,y),D(x,y),E(x,y))=</math>
|-
| style="text-align: center;" | <math> A(x,y)\frac{\partial ^{2}u}{\partial x^{2}}+B(x,y)\frac{\partial ^{2}u}{\partial y^{2}}+C(x,y)\frac{\partial u}{\partial x}+D(x,y)\frac{\partial u}{\partial y}+E(x,y)u </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |(2)
|}

at the grid point <math display="inline">p_0</math> employing the combination

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L_0:=\sum _{l=0}^k\Gamma _l(p_0,p_l,A(x_0,y_0),B(x_0,y_0),C(x_0,y_0),D(x_0,y_0),E(x_0,y_0)) u(p_l), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

must satisfy the condition

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau (p_0):=[L(u,A(x,y),B(x,y),C(x,y),D(x,y),E(x,y))]_{(x_0,y_0)}-L_0\to 0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

as <math display="inline">p_1,...p_k\to p_0</math> <span id='citeF-23'></span>[[#cite-23|[23]]].

The coefficients <math display="inline">\Gamma _l</math> may change from one inner grid point to another in an irregular region, since they depend on the positions of the nodes <math display="inline">p_i</math> relative to <math display="inline">p_0</math>. Nevertheless, in the case of nearly rectangular grids and the Laplace equation, since the coefficients <math display="inline">\Gamma _l</math> are invariant under rotations and translations, essentially one set of coefficients is required.

For the sake of brevity, in the following paragraphs the coefficient

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \Gamma _l(p_0,p_l,A(x_l,y_l),B(x_l,y_l),C(x_l,y_l),D(x_l,y_l),E(x_l,y_l))</math>
|}
|}

will be denoted as <math display="inline">\Gamma _l</math>.

Let <math display="inline">\Delta x_l</math> and <math display="inline">\Delta y_l</math> be the <math display="inline">x</math> and <math display="inline">y</math> components of <math display="inline">p_l-p_0</math>, respectively. Thus, the local truncation error <math display="inline">\tau (p_0)</math>  (up to second-order) yields

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau (p_0)  =  </math>
|-
| style="text-align: center;" | <math>   \left(E(x_0,y_0)- \sum _{i=0}^k\Gamma _i \right)u(p_0) + \left(C(x_0,y_0)-\sum _{i=1}^k\Gamma _i\Delta x_i\right)\frac{\partial u}{\partial x}(p_0)  +  </math>
|-
| style="text-align: center;" | <math>   \left(D(x_0,y_0)-\sum _{i=1}^k\Gamma _i\Delta y_i\right)\frac{\partial u}{\partial y}(p_0)+\left(A(x_0,y_0)-\sum _{i=1}^k\frac{\Gamma _i (\Delta x_i)^2}{2}\right)\frac{\partial ^{2}u}{\partial x^{2}}(p_0)+</math>
|-
| style="text-align: center;" | <math>   \left(-\sum _{i=1}^k\Gamma _i\Delta x_i\Delta y_i\right)\frac{\partial ^{2}u}{\partial x \partial y}(p_0)+\left(B(x_0,y_0)-\sum _{i=1}^k\frac{\Gamma _i (\Delta y_i)^2}{2}\right)\frac{\partial ^{2}u}{\partial y^{2}}(p_0)+</math>
|-
| style="text-align: center;" | <math>   \mathcal{O}\left(\max \{ \Delta x_i, \Delta y_i\} \right)^3.   </math>
|}
|}

Therefore, assuming smoothness of <math display="inline">u</math>, to get second-order local consistency the coefficients of the derivatives of <math display="inline">u</math> at the right-hand side of the expansion must satisfy the linear system

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left( \begin{array}{rrrrr}1&  &1 &... & 1\\ 0&  &\Delta x_1 &... &\Delta x_q\\ 0&  &\Delta y_1 &... &\Delta y_q\\ 0&  &(\Delta x_1)^2 &... &(\Delta x_q)^2\\ 0&  &\Delta x_1\Delta y_1 &... &\Delta y_q\Delta z_q\\ 0&  &(\Delta y_1)^2 &... &(\Delta y_q)^2\\ \end{array} \right) \left( \begin{array}{c}\Gamma _0\\ \Gamma _1\\ \Gamma _2\\ \cdot\\ \cdot\\ \cdot\\ \Gamma _q\\ \end{array} \right) = \left( \begin{array}{c}E(x_0,y_0)\\ C(x_0,y_0)\\ D(x_0,y_0)\\ 2A(x_0,y_0)\\ 0\\ 2B(x_0,y_0)\\ \end{array} \right). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

One must note that, in general, this system is not well-determined.  This poses the calculation of the coefficients in terms of a well studied algebraic problem which can be solved robustly.  Hence, the main question is how to produce a solution of equation ([[#eq-5|5]]) which provides a consistent second-order scheme. Several alternatives are possible, for example, in <span id='citeF-10'></span>[[#cite-10|[10]]], an efficient heuristic calculation to solve numerically the diffusion equation was proposed. It was based on an unconstrained optimization problem for every inner grid point. To apply it, first, we separate the first equation of the matrix system ([[#eq-5|5]])

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum _{i=0}^q \Gamma _i=E(x_0,y_0) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

and then we solve the system

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>M\Gamma =\beta , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where
<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>M=\left( \begin{array}{clcl}\Delta x_1 &... &\Delta x_q\\ \Delta y_1 &... &\Delta y_q\\ (\Delta x_1)^2 &... &(\Delta x_q)^2\\ \Delta x_1\Delta y_1 &... &\Delta x_q\Delta y_q\\ (\Delta y_1)^2 &... &(\Delta y_q)^2\\ \end{array} \right), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Gamma =\left( \begin{array}{c}\Gamma _1\\ \Gamma _2\\ .\\ .\\ .\\ \Gamma _q\\ \end{array} \right),\quad  \beta = \left( \begin{array}{c}C(x_0,y_0)\\ D(x_0,y_0)\\ 2A(x_0,y_0)\\ 0\\ 2B(x_0,y_0)\\ \end{array} \right). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

Next, <math display="inline">\Gamma _0</math> is obtained from equation ([[#eq-6|6]]).

Equation  ([[#eq-7|7]]) defines the generalized finite differences in general grids. To solve it is straightforward when the grid is rectangular or nearly rectangular, and a sufficient number of nodes or symmetries are available. 

In irregular grids or clouds of points,  which lack symmetry, the system defined by <math display="inline">M</math> is underdetermined if <math display="inline">q>6</math>, which means that, in that case, at least seven neighbor nodes around <math display="inline">p_0</math> are required to get local consistency. 

On the other hand, in a rather regular structured -or mapping- grid, symmetries reduce the number of required neighbors. An important consequence is that the symmetric location of nodes around the central node avoids rank deficiency of <math display="inline">M</math> in the system ([[#eq-7|7]]).

==3. Grids and stencils==

The main issue to consider when solving the Motz problem is the singularity in the partial derivative <math display="inline">\frac{\partial u}{\partial x}</math> at the origin of coordinates. The usual approach to produce an accurate approximation at that point is to refine the grid around it. However, this causes not only an increase in the computational cost, but often the results are not satisfactory due to the high gradient values of the solution on the positive <math display="inline">x</math>-axis; even with robust techniques as those based on radial basis functions <span id='citeF-6'></span>[[#cite-6|[6]]], a special treatment is required near the origin. Some other approaches, like that of Hong, require a rather sophisticated process of patch mapping to produce accurate results <span id='citeF-9'></span>[[#cite-9|[9]]].

When finite differences are considered, it is important to emphasize that the method has been successfully applied to problems with strong boundary layers, for instance, by Ivanenko  <span id='citeF-24'></span>[[#cite-24|[24]]] and Pereyra <span id='citeF-25'></span>[[#cite-25|[25]]], when the equidistribution principle is addressed. In other words, the key is to produce a suitably adapted grid through a simple procedure. In this paper, to derive the rule that controls the refinement around the origin, the clue is given by the shape of the significative terms in the closed-form solution series of equation ([[#eq-1|1]]). 

Bearing in mind the equidistribution principle as discussed by Ivanenko,  since the graph of the solution <math display="inline">u(x)</math>  is  monotonic along the positive <math display="inline">x</math>-axis, and resembles that of  a boundary layer, there is a mapping <math display="inline">\xi \mapsto x=\xi ^\alpha </math> such that the expression

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{\partial u}{\partial \xi }=\frac{\partial u}{\partial x}\frac{\partial x}{\partial \xi } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

is constant. According to eq. ([[#eq-1|1]]), each term in the solution along such axis equals <math display="inline">x^{(i+\frac{1}{2})}</math>, <math display="inline">i=0,1,...</math>. Therefore, the constant right-hand side of ([[#eq-10|10]]) takes the form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{\partial u}{\partial x}\frac{\partial x}{\partial \xi }=\left[\frac{\partial x}{\partial \xi }\right]\left[\left(i+\frac{i}{2}\right)x^{i-\frac{1}{2}}\right],\qquad i=0,1,...</math>
|}
|}

This implies that the partial derivative <math display="inline">\frac{\partial x}{\partial \xi }</math> must have the form <math display="inline">x^{(\frac{1}{2}-i)}</math> and, in consequence, the grid refinement is of the form

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\xi \mapsto x=\xi ^{\frac{2}{2i+1}},\qquad i=0,1,...</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

According to the estimations presented by Lu <span id='citeF-1'></span>[[#cite-1|[1]]], the first five coefficients in the equation ([[#eq-1|1]]) are the  significative ones: Those terms correspond to <math display="inline">i=0,1,...,5</math>. Since the value <math display="inline">i>1</math> in equation ([[#eq-11|11]]) represents an increase in the distance between the grid nodes around the origin, only the value <math display="inline">i=0</math>, which corresponds to <math display="inline">\alpha=2</math>, is considered.

Thus, the mesh adaptation procedure is as follows:
<ol>
<li>The grid for the whole domain of the problem is generated in the first step as a regular rectangular grid. </li>
<li>Every cell is subdivided into two triangles in the second step to define six neighbours in the stencil for every inner grid node. </li>
<li>Finally, the grid is adapted by modifying the grid coordinates following the mapping given by equation([[#eq-11|11]]). </li>
</ol>


In the numerical tests, the interval <math display="inline">[1.5,2.5]</math> was considered for the exponent <math display="inline">\alpha </math>. For comparison purposes, the mapping was applied to the <math display="inline">x</math>-axis (Type I grids), to the <math display="inline">y</math>-axis (Type II grids), and both axes (Type III). An example of the grids corresponding to the different distributions is shown in Figure [[#img-2|2]].

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto%;"
|-
|style="padding:10px;"|[[Image:Draft_Dominguez Mota_636399092-quadx.png|180px|]]
|style="padding:10px;"|[[Image:Draft_Dominguez Mota_636399092-quady.png|180px|]]
|style="padding:10px;"|[[Image:Draft_Dominguez Mota_636399092-quad.png|180px|Examples of the kinds of grids used in the tests ]]
|- style="text-align: center; font-size: 75%;"
| colspan="3" style="padding:10px;"| '''Figure 2'''. Examples of the kinds of grids used in the tests 
|}


The three kinds of grids define a regular location of neighbors at the inner nodes: star-shaped subgrids  with 6 neighbours that can be used to approximate the derivatives in the left-hand side of equation ([[#eq-6|6]]); an example of such subgrid can be seen in Figure [[#img-3|3]]. The Laplacian operator is approximated at the point <math display="inline">p_0</math>, and the corresponding neighbours are <math display="inline">p_1</math>, <math display="inline">p_2</math>, <math display="inline">p_3</math>, <math display="inline">p_4</math>, <math display="inline">p_5</math>, and <math display="inline">p_6</math>. The values <math display="inline">h_1</math>, <math display="inline">h_2</math>, <math display="inline">k_1</math>, and <math display="inline">k_2</math> are the step size lengths in the stencil. Bearing in mind equation ([[#eq-5|5]]), this poses a problem with 6 equations and <math display="inline">k=7</math> unknowns; since the values <math display="inline">h_1</math>, <math display="inline">h_2</math>, <math display="inline">k_1</math>, and <math display="inline">k_2</math> are non zero, after preconditioning of equation ([[#eq-6|6]]) by scaling, the rank of <math display="inline">M</math> is equal to five,  which allows a robust calculation of the coefficients and leaves one degree of freedom available. 

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
|-
|style="padding:10px;"| [[Image:Draft_Dominguez Mota_636399092-inner.png|498px|Selection of neighbours for an inner grid node]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"|  '''Figure 3'''. Selection of neighbours for an inner grid node
|}


Regarding the boundary condition, a straightforward discretization using  the GFDM method was addressed in <span id='citeF-26'></span>[[#cite-26|[26]]]: at a boundary node of a structured grid, equation ([[#eq-7|7]]) is evaluated at its five neighbors, taking the right-hand side <math display="inline">\beta </math> to define the directional derivative. However, since the neighbors are often not balanced around the boundary node, such selection might produce first approximations which could cause loss of order in the Motz problem with the grids considered here.

In this paper, we define an ''ad hoc'' ghost point to avoid loss of accuracy. The boundary condition

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \frac{\partial u}{\partial n}(p_0)=l_x\frac{\partial u}{\partial x}+l_y\frac{\partial u}{\partial y} </math>
|}
|}

where <math display="inline">n=(l_x,l_y)</math> is the outer normal vector at <math display="inline">p_0</math>, is approximated using the scheme of equation ([[#eq-5|5]]) with <math display="inline">A=B=E=0</math> and <math display="inline">C=l_x</math>, <math display="inline">D=l_y</math> using 7 points: a ghost point <math display="inline">p_g</math>, and the nodes <math display="inline">p_0</math>, <math display="inline">p_1</math>, <math display="inline">p_2</math>, <math display="inline">p_3</math>, <math display="inline">p_4</math> and <math display="inline">p_5</math> (Figure [[#img-4|4]]), which yields the equation ([[#eq-12|12]])

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>l_x\frac{\partial u}{\partial x}+l_y\frac{\partial u}{\partial y}\approx \hat{\Gamma }_g u(p_g)+\sum _{l=0}^k\hat{\Gamma }_l u(p_l),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

or, in other words,

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u(p_g)=\frac{1}{\hat{\Gamma }_g}\frac{\partial u}{\partial n}(p_0)-\sum _{l=0}^k\frac{\hat{\Gamma }_l}{{\hat{\Gamma }_g}} u(p_l). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

For convenience, the ghost point <math display="inline">p_g</math> is calculated in such a way that <math display="inline">p_0</math> is the midpoint of <math display="inline">p_g</math> and <math display="inline">p_c</math>, the latter being the centroid of the nodes <math display="inline">p_0</math>, <math display="inline">p_1</math>, <math display="inline">p_2</math>, <math display="inline">p_3</math>, <math display="inline">p_4</math>, and <math display="inline">p_5</math> (Figure [[#img-4|4]]).

The value of <math display="inline">u(p_g)</math> given by the equation ([[#eq-13|13]]) is substituted in the approximation to the Laplacian (equation ([[#eq-7|7]]) at <math display="inline">p_0</math> using the same 7 points <math display="inline">p_g</math>, <math display="inline">p_0</math>, <math display="inline">p_1</math>, <math display="inline">p_2</math>, <math display="inline">p_3</math>, <math display="inline">p_4</math>, and <math display="inline">p_5</math>, which eliminates the value of <math display="inline">u(p_g)</math>.

Te assembled system for the whole grid becomes sparse as the number of grid points per side increases; it was solved using a sparse lu decomposition to produces the approximation to the solution of the problem.  

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 45%;"
|-
|style="padding:10px;"| [[File:Chavez_et_al_2020a_2003_498px-Draft_Dominguez_Mota_636399092-ghost.png]]
<!--style="padding:10px;"| [[Image:Draft_Dominguez Mota_636399092-ghost.png|498px|Selection of ghost point at a boundary node of a general structured grid with Neumann condition]]-->
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"|  '''Figure 4'''. Selection of ghost point at a boundary node of a general structured grid with Neumann condition
|}

==4 Numerical tests==

Several tests were performed  with <math display="inline">m</math> and <math display="inline">2m-1</math> nodes per side on the <math display="inline">y</math> and <math display="inline">x</math> directions, respectively, varying <math display="inline">m</math> from 31 to 351. Grids of the three types in the previous section were generated, and the Motz problem was solved using the proposed scheme for every grid size and distribution. An example of the finite difference discretization described in the previous section is sketched in Figure [[#img-5|5]].

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
|-
|style="padding:10px;"| [[Image:Draft_Dominguez Mota_636399092-malla2.png|396px|Grid and calculated solution]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"|  '''Figure 5'''. Grid and calculated solution
|}


To assess the accuracy of the numerical solution, since  the approximation was calculated using structured grids, the quadratic error

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\epsilon =\sqrt{\sum _{i,j}\left[u_{i,j}-u(p_{i,j}) \right]^2 A_{i,j}}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

was used, where <math display="inline">u_{i,j}</math> is the solution calculated using the difference scheme at the point <math display="inline">p_{i,j}</math>, <math display="inline">u</math> is the exact solution (equation [[#eq-1|(1)]]) at the same point, and <math display="inline">A_{i,j}</math> is the area of the rhombus defined by the nodes <math display="inline">p_{i+1,j}, p_{i-1,j}, p_{i,j+1}, p_{i,j-1}</math>.

The measure of the roughness of the grid is given by the meh norm <math display="inline">h</math>; the log graph of the  mean square errors versus the grid norm for the type I and II grids are shown in Figures [[#img-6|6]] and [[#img-7|7]].  One can note that there is a clear convergence tendency on all the tests despite the boundary singularity. But, when only a nonuniform distribution either along the <math display="inline">x</math> or the <math display="inline">y</math>-axis is used,  the second-order approximation is lost, in a similar way as reported by Bernal using radial basis functions. The best results are obtained with the grids which are non-uniform on <math display="inline">x</math> and <math display="inline">y</math>, as can be seen in Figure [[#img-7|7]]. 

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
|-
|style="padding:10px;"| [[Image:Draft_Dominguez Mota_636399092-figErrorX20.png|300px|]]
|style="padding:10px;"| [[Image:Draft_Dominguez Mota_636399092-figErrorY20.png|300px|Log values of the quadratic errors, nonuniform node distribution on x or y]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding:10px;"|  '''Figure 6'''. Log values of the quadratic errors, nonuniform node distribution on <math>x</math> or <math>y</math>
|}


<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
|-
|style="padding:10px;"|[[Image:Draft_Dominguez Mota_636399092-figErrorXY20.png|500px|Log values of the quadratic errors, nonuniform node distribution on x and y]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figure 7'''. Log values of the quadratic errors, nonuniform node distribution on <math>x</math> and <math>y</math>
|}


The values of the errors, grid norms, and estimated convergence order are shown in Tables  [[#table-1|1]] to [[#table-3|3]]. The values in the last table show that, indeed, quadratic convergence is obtained when the same nonuniform node distribution along  both axes is used, even in the presence of the boundary singularity.

<div class="center" style="font-size: 75%;">'''Table 1'''. Quadratic errors and approximation order for type I grids</div>

<div id='table-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! <math display="inline">\alpha </math> !!  <math>m</math> !! <math>n</math> !!  <math>h</math> !! <math>\epsilon </math> !!  Order 
|- style="text-align:center"
|  1.00 
| 31 
| 61 
| 3.3333E-02 
| 6.4941E-03 
| &#8211;&#8211; 
|-style="text-align:center"
| 1.00 
| 71 
| 141 
| 1.4286E-02 
| 2.7556E-03 
| 1.01 
|-style="text-align:center"
| 1.00 
| 111 
| 221 
| 9.0909E-03 
| 1.7489E-03 
| 1.01 
|-style="text-align:center"
| 1.00 
| 151 
| 301 
| 6.6667E-03 
| 1.2810E-03 
| 1.00 
|-style="text-align:center"
| 1.00 
| 191 
| 381 
| 5.2632E-03 
| 1.0106E-03 
| 1.00 
|-style="text-align:center"
| 1.00 
| 231 
| 461 
| 4.3478E-03 
| 8.3446E-04 
| 1.00 
|-style="text-align:center"
| 1.00 
| 271 
| 541 
| 3.7037E-03 
| 7.1060E-04 
| 1.00 
|-style="text-align:center"
| 1.00 
| 311 
| 621 
| 3.2258E-03 
| 6.1877E-04 
| 1.00 
|-style="text-align:center"
| 1.25 
| 31 
| 61 
| 4.1492E-02 
| 6.4057E-03 
| &#8211;&#8211; 
|-style="text-align:center"
| 1.25 
| 71 
| 141 
| 1.7825E-02 
| 2.7333E-03 
| 1.01 
|-style="text-align:center"
| 1.25 
| 111 
| 221 
| 1.1351E-02 
| 1.7503E-03 
| 0.99 
|-style="text-align:center"
| 1.25 
| 151 
| 301 
| 8.3264E-03 
| 1.2920E-03 
| 0.98 
|-style="text-align:center"
| 1.25 
| 191 
| 381 
| 6.5746E-03 
| 1.0260E-03 
| 0.98 
|-style="text-align:center"
| 1.25 
| 231 
| 461 
| 5.4318E-03 
| 8.5207E-04 
| 0.97 
|-style="text-align:center"
| 1.25 
| 271 
| 541 
| 4.6275E-03 
| 7.2925E-04 
| 0.97 
|- style="text-align:center"
| 1.25 
| 311 
| 621 
| 4.0306E-03 
| 6.3784E-04 
| 0.97 
|}


<div class="center" style="font-size: 75%;">'''Table 2'''. Quadratic errors and approximation order for type II grids</div>

<div id='table-2'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! <math display="inline">\alpha </math> !!  <math>m</math> !! <math>n</math> !!  <math>h</math> !! <math>\epsilon </math> !!  Order 
|- style="text-align:center"
|  1.00 
| 31 
| 61 
| 0.033333 
| 0.0064941 
| &#8211;&#8211; 
|-style="text-align:center"
| 1.00 
| 71 
| 141 
| 0.014286 
| 0.0027556 
| 1.01 
|-style="text-align:center"
| 1.00 
| 111 
| 221 
| 0.0090909 
| 0.0017489 
| 1.01 
|-style="text-align:center"
| 1.00 
| 151 
| 301 
| 0.0066667 
| 0.001281 
| 1.00 
|-style="text-align:center"
| 1.00 
| 191 
| 381 
| 0.0052632 
| 0.0010106 
| 1.00 
|-style="text-align:center"
| 1.00 
| 231 
| 461 
| 0.0043478 
| 0.00083446 
| 1.00 
|-style="text-align:center"
| 1.00 
| 271 
| 541 
| 0.0037037 
| 0.0007106 
| 1.00 
|-style="text-align:center"
| 1.00 
| 311 
| 621 
| 0.0032258 
| 0.00061877 
| 1.00 
|-style="text-align:center"
| 1.25 
| 31 
| 61 
| 0.041492 
| 0.0060306 
| &#8211;&#8211; 
|-style="text-align:center"
| 1.25 
| 71 
| 141 
| 0.017825 
| 0.0026189 
| 0.99 
|-style="text-align:center"
| 1.25 
| 111 
| 221 
| 0.011351 
| 0.0016827 
| 0.98 
|-style="text-align:center"
| 1.25 
| 151 
| 301 
| 0.0083264 
| 0.001243 
| 0.98 
|-style="text-align:center"
| 1.25 
| 191 
| 381 
| 0.0065746 
| 0.00098714 
| 0.98 
|-style="text-align:center"
| 1.25 
| 231 
| 461 
| 0.0054318 
| 0.00081951 
| 0.97 
|-style="text-align:center"
| 1.25 
| 271 
| 541 
| 0.0046275 
| 0.00070109 
| 0.97 
|- style="text-align:center"
| 1.25 
| 311 
| 621 
| 0.0040306 
| 0.00061294 
| 0.97 
|}


<div class="center" style="font-size: 75%;">'''Table 3'''. Quadratic errors and approximation order for type III grids</div>

<div id='table-2'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! <math display="inline">\alpha </math> !!  <math>m</math> !! <math>n</math> !!  <math>h</math> !! <math>\epsilon </math> !!  Order 
|- style="text-align:center"
|  2.00 
| 71 
| 141 
| 2.8367E-02 
| 9.3949E-04 
| 1.77 
|-style="text-align:center"
| 2.00 
| 111 
| 221 
| 1.8099E-02 
| 4.1865E-04 
| 1.80 
|-style="text-align:center"
| 2.00 
| 151 
| 301 
| 1.3289E-02 
| 2.3926E-04 
| 1.81 
|-style="text-align:center"
| 2.00 
| 191 
| 381 
| 1.0499E-02 
| 1.5584E-04 
| 1.82 
|-style="text-align:center"
| 2.00 
| 231 
| 461 
| 8.6767E-03 
| 1.1006E-04 
| 1.82 
|-style="text-align:center"
| 2.00 
| 271 
| 541 
| 7.3937E-03 
| 8.2131E-05 
| 1.83 
|-style="text-align:center"
| 2.00 
| 311 
| 621 
| 6.4412E-03 
| 6.3786E-05 
| 1.83 
|-style="text-align:center"
| 2.25 
| 31 
| 61 
| 7.3442E-02 
| 4.5265E-03 
| &#8211;&#8211; 
|-style="text-align:center"
| 2.25 
| 71 
| 141 
| 3.1856E-02 
| 9.2020E-04 
| 1.91 
|-style="text-align:center"
| 2.25 
| 111 
| 221 
| 2.0338E-02 
| 3.8741E-04 
| 1.93 
|-style="text-align:center"
| 2.25 
| 151 
| 301 
| 1.4938E-02 
| 2.1315E-04 
| 1.94 
|-style="text-align:center"
| 2.25 
| 191 
| 381 
| 1.1803E-02 
| 1.3496E-04 
| 1.94 
|-style="text-align:center"
| 2.25 
| 231 
| 461 
| 9.7560E-03 
| 9.3190E-05 
| 1.94 
|-style="text-align:center"
| 2.25 
| 271 
| 541 
| 8.3140E-03 
| 6.8259E-05 
| 1.95 
|-style="text-align:center"
| 2.25 
| 311 
| 621 
| 7.2434E-03 
| 5.2179E-05 
| 1.95 
|-style="text-align:center"
| 2.50 
| 31 
| 61 
| 8.1262E-02 
| 4.9361E-03 
| &#8211;&#8211; 
|-style="text-align:center"
| 2.50 
| 71 
| 141 
| 3.5333E-02 
| 9.2007E-04 
| 2.02 
|-style="text-align:center"
| 2.50 
| 111 
| 221 
| 2.2573E-02 
| 3.7194E-04 
| 2.02 
|-style="text-align:center"
| 2.50 
| 151 
| 301 
| 1.6583E-02 
| 1.9951E-04 
| 2.02 
|-style="text-align:center"
| 2.50 
| 191 
| 381 
| 1.3106E-02 
| 1.2407E-04 
| 2.02 
|-style="text-align:center"
| 2.50 
| 231 
| 461 
| 1.0834E-02 
| 8.4506E-05 
| 2.02 
|-style="text-align:center"
| 2.50 
| 271 
| 541 
| 9.2336E-03 
| 6.1224E-05 
| 2.02 
|-style="text-align:center"
| 2.50 
| 311 
| 621 
| 8.0450E-03 
| 4.6380E-05 
| 2.02 
|-style="text-align:center"
| 2.75 
| 31 
| 61 
| 8.9015E-02 
| 5.3606E-03 
| &#8211;&#8211; 
|-style="text-align:center"
| 2.75 
| 71 
| 141 
| 3.8796E-02 
| 9.3491E-04 
| 2.10 
|-style="text-align:center"
| 2.75 
| 111 
| 221 
| 2.4802E-02 
| 3.6799E-04 
| 2.08 
|-style="text-align:center"
| 2.75 
| 151 
| 301 
| 1.8227E-02 
| 1.9446E-04 
| 2.07 
|-style="text-align:center"
| 2.75 
| 191 
| 381 
| 1.4407E-02 
| 1.1975E-04 
| 2.06 
|-style="text-align:center"
| 2.75 
| 231 
| 461 
| 1.1911E-02 
| 8.0998E-05 
| 2.06 
|-style="text-align:center"
| 2.75 
| 271 
| 541 
| 1.0152E-02 
| 5.8375E-05 
| 2.05 
|-style="text-align:center"
| 2.75 
| 311 
| 621 
| 8.8459E-03 
| 4.4042E-05 
| 2.05 
|-style="text-align:center"
| 3.00 
| 31 
| 61 
| 9.6704E-02 
| 5.7955E-03 
| &#8211;&#8211; 
|-style="text-align:center"
| 3.00 
| 71 
| 141 
| 4.2248E-02 
| 9.6029E-04 
| 2.17 
|-style="text-align:center"
| 3.00 
| 111 
| 221 
| 2.7026E-02 
| 3.7142E-04 
| 2.13 
|-style="text-align:center"
| 3.00 
| 151 
| 301 
| 1.9867E-02 
| 1.9454E-04 
| 2.10 
|-style="text-align:center"
| 3.00 
| 191 
| 381 
| 1.5707E-02 
| 1.1915E-04 
| 2.09 
|-style="text-align:center"
| 3.00 
| 231 
| 461 
| 1.2987E-02 
| 8.0299E-05 
| 2.08 
|-style="text-align:center"
| 3.00 
| 271 
| 541 
| 1.1070E-02 
| 5.7715E-05 
| 2.07 
|- style="text-align:center"
| 3.00 
| 311 
| 621 
| 9.6462E-03 
| 4.3453E-05 
| 2.06 
|}

==5. Results and discussion==

Several assertions can be drawn from the numerical results. First, despite the presence of the boundary singularity, the results are accurate enough. This is certainly a strength of the proposed finite difference scheme, which is rather simple. Since the main differences between the scheme discussed in this paper and the standard five-point scheme for the Laplacian are the addition of two extra “diagonal” nodes in the inner grid nodes and the use of a suitable grid spacing, this suggests that the proposed stencils are suitable to produce accurate results.

The convergence orders in table [[#table-3|3]] must be highlighted. Having in mind the presence of the singularity and the simplicity of the addressed difference stencils, this is indeed a remarkable result. Since the  boundary node redistribution was defined by a very simple rule,  this suggests that it is unnecessary to consider another kind of grids or more complex refinements near the singularity to solve the Motz problem, although, it must be acknowledged that the mesh adapting strategy was proposed ''ad hoc'' for it. Naturally, the numerical treatment of different singularities in other problems require a different application of the equidistribution principle.

==6. Conclusions==

Indeed, generalized finite differences can be used to solve the Motz problem numerically. The numerical results show that generalized finite differences  produce accurate solutions to a complicated benchmark problem, and the most noteworthy fact is that this can be done even with very simple structured grids. This is an important feature of the GFDM method, considering the difficulty of the problem. Thinking about the solution to several linear and nonlinear problems with this method which appear in the literature <span id='citeF-27'></span><span id='citeF-10'></span><span id='citeF-22'></span><span id='citeF-13'></span>[[#cite-27|[27,10,22,13]]], it is certainly a reliable numerical tool. 

The results are very encouraging since boundary configurations similar to those of Motz problem appear in several relevant applications, for instance, in interfaces between different media in infiltration problems, or advection-diffusion problems with sources and sinks along the boundary. Both problems are currently under investigation and advances in this sense will be reported in a future paper.

==Acknowledgments==

We want to thank  CIC-UMSNH  9.16,  and Aula-CIMNE Morelia for the financial support for this work.

==References==

<div id="cite-1"></div>
[[#citeF-1|[1]]] Lu T. T., Hu H. Y., Li Z. C. Highly accurate solutions of Motz's and the cracked beam problems. Engineering Analysis with Boundary Elements, 28:1387–1403, 2004.

<div id="cite-2"></div>
[[#citeF-2|[2]]] Li Z. C., Mathon R., Sermer, P. Boundary methods for solving elliptic problems with singularities and interfaces. SIAM Journal of Numerical Analysis, 24: 487–498, 1987.

<div id="cite-3"></div>
[[#citeF-3|[3]]] Li Z. C. An approach for combining the Ritz-Galerkin and Finite-Element Methods. Journal of Approximation Theory, 39:132-152, 1983.

<div id="cite-4"></div>
[[#citeF-4|[4]]] Arad M., Yosibash Z., Ben-Dor G.,  Yakhot, A. Computing flux intensity factors by a boundary method for elliptic equations with singularities. Communications in Numerical Methods in Engineering, 14: 657–670, 1998.

<div id="cite-5"></div>
[[#citeF-5|[5]]] Georgiou G. C., Olson L., Smyrlis Y. S. A singular function boundary integral method for the Laplace equation. Communications in Numerical Methods in Engineering, 12:127-134, 1996.

<div id="cite-6"></div>
[[#citeF-6|[6]]] Bernal F., Kindelan, M. Radial basis function solution of the Motz problem. Engineering Computations, 27: 606–620, 2010.

<div id="cite-7"></div>
[[#citeF-7|[7]]] Kansa E. J. Multiquadrics, a scattered data approximation scheme with applications to computational fluid dynamics II: solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers and Mathematics with Applications, 19:147–161, 1990.

<div id="cite-8"></div>
[[#citeF-8|[8]]] Shu C., Ding H., Yeo K. S. Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering, 192: 941–954, 2003.

<div id="cite-9"></div>
[[#citeF-9|[9]]] Hong W.T. Enriched Meshfree Method for an Accurate Numerical Solution of the Motz Problem. Advances in Mathematical Physics, 2016:1-12, 2016.

<div id="cite-10"></div>
[[#citeF-10|[10]]] Domínguez-Mota F. J., Mendoza-Armenta S., Tinoco-Guerrero G., Tinoco-Ruiz J. G. Finite difference schemes satisfying an optimality condition for the unsteady heat equation. Mathematics and Computers in Simulation, 106:76–83, 2014.

<div id="cite-11"></div>
[[#citeF-11|[11]]] Chávez-Negrete C., Domínguez-Mota F.J., Santana-Quinteros D. Numerical solution of Richards' equation of water flow by generalized finite differences. Computers and Geotechnics, 101:168–175, 2018.

<div id="cite-12"></div>
[[#citeF-12|[12]]] Tinoco-Guerrero G., Domínguez-Mota F. J., Gaona-Arias A., Ruiz-Zavala M. L., Tinoco-Ruiz J. G. A stability analysis for a generalized finite-difference scheme applied to the pure advection equation. Mathematics and Computers in Simulation, 147:293–300, 2018.

<div id="cite-13"></div>
[[#citeF-13|[13]]] Benito J. J., Ureña F., Gavete L. Solving parabolic and hyperbolic equations by the generalized finite difference method. Journal of Computational and Applied Mathematics, 209:208–233, 2007.

<div id="cite-14"></div>
[[#citeF-14|[14]]] Benito J. J., Ureña F., Gavete, L. Application of the generalized finite difference method to solve the advection-diffusion equation. Journal of Computational and Applied Mathematics, 235:1849–1855, 2011.

<div id="cite-15"></div>
[[#citeF-15|[15]]] Benito J. J., Ureña F., Salete E., Gavete, L. A note on the application of the generalized finite difference method to seismic wave propagation in 2D. Journal of Computational and Applied Mathematics, 236:3016–3025, 2012.

<div id="cite-16"></div>
[[#citeF-16|[16]]] Gavete L., Ureña F., Benito J.J., García A., Ureña M., Salete E. Solving second order non-linear elliptic partial differential equations using generalized finite difference method. Journal of Computational and Applied Mathematics, 318:378–387, 2017.

<div id="cite-17"></div>
[[#citeF-17|[17]]] Fan C.M., Li P.W. Generalized finite difference method for solving two-dimensional Burgers' equations. Procedia Engineering, 79: 55–60, 2014.

<div id="cite-18"></div>
[[#citeF-18|[18]]] Chan H.F., Fan C.M., Kuo C.W. Generalized finite difference method for solving two-dimensional non-linear obstacle problem.  Engineering Analysis with Boundary Elements, 37:1189–1196, 2013.

<div id="cite-19"></div>
[[#citeF-19|[19]]] Li P.W., Fan C.M. Generalized finite difference method for two-dimensional shallow water equations. Engineering Analysis with Boundary Elements, 80:58–71, 2017.

<div id="cite-20"></div>
[[#citeF-20|[20]]] Li P.W. , Chen W., Fu Z.J., Fan C.M. Generalized finite difference method for solving the double-diffusive natural convection in fluid-saturated porous media. Engineering Analysis with Boundary Elements 95: 175–186, 2018.

<div id="cite-21"></div>
[[#citeF-21|[21]]] Liszka T., Orkisz J. The finite difference method at arbitrary irregular grids and its application in applied mechanics. Computers and Structures, 11:83–95, 1980.

<div id="cite-22"></div>
[[#citeF-22|[22]]] Benito J. J., Ureña F., Gavete L., Alonso B. A posteriori error estimator and indicator in generalized finite differences. Application to improve the approximated solution of elliptic PDEs. International Journal of Computer Mathematics, 85:359–370, 2008.

<div id="cite-23"></div>
[[#citeF-23|[23]]] Celia, M., Gray, W. Numerical Methods for Differential Equations: fundamental concepts for scientific and engineering applications. Prentice-Hall, 1992.

<div id="cite-24"></div>
[[#citeF-24|[24]]] Ivanenko, S. In Selected chapters on grid generation and applications, B. Azarenok (Ed.), Russian Academy of Sciences, 2009.

<div id="cite-25"></div>
[[#citeF-25|[25]]] Pereyra V., Sewell G. A singular function boundary integral method for the Laplace equation. Numer. Math., 23:261-268, 1975.

<div id="cite-26"></div>
[[#citeF-26|[26]]] Domínguez-Mota F.J., Mendoza S., Tinoco-Ruiz J.G, Tinoco-Guerrero G. Numerical Solution of Poisson-like Equations with Robin Boundary Conditions Using a Finite Difference Scheme Dfined by an Optimality Condition. Proceedings of MASCOT 2011, 17, 101-110, 2011.

<div id="cite-27"></div>
[[#citeF-27|[27]]] Domínguez-Mota F. J., Equihua M., Mendoza S.,Tinoco-Ruiz J. G. Solución de ecuaciones diferenciales elípticas en regiones planas irregulares usando mallas convexas generadas por métodos variacionales empleando elementos finitos. Rev. Int. Mét. Num. Cálc. Dis. Ing., 26:187–194, 2010.