A Matlab software for approximate solution of 2D elliptic problems by means of the meshless Monte Carlo random walk method

This paper is devoted to the development of an innovative Matlab software, dedicated to the numerical analysis of two-dimensional elliptic problems, by means of the probabilistic approach. This approach combines features of the Monte Carlo random walk method with discretization and approximation techniques, typical for meshless methods. It allows for determination of an approximate solution of elliptic equations at the specified point (or group of points), without a necessity to generate large system of equations for the entire problem domain. While the procedure is simple and fast, the final solution may suffer from both stochastic and discretization errors. The attached Matlab software is based on several original author’s concepts. It permits the use of arbitrarily irregular clouds of nodes, non-homogeneous right-hand side functions, mixed type of boundary conditions as well as variable material coefficients (of anisotropic materials). The paper is illustrated with results of analysis of selected elliptic problems, obtained by means of this software.

electrical potential or filtration through porous media). Two different types of mathematical models may be investigated, depending on the problem variables' nature, namely deterministic and probabilistic ones. In deterministic models, a particular state (either continuous or discrete) is uniquely assigned to input data and no element of randomness occurs. In other words, all input parameters are represented by their mean values and one fixed data set corresponds to one unique solution. In most cases, the specified numerical solution approach has to be incorporated, as the analytical solution is limited to problems with simple geometry and uniformly distributed subjected load. Therefore, such problems ascribed by deterministic models are usually analyzed by means of the standard numerical approaches (hard computing methods, [22]), which are based upon appropriate domain discretization (partition into nodes and/or elements) and unknown function approximation, built upon those geometrical quantities (e.g., finite element method (FEM, [37]), finite difference method (FDM, [31,32]), meshless methods (MM, [29,31,33,35,45]), or boundary element method (BEM)).
On the other hand, the probabilistic class of mathematical models specifies governing differential equations that combine one or more random variables (with assigned probability distribution function) with other non-random variables. Regardless of the applied solution approach, taking advantage from the problem's stochastic nature, one deals with the analysis of multiple standard problems, varying with input parameters' values. These auxiliary problems are usually investigated by means of the standard aforementioned deterministic methods. As a consequence, the entire family of solutions is obtained, out of which the optimal solution is selected, using appropriate averaging criteria. Among the most commonly applied probabilistic approaches (soft computing methods), which rule the principles of a random selection of the input data and averaging techniques, one may distinguish the following: -Fuzzy sets (FS, [36]) analysis, in which the fuzzy models are applied to selection of input data, with assigned membership function, characterizing data randomness. -Genetic and evolutionary algorithms (GA and EA, [24]), inspired by the biological evolution and dedicated to non-gradient analysis of multidimensional optimization problems. -Artificial neural networks (ANN, [39]), inspired by biological nervous system, tailored for analysis of variety of problems by means of a-priori prepared patterns. -Monte Carlo stochastic approach, in which series of simulations (trials), representing the analyzed problem, with randomly selected input values, are performed. Among these trials, the specified number of properly defined successes is achieved. The ratio between the number of success trials to the number of all trials, scaled by dimensional quantity (e.g., area, function value) allows for the estimation of the unknown solution, providing the number of trials is large enough.
The first mature formulation as well as technical applications of Monte Carlo (MC) concept were derived by S. Ulam and J. von Neumann in their pioneering research [1,21], devoted to the neutron transport in an atomic bomb framework. It was followed by a variety of works delivered by many other researchers, who applied this simple concept to analysis of algebraic and differential problems, like eigenvalues estimation [3], solution of linear systems of equations [2,4,5], numerical integration in multidimensional spaces [13,26] as well as numerical solution of Laplace differential equations [6,7,9,15,17,20,23,27,28,30], at the selected internal point of the problem domain. This application is discussed in more detailed manner in this paper.
The main idea is based upon the classic MC concept combined with the so-called fixed random walk (RW, [6,7,9,42,53]) technique. This technique allows for a random selection of the path consisting of the horizontal and vertical walks. Those walks are performed on the regular grid of lines, starting from the given internal point (thus forming a Markov chain). This path has to be terminated at the domain boundary and the entire procedure is repeated as long as all random walks are completed. Therefore, the net of connections is generated, linking the selected internal point of interest with all points located on the boundary. The total sum of all numbers of boundary indications (boundary hits), scaled by their boundary values, estimates the solution of the Laplace equation at this particular point. It was proved [7], that this estimation corresponds to the finite difference solution, obtained on the same regular grid of nodes. This simple solution approach was further developed by many investigators in subsequent years. Apart from the standard RW, in which both step size and move directions are fixed and preassigned, two main random walk types were derived, namely semi floating random walk (step size is fixed, though the move direction is not limited, [47,52]) and full floating random walk (step size is not preassigned and changes at each step, whereas the move direction is not limited, [8,42,47,52]). Moreover, a continuous random walk procedure, minimizing the solution error has been proposed [17] along with the self-adaptive, grid-free algorithm, with improved solution smoothness and its application to diffusion equations [20]. State-of-the-art as well as recent developments (up to date) may be found in [23]. Further reading may include the following subjects: effective reduction of random walk steps' number [28], weighted version of MC method [30], MC treatment of complex 2D geometries [38]. Within the last several years, investigated were, for instance, reduction of MC error by new probabilistic sampling [41], development of stochastic limit theory and governing equations of continuous time random walk [44,50], improvement of the simulation efficiency for radiative transfer problems with strongly frequency-dependent opacities [51], as well as the analysis of heat conduction equations [52].
In this research, special emphasis is laid upon the development of the meshless random walk technique, combined with the Monte Carlo approach, towards analysis of wider class of 2D elliptic problems. While the original MC/RW technique (Section 2) may be applied to Laplace equation with regular grid of lines, essential boundary conditions as well as homogeneous right-hand side function only, the proposed meshless version of MC/RW technique allows for analysis of 2D elliptic equation in more general form, including the following: -More complex geometries, discretized by means of arbitrarily irregular clouds of nodes, without any imposed structure (therefore neither mesh size or move directions are fixed) -Non-homogeneous right-hand side functions (typical for equations of Poisson type) -Mixed type of boundary conditions (both essential and natural types may be taken into the account) -Non-constant material coefficients Similar stochastic concepts were successfully applied in variety of methods based upon the reduction of the solution of multidimensional problems for partial differential equations to integration of stochastic equations [43,49]. Especially, the well-known Feynman-Kac formula [14,34] is worth mentioning. It may be interpreted as a method for evaluating functional integrals of a specified continuous form. However, its practical application is based upon the stochastic approximation of these integrals [49], which incorporates both Monte Carlo approach and random walk procedure. Although the entire numerical procedure is similar to the one presented here, the mathematical formulation of the considered problem is more complex, since it requires determination of the explicit closed form solution formula. Moreover, the right-hand side function of the differential equation has to be integrated, which limits the number of potential applications to cases with smooth solution only (i.e., no concentrated loads). The meshless MC/RW, developed in [53] and implemented in Matlab package in this paper, uses the original deterministic differential equation as well as the computational framework which is typical for elementbased and difference methods (i.e., existence of nodes (though with no imposed structure), approximation schemes, solution smoothing). However, those aspects are incorporated into the pure stochastic approach.
All innovative ideas and concepts, proposed and introduced in [53], are briefly presented and discussed here (Section 3). They are based upon the incorporation of selected discretization and approximation schemes, typical for meshless methods (MM). Especially, features of the moving weighted least squares approximation (MWLS) are applied in order to reformulate classic random walk principles. However, the main part of this paper is devoted to the presentation of the innovative Matlab software (Section 4), based upon those concepts. It may work in several various modes, for instance, it allows either for a determination of an approximated solution of an elliptic equation at selected point for fixed number of random walks, or for an examination of a solution convergence at selected point(s) with number of random walks increasing. Finally, determination of the solution at all nodes is possible, with additional smoothing of rough MC results. Moreover, for illustration purposes, full graphical interface has been designed, including visualization of construction of single random walks, solution convergence of N-type (i.e., with respect to the number of random walks), graphs of the final MC/RW solution and its first derivatives as well as the simplified FD versus MC/RW error analysis [40]. By default, a user works with rectangular, circular or triangular domain and irregular meshes/clouds of nodes, though any other discretization may be loaded from the disc file. Several examples of clouds of nodes corresponding to more complex geometries are attached as well. Results of their analysis, conducted by means of the attached Matlab software, are presented in Section 5. The paper is briefly concluded and directions of future work are mentioned.

