Modified Moving Least Squares Method for Two-Dimensional Linear and Nonlinear Systems of Integral Equations

Abstract

This work aims at focusing on modifying the moving least squares (MMLS) methods for solving two-dimensional linear and nonlinear systems of integral equations and system of differential equations. The modified shape function is our main aim, so for computing the shape function based on the moving least squares method (MLS), an efficient algorithm is presented. In this modification, additional terms is proposed to impose based on the coefficients of the polynomial base functions on the quadratic base functions (m = 2). So, the MMLS method is developed for solving the systems of two-dimensional linear and nonlinear integral equations at irregularly distributed nodes. This approach prevents the singular moment matrix in the context of MLS based on meshfree methods. Also, determining the best radius of the support domain of a field node is an open problem for MLS-based methods. Therefore, the next important thing is that the MMLS algorithm can automatically find the best neighborhood radius for each node. Then, numerical examples are presented that determine the main motivators for doing this so. These examples enable us to make comparisons of two methods: MMLS and classical MLS.

1. Introduction

In mathematics, there are many functional equations of the description of a real system in the natural sciences (such as physics, biology, Earth science, meteorology) and disciplines of engineering. For instance, we can point to some mathematical model from physics that describe heat as a partial differential equation and the inverse problem of it’s as integro-differential equations. Also, another example in nature is Laplace’s equation which corresponds to the construction of potential for a vector field whose effect is known at the boundary of Domain alone. Especially, the integral equations have wide applicability which has been cited in [1, 2, 3, 4].However, there are many significant analytical methods for solving integral equations but most of them especially in nonlinear cases, finding an analytical representation of the solution is so difficult, therefore, it is required to obtain approximate solutions. The interested reader can find several numerical methods for approximating the solution of these problems in [5, 6, 7, 8, 9, 10, 11, 12, 13, 14] and the references therein.

Moreover, there are various numerical and analytical methods have been used to estimate the solution of integrodifferential equations or Abels integral equations [12, 15, 16, 17, 18]. Recently the meshless based methods, particularly Moving Least Squares (MLS) method, for a solution of partial differential equations and ordinary differential equations have been paid attention. Using this approach some new methods such as meshless local boundary integral equation method [19], Boundary Node Method (BNM) [20], moving least square reproducing polynomial meshless method [21] and other relative methods are constructed. The new class of meshless methods has been developed which only relied on a set of nodes without the need for an additional mesh in the solution of a one-dimensional system of integral equations [22].

A local approximation of unknown function presented in the MLS method give us to possible choose the compact support domain for each data point as a sphere or a parallelogram box centered on a point [23, 24]. So each data point has an associated with the size of its compact support domain as dilatation parameter. Therefore the number of data point and dilatation parameter are direct effects on the MLS, Also by increasing the degree of the polynomial base function for complex data distributions give a more validated fashion. Nevertheless, in this case, it becomes more difficult to ensure the independence of the shape functions, and the least-squares minimization problem becomes ill-posed.

The common solution for increased the number of admissible node distribution is increasing the size of the support domains (a valid node distribution is referred to as an œadmissible node distribution [23]). There have been several proposed for choosing the radius of support domain [25], but one of the efficient suggestion was raised by Chen shen [26]. The author in [27] has introduced a new algorithm for selecting the suitable radius of the domain of influence. Also in [28], presented a modified MLS(MMLS) approximation on the shape function generation algorithm with additional terms based on the coefficients of the polynomial basis functions. It is an efficient method which has been proposed for handling a singular moment matrix in the MLS based methods. The advantage of this method compared to methods based on mesh such as a finite element or finite volume is this the domain of the problem is not important because this approximation method is based on a set of scattered points instead of domain elements for interpolation or approximation. So the geometry of the domain does not interfere in the MLS.

2. Methodology

2.1 Introduction of the MLS approximation

The Moving Least Square (MLS) method is a feasible numerical approximation method that is an extension of the least squares method, also it is the component of the class of meshless schemes that have a highly accurate approximation. The MlS approximation method is a popular method used in the many meshless methods [12, 19, 21, 22, 29, 30]. In many procedures used to construct the MLS shape function is used support-domain concept. The support domain of the shape function is a set of nodes in the problem domain that just those points directly contributes to the construction of the shape function, so the MLS shape function is locally supported. According to the classical least squares method, an optimization problem should be solved as follows

$$\mathit{min}\left(\sum _{j=1}^m{\left(u^h\left({\mathbf{x}}_j\right)-u\left({\mathbf{x}}_j\right)\right)}^2\right)$$

Where Ideally the approximation function $u^h\left(x\right)$ should match the function $u\left(\mathbf{x}\right)$. Therefore, in the MLS approach, a weight optimization problem will be solved which is dependent on nodal points. We start, with the basic idea of taking a set of the nodal points in $\mathrm{\Omega }$ so that $\mathrm{\Omega }\subseteq {\mathbb{R}}^d.$ Also ${\mathrm{\Omega }}_{\mathbf{x}}\subseteq \mathrm{\Omega }$ is neighboring nodes of point $\mathbf{x}$ and finding an approximation function with $m$ basis functions, in a system with $n$ equations as

$$T\left(U\right)=F$$

where $T$ consists of linear and nonlinear operators and $U=\left(u_1,u_2,\dots ,u_n\right)$ is the unknown vector of functions, also $F=\left(f_1,f_2,\dots ,f_n\right)$ is the known vector of functions.

So for the approximation of any of the $u_i,i=1,2,\dots ,n$ in ${\mathrm{\Omega }}_{\mathbf{x}},$$\forall x\in {\mathrm{\Omega }}_{\mathbf{x}},$$u_i^h\left(\mathbf{x}\right)$ can be defined as

$$u_i^h\left(\mathbf{x}\right)=\sum _{j=1}^m{a}_j\left(\mathbf{x}\right)p_j\left(\mathbf{x}\right)=P^T\left(\mathbf{x}\right)a\left(\mathbf{x}\right).$$E1

Let $\mathbf{P}=\left\{p_1,p_2,\dots p_m\right\}$ a set of polynomial of degree at most $m,m\in \mathbb{N}.$ Let $\mathbf{a}\left(\mathbf{x}\right)$ is a vector containing unknown coefficients $a_j\left(\mathbf{x}\right),j=1,2,\dots m$ dependent on the intrest point position. Also $m$ unknown functions $\mathbf{a}\left(\mathbf{x}\right)=\left(a_1\left(\mathbf{x}\right),a_2\left(x\right),\dots a_m\left(\mathbf{x}\right)\right)$ are chosen such that:

$$J\left(\mathbf{x}\right)=\sum _{j=1}^m{\left({\mathbf{P}}^T\left({\mathbf{x}}_j\right)\mathbf{a}\left(\mathbf{x}\right)-u_i\left({\mathbf{x}}_j\right)\right)}^2{w}_i\left(\mathbf{x}\right)={\left[P.\mathbf{a}-{\mathbf{u}}_i\right]}^T.W.\left[P.\mathbf{a}-{\mathbf{u}}_i\right],$$E2

is minimized. Note that the weight function $w_i\left(\mathbf{x}\right)$ is associated with node $j$. As we know, each redial basis function that define in [31] can be used as weight function, we can define $w_j\left(r\right)=\varphi \left(\frac{r}{\delta }\right)$ where $r={‖\mathbf{x}-{\mathbf{x}}_i‖}_2$ (the Euclidean distance between $\mathbf{x}$ and ${\mathbf{x}}_j)$ and $\varphi :{\mathbb{R}}^d\to \mathbb{R}$ is a nonnegative function with compact support. In this chapter, we will use following weight functions and will compare them to each other, corresponding to the node $j$, in the numerical examples.

a. Guass weight function

