Meshless Analysis of Nonlocal Boundary Value Problems in Anisotropic and Inhomogeneous Media

In this work, meshless methods based on a radial basis function (RBF) are applied for the solution of two-dimensional steady-state heat conduction problems with nonlocal multi-point boundary conditions (NMBC). These meshless procedures are based on the multiquadric (MQ) RBF and its modified version (i.e., integrated MQ RBF). The meshless method is extended to the NMBC and compared with the standard collocation method (i.e., Kansa’s method). In extended methods, the interior and the boundary solutions are approximated with a sum of RBF series, while in Kansa’s method, a single series of RBF is considered. Three different sorts of solution domains are considered in which Dirichlet or Neumann boundary conditions are specified on some part of the boundary while the unknown function values of the remaining portion of the boundary are related to a discrete set of interior points. The influences of NMBC on the accuracy and condition number of the system matrix associated with the proposed methods are investigated. The sensitivity of the shape parameter is also analyzed in the proposed methods. The performance of the proposed approaches in terms of accuracy and efficiency is confirmed on the benchmark problems.


Introduction
The partial differential equations with nonlocal boundary conditions that emerged in the literature have significant applications in the fields of engineering, astrophysics, and biology. Some of the first nonlocal equations that appeared in the literature are encountered in the field of phase transition and are related to theories due to Van der Waals, Ginzburg, Landau, and Cahn & Hilliard [1]. For instance, such models with nonlocal spatial terms are encountered in the Ohmic heating production [2], in the shear banding formation in metals being deformed under high strain rates [3,4], in the theory of gravitational equilibrium of polytropic stars [5], in the investigation of the fully turbulent behavior of real flows, using invariant measures for the Euler equation [6], in population dynamics [7], in modeling aggregation of cells via interaction with a chemical substance (chemotaxis) [8], and in image processing [9]. Therefore, nonlocality is not a technical obstruction to scientific research but it provides the essence of what happens in reality, and for that purpose, its mathematical study can provide very useful predictions in many areas of applications. In general, nonlocal models provide more accurate predictions compared to their local counterparts since they use all the available information regarding the evolution of the inspected process.
Nonlocal boundary conditions have significant applications in engineering and biological sciences. These nonlocal boundary conditions arise in mathematical models of various biological, physical, and chemical processes. In addition, nonlocal boundary conditions also have wide applications in heat conduction processes, control theory, chemical diffusion, population dynamics, blood flow problems, underground water flow, and chemical engineering. Related references and examples of mathematical models involving nonlocal boundary conditions can be found in [10][11][12][13][14][15]. Nonlocal boundary value problems arise when the boundary values of the solution are not known on the boundaries. In such circumstances, the solution at the boundaries is connected with the solution inside the given domain. This type of formulation is placed in a separate class known as nonlocal boundary value problems.
It is sometimes better to impose nonlocal boundary conditions since the measurements needed by a nonlocal boundary condition may be more precise than the measurement given by a local condition. For example, consider the second order differential equation [12] u (t) = g(t, u(t), u (t)), (1) subject to local boundary conditions Equation (2) is called a mixed type boundary condition because on t = 0, u is prescribed and the derivative is prescribed on t = 1. Since we know that, in the process of numerical computation and a scientific experiment, the value of u (1) is more difficult to determine than that of u(ξ)−u (1) ξ−1 where (0 < ξ < 1); therefore, by replacing u (1) in (2) with u(ξ)−u (1) ξ−1 we have nonlocal boundary conditions: Equation (1) along with (3) is called a nonlocal boundary value problem. Equation (3) is called two-point boundary condition which is a special type of multi-point boundary condition.
Meshless methods for the solutions of differential equations are being successfully applied in various fields of science and engineering; therefore, they are more effective and suitable for complex domains as compared to mesh-based numerical methods. Since a set of independent points is required for meshless methods, they are often better suited for complex geometries. The computational costs associated with mesh generation in mesh-based methods are high [16] as compared to meshless methods. In comparison to meshless methods, other traditional numerical methods such as finite difference methods, finite element methods, and finite volume methods are usually limited to problems involving two or three spacial variables (space dimensions).
Literature suggests that multiquadric (MQ) RBF is considered the best in accuracy among different types of RBFs; however, the shape of the function is controlled by shape parameter . The effects the accuracy of the method where its optimal value can lead to good approximations. To determine the optimal value of the is an open problem and there is no such mathematical theory developed for optimal value of . Another drawback is the high condition number of the interpolation matrix. The ill-conditional system matrix arises in the case of globally supported RBFs, as the separation distance 1 2 max i =j x i − x j 2 , i, j = 1, 2, 3, . . . , N, between data points decreases for fixed values of the shape parameter [17]. Many prominent researchers have used a different value of the shape parameter in the hope of accurately handling such problems. Both constant and variable shape parameters are recommended. Hardy introduced the following shape parameter formula in [18]: where N is number of collocation points and d j is the distance from the jth center to its nearest neighbor. Later on, a smallest circle enclosing all data points is considered and the shape parameters ( = 1.25 √ N × diameter of this circle) were recommended in [19]. The shape parameter 2 √ N was used for solution of nonlinear PDEs in [20]. In [21], two optimal values of the shape parameter 1.03 and 1.42 were recommended for MQ RBF, whereas for Gaussian basis, a range 0.003-0.03 was given, for solving cantilever beam and perforated strip plate problems. In [11,22,23], it is concluded that splitting of shape parameter into x-and y-axes give good results, i.e., x = h x N 2 x and y = h y N 2 y , where h x and h y are positive constants, but the condition number of the collocation matrix is high. Different selection procedures of shape parameter can also be found in [24][25][26][27][28][29][30][31][32].
In this paper, a meshless discretization technique based on MQ RBF and its modified form, i.e., integrated MQ RBF, is applied for the numerical solution of elliptic PDEs with nonlocal boundary conditions (NBC). In the last few decades, the RBF-based collocation technique for the solution of nonlocal boundary value problems has been a vital research area in many branches of science and hence was successfully applied to solve such problems numerically [33][34][35]. A one-dimensional hyperbolic equation with an integral boundary condition was studied in [36]. The two-dimensional diffusion equation subject to a nonlocal condition involving a double integral along with Neumann's boundary conditions was formulated in a rectangular region. Such a problem is solved numerically using a meshless collocation method [37]. In paper [38], the RBF-based collocation technique was also applied to investigate the influences of nonlocal conditions on the optimal selection of shape parameter , ill-conditioning, and accuracy of the method. A nonlocal, two-dimensional Poisson equation was solved with two-point and integral boundary conditions; it was formulated on one side of the rectangular region. Later, in [11], a comparison of the RBF-based collocation technique with the Haar wavelet-based collocation technique was studied for the same problem. Due to the splitting of the shape parameter into x and y directions, its distribution effect on the accuracy of the problem at hand was also investigated. The influences of shape parameter and the distribution of nodes on the accuracy of the RBF-based collocation technique were investigated by considering the multidimensional linear elliptic equation with integral conditions [39]. Recently, the RBF-based meshless method has been applied to elliptic PDEs with nonlocal multi-point boundary conditions (NMBC) [40].
It is known that the meshless method can be successfully applied to solve such nonlocal problems. However, it was shown (e.g., see paper [11,38,40]) that due to the nonlocal parameter γ, the condition number κ of the system matrix involved and accuracy of the method are affected. It is shown that the RBF-based collocation technique gives quite good results for the optimal value of shape parameter [38,40]. On the other hand, for the optimal value of the shape parameter, the accuracy and conditioning are also slightly affected by NBC. This paper aims to investigate the effect of this parameter on the property of the RBF-based method using integrated MQ RBF.
Nowadays, integrated RBF(s) are successfully implemented in solving miscellaneous differential equations. An MQ RBF integrated twice with a fixed value of shape parameter was used by [52] to solve one-dimensional problems. Elliptic problems were solved by [53] through a multi-domain integrated RBF collocation method. Higher-order ODEs and PDEs were solved by an integrated RBF method in the papers [54,55]. Further detail about integrated RBFs can be found in [17,[56][57][58].
We have implemented four different methods based on radial basis functions, i.e., the Kansa method (multi-quadric is a basis function) and the Kansa method (integrated-multi quadric is a basis function), and we extended the method proposed in [59] (with these two different basis functions) to the complicated nonlocal boundary value problems. For the first time, these methods were applied to the nonlocal Neumann boundary value problems accurately, which can be very complicated in the case of finite difference and finite element methods. An accurate numerical approximation for the unknown boundary function has been obtained through these techniques.
The paper is distributed into five sections, including the present introductory Section 1. In Section 2, we discuss the governing equations along with different types of boundary conditions. In Section 3, the numerical procedure is described. Numerical experiments are given in Section 4. Some conclusions are drawn in Section 5.

Governing Equation
A general form of the diffusion equation is given by [59] The above model is subject to mixed boundary conditions: • Classical Dirichlet boundary conditions: • Nonlocal multipoint boundary conditions: where Ω ⊂ Ω are set of discrete distinct scattered points inside the solution domain and n is the number of these scattered points. Throughout the paper we will take n = 10; if not, it will be mentioned. It is also assumed that the distribution of the pointsv l on the solution domain will be fixed. The given functions are k 1,1 , k 1,2 , k 2,2 , f , g, h, and the parameter γ. Ω is the interior portion of the solution domain and Γ = Γ 1 ∪ Γ 2 are the boundary parts of the solution domain such that Γ 1 ∩ Γ 2 = ∅ and Γ 2 = ∅. The NBC defined on the boundary Γ 2 is reduced to Dirichlet boundary condition when γ = 0.

The Numerical Scheme
Let us divide the solution domain into three disjoint sets. The set of inner distinct scattered nodes: and the set of distinct scattered boundary points: such thatΩ ⊂ Ω,Γ 1 ⊂ Γ 1 , andΓ 2 ⊂ Γ 2 . N is the total number of data points inside the solution domain and M = M 1 + M 2 is the total number of boundary points on the boundary Γ. We also assume that the pointsv are not included inΩ, i.e., Ω ⊂ Ω, but Ω Ω . Furthermore, throughout the paper a given RBF inside the solution domain is denoted by where the centersv i = (x i , y i ), i = 1, 2, . . . , N are distributed inside the solution domain. Similarly, if the centersv k = (x k , y k ), k = 1, 2, . . . , M are distributed on the boundary Γ, then the RBF on the boundary is given by Recently, an RBF-based method was used to solve two-dimensional, steady-state PDEs with three different types of boundary conditions (Dirichlet, Neumann, and Robin boundary conditions) in anisotropic and inhomogeneous media [59]. In [59], an approximate solution was sought in the form where basis u g (v) is a smooth function which satisfies the boundary conditions, and ψ i (v) is a sum of the given RBF. ω i (v), called the correcting function, is chosen such that: It is to be noted that in paper [59], a trigonometric basis, a polynomial basis, and an RBF basis were used to find the functions u g (v) and the correcting function ω i (v). However, we use here just MQ RBF and integrated MQ RBF to find these functions.
This algorithm strictly divides the approximation of boundary conditions and the approximation of the differential equation inside the solution domain in cases of Dirichlet, Neumann, and Robin boundary conditions (for more detail see [59]). However, in the case of NBC the values of the solution u on the boundary part Γ 2 are related to the values at the interior points. Thus, boundary conditions cannot be approximated separately because the values of the solution u at the boundary Γ 2 are unknown.
Equation (12) can be rewritten in the form where The correcting function ω i (v) is determined from condition (13). The unknown parameters β k , k = 1, 2, 3, . . . , M, and λ i , i = 1, 2, 3, . . . , N are determined by substituting (14) in the governing Equation (4) and collocating at the interior and boundary points of the solution domain; i.e., The above system (15) is solved for the unknowns β k , k = 1, 2, 3, . . . , M, and λ i , i = 1, 2, 3, . . . , N. The approximate solution u is obtained by substituting these unknowns into (14). The numerical implementation of this method contains two stages. At first stage, the correcting function ω i (v) is approximated. At the second stage, the correcting function ω i (v) is used to solve system (15). When φ(v) is MQ RBF, we named the method RCM1, and when it integrates MQ RBF, the method is named IRCM1.
The standard collocation method known as Kansa's approach is explained here. In Kansa's approach, the approximate solution is sought in the form where the centers are distributed inside the solution domain and on the boundary, i.e.,v i ∈Ω ∪Γ 1 ∪Γ 2 . By substituting (16) in (4), (5), and (6), we have This system (17) of linear equations is solved for λ i , (i = 1, 2, 3, . . . , N + M) to get the approximate solution (16) of (4). In this procedure when φ(v) is MQ RBF, the method is named RCM2, and when it integrates MQ RBF, the method is named IRCM2.

Integrated MQ RBF
It is shown in the numerical section that integrated RBF methods give quite a good result as compared to standard non-integrated RBF methods. In order to achieve new basis for numerical solution of differential equations, RBF is integrated six times with respect to r. In the present paper, we will use MQ RBF as where r is the Euclidean norm.
Six-time integration of MQ RBF (18) with respect to r, gives the basis functions [17] The parameter o in (18) is a shape parameter but the symbol is changed to distinguish it from the shape parameter in (19). In our study o will represent the shape parameter for MQ and will be denoted the shape parameter for IMQ6. The MQ RBF can be integrated easily using a computer algebra system. In paper [17], the first four members of the integrated RBF family based on MQ RBF were listed. The benefit of IMQ6 is that it is not very sensitive to the shape parameter as compared to MQ. The IMQ6 basis function gives quite good results for a wide range of , but to find optimal is still challenge in this case as well. In case of IMQ6 basis, the value of is taken = 80 in our study. However, this value is not the optimal value of the shape parameter. In the case of MQ basis, we will take the value of o = 2; otherwise it will be mentioned during the investigation.

Numerical Experiments
In this section, the numerical results and the effect of nonlocal boundary condition (6) are demonstrated for (4). We discuss here different types of RBF-based collocation methods, i.e., RCM1, IRCM1, RCM2, and IRCM2. The effects of the shape parameters on the accuracy and κ of the system matrices of the proposed methods are also investigated. This section shows that the proposed methods IRCM1 and IRCM2 are less sensitive to the selection of shape parameters as compared to RCM1 and RCM2.
Various error measures are used to estimate the accuracy of the numerical method. Among them, one is L rms error norm which is defined as where u(v) and u(v) represent the exact and the approximate solutions respectively. One can also check the numerical stability of the method by κ of the system matrix involved. A well conditioned system matrix has a small value of κ, while a system with larger value of κ indicates that it is ill-conditioned. Problem 1. The partial differential equation is given in (4) with constant coefficients k 1,1 (v) = −1 and k 1,2 (v) = 0, and k 2,2 (v) = −1 reduces to a Poisson equation and the one studied in reference [40] with the same nonlocal multipoint boundary condition. We consider the same problem with the same analytical solution The domain is a unit rectangle having a cut of the radius R = 0.5, which is depicted in Figure 1 (left). The functions f , g, and h are calculated according to the exact solution (21).
The representation of the solution domain and distribution of the nodes are depicted in Figure 1. In Figure 1 (right), the distinct scattered points marked with black dot represent an interior portion of the solution domain. The regularly distributed points marked with red dots represent the boundary parts Γ 2 , and the regularly distributed points marked with green dots represent the boundary parts Γ 1 . The scattered points marked with red asterisk are the points of the set Ω . In case of Test Problem (1), the effects of the shape parameters ( or o ), the parameter γ, and the nodes are shown in Figures 2 and 3 on the proposed methods using IMQ6 and MQ RBF as basis functions. From these figures, we can see that the methods based on IMQ6 RBF are better in accuracy and conditioning as compared to the methods based on MQ RBF. However, when the shape parameter is large, κ of the methods based on MQ RBF is small but the accuracy of the methods is low. The methods RCM1 and RCM2 also give quite good results for some value of o (see Figure 3). The proposed methods IRCM1 and IRCM2 give reasonable accuracy for a wide range of (see Figure 2). Figure 3 (right) shows that the RCM2 is more accurate than the RCM1 for dense grid and o = 2. However, κ of RCM2 is worse than RCM1. Figure 2 indicates that the IRCM1 is slightly better at conditioning than IRCM2. However, IRCM2 are much more accurate for a small value of . When the value of increases, the error norms of the methods get closer to each other (see Figure 2 (left)). The conditioning performances of both methods IRCM1 and IRCM2 become better, up to certain limit when the value of increases. After this limit, it does not much affect the κ of the coefficient matrix. The behavior of IMQ6 RBF is the same as the MQ RBF, which gives quite good results for a small value of shape parameter along with high κ of the coefficient matrix. Figures 2 (left) and 3 (left) also indicate that accuracy of the IRCM1 and IRCM2 is better for a very wide range of the value of as compared to the RCM1 and RCM2.
The effect of the collocation nodes on the accuracy of the proposed methods based on IMQ6 basis is also investigated. The value of the shape parameter is considered = 80. The methods IRCM1 and IRCM2 give quite good results. As the nodal points increase, we see good agreement between solutions of IRCM1 and IRCM2 and the exact one. However, κ increases for dense nodes while keeping the fixed value of the shape parameter = 80.
The NBC (6) is defined comparatively on a small part of the boundary. However, from Figures 2 and 3 (middle), we see that the NBC has a fairly strong negative influence on the κ of the coefficient matrices corresponding to the methods RCM1, RCM2, IRCM1, and IRCM2. However, this influence is weak in case of IRCM1 and IRCM2 as compared to the methods RCM1 and RCM2. From Figure 2 (middle), we can observe that both the methods IRCM1 and IRCM2 give the same lowest and highest values of error norm. The highest and lowest error norms are 6.222 × 10 −7 and 6.184 × 10 −7 corresponding to γ = 2 and γ = 0. The influence of the NBC (6) is also weak on the κ of the coefficient matrices corresponding to the methods IRCM1 and IRCM2 as compared to the methods RCM1 and RCM2. When o = 2 and = 80, the methods RCM1 and RCM2 are the most ill-conditioned, and the absolute value of γ increases as compared to the methods IRCM1 and IRCM2.  (4) with variables coefficients k 1,1 (v) = 5(1 + x 2 + y 2 ), k 1,2 (v) = 1 + xy 2 and k 2,2 (v) = 3(1 + x 2 + y 2 ). The functions f , g, and h are chosen such that the exact solution gets the form [59]: u(v) = cos(y) exp(x − y).

Problem 2. Consider the equation given in
Here in this example, we have considered an irregular-shaped domain defined by where p(θ) is the radius of the domain which depends on the polar angle θ. The representation of the domain for Test Problem 2 is depicted in Figure 4. In Figure 4 (left), the boundary parts Γ 2 (0 ≤ θ ≤ 2π/3) and Γ 1 (2π/3 < θ < 2π) are represented by red dashed and green solid lines respectively. The influences of the shape parameters, the parameter γ, and nodal points on the accuracy and conditioning of the methods are shown in Figures 5 and 6. The accuracy of IRCM2 is good as compared to the method IRCM1 for different values of shape parameter , as shown in Figure 5 (left). Both the methods are most poorly ill-conditioned for a small value of shape parameter , while the method IRCM2 gives high accuracy for a smaller value of than IRMC1. It can be observed from these figures that the accuracy of the methods and κ of the coefficient matrices for a fixed value of increases when the number of nodes increases. In contrast to the methods based on MQ RBF, the κ of the coefficient matrices of IRCM1 and IRCM2 is small. This is because of which is large and is equal to 80. The system matrices of the proposed methods are ill-conditioned for a small value of . In addition, the effect of the parameter γ on the methods IRCM1 and IRCM2 is also small. The method IRCM2 gives somewhat accurate results for all values of γ, even though the κ of the coefficient matrix of IRCM1 is high for some values of γ as compared to the method IRCM2 (see Figure 5 (middle)).
For dense nodes, RCM2 gives quite good results for o = 2 as compared to RCM1, as shown in Figure 6. However, the results obtained from IRCM1 and IRCM2 and the conditioning of the methods are much better there. Again in this example, the NBC (6) has a strong negative influence on the conditioning of the methods based on MQ RBF. From Figure 6, we see a fluctuation in κ of the coefficient matrices associated with the methods RCM1 and RCM2 and in error norms as the value of γ changes. This fluctuation in error norms occurs due to the small value of o . Problem 3. Reconsider the differential equation given in (4) with constant coefficients k 1,1 (v) = 3, k 1,2 (v) = 1.5, and k 2,2 (v) = 2. The functions f , g, and h are chosen such that the exact solution gets the form [59] The representation of the domain for Test Problem 3 is depicted in Figure 7. The curve Γ 2 ∪ Γ 1 is called Cassini curve and is defined by parametric equation [59]: and rest of the terms have the same meanings as defined in Test Problem 2. We consider the Neumann boundary condition on the boundary part Γ 1 . The normal component of the heat flux on the boundary is of the form: where q(v) is some given function. n x and n y are the components of unit outward normal vectors. Since it was shown in the paper [59] that an RBF approximation of the mixed connected boundary conditions, i.e., Dirichlet boundary condition, is specified on some part of the boundary and Neumann condition, it is defined on the rest of boundary, and is too low in accuracy. If we specified Neumann boundary condition on boundary part Γ 1 , then the accuracy of the boundary condition would be too low. Since we used RBF approximation (i.e., using MQ or IMQ6 RBF) to approximate the mixed connected boundary condition, i.e., the multi-point and Neumann boundary condition, and the solution should be inside the domain. Hence, obviously, it would also affect the approximate solution inside the solution domain. To avoid this problem, we used the same procedure as was done above for the Dirichlet boundary condition but replaced the system of equation obtained from boundary part The influences of shape parameters ( or o ), the parameter γ, and the number of nodes for Test problem 3 are shown in Figures 8 and 9. It can be seen from Figures 8 and 9 (right) that the proposed methods give good accuracy for dense grid whether the RBF basis is MQ or IMQ6. However, the κ of the coefficient matrices associated with the methods RCM1 and RCM2 are high. From Figure 8 (top middle), we see that when γ = 0 the error norms of the methods IRCM1 and IRCM2 are 1.175 × 10 −6 and 1.226 × 10 −6 respectively. Now we check the influence of the parameter γ on the methods IRCM1 and IRCM2 when the shape parameter is small. The value of the shape parameter = 2 is taken for each problem. We see from Figure 10 that IRCM2 give quite a good accuracy for a small value of shape parameter in all Test Problems 1-3. The accuracy is good when the value of γ is zero or near to zero. It can be observed from this figure that the κ and the accuracies of the proposed methods deteriorate when the absolute value of γ increases. In contrast, when the value of shape parameter is high (i.e., = 80) the κ deteriorates but the accuracy does not affect much when |γ| is increased (see middle of Figures 2,5 and 8). The absolute error of IRCM1 and IRCM2 of the unknown solution on the boundary part Γ 2 is shown in Figure 11. From this figure, we observed that both the methods IRCM1 and IRCM2 give quite good accuracy. Accuracy in terms of error norm and conditioning of the coefficient matrices corresponding to the proposed methods is shown in Figure 12. The values of the parameters are = 150, o = 2, and γ = 20. We observed from this figure that the proposed methods based on integrated MQ RBF performed well both in terms of accuracy and conditioning. We have also seen that different types of boundary conditions (i.e., Dirichlet and Neumann conditions) do not much affect the conditioning and accuracy of the proposed methods IRCM1 and IRCM2.

Conclusions
The collocation methods based on MQ and IMQ6 were set up for the numerical solution of elliptic NMBC with variable coefficients. The results of the propsed extended methods were compared with Kansa's approach. The method based on IMQ6 RBF is well conditioned as compared to the collocation method based on MQ RBF (see Figure 11) for a large value of the shape parameter . However, the methods IRCM1, IRCM2, RCM1, and RCM2 became ill-conditioned for small values of shape parameter (see Figures 2,3,5,6,8,and 9) but the results are reasonably accurate. This approach can be extended to a higher dimension easily. In addition, the proposed approach can be applied to time-dependent PDEs with nonlocal multi-point boundary conditions. suggestions which are fully incorporated in the revised version. As a result, the quality of the paper has been improved considerably. We appreciate their time and effort. Moreover, The first author appreciates the support provided by CECOS University of IT and Emerging Sciences Peshawar Pakistan.

Conflicts of Interest:
The authors declare no conflict of interest. Kansa's radial basis function based on mulitquadric IRCM2

Nomenclature
Kansa's radial basis function based on integrated Mulitquadric RCM1 radial basis function based on mulitquadric IRCM1 radial basis function based on integrated Mulitquadric