Standard Monte Carlo random walk technique
The following 2D Laplace equation is considered as follows: with essential boundary conditions as follows: where = {(x, y)} ∈ R 2 is the problem domain, ∂ -its boundary, F = F (x, y) : → R C 2 is the unknown scalar function, with given valuesF (x, y) at every boundary point. Let us assume that the grid h ⊂ of mutually perpendicular, n x horizontal and n y vertical lines, was generated, forming a mesh of n = n x × n y regularly spaced points (Fig. 1). Function values F i,j at those points (i = 1, ..., n x , j = 1, ..., n y ) constitute the set of unknowns, while values at all boundary nodes (located on ∂ h ⊂ ∂ ) are known. Regardless of whether the solution is sought at all grid points or at selected point (or group of points) only, it is required to combine all known (boundary) and unknown (internal) function values into one system of algebraic equations, providing the standard deterministic method is applied. Its generation and solution may be time-consuming, especially in case when h is very dense (number of nodes n → ∞ or mesh size h → 0, e.g., due to global mesh refinement).
Instead, the following stochastic procedure may be adopted, in case only one function value (e.g., at domain center point) is to be determined in relatively fast manner and with reasonable accuracy: 1. Initiate the subsequent random walk starting from the given internal node x i , y j , randomly selecting one of the fourth equally possible direction senses (north, east, south, or west). Move to the closest neighboring node, located in the selected direction ((i, j + 1), (i − 1, j), (i, j − 1) or (i + 1, j), respectively ( Fig. 1)), and repeat the random walk from this new node as long as the first boundary node (x r , y s ) is reached. Evaluate the function value at this boundary node as follows:F r,s =F (x r , y s ) (3) and modify the number of its indications N (e) r,s by one (we initialize all N (e) r,s with zeros). 2. Return to the node of interest x i , y j . Repeat the entire procedure described above until all N random walks are performed (N is the assumed total number of random walks). 3. Evaluate estimation of function value F i,j , according to the Monte Carlo concept (in which reaching the boundary node is treated as a success trial), by means of the simple formula as follows: where n (e) b denotes the number of all boundary nodes with essential boundary conditions (i.e., all boundary nodes in this case).
It has been proved [7] that this stochastic formula (4) corresponds to the finite difference solution, obtained for the considered regular mesh of nodes h , namely The first equation in (5) is the difference equation, generated by means of the nodal collocation technique [46]. It is based upon the standard difference (FD) operator of Laplace type [31], replacing the differential operator from (1) at the central node (i, j ) ( Fig. 1), forming a configuration of five nodes, regularly spaced, called a FD star or stencil [31,45]. It should be noted here that the difference coefficients assigned to four external nodal values (all equal to 1 4 ) are the same as the probabilities of direction sense selection. On the other hand, a function value at the central node in (5) corresponds to the expected value defining the result of a random walk r,sF r,s P r,s i,j . Here, P r,s i,j is the probability of random walk termination at boundary point (r, s), starting from the internal point (i, j ). The number of random walks N has the crucial influence on the accuracy of estimation (4). Generally, the Monte Carlo solution error may be upper-bounded by the non-linear function of N, namely However, this formula is valid for the simplest cases only, in which there is no curse of dimensionality (i.e., mesh size h has no influence on the solution estimation). In more complex problems, discussion on the optimal selection of N (as well as the relation between N and h) may be found in [20,23,28]. Although the standard Monte Carlo method with random walk technique (MC/RW) exhibits numerous advantages, for instance, an existence of an explicit formula for an approximate solution, algorithm simplicity or low computational cost, serious drawbacks may be noted either. First of all, the formula (4) works for regular meshes only, with equal and fixed mesh spacing (neither curvilinear edges nor local mesh refinements may be introduced). Moreover, it may be applied to Laplace equations only (no heterogeneity in right-hand side function), with constant material coefficients. Furthermore, boundary conditions (2) may be of essential type only. This paper, accompanied with the attached Matlab software, aims to extend the standard MC/RW method towards issues listed above. However, several aspects of the approach have to be reformulated and extended, including the following: -Selection of potential random walk directions for the case of arbitrarily irregular cloud of nodes -Determination of probabilities of a next move, depending on the equation type -Significant reconstruction of the final MC/RW formula, taking all information concerning the analyzed problem into account in a stochastic manner