$$w\left(r\right)=\left\{\begin{array}{ll}\frac{exp\left(\frac{-r^2}{c^2}\right)-exp\left(\frac{-{\delta }^2}{c^2}\right)}{1-exp\left(\frac{-{\delta }^2}{c^2}\right)}& 0\le r\le \delta \\ 0& \mathit{\text{elsewhere}}.\end{array}\right.$$E3

b. RBF weight function

$$w\left(r\right)=\left\{\begin{array}{ll}{\left(1-r\right)}^6\left(6+36r+82r^2+72r^3+30r^4+5r^5\right)& 0\le r\le \delta \\ 0& \mathit{\text{elsewhere}}.\end{array}\right.$$E4

c. Spline weight function

$$w\left(r\right)=\left\{\begin{array}{ll}1-6{\left(\frac{r}{\delta }\right)}^2+8{\left(\frac{r}{\delta }\right)}^3-3{\left(\frac{r}{\delta }\right)}^4& 0\le r\le \delta \\ 0& \mathit{\text{elsewhere}}.\end{array}\right.$$E5

Where $c$ is constant and is called shape parameter. Also $\delta$ is the size of support domain.

$N$ is the number of nodes in ${\mathrm{\Omega }}_{\mathbf{x}}$ with $w_i\left(x\right)>0,$ the matrices $P$ and $W$ are defined as

$$P={\left[{\mathbf{p}}^T\left({\mathbf{x}}_1\right),{\mathbf{p}}^T\left({\mathbf{x}}_2\right),\dots {\mathbf{p}}^T\left({\mathbf{x}}_N\right)\right]}_{N\times \left(m+1\right)}^T$$E6

$$W=\mathit{diag}(\left(w_i\left(\mathbf{x}\right)\right),i=1,2,\dots ,N$$E7

and

$${\mathbf{u}}^h=\left[u_1^h,u_2^h,\dots u_n^h\right].$$E8

It is important to note that $u_i^T,i=1,2,\dots n,$ in (2) and (8) are the artificial nodal values, and not the nodal values of the unknown trial function $u^h\left(\mathbf{x}\right)$ in general. With respect to $\mathbf{a}\left(\mathbf{x}\right)$ and $u_i^T$ will be obtained,

$$A\left(\mathbf{x}\right)a\left(\mathbf{x}\right)=B\left(\mathbf{x}\right){\mathbf{u}}_i,$$E9

where the matrices $A\left(\mathbf{x}\right)$ and $B\left(\mathbf{x}\right)$ are defined by:

$$B\left(\mathbf{x}\right)=\left[w_1\mathbf{p}\left({\mathbf{x}}_1\right),w_2\mathbf{p}\left({\mathbf{x}}_2\right),\dots ,w_N\mathbf{p}\left({\mathbf{x}}_N\right)\right]$$E10

$$A\left(\mathbf{x}\right)=\sum _{i=1}^N{w}_i\left(\mathbf{x}\right){\mathbf{p}}^T\left({\mathbf{x}}_i\right)\mathbf{p}\left({\mathbf{x}}_i\right)={\mathbf{p}}^T\left(\mathbf{x}\right)w\left(\mathbf{x}\right)\mathbf{p}\left(\mathbf{x}\right).$$E11

The matrix $A\left(\mathbf{x}\right)$ in (11) is non-singular when the rank of matrix $P\left(\mathbf{x}\right)$ equals to $m$ and vice versa. In such a case, the MLS approximation is well-defined. With computing $\mathbf{a}\left(\mathbf{x}\right),u_i^h$ can be obtained as follows:

$$u_i^h\left(\mathbf{x}\right)=\sum _{j=1}^N{\varphi }_j\left(\mathbf{x}\right)u_i\left({\mathbf{x}}_j\right)={\phi }^T.{\mathbf{u}}_i$$E12

${\varphi }_j\left(\mathbf{x}\right)$ is called the shape function of the MLS approximation corresponding to the nodal point ${\mathbf{x}}_j,$ where

$$\phi \left(\mathbf{x}\right)={\mathbf{p}}^T\left(\mathbf{x}\right)A^{-1}\left(\mathbf{x}\right)B\left(\mathbf{x}\right)$$E13

Also with use the weight function, matrix $A,B$ are computed and then ${\varphi }_i\left(x\right)$ is determined from (13), If, further, $\varphi$ is sufficiently smooth, derivatives of $U$ are usually approximated by derivatives of $U^h,$

$$D^{\alpha }{u}_i\mathbf{x}\approx D^{\alpha }{u}_i^h\left(\mathbf{x}\right)=\sum _{j=1}^N{D}^{\alpha }{a}_j\left(\mathbf{x}\right)u_i\left({\mathbf{x}}_j\right),x\in \mathrm{\Omega }$$E14

2.2 Modify algorithm of MLS shape function

In the MLS approximation method, a local evaluation of the approximating unknown function is desired, and therefore for any nodal points the compact support domain is chosen as a sphere or a parallelogram box centered on the point [23, 29, 32]. This finding which the support domains contain what points. Each data point has a connected dilatation parameter $\lambda$ which is given ${\delta }_i=\lambda h_i$. Also, ${\delta }_i$ is the size of compact support domain in a node point ${\mathbf{x}}_i$.

Also, the necessary condition for that the moment matrix A be nonsingular is that for any point $x_i\in \mathrm{\Omega },i=1,2,\dots ,N$, [31].

$$\aleph \left(\left\{j|x_j\in {\mathrm{\Omega }}_{{\delta }_i}\right\}\right)\ge m,j=1,2,\dots ,N$$

So the dilatation parameters $\lambda$ determine the number of points of support domain, Also these points participate in the production of the shape function Therefore, $\lambda$ is very important and should be chosencorrectly so that the moment matrix $A$ is nonsingular.

In the remainder of this section, we introduce the new algorithm, with the aim of avoiding the singularity of the matrix A by choosing the correct $\lambda$ parameter by the algorithm.

Algorithm 1

Require:$X=\left\{x_i:i=1,2,\dots ,N\right\}$- Coordinates of points whose MLS shape function to be evaluated.

1: procedure Matrix A

2: $${\lambda }_{\mathit{new}}\leftarrow \lambda$$

3: $\alpha \leftarrow 0.01$ (This value selected experimentally.)

4: $\delta ={\lambda }_{\mathit{new}}\times h$ (h: the fill distance is defined to be $h={sup}_{x\in \mathrm{\Omega }}{min}_{1\le j\le N}\parallel x-x_i{\parallel }_2$)

5: Loop

6: set $I\left(x\right)=\left\{j\in \left\{1,2,\dots ,N\right\},\parallel x-x_i{\parallel }_2\le \delta \right\}$ (Using set of indices $I$, by localization at a fixed point $x$)

7: for$j\in I\left(x\right)$do

8: for$i=1:N$do

9: Compute $w_i$ for any $x_j\in {\mathrm{\Omega }}_i$

10: $$A=A+w_i{p}_i^2$$

11: end

12: end

13: if$\mathit{\text{cond}}\left(A\right)\ge \frac{1}{\epsilon }$then

14: {

15: $${\lambda }_{\mathit{new}}={\lambda }_{\mathit{new}}+\mathit{\alpha \lambda }$$

16: $$\delta ={\lambda }_{\mathit{new}}\times h$$

17: else

18: gotoend

19: }

20: if${\delta }_i\le \parallel X_{\mathrm{\Omega }}{\parallel }_2$then

21: gotoLoop

22: end

In Algorithm 1.

$X$: is a set containing N scattered points which are called centers or data site and I(x) is the Index of points which MLS shape function is evaluated.

$\alpha$: is a small positive number that is selected experimentally.

Then in every node points, matrix A is computed.

By running the algorithm the new value is assigned to $\lambda$, this value is related to the condition number of matrix A and its amount will increase. Therefore, the size of the support domain is increased and then the matrix A with new nodal points in the support domain is reproduced. This loop is repeated until $\frac{1}{\mathit{\text{cond}}\left(A\right)}\ge \epsilon$.

The main idea of the moving least squares approximation is that for every point x can solve a locally weighted least squares problem [30], it is a local approximation, thus the additional condition to stop the loop is the size of the local support domain, the value of $\lambda$ should be well enough to pave the local approximation, Line 20 is said to satisfy this condition.

2.3 Modified MLS approximation method

One of the common problems in Classic MLS method is the singularity of the moment matrix A in irregularity chosen nodal points. To avoid the nodal configurations which lead to a singular moment matrix, the usual solution is to enlarge the support domains of any nodal point. But it is not an appropriate solution, in [31] to tackle such problems is proposed a modified Moving least squares(MMLS)approximation method. This modifies allows, base functions in $m\ge 2$ to be used with the same size of the support domain as linear base functions $\left(m=1\right).$ We should note that,impose additional terms based on the coefficients of the polynomial base functions is the main view of the modified technique. As we know, in the basis function $\mathbf{p}\left(\mathbf{x}\right)$ is

$$\mathbf{p}\left(\mathbf{x}\right)={\left[1,x,x^2,\dots ,x^m\right]}^T$$E15

where $\mathbf{x}\in \mathbb{R}$, Also the correspond coefficients $a_j,$ that should be determined are [24]:

$$\mathbf{a}\left(\mathbf{x}\right)={\left[a_1,a_x,a_{x^2},\dots ,a_{x^m}\right]}^T$$E16

For obtaining these coefficients, the functional (2) rewrite as follows:

$$\overline{J}\left(\mathbf{x}\right)=\sum _{j=1}^m{\left({\mathbf{P}}^T\left({\mathbf{x}}_j\right)\mathbf{a}\left(\mathbf{x}\right)-u_i\left({\mathbf{x}}_j\right)\right)}^2{w}_i\left(\mathbf{x}\right)+\sum _{\nu =1}^{m-2}{\overline{w}}_{\nu }\left(x\right){\overline{\mathbf{a}}}_{\nu }^2\left(\mathbf{x}\right),i=1,2,\dots ,n$$E17

Where $\overline{w}$ is a vector of positive weights for the additional constraints, also $\overline{\mathbf{a}}={\left[a_{x^2},a_{x^3},\dots ,a_{x^m}\right]}^T$ is the modified matrix.

The matrix form of (17) is as follows:

$$\overline{J}\left(\mathbf{x}\right)={\left[P.\mathbf{a}-{\mathbf{u}}_i\right]}^T.W.\left[P.\mathbf{a}-{\mathbf{u}}_i\right]+{\mathbf{a}}^{T}H\mathbf{a},i=1,2,\dots ,n$$E18

where $H$ is as,

$$H=\left[\begin{array}{ll}{O}_{2,2}& O_{m-2,m-2}\\ O_{2,2}& \mathit{diag}\left(\overline{w}\right)\end{array}\right],$$E19

where, $O_{i,j}$ is the null matrix. By minimizing the functional (18), the coefficients $a\left(\mathbf{x}\right)$ will be obtained. So we have

$$\overline{A}\left(\mathbf{x}\right)\mathbf{a}\left(x\right)=B\left(\mathbf{x}\right){\mathbf{u}}_i,$$E20

where

$$\overline{A}=P^T.W.P+H$$E21

And the matrics $B\left(\mathbf{x}\right)$ is determined as the same before. So we have

$${\phi }_m\left(\mathbf{x}\right)=\mathbf{a}\left(\mathbf{x}\right)=p^T\left(\mathbf{x}\right){\overline{A}}^{-1}\left(\mathbf{x}\right)B\left(\mathbf{x}\right)$$E22

where ${\phi }_m\left(\mathbf{x}\right)$ is the shape function of the MMLS approximation method.

3. Stiff systems of ordinary differential equations

In this section, we use MLS approximation method for numerical solution of the Stiff system of ordinary differential equations so consider the following differential equation

$$A\left(U\right)-F\left(\mathbf{x}\right)=0,U\left(0\right)=U_0,\mathbf{x}\in \mathrm{\Omega }$$E23

with boundary conditions,

$$B\left(U,\frac{\partial U}{\partial \mathbf{x}}\right)=0,\mathbf{x}\in \mathrm{\partial \Omega }.$$

where A is a general differential operator, $U_0$ is an initial approximation of (23), $F\left(\mathbf{x}\right)$ is a vector of known analytical functions on the domain $\mathrm{\Omega }$ and $\mathrm{\partial \Omega }$ is the boundary of $\mathrm{\Omega }$. The operator can be divided by $A=L+N,$ where $L$ is the linear part, and $N$ is the nonlinear part of its. So (23) can be, rewritten as follows;

$$L\left(U\right)+N\left(U\right)-F\left(\mathbf{x}\right)=0$$E24

We assume that $\mathbf{a}=\left\{a_1,a_2,\dots ,a_m\right\}$ are the MLS shape functions so in order to solve system (24), $N$ nodal points $x_i$ are selected on the $\mathrm{\Omega }$, which $\left\{x_i|i=1,2,\dots ,N\right\}$ is q-unisolvent. The distribution of nodes could be selected regularly or randomly. Then instead of $u_j$ from $U$, we can replace $u_j^h$ from (13). So we have

$$u_j^h\left(\mathbf{x}\right)=\sum _{i=1}^N{a}_i\left(\mathbf{x}\right)u_j\left({\mathbf{x}}_i\right)$$E25

where $j=1,2,\dots ,n$ is the number of unknown functions. we estimate the unknown functions $u_i$ by (25), so the system (24) becomes the following form

$$\begin{array}{l}L\left(\sum _{i=1}^N{a}_i\left(\mathbf{x}\right)u_1\left({\mathbf{x}}_i\right),\sum _{i=1}^N{a}_i\left(\mathbf{x}\right)u_2\left({\mathbf{x}}_i\right),\dots ,\sum _{i=1}^N{a}_i\left(\mathbf{x}\right)u_n\left({\mathbf{x}}_i\right)\right)+\\ N\left(\sum _{i=1}^N{a}_i\left(\mathbf{x}\right)u_1\left({\mathbf{x}}_i\right),\sum _{i=1}^N{a}_i\left(\mathbf{x}\right)u_2\left({\mathbf{x}}_i\right),\dots ,\sum _{i=1}^N{a}_i\left(\mathbf{x}\right)u_n\left({\mathbf{x}}_i\right)\right)=\left(f_1\left(\mathbf{x}\right),f_2\left(\mathbf{x}\right),\dots ,f_n\left(\mathbf{x}\right)\right)+\mathbf{r}\left(\mathbf{x}\right).\end{array}$$E26

where $\mathbf{r}\left(\mathbf{x}\right)$ is residual error function which vanishes to zero in collocation points thus by installing the collocation points ${\mathbf{y}}_r;r=1,2,\dots ,N$, so

$$\begin{array}{l}L\left(\sum _{i=1}^N{a}_i\left({\mathbf{y}}_r\right)u_1\left({\mathbf{x}}_i\right),\sum _{i=1}^N{a}_i\left({\mathbf{y}}_r\right)u_2\left({\mathbf{x}}_i\right),\dots ,\sum _{i=1}^N{a}_i\left({\mathbf{y}}_r\right)u_n\left({\mathbf{x}}_i\right)\right)+\\ N\left(\sum _{i=1}^N{a}_i\left({\mathbf{y}}_r\right)u_1\left({\mathbf{x}}_i\right),\sum _{i=1}^N{a}_i\left({\mathbf{y}}_r\right)u_2\left({\mathbf{x}}_i\right),\dots ,\sum _{i=1}^N{a}_i\left({\mathbf{y}}_r\right)u_n\left({\mathbf{x}}_i\right)\right)=\\ \sum _{i=1}^{N}L\left(a_i\left({\mathbf{y}}_r\right)\right)u_1\left({\mathbf{x}}_i\right),\sum _{i=1}^{N}L\left(a_i\left({\mathbf{y}}_r\right)\right)u_2\left({\mathbf{x}}_i\right),\dots ,\sum _{i=1}^{N}L\left(a_i\left({\mathbf{y}}_r\right)\right)u_n\left({\mathbf{x}}_i\right))+\\ N\left(\sum _{i=1}^N{a}_i\left({\mathbf{y}}_r\right)u_1\left(x_i\right),\sum _{i=1}^N{a}_i\left({\mathbf{y}}_r\right)u_2\left({\mathbf{x}}_i\right),\dots ,\sum _{i=1}^N{a}_i\left({\mathbf{y}}_r\right)u_n\left({\mathbf{x}}_i\right)\right)=\\ \left(f_1\left({\mathbf{y}}_r\right),f_2\left(y_r\right),\dots ,f_n\left({\mathbf{y}}_r\right)\right)\end{array}$$E27

therefore

$$\mathit{CU}=\left[\begin{array}{cccc}L\left(a_1\left(y_1\right)\right)& L\left(a_2\left(y_1\right)\right)& \dots & L\left(a_N\left(y_1\right)\right)\\ L\left(a_1\left(y_2\right)\right)& L\left(a_2\left(y_2\right)\right)& \dots & L\left(a_N\left(y_2\right)\right)\\ ⋮& & & \\ L\left(a_1\left(y_N\right)\right)& L\left(a_2\left(y_N\right)\right)& \dots & L\left(a_N\left(y_N\right)\right)\end{array}\right]\left[\begin{array}{cccc}{u}_1\left(x_1\right)& u_2\left(x_1\right)& \dots & u_n\left(x_1\right)\\ u_1\left(x_2\right)& u_2\left(x_2\right)& \dots & u_n\left(x_2\right)\\ ⋮& & & \\ u_1\left(x_N\right)& u_2\left(x_N\right)& \dots & u_n\left(x_N\right)\end{array}\right]$$E28

And the matrix form of (27) as follows

$$C_{N\times N}{U}_{N\times n}+N_{N\times n}\left(\mathbf{a},U\right)=F_{N\times n}\left(y_r\right)$$E29

where

$$\begin{array}{c}{C}_i={\left[L\left(a_1\left({\mathbf{y}}_r\right)\right),\dots ,L\left(a_N\left({\mathbf{y}}_r\right)\right)\right]}_{i=1}^n\\ U_i={\left[{\left(u_i\left(x_1\right),u_i\left(x_2\right),\dots ,u_i\left(x_N\right)\right)}^T\right]}_{i=1}^n\\ F\left({\mathbf{y}}_r\right)={\left({\left[{\left(f_1\left({\mathbf{y}}_r\right)\right)}_{r=1}^N\right]}^T,{\left[{\left(f_2\left({\mathbf{y}}_r\right)\right)}_{r=1}^N\right]}^T,\dots \left[{\left(f_n\left({\mathbf{y}}_r\right)\right)}_{r=1}^N\right]\right)}^T.\end{array}$$E30

by imposing the initial conditions at $t=0$, we have

$$\left(\sum _{i=1}^N{a}_i\left(0\right)u_1\left(t_i\right),\sum _{i=1}^N{a}_i\left(0\right)u_2\left(t_i\right),\dots ,\sum _{i=1}^N{a}_i\left(0\right)u_n\left(t_i\right)\right)=U_0$$E31

and Solving the non-linear system (29) and (31), lead to finding quantities $u_j\left(x_i\right)$. Then the values of $u_j\left(x\right)$ at any point $x\in \mathrm{\Omega }$, can be approximated by Eq. (25) as:

$$u_j\left(\mathbf{x}\right)\simeq \sum _{i=1}^N{a}_i\left(\mathbf{x}\right)u_j\left({\mathbf{x}}_i\right)$$

Remark

Note that, for simplicity, the modification derivation is made for bivariate data, but can be easily extended to higher dimensions. Also, for implementation, the modified moving least squares approximation method it is sufficient to use ${\phi }_m$ instead of $\phi .$

3.1 Error analysis

The convergence analysis of the method in matrix norm has been investigated for the systems of one and two-dimensional Fredholm integral equations by authors of [22]. It should be noted that The convergence analysis of the method for implementation classic moving least squares approximation method on a system of integral equations has been discussed and it should be investigated for modified Mls method. we can use the results for this type of equations.

So in continuation of this section, the error estimations for the system of differential equations is developed. In [26], has obtained error estimates for moving least square approximations in the one-dimensional case. Also in [33], is developed for functional in n-dimensional and was used the error estimates to prove an error estimate in Galerkin coercive problems. In this work, have improved error estimate for the systems of stiff ordinary differential equations.

Given $\delta >0$ let $W_{\delta }\ge 0$ be a function such that $\mathit{supp}\left(w_{\delta }\right)\subset \overline{B_{\delta }\left(0\right)}=\left\{z\Vert z|\le \delta \right\}$ and $X_{\delta }=\left\{x_1,x_2,\dots ,x_n\right\},$$n=n\left(\delta \right)$, a set of points in $\mathrm{\Omega }\subset \mathbb{R}$ an open interval and let $U=\left(u_1,u_2,\dots ,u_N\right)$ be the unknown function such that $u_{i1},u_{i2},\dots ,u_{\mathit{\text{in}}}$ be the values of the function $u_i$ in those points, i.e., $u_{i,j}=u_i\left(x_j\right),i=1,\dots ,N,j=1,\dots ,n$. A class of functions $W={\left\{{\omega }_j\right\}}_{j=1}^N$ is called a partition of unity subordinated to the open covering $I_N$ if it possesses the following properties:

• $W_j\in C_s^0,s>0\mathit{\text{or}}s=\infty ,$

• $sup\left({\omega }_j\right)\subseteq {\overline{\mathrm{\Lambda }}}_j,$

• ${\omega }_j\left(x\right)>0,x\in {\mathrm{\Lambda }}_j,$

• $\sum _{i=1}^N{\omega }_j=1 \mathit{\text{for\hspace{0.17em}\hspace{0.17em}every}}x\in \overline{\mathrm{\Omega }}$

There is no unique way to build a partition of unity as defined above [34].

As we know, the MLS approximation is well defined if we have a unique solution at every point $x\in \overline{\mathrm{\Omega }}$. for minimization problem. And non-singularity of matrix $A\left(x\right)$, ensuring it is. In [33] the error estimate was obtained with the following assumption about the system of nodes and weight functions $\left\{{\mathrm{\Theta }}_N,W_N\right\}:$

We define

$$⟨u,v⟩=\sum _{j=1}^{n}w\left(x-x_j\right)u\left(x_j\right)v\left(x_j\right)$$

then

$$\parallel u{\parallel }_x^2=\sum _{j=1}^{n}w\left(x-x_j\right)u{\left(x_j\right)}^2$$

Also for vector of unknown functions, we define

$$\parallel U{\parallel }_{\infty }=\mathit{max}\left\{{\left|u_i\right|}_x,i=1,2,\dots ,N\right\}$$

are the discrete norm on the polynomial space ${\mathbb{P}}_m^1$ if the weight function $w$ satisfy the following properties.

a. For each $x\in \mathrm{\Omega }$, $w\left(x-x_j\right)>0$ at least for $\left(m+1\right)$ indices $j$.

b. For any $x\in \mathrm{\Omega }$, the moment matrix $A\left(x\right)=w\left(x\right)P^T$ is nonsingular.

Definition 3.1.Given$\mathit{x}\in \overline{\mathrm{\Omega }}$, the set$\mathit{ST}\left(\mathit{x}\right)=\left\{j:{\omega }_j\ne 0\right\}$will be called the star of$x$.

Theorem 3.1.[34,35] A necessary condition for the satisfaction of Property b is that for any$\mathit{x}\in \overline{\mathrm{\Omega }}$

$$n=\mathit{card}\left(\mathit{ST}\left(\mathit{x}\right)\right)\ge \mathit{card}\left({\mathbb{P}}_m\right)=m+1$$

For a sample point $\mathbf{c}\in \overline{\mathrm{\Omega }}$, if $\mathit{ST}\left(\mathbf{c}\right)=\left\{j_1,\dots j_k\right\}$, the mesh-size of the star $\mathit{ST}\left(\mathbf{c}\right)$ defined by the number is $h\left(\mathit{ST}\left(\mathbf{c}\right)\right)=max\left\{h_{j1},\dots h_{\mathit{jk}}\right\}.$

Assumptions. Consider the following global assumptions about parameters. There exist

$\left(a_1\right)$ An over bound of the overlap of clouds:

$$E={sup}_{c\in \overline{\mathrm{\Omega }}}\left\{\mathit{card}\left(\mathit{ST}\left(c\right)\right)\right\}.$$

$\left(a_2\right)$ Upper bounds of the condition number:

$${\mathit{CB}}_q={sup}_{c\in \overline{\mathrm{\Omega }}}\left\{{\mathit{CN}}_q\left(\mathit{ST}\left(c\right)\right),q=1,2\right\}.$$

where the numbers ${\mathit{CN}}_q\left(\mathit{ST}\left(\mathbf{c}\right)\right)$ are computable measures of the quality of the star $\mathit{ST}\left(c\right)$ which defined in Theorem 7 of [19].

$\left(a_3\right)$ An upper bound of the mesh-size of stars:

$$R={sup}_{c\in \overline{\mathrm{\Omega }}}\left(\mathit{hST}\left(c\right)\right).$$

$\left(a_4\right)$ An uniform bound of the derivatives of $\left\{w_j\right\}.$ That is the constant $G_q>0,q=1,2$, such that

$${‖D^{\mu }{W}_j‖}_{L_{\infty }}\le \frac{G_q}{R^{\mid \mu \mid }}1<\mu <q,$$

$\left(a_5\right)$ There exist the number $\gamma >0$ such that any two points $\mathbf{x}\mathbf{,}\mathbf{y}\in \mathrm{\Omega }$ can be joined by a rectifiable curve $\mathrm{\Gamma }$ in $\mathrm{\Omega }$ with length $\mid \mathrm{\Gamma }\mid \le \gamma ‖\mathbf{x}‐\mathbf{y}‖.$

Assuming all these conditions, Zuppa [34] proved.

Lemma 3.1.$U=\left(u_1,u_2,\dots u_n\right)$such that$u_i\in C^{m+1}\left(\overline{\mathrm{\Omega }}\right)$and${‖U‖}_{\infty }=u_k,1<k<n$, There exist constants$C_q,q=1 \mathit{\text{or}} 2,$

$$C_1=C_1\left(\gamma ,d,E,G_1,{\mathit{CB}}_1\right),$$

$$C_2=C_1\left(\gamma ,d,E,G_2,{\mathit{CB}}_1,{\mathit{CB}}_2\right),$$

then

$${‖D^{\mu }U-D^{\mu }{U}^h‖}_{\infty }<C_q{R}^{q+1-\mid \mu \mid }\parallel u_k^{\left(m+1\right)}{\parallel }_{L^{\infty }\left(\mathrm{\Omega }\right)}0<\mu <q$$E32

Proof: see [36].

3.2 System of ODE

If in (24) the non-linear operator $N$ be zero, we have

$$L\left(U\right)=\left(f_1,f_2,\dots ,f_n\right)$$E33

where $U$ is the vector of unknown function and $L$ is a matrix of derivative operators,

$$L\left(U\left(.\right)\right)=\sum _{\varsigma =1}^n{\lambda }_{\varsigma }\frac{{\partial }^{\varsigma }}{{\left(.\right)}^{\varsigma }}U\left(.\right).$$E34

And from (25), we define

$$U^h\left(t\right)=\sum _{i=1}^N{a}_i\left(t\right)U\left(t_i\right)$$

where ${\left(a_i\right)}_{i=1}^N$ are the MLs shape functions defined on the interval $\left[0,1\right]$ satisfying the homogeneous counterparts of the boundary conditions in (23). Also if the weight function $w$ possesses $k$ continuous derivatives then the shape functions $a_j$ is also in $C^k$ [33]. By the collocation method, is obtained an approximate solution $U^h\left(t\right)$. And demand that

$$L^h\left(U\left(.\right)\right)=\sum _{\varsigma =0}^n{\lambda }_{\varsigma }\frac{{\partial }^{\varsigma }}{{\left(.\right)}^{\varsigma }}{U}^h\left(.\right)$$E35

where ($\lambda =0 \mathit{\text{or}} 1$). It is assumed that in the system of ODE derivative of order at most $n=2$. Each of the basis functions $a_i$ has compact support contained in $\left(0,1\right)$ then the matrix $C$ in (30) is a bounded matrix. If $\delta$ be fixed, independent of $N$, then the resulting system of linear equations can be solved in $O\left(N\right)$ arithmetic operations.

Lemma 3.2.Let$U=\left(u_1,u_2,\dots u_n\right)$and$F=\left(f_1,f_2,\dots f_n\right)$so that$u_i\in C^{m+1}\left(\overline{\mathrm{\Omega }}\right)$$m\ge 1$and$\parallel u_i{\parallel }_{\infty }=u_k,k\in \left\{1,2,\dots ,n\right\}$where$\mathrm{\Omega }$be a closed, bounded set in$R$. Assume the quadrature scheme is convergent for all continuous functions on$\mathrm{\Omega }.$Further, assume that the stiff system of ODE(23)with the fixd initial condition is uniquely solvable for given$f_i\in C\left(\mathrm{\Omega }\right).$Moreover take a suitable approximation$U^h$of$U$Then for all sufficiently large$n$, the approximate matrix$L$for linearly case exist and are uniformly bounded,$\mid L\mid \le M$with a suitable constant$M<\infty .$For the equations$L\left(U\right)=F$and$L^h\left(U\right)=F$we have

$$E_t=\parallel L\left(U\left(t\right)\right)-L^h\left(U\left(t\right)\right){\parallel }_{\infty }$$

so that

$$\parallel E_t{\parallel }_{\infty }\le {\mathrm{C}}_{q}K\left(\lambda ,\varsigma \right)R^{m+1-\mu }\parallel u_k^{\left(m+1\right)}{\parallel }_{L_{\infty }}.$$

Proof. we have

$$\parallel L\left(U\left(t\right)\right)-L^h\left(U\left(t\right)\right){\parallel }_{\infty }=\parallel \sum _{\varsigma =0}^n{\lambda }_{\varsigma }\frac{{\partial }^{\varsigma }}{t^{\varsigma }}U\left(t\right)-\sum _{\varsigma =0}^n{\lambda }_{\varsigma }\frac{{\partial }^{\varsigma }}{t^{\varsigma }}{U}^h\left(t\right){\parallel }_{\infty }$$

so due to the lemma (36),

$$\begin{array}{l}\parallel L\left(U\left(t\right)\right)-L^h\left(U\left(t\right)\right){\parallel }_{\infty }\le \sum _{\varsigma =0}^n\mid {\lambda }_{\varsigma }\mid \parallel \frac{{\partial }^{\varsigma }}{t^{\varsigma }}U\left(t\right)-\frac{{\partial }^{\varsigma }}{t^{\varsigma }}{U}^h\left(t\right){\parallel }_{\infty }\\ \le \underset{i}{max}\sum _{\varsigma =0}^n\mid {\lambda }_{\varsigma }\Vert \frac{{\partial }^{\varsigma }}{t^{\varsigma }}{u}_i\left(t\right)-\frac{{\partial }^{\varsigma }}{t^{\varsigma }}{u}_i^h\left(t\right)\mid \\ \le \sum _{\varsigma =0}^n{\mathrm{C}}_q\mid {\lambda }_{\varsigma }\mid \parallel u_k^{\left(m+1\right)}{\parallel }_{L_{\infty }}{R}^{m+1-\varsigma }\end{array}$$

where should be $m\ge \varsigma$ so,

$$\sum _{\varsigma =0}^n\mid {\lambda }_{\varsigma }\mid R^{m+1-\varsigma }\le K\left(\lambda ,\varsigma \right)R^{m+1-\mu }$$

where $\mu$ is the highest order derivative And $K\left(\lambda ,\varsigma \right)=\sum _{\varsigma =0}^n\mid {\lambda }_{\varsigma }\mid$, so demanded that

$$\parallel E_t{\parallel }_{\infty }\le {\mathrm{C}}_{q}K\left(\lambda ,\varsigma \right)R^{m+1-\mu }\parallel u_k^{\left(m+1\right)}{\parallel }_{L_{\infty }}.$$

It should be noted that in the nonlinear system the upper bound of error depends on the nonlinear operator.

4. Two-dimensional linear systems of integral equations

4.1 Fredholm type

In this section, we will provide an approximation solution of the 2-D linear system of Fredholm integral equations by the MLS method. The matrix form of this system could be considered as

$$\mathbf{U}\left(x,y\right)=\mathbf{F}\left(x,y\right)+\underset{\mathrm{\Omega }}{\int }K\left(x,y,\theta ,s\right)\mathbf{U}\left(\theta ,s\right)\mathit{\text{dθds}},\left(x,y\right)\in \mathrm{\Omega },$$E36

where $\mathrm{\Omega }=\left[a,b\right]\times \left[c,d\right]$ as $\mathrm{\Omega }\subset {\mathbb{R}}^2$, Also $K\left(x,y,\theta ,s\right)=\left[{\kappa }_{\mathit{ij}}\left(x,y,\theta ,s\right)\right], i,j=1,2,\dots ,n$ is the matrix of kernels, $\mathbf{U}\left(x,y\right)={\left(u_1\left(x,y\right),u_2\left(x,y\right),\dots u_n\left(x,y\right)\right)}^T$ is the vector of unknown function and $\mathbf{F}\left(x,y\right)={\left(f_1\left(x,y\right),f_2\left(x,y\right),\dots f_n\left(x,y\right)\right)}^T$ is the vector of known functions.

In addition, is took that two cases for the domain, the rectangular shape, and nonrectangular one and three cases relative to the geometry of the nonrectangular domain are considered where can be transformed into the rectangular shape [35].

The first one is $\mathrm{\Omega }=\left\{\left(\theta ,s\right)\in {\mathbb{R}}^2:a\le s\le b,g_1\left(s\right)\le \theta \le g_2\left(s\right)\right\}$ where $g_1\left(s\right)$ and $g_2\left(s\right)$ are continues functions of $s$, the second one can be consider as $\mathrm{\Omega }=\left\{\left(\theta ,s\right)\in {\mathbb{R}}^2:c\le \theta \le d,g_1\left(\theta \right)\le s\le g_2\left(\theta \right)\right\}$ where $g_1\left(\theta \right)$ and $g_2\left(\theta \right)$ are continues functions of $\theta$, Also the last one is a domain which is neither of the first nor the second kinds but could be separated to finite numbers of first or second sub-domains, it is labeled as a domain of third kind. Without loss of generality, the first kind domain can be assumed as

$$\mathrm{\Omega }=\left\{\left(\theta ,s\right)\in {\mathbb{R}}^2:-1\le s\le 1,g_1\left(s\right)\le \theta \le g_2\left(s\right)\right\}$$E37

by the following linear transformation

$$\theta \left(t,s\right)=\frac{g_2\left(\theta \right)-g_1\left(\theta \right)}{2}t+\frac{g_2\left(\theta \right)-g_1\left(\theta \right)}{2},$$E38

the interval$\left[g_1\left(\theta \right),g_2\left(\theta \right)\right]$ is converted to the fixed interval $\left[-1,1\right]$, so we have

$$\mathbf{U}\left(x,y\right)=\mathbf{F}\left(x,y\right)+\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\mathrm{K}\left(x,y,t,s\right)\mathbf{U}\left(t,s\right)\mathit{\text{dtds}}, \mathit{\text{such\hspace{0.17em}that}}\left(x,y\right)\in \left[-1,1\right]\times \left[-1,1\right]$$E39

where

$$\mathrm{K}\left(x,y,t,s\right)=\frac{g_2\left(\theta \right)-g_1\left(\theta \right)}{2}\mathrm{K}\left(x,y,\theta ,s\right)$$E40

Also, the second kind is straight similarly by commuting the variables and the third kind can be separated to finite numbers of sub-domains of the first or second kinds, so the method can be applied in each sub-domain as described earlier.

So, for the numerical integration $\mathrm{\Omega }={\bigcup }_{l=1}^L{\mathrm{\Omega }}_l$ and ${\mathrm{\Omega }}_l\bigcap {\mathrm{\Omega }}_k\ne \varnothing$, $1\le k,l\le L$

$$\underset{\mathrm{\Omega }}{\int }g\left(s\right)\mathit{ds}=\sum _{l=1}^L\underset{{\mathrm{\Omega }}_l}{\int }g\left(s\right)\mathit{ds}$$E41

Here, the MLS method is applied for the general case where the domain is $\left[a,b\right]\times \left[c,d\right]$.

To apply the method, as described in section 2.1, instead of $u_i$ from $U$, we can replace ${\mathit{u}}_{\mathit{i}}^{\mathit{h}}$ from (12). So we have

$$U^h\left(\mathbf{x}\right)={\left(u_1^h\left(\mathbf{x}\right),u_2^h\left(\mathbf{x}\right),\dots ,u_n^h\left(\mathbf{x}\right)\right)}^T$$E42

Also, obviously from (12)

$$u_j^h\left(x,y\right)=\sum _{i=1}^N{\varphi }_i\left(x,y\right)u_j\left(x_i,y_i\right)$$E43

in this section, is assumed that the domain has rectangular shape, so system (36) becomes as follows

$${\left(u_1^h\left(x,y\right),u_2^h\left(x,y\right),\dots u_n^h\left(x,y\right)\right)}^T=f\left(x,y\right)+\underset{c}{\overset{d}{\int }}\underset{a}{\overset{b}{\int }}\left[k_{\mathit{ij}}\left(x,y,t,s\right)\right]{\left(u_1^h\left(t,s\right),u_2^h\left(t,s\right),\dots u_n^h\left(t,s\right)\right)}^T\mathit{\text{dtds}}.$$E44

By substituting (42) in (44), and it holds at points $\left(x_r,y_r\right),r=1,2,\dots ,N$ we have

$$\begin{array}{l}\mathbf{f}\left(x_r,y_r\right)={\left(\sum _{i=1}^N{\varphi }_i\left(x_r,y_r\right)u_1\left(x_i,y_i\right),\sum _{i=1}^N{\varphi }_i\left(x_r,y_r\right)u_2\left(x_i,y_i\right),\dots \sum _{i=1}^N{\varphi }_i\left(x_r,y_r\right)u_n\left(x_i,y_i\right)\right)}^T\\ -\underset{c}{\overset{d}{\int }}\underset{a}{\overset{b}{\int }}\left[k_{\mathit{ij}}\left(x_r,y_r,t,s\right)\right]{\left(\sum _{i=1}^N{\varphi }_i\left(t,s\right)u_1\left(x_i,y_i\right),\dots ,\sum _{i=1}^N{\varphi }_i\left(t,s\right)u_n\left(x_i,y_i\right)\right)}^T\mathit{\text{dtds}}.\end{array}$$E45

We consider the $m_1$-point numerical integration scheme over $\mathrm{\Omega }$ relative to the coefficients $\left({\tau }_k,{\varsigma }_p\right)$ and weights ${\omega }_k$ and ${\omega }_p$ for solving integrals in (45), i.e.,

$${\left({\mathcal{Ϝ}}_{\mathit{N}}\right)}_j{u}_j\left(x,y\right)=\sum _{p=1}^{m_1}\sum _{k=1}^{m_1}\left(\sum _{i=1}^N{k}_{\mathit{ji}}\left(x_r,y_r,{\tau }_k,{\varsigma }_k\right){\varphi }_i\left({\tau }_k,{\varsigma }_k\right){\omega }_k{\omega }_p\right),\left(x,y\right)\in \mathrm{\Omega },u_i\in \left(-\infty ,\infty \right)$$E46

Applying the numerical integration rule (46) yields

$$\begin{array}{l}\mathbf{f}\left(x_r,y_r\right)=(\sum _{i=1}^N\left({\varphi }_i\left(x_r,y_r\right)-\sum _{p=1}^{m_1}\sum _{k=1}^{m_1}\left(\sum _{j=1}^N{k}_{j1}\left(x_r,y_r,{\tau }_k,{\varsigma }_k\right){\varphi }_i\left({\tau }_k,{\varsigma }_k\right){\omega }_k\mathit{\omega p}\right)\right)u_{1i},\\ \sum _{i=1}^N\left({\varphi }_i\left({\tau }_k,{\varsigma }_k\right)-\sum _{p=1}^{m_1}\sum _{k=1}^{m_1}\left(\sum _{j=1}^N{k}_{j2}\left(x_r,y_r,{\tau }_k,{\varsigma }_k\right){\varphi }_i\left({\tau }_k,{\varsigma }_k\right){\omega }_k\mathit{\omega p}\right)\right)u_{2i},\\ \dots \sum _{i=1}^N\left({\varphi }_i\left({\tau }_k,{\varsigma }_k\right)-\sum _{p=1}^{m_1}\sum _{k=1}^{m_1}\left(\sum _{j=1}^N{k}_{\mathit{jn}}\left(x_r,y_r,{\tau }_k,{\varsigma }_k\right){\varphi }_i\left({\tau }_k,{\varsigma }_k\right){\omega }_k\mathit{\omega p}\right)\right)u_{\mathit{ni}}){}^T,r=1,2,\dots N\end{array}$$

where ${\left(u_j\right)}_i$ are the approximate quantities of $u_j$ when we use a quadrature rule instead of the exact integral. Now if we set $F_l$, $l=1,2,\dots n$ as a $N$ by $N$ matrices defined by:

$${\left(F_l\right)}_{i,j}={\varphi }_i\left(x_r,y_r\right)-\sum _{p=1}^{m_1}\sum _{k=1}^{m_1}\left(\sum _{j=1}^N{k}_{\mathit{jl}}\left(x_r,y_r,{\tau }_k,{\varsigma }_p\right){\varphi }_i\left({\tau }_k,{\varsigma }_p\right){\omega }_k\right){\omega }_p$$E47

So, the moment matrix $\mathrm{F}$ is defined by (47) as follows

$$\mathrm{F}={\left[F_1,F_2,\dots F_n\right]}_{\mathit{nN}\times \mathit{nN}}$$E48

And

$$\begin{array}{c}U={\left[{\left(u_{11},u_{12},\dots u_{1N}\right)}^T,{\left(u_{21},u_{22},\dots u_{2N}\right)}^T,\dots {\left(u_{n1},u_{n2},\dots u_{\mathit{nN}}\right)}^T\right]}^T\\ \mathbf{f}\left(x_r,y_r\right)={\left({\left[{\left(f_1\left(x_r,y_r\right)\right)}_{r=1}^N\right]}^T,{\left[{\left(f_2\left(x_r,y_r\right)\right)}_{r=1}^N\right]}^T,\dots \left[{\left(f_n\left(x_r,y_r\right)\right)}_{r=1}^N\right]\right)}^T.\end{array}$$E49

So by solving the following linear system of equations with a proper numerical method such as Gauss elimination method or etc. leads to quantities, $u_{\mathit{ji}}$.

$$\mathrm{F}U=\mathbf{f}$$E50

Therefore any $u_j\left(x,y\right)$ at any arbitrary point $\left(x,y\right)$ from the domain of the problem, can be approximated by Eq. (43) as

$$u_j\left(x,y\right)\approx \sum _{i=1}^N{\varphi }_i\left(x,y\right)u_{\mathit{ji}}\left(x_i,y_i\right)$$E51

4.2 Volterra type

Implementation of the proposed method on the Volterra integral equations is very simple and effective. In this case, the domain under study is as $\mathrm{\Omega }=\left[a,x\right]\times \left[c,y\right]$ such that $0\le x\le 1,0\le y\le 1$ and $a,c$ are constant, so a Volterra system type of integral equations can be consider as

$$\mathbf{U}\left(x,y\right)=\mathbf{F}\left(x,y\right)+\underset{\mathrm{\Omega }}{\int }\mathrm{K}\left(x,y,t,s\right)\mathbf{U}\left(t,s\right)\mathit{\text{dtds}},\left(x,y\right)\in \mathrm{\Omega },$$E52

like the Fredholm type, it is the matrix form of a system, so we have

$$\begin{array}{c}\mathbf{U}\left(x,y\right)={\left(u_1\left(x,y\right),u_2\left(x,y\right),\dots u_n\left(x,y\right)\right)}^T,\mathit{\text{the\hspace{0.17em}\hspace{0.17em}vector\hspace{0.17em}of\hspace{0.17em}unknown\hspace{0.17em}functions}}\\ \mathbf{F}\left(x,y\right)={\left(f_1\left(x,y\right),f_2\left(x,y\right),\dots f_n\left(x,y\right)\right)}^T,\mathit{\text{the\hspace{0.17em}\hspace{0.17em}vector\hspace{0.17em}of\hspace{0.17em}known\hspace{0.17em}functions}}\\ \mathrm{K}\left(x,y,t,s\right)=\left[{\kappa }_{\mathit{ij}}\left(x,y,t,s\right)\right]i,j=1,2,\dots ,n\mathit{\text{the\hspace{0.17em}matrix\hspace{0.17em}of\hspace{0.17em}kernels}}.\end{array}$$E53

By the following transformation the interval $\left[a,x\right]$ and $\left[c,y\right]$ can be transferred to a fixed interval $\left[a,b\right]$ and $\left[c,d\right]$,

$$t\left(x,\theta \right)=\frac{x-a}{b-a}\theta +\frac{b-x}{b-a}a.$$E54

$$s\left(y,\xi \right)=\frac{y-c}{d-c}\xi +\frac{d-y}{d-c}c.$$E55

Then instead of $u_i$ from $U,$ we can replace $u_i^h$ from (12). So we have

$$U^h\left(\mathbf{x}\right)={\left(u_1^h\left(\mathbf{x}\right),u_2^h\left(\mathbf{x}\right),\dots u_n^h\left(\mathbf{x}\right)\right)}^T$$E56

where

$$u_i^h\left(\mathbf{x}\right)=\sum _{j=1}^N{\varphi }_j\left(\mathbf{x}\right)u_i\left({\mathbf{x}}_j\right)$$E57

where $\mathbf{x}=\left(x,y\right)\in \left[a,b\right]\times \left[c,d\right],$ thus, system (52) becomes

$${\left(u_1^h\left(\mathbf{x}\right),u_2^h\left(\mathbf{x}\right),\dots u_n^h\left(\mathbf{x}\right)\right)}^T=F\left(\mathbf{x}\right)+\underset{a}{\overset{\mathbf{x}}{\int }}\underset{c}{\overset{\mathbf{y}}{\int }}\left[{\kappa }_{\mathit{ij}}\left(x,y,t,s\right)\right].{\left(u_1^h\left(x,y\right),u_2^h\left(x,y\right),\dots u_n^h\left(x,y\right)\right)}^T\mathit{\text{dtds}}.$$E58

Therefore from (54) and (55), the system (58) takes the following form

$${\left(u_1^h\left(\mathbf{x}\right),u_2^h\left(\mathbf{x}\right),\dots u_n^h\left(\mathbf{x}\right)\right)}^T=F\left(\mathbf{x}\right)+\underset{a}{\overset{\mathbf{b}}{\int }}\underset{c}{\overset{\mathbf{d}}{\int }}\left[{\kappa }_{\mathit{ij}}\left(x,y,t,s\right)\right].{\left(u_1^h\left(x,y\right),u_2^h\left(x,y\right),\dots u_n^h\left(x,y\right)\right)}^T\mathit{\text{dθdξ}}.$$E59

Where

$$K\left(.,.,.,.\right)=\frac{x-a}{b-a}\frac{y-c}{d-c}\mathrm{K}\left(.,.,.,.\right),$$E60

Using techniques employed in the Fredholm case yields the same final linear system where

$${\left(F_l\right)}_{i,j}={\varphi }_i\left(x_r,y_r\right)-\sum _{p=1}^{m_1}\sum _{k=1}^{m_1}\left(\sum _{j=1}^N{k}_{\mathit{jl}}\left(x_r,y_r,{\tau }_k,{\varsigma }_p\right){\varphi }_i\left({\tau }_k,{\varsigma }_p\right){\omega }_k\right){\omega }_p$$E61

where $l=1,2,\dots ,n$.

5. Nonlinear systems of two-dimensional integral equation

5.1 Fredholm type

In the nonlinear system, the unknown function cannot be written as a linear combination of the unknown variables or functions that appear in them, so the matrix form of Fredholm integral equations defined as the following form [27].

$$\mathbf{U}\left(x,y\right)=\mathbf{F}\left(x,y\right)+\underset{\mathrm{\Omega }}{\int }\mathbf{K}\left(x,y,\theta ,s,\mathbf{U}\left(\theta ,s\right)\right)\mathit{\text{dθds}},\left(x,y\right)\in \mathrm{\Omega },$$E62

Where $\mathbf{U}\left(x,y\right),$$\mathrm{K}$ and $\mathbf{F}$ are defined as,

$$\begin{array}{cc}\mathbf{U}\left(x,y\right)& ={\left(u_1\left(x,y\right),u_2\left(x,y\right),\dots ,u_n\left(x,y\right)\right)}^T\\ \mathbf{K}\left(x,y,\theta ,s,\mathbf{U}\left(\theta ,s\right)\right)& =\left[k_{\mathit{ij}}\left(x,y,\theta ,s,\mathbf{U}\left(\theta ,s\right)\right)\right],i,j=1,2,\dots ,n\\ \mathbf{F}& ={\left(f_1,f_2,\dots ,f_n\right)}^T\end{array}$$

As mentioned above, we assume that $\mathrm{\Omega }=\left[a,b\right]\times \left[c,d\right]$.

To apply the aproximation MLS method, we estimate the unknown functions $u_i$ by (12), so the system (62) becomes the following form

$${\left(u_1^h\left(x,y\right),u_2^h\left(x,y\right),\dots u_n^h\left(x,y\right)\right)}^T=f\left(x,y\right)+\underset{c}{\overset{d}{\int }}\underset{a}{\overset{b}{\int }}\left[k_{\mathit{ij}}\left(x,y,t,s,{\mathbf{U}}^h\left(t,s\right)\right)\right]\mathit{\text{dtds}}$$E63

We consider the $m_1$-point numerical integration scheme over the domain under study relative to the coefficients $\left({\tau }_k,{\varsigma }_p\right)$ and weights ${\omega }_k$ and ${\omega }_p$ for solving tow-dimentional integrals in (63), i.e.,

$${\left({\mathcal{Ϝ}}_N\right)}_i{u}_i\left(x,y\right)=\sum _{p=1}^{m_1}\sum _{k=1}^{m_1}{k}_{\mathit{ji}}\left(x,y,{\tau }_k,{\varsigma }_p,\sum _{j=1}^N{\varphi }_j\left({\tau }_k,{\varsigma }_p\right)u_i\left(x,y\right){\omega }_k{\omega }_p\right),\left(x,y\right)\in \mathrm{\Omega },u_i\in \left(-\infty ,\infty \right)$$E64

Applying the numerical integration rule (64) in (63) yields

$${\left(u_1^h\left(x_r,y_r\right),u_2^h\left(x_r,y_r\right),\dots ,u_n^h\left(x_r,y_r\right)\right)}^T=f\left(x_r,y_r\right)+\sum _{p=1}^{m_1}\sum _{k=1}^{m_1}\left[k_{\mathit{ij}}\left(x_r,y_r,{\tau }_k,{\varsigma }_p,{\mathbf{U}}^h\left({\tau }_k,{\varsigma }_p\right)\right)\right]\mathit{\text{dtds}}$$E65

Finding unknowns ${\mathbf{U}}^h$ by solving the nonlinear system of algebraic Eq. (65) yields the following approximate solution at any point $\left(x,t\right)\in \mathrm{\Omega }.$

$$u_j\left(x,y\right)\approx \sum _{i=1}^N{\varphi }_i\left(x,y\right)u_{\mathit{ji}}\left(x_i,y_i\right)$$E66

5.2 Volterra type

Two-dimensional nonlinear system of Volterra integral equations can be considered as the following form

$$\mathbf{U}\left(x,y\right)=\mathbf{F}\left(x,y\right)+\underset{a}{\overset{x}{\int }}\underset{c}{\overset{y}{\int }}\mathrm{K}\left(x,y,\theta ,s,\mathbf{U}\left(\theta ,s\right)\right)\mathit{\text{dθds}},\left(x,y\right)\in \mathrm{\Omega },$$E67

where $\mathrm{K}$, $\mathbf{F}$ are known function and $\mathbf{U}$ the vector of unknown functions are defined in (63) [27]. In order to apply the MLS approximation method, as same as the linear type, the interval $\left[a,x\right]$ and $\left[c,y\right]$ transferred to a fixed interval $\left[a,b\right]$ and $\left[c,d\right]$. Then $u_i^h,i=1,2,\dots ,n$ instead of $u_i$ in $U=(u_1,u_2,\dots ,u_n$ from (12) is replaced. So the nonlinear system (67) is converted to

$${\left(u_1^h\left(x,y\right),u_2^h\left(x,y\right),\dots ,u_n^h\left(x,y\right)\right)}^T=f\left(x,y\right)+\underset{c}{\overset{d}{\int }}\underset{a}{\overset{b}{\int }}\overline{K}\left(x,y,t,s,{\mathbf{U}}^h\left(t,s\right)\right)\mathit{\text{dtds}}$$E68

where

$$\overline{K}\left(x,y,t,s,{\mathbf{U}}^h\left(t,s\right)\right)=\frac{\mathbf{x}-a}{b-a}\frac{\mathbf{y}-c}{d-c}K\left(x,y,t,s,{\mathbf{U}}^h\left(t,s\right)\right),$$E69

Using the numerical integration technique (64) which applied in the Fredholm case yields the same final nonlinear system (65), so the approximation solution of $\mathbf{U}$ would be found by solving this system of equations.

6. Examples

In this section, the proposed method can be applied to the system of 2-dimensional linear and nonlinear integral equations [37] and the system of differential equations. Also, the results of the examples illustrate the effectiveness of the proposed method Also the relative errors for the collocation nodal points is used.

$$\parallel e_i{\parallel }_{\infty }=\frac{\parallel u_{\mathit{iex}}\left(x,y\right)-u_i^h\left(x,y\right)\parallel }{\parallel u_{\mathit{iex}}\left(x,y\right)\parallel }$$

where $u_i^h$ is the approximate solution of the exact solution $u_{\mathit{iex}}$. Linear and quadratic basis functions are utilized in computations.

6.1 Example 1

As the first example, we consider the following system of nonlinear Fredholm integral equations [27].

$$\begin{array}{c}{u}_1\left(x,y\right)=f_1\left(x,y\right)+\underset{\mathrm{\Omega }}{\int }{u}_1\left(s,t\right)u_2\left(s,t\right)\mathit{\text{dsdt}}\\ u_2\left(x,y\right)=f_2\left(x,y\right)+\underset{\mathrm{\Omega }}{\int }{u}_1\left(s,t\right)u_2\left(s,t\right)+u_2^2\left(s,t\right)\mathit{\text{dsdt}}\end{array}$$

where $\mathrm{\Omega }=\left[0,1\right]\times \left[0,1\right].$ The exact solutions are $U\left(x,y\right)=\left(x+y,x\right)$ and the $F\left(x,y\right)=\left(x+y-\frac{7}{12},x-\frac{11}{12}\right).$ The distribution of randomly nodes is shown in Figure 1. By attention to the irregular nodal points distribution, unsuitable $\delta$ can lead to a singular matrix A. So in this example, the adapted algorithm can tackle such problems. The MLS and MMLS shape functions are computed by using Algorithm 1, so the exact value of the radius of the domain of influence is not important; in fact, it is chosen as an initial value.

Figure 1.

The scatter data of example 1.

The condition numbers of the matrix A is shown in Figure 2 and the determinant of A at sample points $p$ is shown in Figure 3, where the radius of support domains for any nodal points is started from $\delta =0.05$. Note that, there is a different radius of support domain for any node point, it might be increased due to the inappropriate distribution of scattered points by the algorithm. These variations are shown in Table 1 for sample points $\left(x,y\right)$, where $\mathit{\text{Cond}}\left(A\right)$ is the conditions number A, its initial case ($\delta =0.005$) and final case ($\mathit{\text{Newδ}},$) and $N.O.\mathit{\text{iteration}}$ is the number of iteration of the algorithm for determining a suitable radius of support domain.

Figure 2.

The condition numbers of a at a sample point $\mathit{p}$ and $\delta =\mathit{0.05}$ for example 1. Using algorithm 1.

	Sample points		Cond(A)		Result of algorithm
n	$x$	$y$	initial	final	Newδ	N.O.iteration
1	$0.2575$	$0.4733$	$1.1005e-17$	$1.0117e-06$	$0.1297$	$11$
2	$0.2575$	$0.6160$	$0$	$3.1111e-06$	$0.1569$	$13$
3	$0.2575$	$0.9293$	$5.3204e-17$	$5.2445e-07$	$0.0974$	$8$
4	$0.2551$	$0.6160$	$0$	$3.1115e-06$	$0.1569$	$13$
5	$0.6991$	$0.2435$	$2.7166e-17$	$1.4652e-07$	$0.0550$	$2$

Table 1.

Change of the radius and the condition number A at sample points (x,y) using algorithm 1, for example 1.

In computing, $\delta =2r$ where $r=0.05$ and $c=\frac{2}{\sqrt{3}}r.$ Also In MMLS, ${\overline{w}}_{\nu }=0.1,\nu =\mathrm{1,2,3}.$ It should be noted that, these values were also selected experimentally. Relative errors of the MLS method for different Gauss-Legendre number points at $m=1,2$ and $3$ compared in Table 2, also investigating the proposed methods shown that increasing the number of numerical integration points does not guarantee the error decreases. jkjk.

	m = 1 ${‖e‖}_{\infty }$		m = 2 ${‖e‖}_{\infty }$		m = 3 ${‖e‖}_{\infty }$
N	$u_1$	$u_2$	$u_1$	$u_2$	$u_1$	$u_2$
5	$6.28\times {10}^{-4}$	$2.7\times {10}^{-3}$	$3.45\times {10}^{-7}$	$1.15\times {10}^{-6}$	$3.1\times {10}^{-6}$	$2.8\times {10}^{-6}$
10	$1.2\times {10}^{-4}$	$5.9\times {10}^{-4}$	$2.46\times {10}^{-7}$	$8.41\times {10}^{-7}$	$2.1\times {10}^{-7}$	$5.41\times {10}^{-7}$
15	$2.3\times {10}^{-4}$	$5.9\times {10}^{-4}$	$4.47\times {10}^{-8}$	$1.34\times {10}^{-7}$	$4.47\times {10}^{-8}$	$2.15\times {10}^{-7}$
20	$2.3\times {10}^{-4}$	$5.14\times {10}^{-4}$	$6.12\times {10}^{-7}$	$2.46\times {10}^{-6}$	$3.2\times {10}^{-7}$	$1.9\times {10}^{-6}$
30	$3.2\times {10}^{-4}$	$5.9\times {10}^{-4}$	$1.84\times {10}^{-6}$	$6.64\times {10}^{-6}$	$2.34\times {10}^{-6}$	$7.14\times {10}^{-6}$

Table 2.

Maximum relative errors for different points Gauss-Legendre quadrature rule $\delta =\mathit{2}r,$ for example 1.

In Table 3, we can see that the CPU times for solving the nonlinear system (65) are much larger in MMLS method; but, the errors are very smaller (Figure 3).

	MMLS ${‖e‖}_{\infty }$		MLS ${‖e‖}_{\infty }$		CPU times
N	$u_1$	$u_2$	$u_1$	$u_2$	MLS	MMLS
10	$3.78\times {10}^{-10}$	$8.13\times {10}^{-10}$	$2.46\times {10}^{-7}$	$8.41\times {10}^{-7}$	$389.9574$	$2.9776\times {10}^3$
15	$1.35\times {10}^{-10}$	$2.00\times {10}^{-11}$	$4.47\times {10}^{-8}$	$1.34\times {10}^{-7}$	$410.9083$	$2.1109\times {10}^3$
20	$8.63\times {10}^{-11}$	$2.51\times {10}^{-9}$	$6.12\times {10}^{-7}$	$2.46\times {10}^{-6}$	$634.8373$	$3.2115\times {10}^3$
30	$3.99\times {10}^{-10}$	$1.58\times {10}^{-9}$	$1.84\times {10}^{-6}$	$6.64\times {10}^{-6}$	$1.0331\times {10}^3$	$2.5844\times {10}^3$

Table 3.

Compare relative errors and CPU times of MLS and MMLS for different points Gauss-Legendre quadrature rule,$\left(m=\mathit{2}\right),$ for example 1.

Figure 3.

The determinant of a at a sample point $\mathit{p}$ and $\delta =\mathit{0.05}$ for example 1. Using algorithm 1.

6.2 Example 2

Consider the system of linear Fredholm integral equations with [27].

$$K\left(x,y,t,s\right)=\left(\begin{array}{cc}x\left(t+s\right)& -t\\ \mathit{ts}& \left(y+x\right)t\end{array}\right),$$E70

Such that $U\left(x,y\right)=\left(x+y,x\right)$ is the vector of The exact solutions and the vector of unknown function is $F\left(x,y\right)=\left(-\frac{1}{6}\left(x+y\right)+\frac{1}{3},\frac{4}{3}x-\frac{1}{3}+\frac{1}{3}y\right).$ Also the domain of the problem determine by $\mathrm{\Omega }=\left[0,1\right]\times \left[0,1\right].$ In this example, initial value of $r$ as radius of support domain set by $0.05$. Also Algorithm 1 is used for producing shape function at $m=\mathrm{1,2,3}$. In computing, we put ${\overline{w}}_{\nu }=0.1,\nu =\mathrm{1,2,3}.$

Table 4 shows relative errors and CPU times of MLS for different Gauss-Legendre number points at $m=1,3.$ As shown in Table 5, comparing the errors of MMLS and MLS method determines the capability and accuracy of the proposed technique to solve systems of linear Fredholm integral equations. This indicates the advantage of the proposed method over these systems of equations.

	m = 1			m = 3
N	$u_1$	$u_2$	CPU.T.	$u_1$	$u_2$	CPU.T.
5	$2.94\times {10}^{-4}$	$6.02\times {10}^{-4}$	$108.957$	$2.08\times {10}^{-4}$	$2.91\times {10}^{-4}$	$152.1474$
10	$1.79\times {10}^{-4}$	$3.89\times {10}^{-4}$	$159.258$	$1.7\times {10}^{-4}$	$2.37\times {10}^{-4}$	$170.7879$
15	$1.86\times {10}^{-4}$	$3.989\times {10}^{-4}$	$198.135$	$2.47\times {10}^{-4}$	$3.30\times {10}^{-4}$	$252.75424$
20	$1.86\times {10}^{-4}$	$3.986\times {10}^{-4}$	$221.321$	$2.41\times {10}^{-4}$	$3.29\times {10}^{-4}$	$247.0093$
30	$1.85\times {10}^{-4}$	$3.987\times {10}^{-4}$	$308.987$	$2.25\times {10}^{-4}$	$3.37\times {10}^{-4}$	$314.8173$
40	$1.85\times {10}^{-4}$	$3.980\times {10}^{-4}$	$395.125$	$1.87\times {10}^{-4}$	$2.86\times {10}^{-4}$	$402.9594$

Table 4.

Relative errors and CPU times of MLS for different points Gauss-Legendre quadrature rule at $m=\mathit{1},\mathit{3},$ for example 2.

		MLS			MMLS
m	N	$u_{\mathit{1}}$	$u_{\mathit{2}}$	CPU.T.	$u_{\mathit{1}}$	$u_{\mathit{2}}$	CPU.T.
2	5	$1.24\times {10}^{-7}$	$3.89\times {10}^{-6}$	$289.75$	$4.47\times {10}^{-7}$	$6.12\times {10}^{-6}$	$297.7$
	10	$3.16\times {10}^{-7}$	$5.23\times {10}^{-7}$	$189.99$	$2.16\times {10}^{-8}$	$6.26\times {10}^{-8}$	$210.09$
	15	$1.68\times {10}^{-7}$	$4.13\times {10}^{-7}$	$389.95$	$2.02\times {10}^{-8}$	$8.12\times {10}^{-8}$	$390.15$
	20	$1.61\times {10}^{-7}$	$3.88\times {10}^{-7}$	$410.90$	$2.04\times {10}^{-8}$	$5.24\times {10}^{-8}$	$425.87$
	30	$1.45\times {10}^{-7}$	$3.80\times {10}^{-7}$	$604.80$	$1.98\times {10}^{-8}$	$3.25\times {10}^{-8}$	$618.44$
	40	$1.44\times {10}^{-7}$	$3.84\times {10}^{-7}$	$689.9574$	$1.04\times {10}^{-8}$	$6.87\times {10}^{-8}$	$697.04$

Table 5.

Compare relative errors and CPU times of MLS and MMLS for different points Gauss-Legendre quadrature rule at $m=2,$ for example 2.

Comparing the errors of MMLS and MLS method determines the capability and accuracy of the proposed method to solve systems of linear Fredholm integral equations.

6.3 Example 3

The third example that we want to approximate is the system of linear Volterra-Fredholm integral equations with [27].

$$K\left(x,y,t,s\right)=\left(\begin{array}{cc}\left(x+y\right){exp}^{\left(t+s\right)}& \left(x+y\right){exp}^{\left(t+s\right)}\\ 1& -1\end{array}\right),$$E71

The domain is considered as $\mathrm{\Omega }=\left\{\left(x,y\right)\in {\mathbb{R}}^2:0\le x\le 1,0\le y\le \left(1-y\right)\right\}$ so that $y\in \left[0,1\right]$. Also the exact solutions are $\left({exp}^{x+y},{exp}^{x-y}\right).$ It is important to note that the linear transformation used in the experiment is only (55) and from (40) the kernel becomes

$$K\left(.,.,.,.\right)=\frac{y-c}{d-c}\mathrm{K}\left(.,.,.,.\right).$$E72

In this example, the effect of increasing radius of the domain of influence ${\delta }_i$ in MLS method on error has been investigated. Therefore the ${\delta }_i$ was considered as follows

$${\delta }_i=\mathit{\lambda r}i=1,2,\dots ,N$$E73

In this way, useful information is obtained about the performance of the proposed method. By investigating the results in Tables 6 and 7 we found that the relative error in MLS was also related to the radius of the domain of influence (i.e. $\delta =\mathit{\lambda r}$ so that $\lambda =\mathrm{3,5,7}$); however, it cannot be greater than $7.$ For example, the relative errors by choosing $\lambda =10$ (i.e. $\delta =0.05\lambda$) and $10-$point GaussLegendre quadrature rule are ${‖e‖}_{\infty u_1}=3.3\times {10}^{-2}$ and ${‖e‖}_{\infty u_2}=1.8\times {10}^{-2}$ at $m=1.$ Also, Table 8 depicts, the number of points in the numerical integration rule cannot be effective to increase the accuracy of the method.

	$\lambda =3{‖e‖}_{\infty }$		$\lambda =5{‖e‖}_{\infty }$		$\lambda =7{‖e‖}_{\infty }$
r	$u_1$	$u_2$	$u_1$	$u_2$	$u_1$	$u_2$
$0.2$	$9\times {10}^{-3}$	$5\times {10}^{-3}$	$1.9\times {10}^{-3}$	$1\times {10}^{-3}$	$6.02\times {10}^{-4}$	$3.36\times {10}^{-4}$
$0.1$	$1.2\times {10}^{-3}$	$6.87\times {10}^{-4}$	$1.34\times {10}^{-4}$	$7.46\times {10}^{-5}$	$4.57\times {10}^{-6}$	$2.91\times {10}^{-6}$
$0.05$	$1.44\times {10}^{-4}$	$8.14\times {10}^{-5}$	$2.23\times {10}^{-5}$	$1.26\times {10}^{-5}$	$6.27\times {10}^{-6}$	$3.6\times {10}^{-6}$

Table 6.

Maximum relative errors of MLS for $\mathit{10}$ Gauss-Legendre points and different values of $\delta =\mathit{\lambda r}$ at $m=\mathit{2},$ for example 3.

	$\lambda =3{‖e‖}_{\infty }$		$\lambda =5{‖e‖}_{\infty }$		$\lambda =7{‖e‖}_{\infty }$
r	$u_1$	$u_2$	$u_1$	$u_2$	$u_1$	$u_2$
$0.2$	$1.5\times {10}^{-3}$	$8.77\times {10}^{-3}$	$6.56\times {10}^{-5}$	$3.01\times {10}^{-5}$	$1.09\times {10}^{-4}$	$6.19\times {10}^{-5}$
$0.1$	$1.08\times {10}^{-4}$	$6.51\times {10}^{-5}$	$3.47\times {10}^{-5}$	$2.08\times {10}^{-5}$	$1.24\times {10}^{-5}$	$7.22\times {10}^{-6}$
$0.05$	$6.63\times {10}^{-6}$	$3.93\times {10}^{-6}$	$3.9\times {10}^{-6}$	$2.31\times {10}^{-6}$	$2.41\times {10}^{-6}$	$1.38\times {10}^{-6}$

Table 7.

Maximum relative errors of MLS for $\mathit{10}$ Gauss-Legendre points and different values of $\delta =\mathit{\lambda r}$ at $m=\mathrm{3}$ for example 3.

	$\lambda =\mathrm{3}{‖e‖}_{\infty }$		$\lambda =\mathrm{5}{‖e‖}_{\infty }$		$\lambda =\mathrm{7}{‖e‖}_{\infty }$
N	$u_{\mathrm{1}}$	$u_{\mathrm{2}}$	$u_{\mathrm{1}}$	$u_{\mathrm{2}}$	$u_{\mathrm{1}}$	$u_{\mathrm{2}}$
5	$\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{6.33}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{4.52}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{2.21}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{1.59}\times {\mathrm{10}}^{-\mathrm{5}}$
10	$\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{7.4}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{4.17}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{2.89}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{2.96}\times {\mathrm{10}}^{-\mathrm{6}}$	$\mathrm{2.72}\times {\mathrm{10}}^{-\mathrm{6}}$
15	$\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{7.9}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{5.07}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{6.31}\times {\mathrm{10}}^{-\mathrm{6}}$	$\mathrm{5.11}\times {\mathrm{10}}^{-\mathrm{6}}$
20	$\mathrm{2.1}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{8.06}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{5.17}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{3.96}\times {\mathrm{10}}^{-\mathrm{6}}$	$\mathrm{3.69}\times {\mathrm{10}}^{-\mathrm{6}}$
30	$\mathrm{2.1}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{8.16}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{5.17}\times {\mathrm{10}}^{-\mathrm{5}}$	$\mathrm{4.11}\times {\mathrm{10}}^{-\mathrm{6}}$	$\mathrm{3.74}\times {\mathrm{10}}^{-\mathrm{6}}$

Table 8.

Maximum relative errors of MLS for different values of $\delta =\mathit{\lambda r},r=\mathit{0.05}$ and Gauss-Legendre points at $m=\mathit{1},$ using algorithm 1 for example 3.

Then the relative error by MMLS shape function described in section (2.3) were obtained, using $\delta =7r$ and ${\overline{w}}_{\nu }=10$, such that $\nu =\mathrm{1,2,3}$ as weights of additional coefficients for MMLS. We can see that in Table 9 the errors of the system of linear Volterra-Fredholm integral equations are similar in both methods (i.e. MLS and MMLS methods).

	MLS		MMLS		CPU time
r	$u_1$	$u_2$	$u_1$	$u_2$	MLS	MMLS
$0.2$	$1.09\times {10}^{-4}$	$6.19\times {10}^{-5}$	$5.52\times {10}^{-4}$	$3.08\times {10}^{-4}$	$4.0556$	$3.3893$
$0.1$	$1.24\times {10}^{-5}$	$1.08\times {10}^{-6}$	$1.25\times {10}^{-5}$	$7.36\times {10}^{-6}$	$38.0616$	$31.7945$
$0.05$	$1.91\times {10}^{-6}$	$1.085\times {10}^{-6}$	$1.54\times {10}^{-5}$	$8.77\times {10}^{-6}$	$496.1746$	$409.0342$

Table 9.

Compare relative errors and CPU times of MLS and MMLS for ${\overline{w}}_{\nu }=\mathit{10}$ and $\mathit{10}$ Gauss-Legendre points and different values of $\delta =\mathit{7}r$ at $m=\mathit{2}$ for example 3.

6.4 Example 4

Consider the following nonlinear stiff systems of ODEs [38].

$$\left\{\begin{array}{l}{u}_1^{\prime }\left(t\right)=-1002u_1\left(t\right)+1000u_2^2\left(t\right)\\ u_2^{\prime }\left(t\right)=u_1\left(t\right)-u_2\left(t\right)-u_2^2\left(t\right)\end{array}\right.$$

With the initial condition $u_1\left(0\right)=1$ and $u_2\left(0\right)=1.$ The exact solution is

$$\begin{array}{l}{u}_1\left(t\right)=\mathit{exp}\left(-2t\right)\\ u_2\left(t\right)=\mathit{exp}\left(-t\right).\end{array}$$

In this numerical example, two scheme are compared and as explained the main task of the modified method tackle the singularity of the moment matrix. Table 10 presents the maximum relative error by MLS on a set of evaluation points (with $h=0.1\mathit{\text{and}}0.02$) and $\delta =\mathrm{4h}\mathit{\text{and}}3h$. Also in Table 11 MLS and MMLS at different number of nodes for $h=0.004$ and $\delta =\mathrm{5h}\mathit{\text{and}}8h$, were compared (Tables 10–12).

	m = 2,$\delta$ = 4 r		m = 2, $\delta$ = 3 r
r	$u_1$	$u_2$	$u_1$	$u_2$
0.1	$5\times {10}^{-3}$	$4.1\times {10}^{-4}$	$8.85\times {10}^{-4}$	$2.2\times {10}^{-3}$
0.02	$5.8\times {10}^{-2}$	$6.5\times {10}^{-5}$	$5.42\times {10}^{-4}$	$6.52\times {10}^{-5}$

Table 10.

Maximum relative errors by MLS, example 4.

	m = 3, $\delta$ = 5 r		m = 3, $\delta$ = 8 r
Type	$u_1$	$u_2$	$u_1$	$u_2$
MLS	$1.03\times {10}^0$	$0.98\times {10}^1$	$1.01\times {10}^0$	$9.2\times {10}^0$
MMLS	$9.23\times {10}^{-4}$	$9.22\times {10}^{-4}$	$6.89\times {10}^{-4}$	$6.96\times {10}^{-4}$

Table 11.

Maximum relative errors for h = 0.004 by MMLS and MLS, example 4.

	MLS			MMLS
weight type	$u_1$	$u_2$	Cputime	$u_1$	$u_2$	Cputime
Guass	$3.06\times {10}^{-3}$	$9.92\times {10}^{-5}$	$61.1706$	$8.5\times {10}^{-4}$	$6.5\times {10}^{-4}$	$0.5598$
Spline	$5.06\times {10}^{-4}$	$5.3\times {10}^{-4}$	$64.5897$	$1.93\times {10}^{-2}$	$4.23\times {10}^{-4}$	$0.6714$
RBF	$6.407\times {10}^{-4}$	$3.02\times {10}^{-4}$	$59.1790$	$6.9\times {10}^{-3}$	$6.9\times {10}^{-3}$	$0.5768$

Table 12.

Maximum relative errors by MLS $t\in \left[\mathit{0},\mathit{5}\right]$,$h=\mathit{0.004}$,, example 2.

6.5 Example 5

In this example, we consider $U\left(t\right)=\left(\frac{1}{47}\left(95{exp}^{\left(-2t\right)}-48\mathit{exp}\left(-96t\right)\right),\frac{1}{47}\left(48exp\left(-96t\right)-exp\left(-2t\right)\right)\right)$ as the exact solution and $U\left(0,0\right)=\left(1,1\right)$ as the initial conditions for the following system of ODE,

$$\left\{\begin{array}{c}{x}^{\prime }\left(t\right)=-x\left(t\right)+95y\left(t\right)\\ y^{\prime }\left(t\right)=-x\left(t\right)-97y\left(t\right)\end{array}\right.$$

Table 12 presents the maximum relative norm of the errors on a fine set of evaluation points (with $h=0.004$) and $\delta =5h$ for MLS and MMLS at different type of weight functions. As seen in this table, one major advantage of MMLS is that the computational time used by MMLS is less than MLS.

7. Conclusion

In this paper, two meshless techniques called moving least squares and modified Moving least-squares approximation are applied for solving the system of functional equations. Comparing the results obtained by these methods with the results obtained by the exact solution shows that the moving least squares methods are the reliable and accurate methods for solving a system of functional equations. Meshless methods are free from choosing the domain and this makes it suitable to study real-world problems. Also, the modified algorithm has changed the ability to select the support range radius In fact, the user can begin to solve any problem with an arbitrary radius from the domain and the proposed algorithm can correct it during execution.