<figure><img src="image_1.png" /></figure>

HAL
open science

# A boundary meshless method with shape functions computed from the PDE

Djédjé Sylvain Zézé, Michel Potier-Ferry, Noureddine Damil

► To cite this version:

Djédjé Sylvain Zézé, Michel Potier-Ferry, Noureddine Damil. A boundary meshless method with shape functions computed from the PDE. Engineering Analysis with Boundary Elements, 2010, 34 (8), pp.747-754. 10.1016/j.enganabound.2010.03.008. hal-01517951

HAL Id: hal-01517951

https://hal.science/hal-01517951

Submitted on 4 May 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# A boundary meshless method with shape functions computed from the PDE

Djédjé Sylvain Zézé<sup>a</sup>, Michel Potier-Ferry<sup>a</sup>, Noureddine Damil<sup>b</sup>

<sup>a</sup> Laboratoire de Physique et Mécanique des Matériaux, FRE CNRS 3236, Université Paul Verlaine-Metz, Ile du Saulcy, 57045 Metz Cedex 01, France

<sup>b</sup> Laboratoire de Calcul Scientifique en Mécanique, Faculté des Sciences Ben M'Sik, Université Hassan II Mohammeda - Casablanca, B.P. 7955, Sidi Othman, Casablanca, Morocco

This paper presents a new meshless method using high degree polynomial shape functions. These shape functions are approximated solutions of the partial differential equation (PDE) and the discretization concerns only the boundary. If the domain is split into several subdomains, one has also to discretize the interfaces. To get a true meshless integration-free method, the boundary and interface conditions are accounted by collocation procedures. It is well known that a pure collocation technique induces numerical instabilities. That is why the collocation will be coupled with the least-squares method. The numerical technique will be applied to various second order PDE's in 2D domains. Because there is no integration and the number of shape functions does not increase very much with the degree, high degree polynomials can be considered without a huge computational cost. As for instance the p-version of finite elements or some well established meshless methods, the present method permits to get very accurate solutions.

Keywords:
Perturbation technique
Least-squares method
Meshless
High degree polynomials
Boundary-only

## 1. Introduction

This paper presents a discretization technique to solve elliptic partial differential equations. It relies on high degree shape functions that are approximate solutions of the PDE. More precisely one requires that the residual is of order $O(|x-c|^{N+1})$ in the neighborhood of a given center c. This permits to discretize only the boundary conditions. These discretized equations could be obtained by two classical ways, first Galerkin procedure with the drawback of the integration cost, second point-collocation that is integration-less, but can lead to numerical instabilities and ill-conditioned matrices [7]. In this paper point-collocation method has been chosen and associated with a least-squares minimization to overcome the numerical instabilities, as proposed by Zhang et al. [18] and used by several others [19-24]. So the proposed method is characterised by closed-form solutions built by Taylor series, boundary discretization and coupling between collocation and least-squares minimization.

To our best knowledge, such a discretization principle has never been presented in the numerical literature, but of course it has some points in common with the many numerical methods that are not based on low degree polynomials. It can be compared with the p-version of finite element and one can hope that the presented technique permits to recover more or less the same accuracy and adaptivity as the p-version [1-3]. The differences lie in the number of shape functions that is much smaller with the

present method, in the computation cost and in the discretization principle. Our discrete problem is deduced from point-wise equations as in many meshless methods [4-6,8,25], but the present method does not use a priori given shape functions, they are built from a local solving of the PDE. There are at least three well known numerical methods that associate a family of exact solutions and a boundary discretization: the integral equation method [9-11], the method of fundamental solutions [12-14] and the scaled boundary finite element method [15-17]. In these three methods, the reference problem has to be linear with constant coefficients while the present Taylor series method can be extended to generic PDE's.

The paper is organised as follows. In the second part, the instabilities due to pure collocation are pointed out and compared with least-squares collocation. In Part 3, a computational technique is sketched that permits to apply Taylor series to PDE's. Finally in Part 4, various 2D applications are discussed to assess the possibilities of the presented numerical method.

## 2. Boundary collocation versus boundary least-squares collocation

### 2.1. Polynomial shape functions

Let us consider the Dirichlet problem in a 2D domain:

$$ \begin{cases} \Delta u = 0 & \text{in } \Omega \\ u(x) = u^d(x) & \text{on } \partial\Omega \end{cases} \tag{1} $$
The aim is to introduce high degree polynomial shape functions, that are exact solutions of the considered partial differential equations. In the case of the Laplace equation $\Delta u = 0$, there are only two such polynomials of degree $n$, namely $Re(x+iy)^n$ and $Im(x+iy)^n$. Next, let us introduce all the polynomials, whose degree is lower or equal to $p$. The dimension of this vectorial space is $((p+1)(p+2))/2$, but if one limits to the solutions of the Laplace equation, this dimension is reduced to $2p+1$. Note that the limitation to the solutions of the PDE permits to reduce strongly the number of shape functions, for instance 101 polynomials instead of 1326 for a degree equal to 50. With this reduction of the number of functions, one can hope to build a numerical method that remains efficient with a large degree.

## 2.2. Boundary collocation

As for instance in the method of fundamental solutions (MFS), exact solutions of the PDE are used. Hence it is only necessary to discretize the boundary conditions. Hence, the cloud of collocation points is located on the boundary. The simplest technique is to choose as many collocation points $\underline{x}_i$ as shape functions $P_i(\underline{x}_i)$. The unknown is written in the classical form as

$$u(x) = \sum_{i=1}^{2p+1} P_i(x)v_i \quad (2)$$

and the discretized equations are

$$\sum_{i=1}^{2p+1} P_i(\underline{x}_j)\underline{v}_i = u^d(\underline{x}_j), \quad 1 \le j \le 2p+1 \quad (3)$$

Let us apply, this simple boundary collocation to a unit disk $x^2+y^2 \le 1$, and with the boundary data

$$u^d(x,y) = \frac{x-x_0}{(x-x_0)^2+(y-y_0)^2} \quad (4)$$

The exact solution is known: $u_{ex}(x,y) = (x-x_0)/((x-x_0)^2+(y-y_0)^2)$.

First, one chooses a uniformly distributed cloud (see Fig. 2). In Fig. 1, we have plotted the error $|(u(x,0)-u_{ex}(x,0))/u_{ex}(x,0)|$ along the horizontal axis for three values of the degree. In this case, the boundary collocation method converges with the order $p$, see Fig. 1. For instance, for $p=32$, the maximal error is about $10^{-3}$ and the error in the center of the disk is about $10^{-6}$. The same conclusion holds also by looking at the error anywhere in the domain. Unfortunately, this simple collocation technique is not robust and it does not work with an irregular cloud. For instance,

<figure><img src="image_2.png" /><figcaption>Fig. 1. Dirichlet problem in a disk. Pure boundary collocation with a uniform cloud, see Fig. 2. Error along the horizontal axis.</figcaption></figure>

for $p=10$ and the collocation points of Fig. 3, the maximal value of the approximated solutions by this boundary collocation technique is about 1075, instead of 4 for the exact one. It is not surprising that this simple collocation technique does not work. Indeed in the present example of a disk, the boundary value of the polynomial is given by a truncated Fourier series and the coefficients $v_i$ are identical to the Fourier coefficients. In the present technique, one tries to identify the Fourier coefficients from pointwise data and with about two points per period $2\pi/p$, which is not sufficient for a stable estimate. Theoretically, the Fourier coefficients are given by integral formulae:

$$\left\{ \begin{array}{l} \frac{1}{\pi} \int_{0}^{2\pi} u(\theta) \cos n\theta \, d\theta \\ \\ \frac{1}{\pi} \int_{0}^{2\pi} u(\theta) \sin n\theta \, d\theta \end{array} \right. \quad (5)$$

To avoid the numerical evaluation of these integrals (5) that involves many integration points, we shall propose to identify these Fourier series from a number of pointwise data that is larger than $2p+1$ (Fig. 2).

## 2.3. Boundary least-squares collocation

It is proposed to identify the coefficients $v_i$ of the polynomial (2) from M collocation points, M being larger than $2p+1$. The Dirichlet boundary condition will be satisfied in a least-square sense. Such a least-square collocation method has been presented by Zhang et al. [18] in another meshless framework and it has been widely applied. One requires that the coefficients $v_i$ minimize the function

$$J(v_i) = \frac{1}{2} \sum_{j=1}^{M} |u(\underline{x}_j) - u^d(\underline{x}_j)|^2 = \frac{1}{2} \sum_{j=1}^{M} \left| \sum_{i=1}^{2p+1} P_i(\underline{x}_j) v_i - u^d(\underline{x}_j) \right|^2 \quad (6)$$

After few calculations, this minimization of (6) leads to a linear system

$$[K]\{v\} = \{b\} \quad (7)$$

where

$$[K] = \sum_{j=1}^{M} [K^j] \quad \text{with} \quad K_{ik}^{j} = P_i(\underline{x}_j) P_k(\underline{x}_j) \quad (8a)$$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Uniform distribution of collocation points, $x_0=1.2$, $y_0=0.3$.</figcaption></figure>
$$ \{b\} = \sum_{j=1}^{M} \{b^j\} \quad \text{with} \quad b_i^j = u^d(\underline{x}_j)P_i(\underline{x}_j) \qquad (8b) $$

This method has been applied to the same problem as in paragraph 2.2, with the same order $p=10$, and with a non-uniform cloud of 41 collocation points, see Fig. 3. In this case and contrarily to the simple collocation technique, the polynomial solution is close to the exact one, see Fig. 4. There are small oscillations in the polynomial approximation, but they can be removed by increasing the order $p$. The convergence of the method with the order is shown in Fig. 5. One sees that the error decreases strongly more or less as in the p-version of finite elements. For instance the error at the boundary becomes lower than $10^{-6}$ for $p$ larger than 30 and the maximal error is about 0.1%. In Table 1, the p-convergence of the method is presented in the case of uniform and a non-uniform cloud. Hence the coupling of collocation and least square permits to get the same accuracy with the non-uniform cloud as the uniform one. The accuracy is represented by the decimal logarithm of the maximal error $|(u(x,y)-u_{ex}(x,y))/u_{ex}(x,y)|$. Finally let us discuss the influence of the number of boundary collocation points, in the discretization error, see Table 2 in the case of a degree 18. One checks again that the pure collocation ($M=37$) leads to wrong results. From a number of points $M=61$, the error becomes constant, and in this case $p=18$, about equal to $10^{-1.6}$. This means that it is not necessary to choose too many collocation points, $M \approx 4p$ appearing to be a good compromise.

### 3. Algorithm to compute the shape functions

It is not very difficult to build some polynomial shape functions for a given partial differential equation. In the generic case, this construction could be done with the help of automatic differentiation [26], as it is well known in the case of ordinary differential equation [27]. In this paper, we limit ourselves to the case of a linear equation with constant coefficients and a general right hand side $f(x,y)$:

$$ -\Delta u + u = f(x,y) \qquad (9) $$

The principle is to consider the PDE as an ordinary differential equation with respect to $y$, the function and its gradient being given on a straight line, for instance the line $x=0$. The data and

<figure><img src="image_4.png" /><figcaption>Fig. 3. Non-uniform distribution of collocation points, $x_0 = 1.2$, $y_0 = 0.3$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Dirichlet problem in a disk. Boundary least-squares collocation method with a non-uniform cloud, see Fig. 3. Comparison of the exact solution and the proposed one with $p=10$ along the boundary.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. p-convergence for the Dirichlet problem in a disk. Boundary least-squares collocation method with a non-uniform cloud, see Fig. 3.</figcaption></figure>

the unknown are assumed to be polynomials:

$$ u(x,y) = \sum_{i=1}^{p} \sum_{j=1}^{p} u_{ij} x^i y^j \qquad (10) $$

For convenience, each polynomial is split into homogeneous polynomials. A homogeneous polynomial of degree $k$ combines $k+1$ monoms and it is written in a matricial form

$$ \sum_{i=0}^{k} u_{k-i, i} x^i y^{k-i} = \langle X^k \rangle \{u^k\} \qquad (11) $$

where

$$ \langle X^k \rangle = \langle x^k, x^{k-1} y, \dots, y^k \rangle, \quad \langle u^k \rangle = \langle u_{k,0}, u_{k-1,1}, \dots, u_{0,k} \rangle $$

Thus all the polynomials whose degree is lower or equal to $p$ can be written in the form

$$ u(x,y) = \sum_{k=0}^{p} \langle X^k \rangle \{u^k\}, \quad \{u^k\} \in \mathbb{R}^{k+1} \qquad (12) $$

Next, if $u(x,y)$ is a homogeneous polynomial of degree $k+2$, $\Delta u$ is also a homogeneous polynomial of degree $k$. Hence, there exists a matrix $[L^k]$ with $k+1$ rows and $k+3$ columns such that

$$ \{\Delta u^k\} = [L^k]\{u^{k+2}\} \qquad (13) $$
<table><caption>**Table 1**
p-convergence of the least-squares collocation technique.</caption><thead><tr><th>Degree</th><th>NDOF</th><th>Number of collocation points</th><th>Uniform cloud</th><th>Non-uniform cloud</th></tr></thead><tbody><tr><td>6</td><td>13</td><td>26</td><td>-0.4097</td><td>-0.4551</td></tr><tr><td>8</td><td>17</td><td>40</td><td>-0.5782</td><td>-0.6078</td></tr><tr><td>14</td><td>29</td><td>60</td><td>-1.4869</td><td>-1.3818</td></tr><tr><td>18</td><td>37</td><td>80</td><td>-1.5963</td><td>-1.6019</td></tr><tr><td>25</td><td>51</td><td>105</td><td>-2.4078</td><td>-2.3333</td></tr><tr><td>30</td><td>61</td><td>125</td><td>-2.7822</td><td>-2.7824</td></tr></tbody></table>

Laplace equation with Dirichlet data (4).

<table><caption>**Table 2**
How to choose the number of collocation points.</caption><thead><tr><th>Degree</th><th>NDOF</th><th>Number of collocation points</th><th>Maximal error (log<sub>10</sub>)</th></tr></thead><tbody><tr><td rowspan="6">18</td><td rowspan="6">37</td><td>37</td><td>3.9456</td></tr><tr><td>49</td><td>-0.6641</td></tr><tr><td>61</td><td>-1.6076</td></tr><tr><td>73</td><td>-1.6052</td></tr><tr><td>85</td><td>-1.6046</td></tr><tr><td>130</td><td>-1.6041</td></tr><tr><td></td><td>190</td><td>-1.6040</td></tr></tbody></table>

Laplace equation with Dirichlet data (4).

So, by the formula (13), one is able to compute the $k+1$ coefficients of the homogeneous polynomial from the $k+3$ coefficients of the homogeneous polynomial $u(x,y)$. The matrix $[L^k]$ is defined in appendix (26). To establish shape functions that are approximate solutions of (9) one assumes that $u(x,0) = \sum_{i=0}^{p} u_{i,0}x^i$ and $(\partial u/\partial y)(x,0) = \sum_{i=0}^{p-1} u_{i,1}x^i$ are given. If $u(x,y)$ is a homogeneous polynomial of degree $k$, the function and its first derivative along the line $y=0$ is given by the two coefficients:

$$ \{v^k\} = \begin{Bmatrix} u_{k,0} \\ u_{k-1,1} \end{Bmatrix} \in \mathbb{R}^2 \qquad (14) $$

The other coefficients of the homogeneous polynomial are collected in the following vector:

$$ \{w^k\} = \begin{Bmatrix} u_{k-2,2} \\ u_{k-3,3} \\ \vdots \\ u_{0,k} \end{Bmatrix} \in \mathbb{R}^{k-1} \qquad (15) $$

Of course there exist two matrices $[S^{vk}]$ and $[S^{wk}]$ such that

$$ \{u^k\} = [S^{vk}][v^k] + [S^{wk}][w^k] \qquad (16) $$

The principle is to assume that the Taylor coefficients of $\Delta u - u + f$ are zero, up to the degree $p$. With account of (13), this leads to the following equation:

$$ [L^k]\{u^{k+2}\} = \{u^k\} - \{f^k\} \qquad (17) $$

The latter formula can be considered as a recurrence one, because it yields $\{u^{k+2}\}$ as a function of $\{u^k\}$, of the data $\{f^k\}$ and of the first two coefficients $\{v^{k+2}\}$ of $\{u^{k+2}\}$. Indeed by inserting (16) into (17), we get

$$ [L^k][S^{wk+2}][w^{k+2}] + [L^k][S^{vk+2}][v^{k+2}] = \{u^k\} - \{f^k\} \qquad (18) $$

One checks easily that the square matrix $[L^k][S^{wk+2}]$, that is given in the appendix, is invertible. The latter result is consistent because (17) comes from the solution of a second order differential equation with the two initial data $\{v^k\}$. This yields

$2p+2$ shape functions $P_i(x,y)$ such that any polynomial solution of (9) in an asymptotic sense can be written in the form

$$ u(x,y) = P_0(x,y) + \sum_{i=1}^{2p+1} P_i(x,y)v_i \qquad (19) $$

where the numbers $v_i$ are defined by

$$ \{v\} = \begin{Bmatrix} v_1 \\ v_2 \\ \vdots \\ v_{2p+1} \end{Bmatrix} = \begin{Bmatrix} \{v^0\} \\ \{v^1\} \\ \vdots \\ \{v^p\} \end{Bmatrix} = \begin{Bmatrix} u_{0,0} \\ u_{1,0} \\ u_{1,1} \\ \vdots \\ u_{p,0} \\ u_{p,1} \end{Bmatrix} \in \mathbb{R}^{2p+1} \qquad (20) $$

# 4. Numerical applications

In this part, few applications will be discussed to assess the validity of the presented techniques of least-squares collocation with polynomial shape functions. A first goal is to show that the calculation of shape function presented in the previous part can be effective. The second goal is to discuss briefly some applications with a splitting of the domain into several subdomains and a polynomial approximation in each subdomain.

## 4.1. Amoeba-like domain

Let us consider the modified Helmholtz equation

$$ \begin{cases} -\Delta u + u = e^x \cos y - e^{x-y} & \text{in } \Omega \\ u(x,y) = e^{x-y} + e^x \cos y & \text{in } \partial\Omega \end{cases} \qquad (21) $$

where $\Omega$ has the amoeba-like boundary (see Fig. 6).

$$ \begin{cases} (x,y) = R(s)(\cos(s),\sin(s)) \\ R(s) = e^{\sin(s)}\sin^2(2s) + e^{\cos(s)}\cos^2(2s), \quad 0 \le s \le 2\pi \end{cases} \qquad (22) $$

This problem has been solved, in [13], by using the method of the fundamental solution coupled with the "quasi-Monte Carlo method". In that paper, 4000 collocation points were used inside the domain. At the boundary, 35 collocation points corresponding to 35 equations were used in [13] to get a maximal error of about $2 \times 10^{-3}$. Here, 19 polynomials of degree 9 and 31 collocation points were used, to get the accuracy presented in Fig. 7, with a maximal error $< 10^{-4}$. This seems much more efficient than the method of [13]. Likely this good behavior is due to the analytical solution inside the domain and also to the smoothness of the sought function that can be easily approximated by polynomials. This establishes that the considered polynomials are able to solve very accurately a PDE in a complex domain.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Amoeba-like domain and collocation points.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Helmholtz equation (21) in the amoeba-like domain. Degree 9. Error along the horizontal axis y=0.</figcaption></figure>

## 4.2. A first piecewise study

In this part, we consider again the Laplace equation (1) in the unit disk with the Dirichlet data (4) and $x_0=1.2$, $y_0=0.3$. The domain will be split into two parts, as in Fig. 8. A small subdomain has been chosen to account for the rapid variation of the solution close to the singularity $(x_0, y_0)$. The parameters of the discretization are the following ones:

* $N_1(N_2)$: degrees of the shape functions in Zone 1 (Zone 2);

* $M_1(M_2)$: number of collocation points at the external boundary of Zone 1 (Zone 2);

* $I$: number of collocation points along the interface;

* $C_1(C_2)$: center of the Taylor expansion in Zone 1 (Zone 2).

The discretization principle is the same as in the case of a single subdomain, with an account of the continuity condition between

<figure><img src="image_9.png" /><figcaption>Fig. 8. Dirichlet problem in a disk. A first splitting in two subdomains.</figcaption></figure>

the two polynomial approximations at the interface. The conditions to be satisfied are the Dirichlet boundary conditions and the continuity of the unknown and of its normal derivative at the interface. More precisely, we are minimizing the following functional $J(v_1, v_2)$:

$$J(v_1, v_2) = \frac{1}{2} \sum_{j=1}^{M_1} \left| P_{10}(x_j) + \sum_{i=1}^{2N_1+1} P_{1i}(x_j) v_{1i} - u^d(x_j) \right|^2 \\ + \frac{1}{2} \sum_{t=1}^{I} \left| P_{10}(x_t) - P_{20}(x_t) + \sum_{i=1}^{2N_1+1} P_{1i}(x_t) v_{1i} - \sum_{k=1}^{2N_2+1} P_{2k}(x_t) v_{2k} \right|^2 \\ + \frac{1}{2} \sum_{t=1}^{I} \left| \frac{\partial P_{10}}{\partial n}(x_t) - \frac{\partial P_{20}}{\partial n}(x_t) + \sum_{i=1}^{2N_1+1} \frac{\partial P_{1i}}{\partial n}(x_t) v_{1i} - \sum_{k=1}^{2N_2+1} \frac{\partial P_{2k}}{\partial n}(x_t) v_{2k} \right|^2 \\ + \frac{1}{2} \sum_{r=1}^{M_2} \left| P_{20}(x_r) + \sum_{k=1}^{2N_2+1} P_{2k}(x_r) v_{2k} - u^d(x_r) \right|^2 \quad (23)$$

The p-convergence of the problem with the two subdomains is presented in Table 3. Clearly, the process converges and it yields accurate solutions for degrees larger than 15. The accuracy obtained with a degree 25 ($10^{-4.55}$) is better than with a single polynom ($10^{-2.33}$) according to Table 1. Other calculations have been done by keeping a smaller degree $N_2=6$ in the larger zone (zone 2) that is away from the singularity. In this case, the error is never smaller than $10^{-1.5}$. This means that a good convergence seems to require a high degree in the two subdomains. Two other splittings in two subdomains have been considered, where the interface is a straight line, see Figs. 9 and 10. The p-convergence is similar as with the previous case: the error is about $10^{-3}$ in the first case and $10^{-4}$ in the second one. All these results establish the robustness of the method, since the p-convergence is always rapid whatever be the choice of the subdomains.

## 4.3. Helmholtz equation

Now we consider the Helmholtz equation in a rectangular domain, $-5 \le x \le 5$, $0 \le y \le 4$

$$\begin{cases} -\Delta u + u = 0 \text{ in } \Omega \\ u(x,0) = u(x,4) = 0 \\ u(\pm 5, y) = \sin(\frac{\pi}{4}y) \end{cases} \quad (24)$$

The exact solution of (24) is given by

$$u(x,y) = \frac{\cosh\left(x\sqrt{1 + \frac{\pi^2}{16}}\right)}{\cosh\left(5\sqrt{1 + \frac{\pi^2}{16}}\right)} \sin\left(\frac{\pi}{4}y\right) \quad (25)$$
<table><caption>**Table 3**
p-convergence of least-squares collocation technique.</caption><thead><tr><th>Degree N<sub>1</sub>-N<sub>2</sub></th><th>NDOF</th><th>Number of collocation points M<sub>1</sub>-I-M<sub>2</sub></th><th>Error along x axis (log<sub>10</sub>)</th></tr></thead><tbody><tr><td>6-6</td><td>26</td><td>20-20-20</td><td>-0.6033</td></tr><tr><td>10-10</td><td>42</td><td>35-20-35</td><td>-1.7639</td></tr><tr><td>15-15</td><td>82</td><td>45-20-45</td><td>-2.0637</td></tr><tr><td>25-25</td><td>102</td><td>95-40-95</td><td>-4.5566</td></tr><tr><td>30-30</td><td>122</td><td>105-50-105</td><td>-4.2842</td></tr></tbody></table>

Laplace equations with Dirichlet data (4). Two subdomains, see Fig. 8.

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Dirichlet problem in a disk. A second splitting in two subdomains.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Dirichlet problem in a disk. A third splitting in two subdomains.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Helmholtz equation (24) in a rectangle. Exact and approximate solution along the mid-line y=2. A single subdomain, degree 12.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Five subdomains.</figcaption></figure>

This function varies rapidly near by the small sides $x = \pm 5$: so this problem involves boundary layers. First, we have tried to discretize the problem with a single subdomain, with the boundary least-squares collocation method. The numerical results show that the algorithm does not converge with the degree $N$, especially in the boundary layers. For instance the maximal error is $10^{-0.28}$ for $N=3$, $10^{-0.09}$ for $N=5$ and $10^{-0.13}$ for $N=10, 20, 25, 30$. Fig. 11 permits to compare the exact and the approximated solution.

To get a better approximation, several geometrical splittings have been considered. For instance, Fig. 12 represents a decomposition in five subdomains. We have chosen the same degree for all the polynomials. The accuracy of the approximate solution is presented in Table 4. Correct approximate solutions have been obtained especially for small degrees ($N=6$). We did not observe the same p-convergence property as with the previous example. Apparently a small error of about 1% persists close to the matching region between the central zone and the other ones. Likely for the same reason, we have not been able to improve this result by increasing the number of subdomains.

## 5. Conclusion

A new meshless method has been presented and evaluated from simple examples involving second order PDE's. The key point is the introduction of a small number of polynomial shape functions that are the approximate solutions of PDE. The principle is to require that the Taylor coefficients of the equation at a given point vanish. This leads to an algorithm to compute these polynomials. This is a true meshless method, the discrete equations being obtained by collocation or by a least-squares collocation technique. Indeed the pure collocation technique is not robust and leads sometimes to numerical instabilities. Because the equation in the domain is solved very accurately, the collocation points are located only on the boundary and also on interface, if one considers several subdomains. Thus the number of degrees of freedom does not increase too much with
<table><caption>**Table 4**
p-convergence of least-squares collocation technique.</caption><thead><tr><th>Degree N<sub>1</sub>-N<sub>2</sub>-N<sub>3</sub>-N<sub>4</sub> - N<sub>5</sub></th><th>NDOF</th><th>Number of collocation points M<sub>1</sub>-M<sub>2</sub>-M<sub>3</sub></th><th>Maximal error (log<sub>10</sub>)</th></tr></thead><tbody><tr><td>3-3-3-3-3</td><td>35</td><td>4-4-4</td><td>-1.4570</td></tr><tr><td>6-6-6-6-6</td><td>65</td><td>8-8-8</td><td>-2.7181</td></tr><tr><td>9-9-9-9-9</td><td>95</td><td>12-12-12</td><td>-2.4173</td></tr><tr><td>16-16-16-16-16</td><td>165</td><td>20-20-20</td><td>-2.3726</td></tr><tr><td>20-20-20-20-20</td><td>205</td><td>25-25-25</td><td>-2.3726</td></tr></tbody></table>

Helmoltz equation in a rectangle (24). Five subdomains, see Fig. 12.

the degree and this permits one to use very high degrees with a moderate cost. The results establish that the solution becomes very accurate, as for instance with the p-version of finite element method.

In this paper, the applications have been restricted to linear equations with constant coefficients. Nevertheless the computation of the shape functions comes from a Taylor series analysis. Hence it is possible to extend the present analysis to non-linear equations or to linear ones with variable coefficients. Such methods have been proposed recently to solve non-linear ordinary differential equations, see [28-30].

As far as we know, the presented method seems to be completely new. There are well known methods based on exact solutions (BEM, MFS), but they are mainly done for linear problems with constant coefficients, while the present method can be extended naturally to non-linear equations. Because of the high degree of the shape functions, it looks like the p-version of finite elements, but it should be much less expensive because the partial differential equation is solved analytically in an almost exact manner. This first paper on the topics has been focused on a basic algorithm to compute the shape functions, on the assessment of a proper collocation technique and a first check of the convergence property.

## Appendix A. Discretization of the Laplace operator

In this appendix, few formulae are presented to define the matrices that are useful to discretize the Laplace operator, see Section 3

$$[L^k] = [G_x^{k+1}] [G_x^{k+2}] + [G_y^{k+1}] [G_y^{k+2}] \quad (26)$$

where

$$[G_x^k] = \begin{bmatrix} k & 0 & 0 & \cdots & 0 & 0 \\ 0 & k-1 & 0 & \cdots & 0 & 0 \\ 0 & 0 & k-2 & & 0 & 0 \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ 0 & 0 & 0 & & 1 & 0 \end{bmatrix} \\ [G_y^k] = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & k-2 & 0 & 0 \\ 0 & 0 & \cdots & 0 & k-1 & 0 \\ 0 & 0 & \cdots & 0 & 0 & k \end{bmatrix} \quad (27)$$

$$[S^{vk}] = \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \end{bmatrix} = \begin{bmatrix} I_2 \\ 0_{k-1,2} \end{bmatrix} \quad \text{et} \quad [S^{wk}] = \begin{bmatrix} 0 & \cdots & \cdots & 0 \\ 0 & \cdots & \cdots & 0 \\ 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & 1 \end{bmatrix} = \begin{bmatrix} 0_{2,k-1} \\ I_{k-1} \end{bmatrix} \quad (28)$$

The matrix $[L^k] [S^{wk+2}]$ is exactly the same as $[L^k]$ where the two first columns have been eliminated.

## References

[1] Szabo BA, Babuska I. Finite element method. New York: Wiley classic library; 1991.

[2] Babuska I, Szabo BA, Katz LN. The p-version of finite element method. SIAM Journal of Numerical Analysis 1981;18:515-45.

[3] Schwab C, Suri M. The p and hp versions of the finite element method for problem with boundary layers. Mathematics of Computation 1996;65: 1403-29.

[4] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Computational Mechanics 1992;10:307-18.

[5] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid dynamics. I. Surface approximations and partial derivative estimates. Computers and Mathematics with Applications 1990;19:127-45.

[6] Golberg MA, Chen CS, Bowman H. Some recent results and proposals for the use of radial basis functions in the bem. Engineering Analysis with Boundary Elements 1999;23:285-96.

[7] Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multi-quadric radial basis functions: applications to elliptic partial differential equations. Computers and Mathematics with Applications 2000;39:123-37.

[8] Bouhamidi A, Jbilou K. Meshless thin plate spline methods for the modified Helmholtz equation. Computer Methods in Applied Mechanics and Engineering 2008;197:3733-41.

[9] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques. Theory and applications in engineering. Berlin and New York: Springer-Verlag; 1984.

[10] Bonnet M. Boundary integral equation methods for solids and fluids. New York: Wiley Publication; 1999.

[11] Wrobel LC, Aliabadi MH. The Boundary element method. New Jersey: Wiley; 2002.

[12] Golberg MA. The method of fundamental solutions for poisson's equation. Engineering Analysis with Boundary Elements 1995;16:205-13.

[13] Chen CS, Golberg MA, Hon YC. Numerical justification of fundamental solutions and the quasi-Monte Carlo method for poisson-type equations. Engineering Analysis with Boundary Elements 1998;22:61-9.

[14] Chen W, Tanaka M. A meshless, integration-free, and boundary-only rbf technique. Computers and Mathematics with Applications 2002;43:379-91.

[15] Song C, Wolf JP. The scaled boundary finite-element method—alias consistent infinitesimal finite-element cell method—for elastodynamics. Computer Methods in Applied Mechanics and Engineering 1997;147:329-55.

[16] Wolf JP, Song C. The scaled boundary finite-element method—a fundamental solutionless boundary-element method. Computer Methods in Applied Mechanics and Engineering 2001;190:5551-68.

[17] Song C, Bazyar MH. A boundary condition in Pade series for frequency-domain solution of wave propagation in unbounded domains. International Journal for Numerical Methods in Engineering 2007;69:2330-58.

[18] Zhang X, Liu XH, Song KZ, Lu MW. Least-squares collocation meshless method. International Journal for Numerical Methods in Engineering 2001;51:1089-100.

[19] Tan JY, Liu LH, Li BX. Least-squares radial point interpolation collocation meshless method for radiative heat transfer. Journal of Heat Transfer 2007;129:669-73.

[20] Wang QX, Li H, Lam KY. Development of a new meshless-point weighted least-squares (pwls) method for computational mechanics. Computational Mechanics 2005;35:170-81.

[21] Kwon KC, Youn SK. The least-squares meshfree method for rigid-plasticity with frictional contact. International Journal of Solids and Structures 2006;43:7450-81.

[22] Kwon KC, Park SH, Youn SK. The least-squares meshfree method for elasto-plasticity and its application to metal forming analysis. International Journal for Numerical Methods in Engineering 2005;64:751-88.

[23] Lee SH, Yoon YC. Meshfree point collocation method for elasticity and crack problems. International Journal for Numerical Methods in Engineering 2004;61:22-48.
[24] Ling L, Kansa EJ. A least-squares preconditioner for radial basis func-
tions collocation methods. Advances in Computational Mathematics 2005;23:
31–54.

[25] Onate E, Idelshon S, Zienkiewicz OC, Taylor RL. A finite point method in computational mechanics. Application to convective transport and fluid flow. International Journal for Numerical Methods in Engineering 1996;39: 3839–3866.

[26] Griewank A. Evaluating derivatives: principles and techniques of algorithmic differentiation frontiers in applied mathematics, vol. 19. Philadelphia: SIAM; 2000.

[27] Corliss G, Chang YF. Solving ordinary differential equations using Taylor series. ACM Transactions on Mathematical Software 2002;8:114–44.

[28] B. Dehbi, B. Braikat, H. Lahmam, M. Potier-Ferry, N. Damil, Une nouvelle méthode dite des éléments finis d'ordre élevé. In: 1er Congrès Tunisien de Mécanique, Hammamet, 2008.

[29] S. Zézé, M. Potier-Ferry, Calcul d'éléments finis par une technique de perturbation. In: 18ème Congrès Français de Mécanique, Grenoble du 27 au 31 août, <http://documents.irevues.inist.fr/bitstream/2042/16495/1>; 2007.

[30] D.S. Zézé, Calcul de fonctions de forme de haut degré par une technique de perturbation. PhD thesis, Université Paul Verlaine de Metz; 2009.