Meshless Monte Carlo random walk technique
The following 2D elliptic equation is considered as follows: with essential (Dirichlet) and natural (Neumann) boundary conditions subjected to boundary parts ∂ e , ∂ n (∂ = ∂ e ∪ ∂ n ). f ∈ C 0 is a scalar righthand side function, a = a (x, y) ∈ C 0 is a scalar positive material function, whereas q ∈ C 1 is a given scalar flux function, applied in the direction represented by the versor n (x, y) = n x (x, y) , n y (x, y) t , normal to the boundary ∂ n . This equation may be a representative mathematical model for a stationary heat flow analysis [19]. Problem domain may have complex geometry, therefore arbitrarily irregular cloud of nodes h = {(x i , y i ) , i = 1, ..., n} ∈ R 2 , without any a-priori imposed structure (like regular mesh, finite element [11], or mapping restrictions [16], has to be considered. Such discretization format is typical for the wide group of computational solution approaches, namely meshless methods (MM, [29,31,33,35,45]).
Moreover, in MM, the unknown function is approximated in terms of nodes only, by means of various techniques, like partition of unity or moving weighted least squares (MWLS, [18,25,31,45,46]). The MWLS approximation has been successfully applied in one of the oldest meshless methods, meshless finite difference method (MFDM, [10,11,16,31,45,46]). The advantages of using MFDM and other MM may be observed in problems with large deformations, concentrated loads, moving boundary, crack development, elastic-plastic boundary, contact of deformable bodies, fluid free surface as well as in h-adaptive approach [12] since nodes may be shifted, added, or removed without any larger influence on nodes topology. In this research, we incorporate selected useful features of both the MWLS technique and the MFDM in order to extend the standard MC/RW technique towards effective stochastic analysis of (7) with (8) and (9).

Selection of random walk directions
A random selection of four mutually perpendicular direction senses, with equal probabilities, which is natural for regular mesh, does not hold in case of irregular clouds of nodes. New direction selection criteria have to be carried out, taking advantage from the nodes' irregular distribution. Let us consider the determination problem of potential walk directions, starting from the central node (i) towards selected nodes (j (i)), (j (i) = 1, 2, ..., m), as it is depicted in Fig. 2a. Both the total number of nodes (m + 1) in such configuration (denoted as the MFD star, or stencil) and their distribution should be assumed in such a manner that the resulting approximation scheme remains non-singular and non-ill-conditioned. Thus, m is usually larger than it is required from the order of differential operator (e.g., six nodes for 2D Laplace equation). The following criteria may be applied here as follows: 1. The simplest criterion is based upon the distance from (i) to (j ) nodes only (Fig. 2b). Though very simple, it does not guarantee the well-conditioning. Therefore, the optimal nodes selection techniques are based rather on the topology information than distance between nodes only. 2. In 2D cross criterion [10,16,31,45] (Fig. 2c), domain is divided into four zones.
Each of four semi-axes is assigned to one of these zones. A specified number of nodes (usually 2), closest to the central node (point) is taken from every zone separately; thus, the number of nodes in the MFD star is always constant. 3. In more complex Voronoi neighbous criterion [31,45] (Fig. 2d), only those nodes are selected to the MFD star which are the Voronoi neighbors (product of domain partitioning into a set of Voronoi polygons; each polygon is assigned to a particular node). Voronoi neighbors are those polygons which have common side (strong neighbors) or common vertex (weak neighbors). It should be stressed that the Voronoi neighbors criterion does not assure the same number of nodes in every star. Moreover, the number of nodes is variable and may not be sufficient to build full FD operator of the specified order. However, such star may be completed by additional nodes to preserve the required approximation order as well as good conditioning (e.g., the last selected ninth node in Fig. 2d, with attached dashed ray, results from the assumed distance criterion). Regardless of the applied selection technique, the final MFD star is ascribed by the topological (m + 1) × 4 matrix S, namely containing local and global nodes' numbers as well as (x, y) coordinates of selected nodes corresponding to potential move directions, respectively, sorted in ascending order.

Determination of direction selection probabilities
Once m potential directions are assigned to the node (i), a procedure continues with the random selection of one of those directions. Respectively to the standard MC/RW method, probabilities of direction selection may be derived from the meshless FD operator, replacing the left-hand side differential operator from (7).
Moreover, the selection probability P (i) j , attached to the j (i)th star node/direction, should be proportional to the distance from node (i) to node (j (i)). Therefore, appropriate approximation scheme has to be determined. We use here the MWLS technique, which is based upon the Taylor series expansion of all star nodal values with respect to the central star node (i). In such a manner, resulting approximation coefficients have simple interpretation as local derivatives, up to the assumed pth order. Fundamentals and general remarks concerning the MWLS technique may be found in [18,25,31], whereas details of its application for the second-order differential equations are given in [45,46,48]. Comparison of MWLS (with singular weight functions) and MLS (with non-singular weight functions) techniques with other meshless approximation methods may be found in [29,31,33,35]. Here, we focus on the final approximation formula (for p = 2) as follows: which may be obtained after the minimization of the appropriate weighted error func- difference coefficients matrix, D F is the 6 × 1 approximation coefficients vector, and q is the (m + 1) × 1 vector of degrees of freedom: in which the singular weights are applied as follows: with being a relatively small though non-zero (up to the machine precision) number (e.g., ∼ 10 −15 for double precision real numbers). Therefore, interpolation is forced at the central star node, in spite of using larger number of nodes than it is required from the assumed second-order differential operator order in (7) and the first-order differential operator in (9). Once M is determined, values of all derivatives included in D F as well as value of any arbitrary differential operators (up to the sec-ond order) may be composed. By means of collocation technique at subsequent nodes with unknown function values, appropriate difference equations may be generated as follows: for nodes (i) ∈ , (k) ∈ ∂ e , (l) ∈ ∂ n , respectively. Here, m i is the number of directions for the internal node, whereas m b is the number of directions for the boundary node. By rearranging terms in (15) and (17), one obtains the selection probabilities P (i) j and P (l) j of each j th potential direction of the subsequent random walk (Fig. 3), initialized from the ith or lth node, respectively  In other words, the probability of selection of any other direction, not included within the FD star assigned to the ith or lth node, equals zero.

Modified Monte Carlo random walk formula
The principles of a random walk technique remain unmodified. We start from the specified arbitrary node located inside the domain (marked as a green circle in Fig. 3) or on its boundary with natural boundary conditions, where the function F is unknown (marked as a blue circle in Fig. 3). We proceed until the first node located on the boundary with essential boundary conditions is reached (marked as red circle in Fig. 3). Eventually, the final closed-form formula for the stochastic estimation of unknown function value F j , by means of the meshless Monte Carlo method with random walk technique is as follows: It takes into account series of all N random walks as well as all additional righthand side components from (15) and (17). Here,F k =F (x k , y k ),q l =q (x l , y l ), f i = f (x i , y i ), a i = a (x i , y i ), n x(l) = n x (x l , y l ) and n y(l) = n y (x l , y l ). Various difference coefficients M correspond to both internal (M (i) ) as well as boundary nodes (M (l) ). Moreover, n i are the numbers (counters) of kth, lth, and ith nodes' random walk indications, respectively. It should be emphasized here that the simplified error estimation formula (6) no longer holds, since the final MC solution (20) depends on both the number of random walks N and the physical discretization parameter(s) h, included in difference coefficients M. Therefore, new combined stochastic-deterministic-type estimation ε has to be defined, in the following general form as follows: Determination of ε as well as the optimal relation between N and h, will be considered in the following author's papers.

A Matlab software-general information
All original author's concepts and ideas, briefly presented in previous section, were successfully implemented in Matlab. The software is available from Netlib (http://www.netlib.org/numeralgo/) as the na52 package. It has been prepared and tested in Matlab R2014b; however, it does not contain any unusual functions and syntax that might cause any problems while using other Matlab versions. A user is supplied with one archive file (meshlessMCRW 2D.zip) which contains several script and function files. It has to be downloaded and extracted to the destination folder.  terion is not implemented yet, due to the requirement of the constrained Voronoi partitioning-work is in progress), m smooth -number of nodes in the FD stars for smoothing purposes (for mode = 3), g smooth -value of a smoothing parameter g (for mode = 3) rw graph -graphical tracking of random walk paths (0, off; 1, on) (for mode = 1) sol graphs -graphs of the final solutions (0, off; 1, on) (for mode = 3) -MC sol -meshless Monte Carlo with random walk analysis (0, off; 1, on) mesh plot -domain and mesh plotting (0, off; 1, on) font weight -font weight of figures' labels ("normal," "bold") marker size -marker size of plots' points font size -font size of figures' labels graph view -graph view style (2, 2D with colorbars displayed; 3, 3D) Three main modes are described below in more detailed manner.

Determination of the approximated solution at the selected point, for the fixed number of random walks (mode=1)
The simplest program mode allows for the determination of one unknown function value approximation at the selected point of the domain (or its boundary, assuming the natural conditions are subjected). Coordinates of this point (variables x0,y0) may be user-defined or default values may be proceeded. In case a regular mesh is applied (dist amp = 0), this given point may not belong to the mesh itself. In that case, solution is obtained at the mesh node, closest to that point. However, in case of irregular cloud (dist amp > 0), the closest internal node is moved to this point (total number of nodes remains unchanged) or additional boundary node is added (the total number of nodes is increased by 1). The simplified flow chart of this algorithm mode is shown in Fig. 5. Graphical interface produces a mesh/cloud of nodes scheme, with the selected point indicated (Matlab's Fig. 1). Moreover, text results are displayed in Matlab's Command Window. Additionally, if a user sets variable rw graph equal to one, all random walk paths are graphically displayed in Matlab's Fig. 1, in a form of a simple animation.

Determination of the approximated solution at the selected point, for the series of random walks' numbers (mode=2)
The second mode may be treated as a generalization of the first one. Instead of one fixed number of random walks (N), a user defines an initial number of random walks (variable N init), increment (variable dN) as well as number of steps (variable adap N). Compared to the flow chart for the first mode, one additional loop appears, which selects the subsequent numbers of random walks out of the series N init , N init + dN, N init + 2dN, ..., N init + adap N · dN. However, this slightly modified flow chart has been omitted here. As the final result, the convergence graph is generated (Matlab's Fig. 2 Moreover, both reference solutions,F and F FDM , are fixed, since the deterministic methods do not rely on the number of random walks. Additionally, the theoretical convergence curve is plotted for the sake of comparison, according to (6).

Determination of the approximated solution at all nodes, for the fixed number of random walks (mode=3)
Though the proposed approach is tailored to estimate the solution of elliptic equations at selected points, it would be interesting to apply the final formula (20) to all nodes with unknown function value, located inside the domain and on its boundary part with natural boundary conditions. In that case, one has to deal with different levels of the rough Monte Carlo solution accuracy, which may change from node to node, due to stochastic error, even though the same number of random walks is used each time. Therefore, the appropriate smoothing may be required to keep the reasonable balance in the solution error. It may be performed by means of the MWLS approximation (11) as well. However, this time the larger number of nodes in FD stars as well as nonsingular weight functions [31] should be used, with built-in smoothing parameter g as follows: Parameter g is a positive scalar number, selected arbitrarily, by observing the smoothing effect or adaptively, in order to minimize the local curvature of the approximating function. Setting g = 0 in (22) leads to singular weights (14) (no smoothing applied). It should be stressed here, that this additional a-posteriori MWLS approximation may be used for recovery of nodal values of the solution derivatives F x and F y as well.
The complete flow chart of this algorithm mode is presented in Fig. 6. In addition to the standard text results displayed in Matlab's Command Window, a user is supplied with graphs of three solutions, namely true solution, FD solution, and Monte Carlo solution as well as the differences of the above (Matlab's Fig. 3). Moreover, comparison of solutions' first derivatives is presented in Matlab's Fig. 4. All graphs are plotted using the additional background plotting mesh of Delaunay triangles [16,31].

Numerical experiments and software testing
A variety of tests were executed in order to examine both the proposed approach and the attached Matlab software. Selected results are presented and briefly discussed. Illustrations of mode = 1 with rw graph set to 1 are shown in Fig. 7, for various domain types with mixed boundary conditions as well as for various analytical solutions. N = 100 random walks were applied, on very coarse irregular clouds, with small numbers of nodes. All other input parameters have their default values. The following points with unknown solution were examined, namely: (1.3333, 1.3333) (Fig. 7a), (0.5, 0.86603), (Fig. 7b), (1, 1.1547) (Fig. 7c), and (0, 2.5) (Fig. 7d). The final distribution of all non-zero indication numbers (N i ) as well as non-zero nodal values (F k ,q l , and f i ), assigned to appropriate nodes (thus giving the non-zero components to the final MC formula), is shown in subsequent graphs. It may be observed how many times the particular node has been reached during the random walk process. Moreover, the exemplary random walk path is plotted (Fig. 7a). The final rough solution estimation as well as its true and numerical errors are presented in the graphs' labels.
The following tests concern the solution convergence with respect to the number of random walks (mode = 2). The series of 50 increasing numbers of walks, starting from N = 50 and with increment dN = 200, were performed. The same domain and solution types were applied as for mode = 1. Results (solution errors) are presented in Fig. 8. The MC solution error is perfectly upper-bounded by its a-priori estimation, corresponding to (6), for the simplest case only (Fig. 8a), whereas it is no longer valid for other cases. Nevertheless, in each case, very good agreement between the FD and MC/RW solutions may be observed, even though there may be enormous true errors (Fig. 8b), caused by a very coarse discretization, when compared to the exact solution.
Eventually, the full analysis of the elliptic (7) was carried out. Non-constant material function a(x, y) = 1 1 + cos 2 (x + y) π 3 was set. The meshless Monte Carlo random walk stochastic formula (20) was applied to every unknown function value at all nodes located inside the domain and on its boundary part with natural boundary conditions (if applicable). Additional a-posteriori solution smoothing was performed, according to (22). Graphs of the final solution, the solution errors as well as its first derivatives are presented in Figs. 9 and 10 (circular domain, with n = 128 nodes irregularly scattered, with N = 500 and solution #2). For this solution type, FD analysis produces the exact solution, within the second polynomial order assumed, up to the machine precision. However, the Monte Carlo solution approach yields a rough solution estimation only, as it is highly influenced by a stochastic error. Solution and  Clouds of nodes generated for the first three domain types (#1, #2, and #3) may be refined in an arbitrary manner by means of the mesh generator, being a part of the presented software (the function file mesh generation.m). However, the last three, more complex domain types have fixed discretization, loaded from the external disc files, as the generation of discretization for more complex shapes would require more sophisticated software. Furthermore, it should be mentioned that the presented Matlab software may be extended by a user in many potential directions, without any difficulties. The modification of the material function may be done in fun a.m function file. New exact solution requires adding new conditional branch (formulae for solution and its derivatives in explicit manner) in analyt sol.m. In case the analytical solution is not known, appropriate boundary conditions as well as right-hand side functions may be formulated in analyt sol.m, fun f.m and fun q.m. A brand new domain shape may be defined in mesh generation.m or another text file with mesh data should be prepared (number of nodes, nodes coordinates and boundary codes (0, internal node; 1, boundary node with essential boundary conditions; 2, boundary node with essential boundary conditions) as well as the number of Delaunay triangles and a list of all triangles' vertices-for potential graphical purposes only).

Final remarks
This research focuses on the modified, meshless Monte Carlo method with random walk (MC/RW) solution approach. It is based on the appropriate discretization and approximation techniques, typical for meshless modeling of boundary value problems. When compared to the standard MC/RW method, it allows for the analysis of wider class of elliptic problems, with non-homogeneous right-hand side functions, variable material coefficients as well as mixed (essential and natural) boundary conditions. Moreover, arbitrarily irregular meshes and unstructured clouds of nodes may be applied. Therefore, analysis of more complex domain shapes is possible.
The proposed approach is based upon the random generation of indications net (random walk technique) starting from the fixed point (node), at which solution estimation is required. Selection probabilities of subsequent walk directions depend on the meshless nodes configuration (equivalent to the FD star) as well as the difference coefficients corresponding to differential operators, obtained by means of the moving weighted least squares approximation. Finally, the Monte Carlo concept is incorporated, according to which the random walk terminates at boundary node located on boundary part with essential boundary conditions (with known solution). Therefore, fast and reasonably accurate solution estimation at node is possible. Moreover, the generation of a large system of equations, combining all unknown function values may be omitted. Full analysis (convenient for potential parallelization) of all internal and boundary nodes with unknown solution has been discussed. Additional postprocessing includes the recovering of solution's derivatives as well as the alternative smoothing of nodal results, which leads to significant reduction of a stochastic error.
All those original author's concepts have been implemented in the attached Matlab software. It may work in three various modes, namely 1. determination of an approximated solution at specified node, 2. solution convergence analysis on the set of random walks' numbers, 3. determination of an approximated solution at all nodes. Graphical interface allows for visualization of the random walk structure and indications' numbers (mode=1), convergence graph (mode=2) as well as graphs of the true, FD and Monte Carlo solutions, corresponding errors and solution's derivatives (mode=3). Six types of domain geometries as well as seven types of the analytical solutions are included. Algorithm's flow charts and the selected results of software usage are presented.
Further research includes, for instance, application of the meshless Monte Carlo random walk solution approach to analysis of the 3D elliptic equations as well as nonlinear and non-stationary thermo-mechanical problems. In case the reference problem has to be solved multiple times (e.g., within an incremental-iterative procedure or an implicit time integration scheme), series of random walks forming indications net, are performed only once. Afterwards, the meshless Monte Carlo formula (20) is applied to each solution increment, allowing for the significant reduction of the computational effort.