Stencil and kernel optimisation for mesh-free very high-order generalised finite difference method

Mathematics and Computers in Simulation 218 (2024) 49–78 Available online 22 November 2023 0378-4754/© 2023 The Author(s). Published by Elsevier B.V. on behalf of International Association for Mathematics and Computers in Simulation (IMACS). This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Mathematics and Computers in Simulation journal homepage: www.elsevier.com/locate/matcom Original articles Stencil and kernel optimisation for mesh-free very high-order generalised finite difference method S. Clain a,∗, J. Figueiredob,c aCenter of mathematics of the FCTUC, Largo D. Dinis, 3000-143 Coimbra, Portugal bCentre of Physics UM - UP, Campus de Gualtar, 4710-057 Braga, Portugal cDepartment of Mathematics, Campus de Azurém, University of Minho, 4080-058 Guimarães, Portugal ARTICLE INFO Keywords: General finite difference Mesh-free Stencil optimisation Very high-order ABSTRACT Generalised Finite Difference Methods and similar mesh-free methods (Pointset method, Multipoint method) are based on three main ingredients: a stencil around the reference node, a polynomial reconstruction and a weighted functional to provide the relations between the derivatives at the reference node and the nodes of the stencil. Very few studies were dedicated to the optimal choice of the stencil together with the other parameters that could reduce the global conditioning of the system and bring more stability and better accuracy. We propose a detailed construction of the very high-order polynomial representation and define a functional that assesses the quality of the reconstruction. We propose and implement several techniques of optimisation and demonstrate the advantages in terms of accuracy and stability. 1. Introduction The Generalised Finite Difference Method (GFDM) dates back to the 70s’ with the articles of Jensen [14], Perrone and Kao [24], Liszka and Orkisz [18] on irregular grids. Such a technique has attracted a lot of researchers and schemes’ designers for several reasons: (1) the simplicity to handle boundary conditions while providing very high-order approximations; (2) the ability to avoid the complex meshing machinery to achieve refinement; (3) the ability to handle the Lagrangian/ALE formulation by moving the nodes. There exists a recent literature on the GFDM method for problem using an Eulerian formulation to take advantage of the ability to handle complex geometries [29,32,33]. Under the umbrella of GFDM, a plethora of mesh-free methods has been proposed and implemented. •The Generalised Finite Difference Method proposed by Liszka & Orkisz [18] and Benito, Gavete, Ureña et al. [1,2,6,30,31], is based on the Taylor expansion to link the partial derivatives with the function at the neighbouring nodes. The least squares minimisation then provides an explicit relation for each derivative with the neighbour values. Such relations are then plug into the physical equations (inner equations or boundary equations) and one obtains a (non)-linear system that only involves the unknown values at the nodes. •The Finite Pointset Method (FPM) proposed by Kuhnert et al. [15,16,28], also uses the Taylor expansion terms, i.e. the partial derivatives, at a given node we connect to the function values at the neighbouring nodes. The principal difference with the GFDM is the inclusion of the physical equations (inner equations and boundary equations) in the weighted least squares procedure. Consequently, the computed partial derivatives do not exactly satisfy the physical equations. Such a technique has also received important developments by Reséndiz-Flores et al. [25–27]. ∗Corresponding author. E-mail address: [email protected] (S. Clain). https://doi.org/10.1016/j.matcom.2023.11.009 Received 28 May 2023; Received in revised form 15 September 2023; Accepted 7 November 2023 Mathematics and Computers in Simulation 218 (2024) 49–78 50 S. Clain and J. Figueiredo •The Finite Point Method of Oñate, Idelsohn et al. [3,4,20–23], derives from a general weighted residual method using a local polynomial representation of the approximation over subsets 𝛺𝑘⊂ 𝛺. Polynomial coefficients are explicitly determined regarding the unknowns located on 𝛺𝑘. Notice that the local representation is the same on the whole domain 𝛺𝑘, contrarily to the GFDM where each representation is related to a node. For that reason, the authors use a Moving Weighted Least Squares method to create a dependence between the polynomial coefficients and the nodes of 𝛺𝑘. The reconstruction is plugged into the weighted residual form of the physical equations that provides the (non)-linear system one has to solve regarding the unknowns at the nodes. •The Multipoint method developed by Jaworska and Orkisz et al. [10–13] is quite different from the previous methodologies, since the authors are dealing with implicit relations between the values and the derivatives. Given a reference node 𝑖and a stencil V𝑖, the derivative at node 𝑖is given by a linear combination of the neighbouring node values together with the derivatives (at the same order) over nodes of V𝑖. Such an approach is close to the compact scheme philosophy, since it combines both the values and the derivatives. As in the Pointset method, the numerical approximations are plugged into the physical equation (inner/core and boundary equation). All the methods are based on three main ingredients: a stencil around the reference node 𝑖, a polynomial reconstruction and a weighted functional to provide the relations between the derivatives at the node 𝑖and the nodes of the stencil. Then arises the question of optimisation of the three ingredients: the choice of the nodes to elaborate the stencil, the scaling parameter ℎ𝑀of the polynomial representation and the scaling parameter ℎ𝑊of the kernel that produces the weights. At the end of the day, one has to solve a linear system of type 𝑀𝑊 𝑀t𝑐=𝑀𝑊 𝛷 where 𝑐represents the coefficients’ vector of the local polynomial reconstruction while 𝛷gathers the values at the neighbouring nodes. The matrix 𝑀=𝑀(ℎ𝑀)gathers the polynomial expressions, while 𝑊=𝑊(ℎ𝑊)is the weights’ diagonal matrix. Condition number of 𝑀𝑊 𝑀tis then a key issue to control the accuracy and stability of the method. There exists several recipes in the literature to improve the computation of the polynomial. Jensen [14] selects the closest nodes from the reference node. Perrone and Kao [24] propose the eight segments criterion, consisting on the selection of a node for each octant in a system of Cartesian axes around the reference node. Liszka and Orkisz [18] introduce the four quadrants technique, which consists in picking the two nearest nodes per quadrant. We also mention the reviews on high-order finite difference method on irregular domains [8] or approximation of non-smooth solutions [9]. We propose an optimisation procedure to achieve the best stencil and coefficients ℎ𝑀and ℎ𝑊that minimise the condition number of matrix 𝑀𝑊 𝑀t. On the one hand, the choice of the stencil is related to a discrete minimisation problem where the nodes are the parameters to optimise. On the other hand, the continuous minimisation regarding parameters ℎ𝑀and ℎ𝑊is achieved by different derivative-free methods. Numerical tests with and without optimisation for different geometrical configurations are carried out to assess the condition number reduction and the accuracy. We restrict the study to the steady state case since we focus on the operator reconstructions in space together with the boundary conditions prescribed on non-polygonal physical domain. Moving boundary domain with Lagrangian formulation is out of the scope of the present paper. The rest of the paper is organised as the following. We first recall some basic facts in Section 2on the polynomial reconstruction and define the material we shall use all along the paper in Section 3. Section 4concerns the optimisation procedures, both for the stencil (discrete minimisation) and the matrix parameters (continuous minimisation). We present an extensive series of benchmarks to assess the advantage of the minimisation procedures and prove the interest in strongly reducing the condition number by some magnitude. We present some examples of simulations with the convection diffusion reaction equation by using the derivatives computed with the optimal reconstruction. 2. Notations The R2space is described by the generic point 𝑥= (𝑥1, 𝑥2)and 𝜙(𝑥) = 𝜙(𝑥1, 𝑥2)stands for a real-valued smooth function. We introduce the multi-index notation 𝛽= (𝛽1, 𝛽2) ∈ N2with |𝛽|=𝛽1+𝛽2and set, 𝜕𝛽𝜙(𝑥) = 𝜕𝛽1 𝑥1𝜕𝛽2 𝑥2𝜙(𝑥1, 𝑥2). Let Abe a finite subset of N2. We define the polynomial of degrees in A, centred at a point 𝑥, by 𝜋(𝑥;ℎ𝑀, 𝑥, 𝑐, A) = ∑ 𝛼∈A 𝑐𝛼(𝑥−𝑥 ℎ𝑀)𝛼 , where 𝑐𝛼are real number coefficients and ℎ𝑀is a scaling parameter. We extend the factorial for multi-index values setting 𝛼!=(𝛼1)! (𝛼2)! and 𝛼!=0if 𝛼1<0or 𝛼2<0. Thanks to this definition, we deduce the relation 𝜕𝛽(𝑥−𝑥 ℎ𝑀)𝛼 =(𝛼−𝛽)! ℎ|𝛽| 𝑀(𝑥−𝑥 ℎ𝑀)𝛼−𝛽 In particular, one has the property 𝜕𝛽(𝑥−𝑥 ℎ𝑀)𝛼 𝑥=𝑥 =𝛽! ℎ|𝛽| 𝑀 𝛿𝛼,𝛽 Mathematics and Computers in Simulation 218 (2024) 49–78 51 S. Clain and J. Figueiredo where 𝛿𝛼,𝛽 = 1 if 𝛼=𝛽, and null otherwise. Let 𝛺be an open bounded domain of R2with boundary 𝜕𝛺. We denote by G= (𝑥𝑖)𝐼 𝑖=1, with 𝑥𝑖= (𝑥𝑖 1, 𝑥𝑖 2), a cloud of 𝐼points of 𝛺. Gcorresponds to the interior nodes of 𝛺while 𝜕Gstands for the nodes on the boundary. For any node 𝑖,V𝑖⊂Gis a stencil of points, usually chosen in the neighbourhood of the node 𝑖. Note that we impose 𝑖∈V𝑖. We further denote by 𝛷= (𝜙𝑖)𝑖∈Gthe vector of the unknowns such that 𝜙𝑖≈𝜙(𝑥𝑖), while we introduce the stencil-vector 𝛷V𝑖= (𝜙𝑗)𝑗∈V𝑖associated to the stencil V𝑖. Similarly, for any 𝛽∈A, we denote 𝛷𝑖,𝛽 ≈𝜕𝛽𝜙(𝑥𝑖)the approximations of the derivatives at node 𝑖. 3. Derivatives’ discretisation The goal is to provide explicit relations for the 𝛽derivative of 𝜙, at node 𝑖, regarding the data over the stencil V𝑖, namely, we seek coefficients 𝑎𝑖,𝛽 𝑗∈Rsuch that 𝜕𝛽𝜙(𝑥𝑖) ≈ ∑ 𝑗∈V𝑖 𝑎𝑖,𝛽 𝑗𝜙𝑗. The structural relations provide the discrete relations where 𝑎𝑖,𝛽 𝑗only depends on the structure of the cloud. To determine the coefficients, we present two distinct methodologies based on the least square method. 3.1. The polynomial approximation Let 𝑖be a node and V𝑖its stencil. We seek a polynomial 𝜋(𝑥;ℎ𝑀, 𝑥𝑖, 𝑐, A)that minimises the loss function 𝐽𝑖(𝑐) ∶= 1 2∑ 𝑗∈V𝑖 𝜔2 𝑖𝑗 (𝜋(𝑥𝑗;ℎ𝑀, 𝑥𝑖, 𝑐, A) − 𝜙𝑗)2 regarding the vector 𝑐∈R|V𝑖|, where 𝜔𝑖𝑗 ≥0stands for the weights between nodes 𝑖and 𝑗. Minimisation provides a vector 𝑐𝑖 solution of the matrix problem [𝑀𝑖𝑊𝑖(𝑀𝑖)t]𝑐𝑖=𝑀𝑖𝑊𝑖𝛷V𝑖(1) with 𝑊𝑖=diag(𝜔2 𝑖𝑗 )the square diagonal |V𝑖|×|V𝑖|matrix of the weights and 𝑀𝑖[𝛼, 𝑗] = 𝑚𝑖,𝛼 𝑗the |A|×|V𝑖|matrix given by 𝑚𝑖,𝛼 𝑗∶= (𝑥𝑗−𝑥𝑖 ℎ𝑀)𝛼 , 𝛼 ∈A, 𝑗 ∈V𝑖.(2) Remark 1. In practice and for implementation purposes, we use a one-to-one index 𝛼∈A→𝓁and a one-to-one local index mapping 𝑗∈V𝑖→𝜄to store the information in 𝑀[𝓁, 𝜄]. We do not refer to the local indexation for simplicity and 𝑀[𝛼, 𝑗]has to be interpreted as 𝑀[𝓁, 𝜄]. The stencil V𝑖is said to be resolvent if 𝑀𝑖has the maximum rank equal to |A|assuming that |V𝑖|≥|A|. Indeed, we have enough information from the neighbouring nodes to completely determine the polynomial coefficients and the square matrix 𝑀𝑖𝑊𝑖(𝑀𝑖)tis non-singular. Noting that, for 𝛽∈A, we have, 𝜕𝛽𝜋(𝑥𝑖;ℎ𝑀, 𝑥𝑖, 𝑐𝑖,A) = 𝑐𝑖 𝛽 𝛽! ℎ|𝛽| 𝑀 . We then obtain the relation between the partial derivatives and the coefficients of the vector 𝑐𝑖. To deliver an explicit relation regarding the value of 𝜙𝑗with 𝑗∈V𝑖, we introduce the diagonal matrix 𝐷=(𝛽! ℎ|𝛽| 𝑀 𝛿𝛽,𝛼)𝛽,𝛼∈A =diag ⎡⎢⎢⎣(𝛽! ℎ|𝛽| 𝑀)𝛽∈A⎤⎥⎥⎦ and define, using relation (1), the |A|×|V𝑖|matrix 𝐴𝑖=𝐷[𝑀𝑖𝑊𝑖(𝑀𝑖)t]−1𝑀𝑖𝑊𝑖. Letting 𝑎𝑖,𝛽 =𝑀𝑖[𝛽, ∶] ∈ R|V𝑖|and vector 𝛷V𝑖the data in the neighbourhood of node 𝑖, we obtain the discretisation of the 𝛽-derivative at node 𝑖given by 𝜙𝑖,𝛽 ∶= 𝑎𝑖,𝛽 𝛷V𝑖=∑ 𝑗∈V𝑖 𝑎𝑖,𝛽 𝑗𝜙𝑗≈𝜕𝛽𝜙(𝑥𝑖),(3) which is exact for all polynomials of degrees in A. Mathematics and Computers in Simulation 218 (2024) 49–78 52 S. Clain and J. Figueiredo 3.2. Lagrangian multipliers formulation We present an alternative approach to produce the relation (3) but with possible extensions that we shall detail in the next section. The idea consists in introducing the relation between the derivatives and the values over the stencil as constraints. To this end, for 𝑖∈Gand 𝛽∈A, we define the functional 𝐸𝑖,𝛽 (𝜙;𝑎) ∶= 𝜕𝛽𝜙(𝑥𝑖) − ∑ 𝑗∈V𝑖 𝑎𝑖,𝛽 𝑗𝜙(𝑥𝑗).(4) To determine the coefficients 𝑎𝑖,𝛽 𝑗, we first prescribe the constraints 𝐸𝑖,𝛽 (𝜙;𝑎)=0for all 𝜙(𝑥) = (𝑥−𝑥𝑖 ℎ𝑀)𝛼,𝛼∈A, that is [𝜕𝛽(𝑥−𝑥𝑖 ℎ𝑀)𝛼]𝑥=𝑥𝑖=∑ 𝑗∈V𝑖 𝑎𝑖,𝛽 𝑗(𝑥𝑗−𝑥𝑖 ℎ𝑀)𝛼 ,∀𝛼∈A. The finite set Acharacterises the space of constraints that controls the accuracy of the approximations. Gathering all the constraints for 𝛼∈A, we get the equality in R|A|, 𝑎𝑖,𝛽 𝑀t=𝑀(𝑎𝑖,𝛽 )t=𝑑𝛽 where 𝑎𝑖,𝛽 is a row vector of R|V𝑖|, matrix 𝑀𝑖is the |A|×|V𝑖|matrix given by (2) and 𝑑𝛽∈R|A|the canonical vector, with zero everywhere except 𝛽!∕ℎ|𝛽| 𝑀for entry 𝛽. Since we have the same relation for every 𝛽∈A, we rewrite the constraints in matrix form 𝐴𝑖(𝑀𝑖)t=𝑀𝑖(𝐴𝑖)t=𝐷=𝐷t, where we collect the row vector 𝑎𝑖,𝛽 in the |A|×|V𝑖|matrix 𝐴𝑖. Since |V𝑖|≥|A|, we do not have uniqueness for the coefficients 𝑎𝑖,𝛽 𝑗(except for the particular case |V𝑖|=|A|), and the number of constraints is not enough to determine the vector 𝑎𝑖,𝛽 . Consequently, we introduce the weighted energy functional 𝐺𝑖(𝑎) ∶= 1 2𝑎(𝑊𝑖)−1𝑎t, 𝑎 ∈R|V𝑖|, and, for each 𝛽∈A, we consider the (constraint) minimisation problem 𝑎𝑖,𝛽 ∶= arg min 𝑎∈R|V𝑖|𝐺𝑖(𝑎),with 𝑀𝑖𝑎t=𝑑𝛽, where 𝑎is a row vector of the coefficients. An important note is: on the contrary to the polynomial approach, we can determine vectors 𝑎𝑖,𝛽 ,𝛽∈A, independently, that is, we do not need to compute the whole matrix 𝐴𝑖but only the derivatives of interest. Minimising 𝐺𝑖under the constraints 𝐸𝑖,𝛽 ((𝑥−𝑥𝑖 ℎ𝑀)𝛼 ;𝑎)= 0 for all 𝛼∈Ausing the Lagrangian multipliers provides the system 𝑎𝑖,𝛽 (𝑊𝑖)−1 +∑ 𝛼∈A 𝜆𝑖,𝛽 𝛼(𝑚𝑖,𝛼)t= 0, with 𝑚𝑖,𝛼 =𝑀[𝛼, ∶]. Setting 𝜆𝑖,𝛽 = (𝜆𝑖,𝛽 𝛼)𝛼∈A, the Lagrangian reads (𝑊𝑖)−1(𝑎𝑖,𝛼)t+𝜆𝑖,𝛽 𝑀𝑖= 0,in R|A|. Gathering the line vectors 𝜆𝑖,𝛽 into the matrix 𝛬𝑖, we obtain the relation 𝐴+𝛬𝑖𝑀𝑖𝑊𝑖= 0 we have to solve together with the constraint 𝐴(𝑀𝑖)t=𝐷. After algebraic operations, we deduce 𝛬𝑖= −𝐷(𝑀𝑖𝑊𝑖(𝑀𝑖)t)−1 and we find, once again, 𝐴𝑖=𝐷(𝑀𝑖𝑊𝑖(𝑀𝑖)t)−1𝑀𝑖𝑊𝑖. Remark 2. As we mentioned above, the method enables to compute the vector 𝑎𝑖,𝛼 corresponding to the derivatives’ discretisation of interest, and, in that way, reduce the computational effort. For example, assume that |A|is large (for example 15 for |𝛼|≤5) and that we only need the discretisation of 𝜕2 𝑥1=𝜕(2,0) and 𝜕2 𝑥2=𝜕(0,2). Therefore, we just determine the two vectors corresponding to these two derivatives. 4. Optimisation The existence of a discrete representation of the derivatives requires that the matrix 𝐾𝑖=𝑀𝑖𝑊𝑖(𝑀𝑖)tis non-singular. We claim that accuracy is strongly related to its condition number, 𝜒𝑖=𝜒(𝐾𝑖), and reducing the condition number shall improve the quality of the approximation of the derivatives (see Section 5.2.1 for the justification). By analysing the construction of 𝑀𝑖and 𝑊𝑖, we identify three ingredients that could help to lower 𝜒𝑖: (1) the stencil choice, (2) the polynomial normalisation parameter, and (3) the weights. We shall consider such ingredients as independent degrees of freedom to compute 𝜒𝑖. Mathematics and Computers in Simulation 218 (2024) 49–78 53 S. Clain and J. Figueiredo To provide the weights, we consider a decreasing function 𝑠→𝜔(𝑠)defined on [0,+∞[ with 𝜔(0) = 1 and 𝜔(+∞) = 0. We then set 𝜔𝑖 𝑗=𝜔(|𝑥𝑗−𝑥𝑖| ℎ𝑊), where the scaling parameter ℎ𝑊controls the kernel range. Consequently, the condition number 𝜒𝑖=𝜒(𝐾𝑖) = 𝜒(V𝑖, ℎ𝑀, ℎ𝑊)depends on V𝑖,ℎ𝑀,ℎ𝑊, which control, respectively, the stencil, matrix 𝑀and the weights 𝑊. We intend to minimise 𝜒regarding these three parameters. The optimisation procedure involves two continuous variables and a discrete one, hence requiring a different approach in function of the nature of the parameter addressed. 4.1. Stencil Size optimisation (SSO) In the first stage of optimisation, we shall deal with parameter V𝑖. The discrete optimisation consists in starting with a large initial stencil of G, usually twice the number of strictly necessary nodes, around the node 𝑖and performing different selections of nodes to reduce both the stencil size (computational effort reduction) and the condition number. Dropping index 𝑖for clarity, we consider the operator V→𝜒(V) = 𝜒(V, ℎ𝑀, ℎ𝑊), where the matrix 𝑀=𝑀(V)has been constructed using the nodes of the stencil V, parameters ℎ𝑀and ℎ𝑊being frozen. The Stencil Size Optimisation (SSO) is based on a brute force-like technique to produce a sequence V(m)such that |V(m+ 1)|=|V(m)|− 1 while avoiding increasing the condition number 𝜒(V(m)). Note that we require that the initial stencil V(0) guarantees that the matrix 𝑀(V)enjoys the maximal rank property. The iterative procedure reads: 1. Set the initial stencil V(0) with maximum rank property; 2. Given the stencil V(m)at stage m, we seek the stencil V(m+ 1) ⊂V(m)that minimises 𝜒by subtracting one node to V(m). In other words, we define the candidate stencil as W= arg min 𝑗∈V(m)𝜒(V(m)⧵{𝑗}). 3. If 𝜒(W)≤(1 + 𝜖)𝜒(V(m)) then we set V(m+ 1) = Wand go back to step (2). 4. Otherwise, the final stencil is given by Vsso =V(m). A strict reduction of the condition number would be quite restrictive, and one can consider a trade-off between stencil size reduction and conditioning. Therefore, we introduce the tolerance factor 𝜖that quantifies the flexibility to let the condition number grow. Notice that 𝜖 > 0is a tolerance, whereas 𝜖 < 0is a penalisation. In this way, we define an operator such that, given an initial stencil V, it provides the optimal stencil Vsso ⊂Vthat minimises 𝜒by using the nodes of Vuniquely. To assess the quality of the optimal stencil, we define the Condition Number Ratio (CNR) and the Stencil Size Ratio (SSR) associated to the stencil Vas CNR(V) = 𝜒(V) 𝜒(Vsso),SSR(V) = |Vsso| |V|≥1. These coefficients assess the gains of the optimisation procedure in terms of condition number and stencil size. 4.2. Scaling Parameters Optimisation (SPO) In addition to the stencil optimisation, parameters ℎ𝑀and ℎ𝑊that control the scaling of matrices 𝑀and 𝑊, respectively, will be optimised to reduce the condition number. We propose to use a greedy algorithm where we successively optimise the stencil and the parameters. To this end, given an initial stencil V(0) = Vand initial parameters ℎ𝑀(0) = ℎ𝑀,ℎ𝑊(0) = ℎ𝑊, we build a sequence (V(m), ℎ𝑀(m), ℎ𝑊(m))such that |V(m+ 1)|=|V(m)|− 1 and 𝜒(V(m+ 1), ℎ𝑀(m+ 1), ℎ𝑊(m+ 1))≤(1 + 𝜖)𝜒(V(m), ℎ𝑀(m), ℎ𝑊(m)), with 𝜖 > 0for a relaxed constraint or 𝜖 < 0for a tighter constraint. The greedy procedure consists in two stages: 1. Given V(m),ℎ𝑀(m),ℎ𝑊(m), we first optimise the discrete problem by determining the new stencil W(m+ 1) = arg min 𝑗∈V(m)𝜒(V(m)⧵{𝑗}, ℎ𝑀(m), ℎ𝑊(m)). 2. With the new stencil frozen, we optimise the continuous parameters under positivity constraints ( ℎ𝑀, ℎ𝑊) = arg min ℎ𝑀>0, ℎ𝑊>0 𝜒(W(m+ 1), ℎ𝑀, ℎ𝑊). Mathematics and Computers in Simulation 218 (2024) 49–78 54 S. Clain and J. Figueiredo Fig. 1. The cloud-like optimisation procedure. Illustration of the (ℎ𝑟 𝑀, ℎ𝑟 𝑊)convergence towards ( ℎ𝑀, ℎ𝑊)for a 7×7points uniform cloud. The panels (from left to right) present the results obtained for three consecutive iterations. The green filled circle corresponds to (ℎ𝑟,1 𝑀, ℎ𝑟,1 𝑊), the orange squares to (ℎ𝑟,𝓁 𝑀, ℎ𝑟,𝓁 𝑊), 𝓁= 2,…,L− 1, and the red triangle to (ℎ𝑟,L 𝑀, ℎ𝑟,L 𝑊). We then update V(m+ 1) = W(m+ 1),ℎ𝑀(m+ 1) =  ℎ𝑀,ℎ𝑊(m+ 1) =  ℎ𝑊if 𝜒(W(m+ 1), ℎ𝑀, ℎ𝑊)≤(1 + 𝜖)𝜒(V(m), ℎ𝑀(m), ℎ𝑊(m)). Otherwise, the procedure ends. The main issue is now the construction and implementation of the minimiser operator O(V, ℎ𝑀, ℎ𝑊) = arg min ℎ𝑀>0, ℎ𝑊>0 𝜒(W, ℎ𝑀, ℎ𝑊). that optimises the two parameters for a given stencil W. 4.2.1. The fminsearch optimisation procedure A first minimisation procedure for O(W, ℎ𝑀, ℎ𝑊)recurs to the MATLAB®function fminsearch. This function finds the minimum of an unconstrained multi-variable function without recurring to numerical or analytic gradients and uses a simplex search method of Lagarias et al. [17]. This kind of approach is adequate for our minimisation problem, since the condition function 𝜒(W, ℎ𝑀, ℎ𝑊)is only piecewise differentiable. 4.2.2. A cloud-like optimisation procedure We propose an alternative optimisation procedure to the fminsearch method. We take advantage of the piecewise differential character of 𝜒(V, ℎ𝑀, ℎ𝑊)by designing a derivative-free procedure to compute the optimal parameters ( ℎ𝑀, ℎ𝑊)(see Fig. 1). Given W(m+1),ℎ𝑀(m),ℎ𝑊(m), we shall build a sub-sequence ℎ𝑟 𝑀=ℎ𝑀(m, 𝑟),ℎ𝑟 𝑊=ℎ𝑊(m, 𝑟)together with 𝛥𝑟 𝑀=𝛥ℎ𝑀(m, 𝑟), 𝛥𝑟 𝑊=𝛥ℎ𝑊(m, 𝑟)and the rectangle 𝑆𝑟= [ℎ𝑟 𝑀−𝛥𝑟 𝑀, ℎ𝑟 𝑀+𝛥𝑟 𝑀]×[ℎ𝑟 𝑊−𝛥𝑟 𝑊, ℎ𝑟 𝑊+𝛥𝑟 𝑊] such that (ℎ𝑟 𝑀, ℎ𝑟 𝑊)→( ℎ𝑀, ℎ𝑊)with 𝛥𝑟 𝑀, 𝛥𝑟 𝑊→0. To this end, we use a cloud of L×Lpoints (Lis an odd, user-defined, integer number) uniformly spanned over the rectangle 𝑆𝑟. We order the points of the cloud (ℎ𝑟,𝓁 𝑀, ℎ𝑟,𝓁 𝑊) ∈ 𝑆𝑟, with index 𝓁= 1,…,L2, such that 𝜒𝑟,𝓁=𝜒(W, ℎ𝑟,𝓁 𝑀, ℎ𝑟,𝓁 𝑊)is a non-decreasing sequence regarding 𝓁. From the sub-sequence, we extract the discrete minimum 𝓁= 1 ℎ𝑟+1 𝑀=ℎ𝑟,1 𝑀, ℎ𝑟+1 𝑊=ℎ𝑟,1 𝑊, while we update the size of the rectangle with the intermediate point 𝓁=L: 𝛥𝑟+1 𝑀=|ℎ𝑟,1 𝑀−ℎ𝑟,L 𝑀|, 𝛥𝑟+1 𝑊=|ℎ𝑟,1 𝑊−ℎ𝑟,L 𝑊|. The procedure ends when both ratios, 𝛥𝑟 𝑀∕ℎ𝑟 𝑀and 𝛥𝑟 𝑊∕ℎ𝑟 𝑊, are smaller than a prescribed value 𝜖𝐶. Remark 3. To initialise the procedure, we take ℎ𝑀(m,0) = ℎ𝑀(m)and ℎ𝑊(m,0) = ℎ𝑊(m)while 𝛥ℎ𝑀(m,0) = 𝜃ℎ𝑀(m), 𝛥ℎ𝑊(m,0) = 𝜃ℎ𝑊(m)with 𝜃∈ ]0,1[ a user constant. Remark 4. To ensure the convergence of the sequence and the positivity of the ℎvalues, we add a threshold correction to ensure that 𝛥𝑟+1 𝑀= [(1 − 𝜃)ℎ𝑟+1 𝑀, 𝜃ℎ𝑟+1 𝑀]. To this end, we perform a correction stage with 𝛥𝑟+1 𝑀∶= max((1 − 𝜃)ℎ𝑟+1 𝑀,min(𝜃ℎ𝑟+1 𝑀, 𝛥𝑟+1 𝑀)). We proceed similarly for ℎ𝑊. Such a correction guarantees the convergence of the procedure, since we automatically reduce the size of the interval of a factor 𝜃. Note that all the benchmarks we present in the numerical section have been carried out with 𝜃= 0.99. Mathematics and Computers in Simulation 218 (2024) 49–78 55 S. Clain and J. Figueiredo Table 1 Approximation errors and convergence orders for the random distribution case with SSO for 𝑝= 2. Approximations for the firstand second-order derivatives are considered for different scales (𝑆) leading to an average converge order denoted ACO (for 𝑆= 1 →1∕16). The values of parameters ℎ𝑀and ℎ𝑊are also presented. 𝑆1 1/2 1/4 1/8 1/16 1/32 1/64 ACO 𝜕𝑥 Error 1.05e−01 2.50e−02 6.17e−03 1.54e−03 3.83e−04 9.58e−05 2.39e−05 Order – 2.1 2.0 2.0 2.0 2.0 2.0 2.0 𝜕𝑦 Error 2.55e−01 6.20e−02 1.54e−02 3.85e−03 9.61e−04 2.40e−04 6.01e−05 Order – 2.0 2.0 2.0 2.0 2.0 2.0 2.0 𝜕𝑥𝑥 Error 6.29e−02 1.45e−02 3.01e−03 4.43e−04 4.62e−04 9.03e−05 6.20e−05 Order – 2.1 2.3 2.8 −0.1 2.4 0.5 1.8 𝜕𝑦𝑦 Error 2.46e−02 6.74e−02 4.75e−02 2.73e−02 1.45e−02 7.48e−03 3.79e−03 Order – −1.5 0.5 0.8 0.9 1.0 1.0 0.2 ℎ𝑀=ℎ𝑊3.22e−01 1.61e−01 8.04e−02 4.02e−02 2.01e−02 1.01e−02 5.03e−03 Table 2 Condition number of matrix 𝑀𝑊 𝑀tand stencil excess size ratio (SESR) for the random distribution case with SSO, before and at the end of the SSO procedure, for different 𝑝values. The values of the condition number and SESR are independent of the scale value 𝑆. 𝑝Condition number SESR (%) Initial Final CNR Initial Final 2 1.7e+01 1.2e+01 1.4 100 17 3 2.2e+02 1.8e+02 1.2 100 10 4 1.1e+03 9.6e+02 1.1 100 7 5 7.4e+03 9.0e+03 0.8 100 10 6 3.3e+04 5.9e+04 0.6 100 0 5. Benchmarks We assess the efficiency of the parameters’ optimisation to produce derivatives with better accuracy and using a smaller stencil. We first deal with the optimisation regarding the stencil and then consider the full optimisation involving the three parameters. Since the optimisation procedure only involves a node and its neighbours, we drop the index 𝑖and use a local indexation where the reference node is 𝑖= 1 and the stencil reads V. Remark. For a better readability of the paper, Tables 28 to 97 are presented in an dedicated appendix. 5.1. Numerical benchmarks (SSO) In a first set of benchmarks, we address the approximation accuracy, convergence order and condition number of 𝑀𝑊 𝑀t provided by the SSO procedure for different geometric distributions of the grid points. In this context, we deal with the approximation of the function 𝑒𝑥+2𝑦at the origin of the coordinate system, (𝑥, 𝑦) = (0,0), using a set of points occupying the domain [−1,1]×[−1,1]. In order to evaluate the convergence order, the original domain is scaled by a factor 𝑆= 1,1∕2,1∕4,1∕8,1∕16,1∕32,and 1∕64. We recall that set Acharacterises the constraints we enforce to provide an expected accuracy. We then use the notation A(𝑝)={𝛼∈N2,|𝛼|≤𝑝}, and by extension, we say that the reconstruction is of 𝑝+ 1th order when using A(𝑝)as the space of constraints. 5.1.1. Random distribution The first benchmark deals with an initial stencil, V, constituted of points randomly distributed around the origin (reference point). We aim at reducing the stencil size while maintaining the condition number within an interval controlled by the tolerance parameter 𝜖. To check the convergence, we produce a sequence of stencils by re-scaling the initial stencil with the scale factor 𝑆, and we report in Table 1 the convergence order for the firstand second-order derivatives for A(2). Similarly, Tables 28–31 present the convergence results obtained for A(𝑝),𝑝= 3,4,5,6respectively. Optimal error convergences are generally achieved for scale factors up to 1∕64, with a 𝑝th-order of convergence for 𝜕𝑥and 𝜕𝑦derivatives and (𝑝−1)th-order for 𝜕𝑥𝑥 and 𝜕𝑦𝑦 derivatives. The convergence order shows an erratic behaviour for 𝑝= 6 when 𝑆= 1∕64 due to loss of accuracy associated to truncation errors that contaminate the results. We note, however, that for 𝑝= 6 the approximation errors for 𝑆= 1∕32 amount to less than 2 × 10−12 and 7 × 10−10 for the firstand second-derivative cases, respectively. We plot in Fig. 12 the final stencil after applying the reduction of the number of points for 𝑝= 2,…,6. We have set 𝜖= 0.2 to alleviate the too restrictive condition on the conditioning leading to a poor reduction of the stencil, and will use this value Mathematics and Computers in Simulation 218 (2024) 49–78 56 S. Clain and J. Figueiredo Table 3 Approximation errors and convergence orders for the uniform distribution case with SSO for 𝑝= 2. Approximations for the firstand second-order derivatives are considered for different scales (𝑆) leading to an average converge order denoted ACO (for 𝑆= 1 →1∕16). The values of parameters ℎ𝑀and ℎ𝑊are also presented. 𝑆1 1/2 1/4 1/8 1/16 1/32 1/64 ACO 𝜕𝑥 Error 4.07e−02 1.05e−02 2.67e−03 6.75e−04 1.70e−04 4.25e−05 1.06e−05 Order – 2.0 2.0 2.0 2.0 2.0 2.0 2.0 𝜕𝑦 Error 1.00e−01 2.59e−02 6.61e−03 1.67e−03 4.20e−04 1.05e−04 2.64e−05 Order – 1.9 2.0 2.0 2.0 2.0 2.0 2.0 𝜕𝑥𝑥 Error 1.21e−01 6.86e−02 3.65e−02 1.88e−02 9.56e−03 4.82e−03 2.42e−03 Order – 0.8 0.9 1.0 1.0 1.0 1.0 0.9 𝜕𝑦𝑦 Error 5.51e−01 3.22e−01 1.74e−01 9.07e−02 4.62e−02 2.33e−02 1.17e−02 Order – 0.8 0.9 0.9 1.0 1.0 1.0 0.9 ℎ𝑀=ℎ𝑊3.61e−01 1.80e−01 9.02e−02 4.51e−02 2.26e−02 1.13e−02 5.64e−03 Table 4 Condition number of matrix 𝑀𝑊 𝑀tand stencil excess size ratio (SESR) for the uniform distribution case with SSO, before and at the end of the SSO procedure, for different 𝑝values. The values of the condition number and SESR are independent of the scale value 𝑆. 𝑝Condition number SESR (%) Initial Final CNR Initial Final 2 2.0e+01 2.9e+01 0.7 100 33 3 5.6e+01 1.3e+02 0.4 100 10 4 3.9e+02 5.0e+02 0.8 100 0 5 3.3e+03 5.5e+03 0.6 100 14 6 9.6e+03 2.7e+04 0.4 100 0 throughout the simulations. The number of points for the initial stencil is the double of the minimal number of points to provide a matrix 𝑀of maximal rank, i.e.,|V|= 12 for 𝑝= 2,|V|= 20 for 𝑝= 3,|V|= 30 for 𝑝= 4,|V|= 42 for 𝑝= 5, and |V|= 56 for 𝑝= 6. This will be the setting for all the simulations performed. We observe a relative uniform distribution of the selected points (red filled circles) regarding the distance and angular configuration. Remark. Notice that the analysis of the location of the selected points in the different panels of Fig. 12, particularly those corresponding to 𝑝= 3,…,6, shows that the SSO methodology manages to avoid the selection of points that are very close to each other by discarding at least one of the points involved. We report in Table 2 the condition number in function of 𝑝before and after the stencil reduction. For 𝑝= 2,3,4, we manage to reduce the condition together with the stencil size (CNR up to 1.4) while for higher values of 𝑝, the stencil size reduction implies a moderate increase in the condition number (CNR down to 0.6). To quantify the stencil size reduction, we define the stencil excess size ratio (SESR) as SESR(V) = |V|−|V𝑚𝑖𝑛| |V𝑚𝑖𝑛|= 2 SSR(V)−1, where |V𝑚𝑖𝑛|is the minimal number of points to provide maximal rank. For example, for 𝑝= 2, the minimal number of points is |V𝑚𝑖𝑛|= 6 and the stencil obtained has 17% more points (i.e.,1). 5.1.2. Uniform distribution The uniform case concerns a uniform distribution of the points, as depicted in Fig. 13. We proceed similarly as we did for the random case and report in Tables 3 (𝑝= 2) and 32–35 (𝑝= 3,…,6) the errors and convergence order for the different values of 𝑝 for the firstand second-order derivatives. The results show evidence of the optimal convergence order of the reconstructions for all the scale factors, except for the 𝑝= 6 and 𝑆= 1∕64 case, where the approximation error for the firstand second-order derivatives stays below 3 × 10−12 and 10−9, respectively. We plot in Fig. 13 the final stencil after size reduction for 𝑝= 2,…,6. We notice the good symmetry of the points’ distribution with a regular distribution around the reference node. Table 4 gives the condition numbers before and after the stencil reduction for the different values of 𝑝. The SESR value is in the range [0%,33%], while the conditioning is well-controlled, only moderately larger than the initial value (CNR ranging from 0.4 to 0.8). Finally, and for the sake of comparison, we present the results obtained for the approximation errors and convergence order, for both firstand second-order derivatives, using the standard (3 points) one-dimension second-order centred differences (cf. Table 36) and the standard (5 points) one-dimension fourth-order centred differences (cf. Table 37) schemes. The convergence order is the one expected in both cases, while the magnitude of the approximation errors is smaller than the one obtained with the SSO procedure Mathematics and Computers in Simulation 218 (2024) 49–78 57 S. Clain and J. Figueiredo Table 5 Approximation errors and convergence orders for the half-plane distribution case with SSO for 𝑝= 2. Approximations for the firstand second-order derivatives are considered for different scales (𝑆) leading to an average converge order denoted ACO (for 𝑆= 1 →1∕16). The values of parameters ℎ𝑀and ℎ𝑊are also presented. 𝑆1 1/2 1/4 1/8 1/16 1/32 1/64 ACO 𝜕𝑥 Error 6.55e−02 1.55e−02 3.77e−03 9.30e−04 2.31e−04 5.75e−05 1.44e−05 Order – 2.1 2.0 2.0 2.0 2.0 2.0 2.0 𝜕𝑦 Error 1.23e−01 2.78e−02 6.63e−03 1.62e−03 4.01e−04 9.98e−05 2.49e−05 Order – 2.1 2.1 2.0 2.0 2.0 2.0 2.1 𝜕𝑥𝑥 Error 6.32e−02 3.10e−02 1.53e−02 7.59e−03 3.78e−03 1.89e−03 9.43e−04 Order – 1.0 1.0 1.0 1.0 1.0 1.0 1.0 𝜕𝑦𝑦 Error 1.15e+00 4.97e−01 2.31e−01 1.12e−01 5.48e−02 2.72e−02 1.35e−02 Order – 1.2 1.1 1.0 1.0 1.0 1.0 1.1 ℎ𝑀=ℎ𝑊3.51e−01 1.75e−01 8.77e−02 4.39e−02 2.19e−02 1.10e−02 5.48e−03 Table 6 Condition number of matrix 𝑀𝑊 𝑀tand stencil excess size ratio (SESR) for the half-plane distribution case with SSO, before and at the end of the SSO procedure, for different 𝑝values. The values of the condition number and SESR are independent of the scale value 𝑆. 𝑝Condition number SESR (%) Initial Final CNR Initial Final 2 3.7e+02 3.0e+02 1.2 100 0 3 1.2e+04 8.8e+03 1.4 100 10 4 1.9e+05 1.3e+05 1.5 100 0 5 4.5e+06 5.8e+06 0.8 100 0 6 4.0e+08 2.5e+08 1.6 100 0 having the same order of convergence. We report that the second-order schemes yield very similar results for the first-order derivative with respect to 𝑦, whereas within a factor 3for both the first-order derivative with respect to 𝑥and the second-order derivative with respect to 𝑦, while within a factor 50 for the second-order derivative with respect to 𝑥. The error ratio for the fourth-order schemes is higher, ranging from 4(first-order derivative with respect to 𝑦) to 135 (second-order derivative with respect to 𝑥). In any case, it should be noted that the fact that the function whose derivatives are being approximated has a separated form for the independent variables favours the one-dimensional character of the standard schemes, helping to justify the performance gap between the two families of schemes. 5.1.3. The half-plane distribution To check the method for a node close to the boundary, we now assume that all points of the stencil belong to the right half-plane, randomly distributed. Such a situation arises when the reference node belongs to the boundary and all the information is located inside the right (or left) domain (see Fig. 14). Following the previous examples, we built an initial stencil with |V|= 12,20,30,42 and 56 points for 𝑝= 2,3,4,5and 6, respectively, and re-scaled the stencil setting 𝑆= 1,1∕2,…,1∕64 to assess the error and the convergence order. The corresponding results are reported in Tables 5 (𝑝= 2) and 38–41 (𝑝= 3,…,6) for the different values of 𝑝. Optimal effective order is obtained for both the firstand second-order derivatives, except for the 𝑝= 5,6cases when 𝑆= 1∕32,1∕64. Again, we observe very low approximation errors associated with the misbehaviour of the convergence. Fig. 14 displays the initial and final stencil for the five reconstruction settings. Uniform distribution of the selected points is again noticeable. Table 6 provides the condition number and the stencil excess size ratio. We manage to reduce the stencil cardinal by a factor two with respect to the initial setting, while the condition number is generally reduced (CNR value below 1 only for 𝑝= 5). On the other hand, we note that the condition number is larger than the one corresponding to the random distribution case on the whole space (roughly, one magnitude higher). Indeed, up-winding naturally results from this one-side configuration and strongly impacts the reconstruction. 5.1.4. The regular convex distribution Boundary conditions may be sensitive to the curvature of the simulation domain frontier, and therefore the assessment of the firstand second-order derivatives is important to check the accuracy of the method for such situations. This paragraph deals with the regular convex case (displayed in Fig. 15), whereas the more sensitive regular concave situation is presented in the next paragraph. The boundary is represented as a piece of circle of unit radius and all the points are situated inside the circle. We re-scale the initial configuration with parameter 𝑆and, therefore, the radius varies with 𝑆(larger curvature when 𝑆decreases). Indeed, if the radius would not change with 𝑆, we would quickly recover the half-plane case since the distance between the points is one or two magnitudes lower than the radius. We report in Tables 7 (𝑝= 2) and 42–45 (𝑝= 3,…,6) the error and convergence order obtained for the different values of 𝑝 with respect to the scale factor 𝑆. Similarly to the former cases, the optimal converge is achieved with a 𝑝th-order for the first-order derivative and (𝑝− 1)th-order for the second-order case as long as 𝑆≥1∕16. We display in Fig. 15 the location of the initial stencil Mathematics and Computers in Simulation 218 (2024) 49–78 64 S. Clain and J. Figueiredo Table 20 Condition number of matrix 𝑀𝑊 𝑀t, condition number ratio (CNR) and scaling parameters (ℎ𝑀, ℎ𝑊)for the convex corner distribution case for different optimisation procedures and scale values (𝑆) when 𝑝= 2. Optimisation 𝑆1 1/2 1/4 1/8 1/16 1/32 1/64 stencil size opt. ℎ𝑀3.35e−01 1.67e−01 8.37e−02 4.19e−02 2.09e−02 1.05e−02 5.23e−03 ℎ𝑊3.35e−01 1.67e−01 8.37e−02 4.19e−02 2.09e−02 1.05e−02 5.23e−03 cond. nb. 2.3e+03 CNR 1.7 cloud opt. ℎ𝑀1.95e−01 9.75e−02 4.87e−02 2.44e−02 1.22e−02 6.09e−03 3.05e−03 ℎ𝑊3.64e−01 1.82e−01 9.11e−02 4.55e−02 2.28e−02 1.14e−02 5.69e−03 cond. nb. 1.4e+03 CNR 3.0 fminsearch opt. ℎ𝑀2.21e−01 1.11e−01 5.53e−02 2.76e−02 1.38e−02 6.91e−03 3.46e−03 ℎ𝑊3.96e−01 1.98e−01 9.91e−02 4.95e−02 2.48e−02 1.24e−02 6.19e−03 cond. nb. 1.6e+03 CNR 2.6 Fig. 3. Generic grid. The open bounded domain 𝛺and its boundary 𝜕𝛺 =𝛤𝐷∪𝛤𝑁∪𝛤𝑅, where 𝛤𝐷,𝛤𝑁and 𝛤𝑅stand, respectively, for the boundary partition where Dirichlet, Neumann and Robin boundary conditions hold. The green filled circles, red squares, blue triangles and orange stars correspond to grid points belonging to 𝛺,𝛤𝐷,𝛤𝑁and 𝛤𝑅, respectively. 6. Generalised finite difference schemes We consider the partial differential equation E(𝜙)=−𝜅1(𝜕(2,0)𝜙+𝜕(0,2)𝜙)+𝜅2𝜕(1,0)𝜙+𝜅3𝜕(0,1)𝜙=𝑓, supplemented with the (generic) Robin boundary condition B(𝜙) = 𝛾𝐷𝜙+𝛾𝑁(𝑛𝑥𝜕(1,0)𝜙+𝑛𝑦𝜕(0,1)𝜙)=𝑔, where 𝜅1,𝜅2,𝜅3and 𝑓are real valued functions defined on 𝛺, while 𝛾𝐷,𝛾𝑁and 𝑔are defined on the boundary 𝜕𝛺. 6.1. Discrete formulation At the discrete level, we substitute the derivatives by the corresponding discretisation relations, and we get for 𝑖∈ G (−𝜅1(𝑥𝑖)(𝑎𝑖,(2,0) +𝑎𝑖,(0,2)) + 𝜅2(𝑥𝑖)𝑎𝑖,(1,0) +𝜅3(𝑥𝑖)𝑎𝑖,(0,1))𝛷V𝑖=𝑓(𝑥𝑖), where 𝑎𝑖,𝛽 is the coefficient of the 𝛽derivative discretisation at node 𝑖. On the other hand, the boundary condition at the grid point 𝑖∈𝜕Greads 𝛾𝐷(𝑥𝑖)𝜙𝑖+𝛾𝑁(𝑥𝑖)(𝑛𝑥(𝑥𝑖)𝑎𝑖,(1,0) +𝑛𝑦(𝑥𝑖)𝑎𝑖,(0,1))𝛷V𝑖=𝑔(𝑥𝑖). Fig. 3 illustrates a generic cloud of points used in the simulations. To compute the weights’ matrix 𝑊, we adopt the Gaussian kernel [7] 𝜔(𝑞, ℎ𝑊) = 1 𝜋ℎ2 𝑊 𝑒−𝑞2, 𝑞 =||𝑥𝑖−𝑥𝑗|| ℎ𝑊 . This function is smooth and considered very stable and accurate [19]. Mathematics and Computers in Simulation 218 (2024) 49–78 65 S. Clain and J. Figueiredo Fig. 4. Exact solution for problem (BVP1). Assembling all the relations for 𝑖∈Gprovides the linear system 𝐺𝛷 =𝐹. Additionally, we use the MATLAB®function “equilibrate” [5], to reduce considerable the condition number of the global system (see further below). It provides a pre-conditioning linear system  𝐺𝛹 = 𝐹we solve with a direct method. To assess the convergence, we compute the 𝐿1and 𝐿∞-errors by 𝐿1-error: 1 𝐼∑ 𝑖∈G|𝜙𝑖−𝜙𝑒𝑥 𝑖|and 𝐿∞-error: max 𝑖∈G|𝜙𝑖−𝜙𝑒𝑥 𝑖|, where (𝜙𝑒𝑥 𝑖)𝑖∈Gand (𝜙𝑖)𝑖∈Gare the exact and the approximated solution, respectively. 6.2. Numerical tests We consider the mixed boundary condition problem ⎧ ⎪ ⎨ ⎪ ⎩ −(𝜕𝑥𝑥𝜙+𝜕𝑦𝑦𝜙)+𝜕𝑥𝜙+ 2𝜕𝑦𝜙= 0,for all (𝑥, 𝑦) ∈ 𝛺, Dirichlet/Neumann conditions, for all (𝑥, 𝑦) ∈ 𝜕𝛺 (BVP1) for which the exact solution is given by 𝜙=𝑒𝑥+2𝑦,for all (𝑥, 𝑦) ∈ 𝛺. As an illustration, we plot the solution over the domain [0,1] × [0,1] in Fig. 4. We also tackle the problem ⎧ ⎪ ⎨ ⎪ ⎩ −(𝜕𝑥𝑥𝜙+𝜕𝑦𝑦𝜙)+𝑥𝜕𝑥𝜙+𝑦𝜕𝑦𝜙=𝑓, for all (𝑥, 𝑦) ∈ 𝛺, Dirichlet conditions, for all (𝑥, 𝑦) ∈ 𝜕𝛺 (BVP2) for which we impose the following exact solution 𝜙= cos(2𝜋𝑥) sin(2𝜋𝑦),for all (𝑥, 𝑦) ∈ 𝛺, by manufacturing the second member, 𝑓, accordingly. 6.2.1. Test case 1 Let 𝛺=] − 1,1[×] − 1,1[ and 𝛤𝐷=𝜕𝛺. We seek a numerical solution of BVP1 with Dirichlet boundary conditions corresponding to set 𝜙=𝑒𝑥+2𝑦, for all (𝑥, 𝑦) ∈ 𝛤𝐷. Numerical simulations have been carried out with four clouds of points using uniform grids (GU-1) with 117,437,1677, and 7565 points. A second set of clouds (GP15-1) is obtained by performing an up to 15% random perturbation on the previous clouds GU-1. At last, another set of clouds (GP30-1) is obtained as an up to 30% random perturbation of clouds GU-1. On the other hand, a set of clouds (GD-1) is obtained with a Delaunay-like triangulation procedure with 209,815,3193, and 12 378 points, using the Gmsh© package. Fig. 5 illustrates the different sets of clouds used in the simulations. Mathematics and Computers in Simulation 218 (2024) 49–78 66 S. Clain and J. Figueiredo Fig. 5. Grids for the test case 1. For 𝛺={(𝑥, 𝑦) ∈ R2∶𝑥, 𝑦 ∈ [−1,1]}the panels illustrate, from left to right, a uniform grid, two grids with perturbation levels up to 15% and 30% with respect to the uniform case, and a grid obtained with a Delaunay-like triangulation procedure. The green filled circles and the red squares correspond, respectively, to interior points and boundary points where Dirichlet conditions are imposed. GU-1 tests. We carry out the simulations using the GU-1 clouds with four optimisation scenarios: no optimisation, stencil size optimisation, cloud-like optimisation and fminsearch optimisation. Errors and convergence orders, as well as for the condition number of 𝐺 and  𝐺are presented in Tables 21 (for 𝑝= 2) and 82–85 (for 𝑝= 3,…,6, respectively). For 𝑝= 2 and 3, we report that sticking to the initial stencil, convergence orders are far from optimal, while the condition number of 𝐺can reach values as high as 2 × 1010. On the other hand, for 𝑝= 4,5and 6, the convergence order of this scheme is close to optimal, but the 𝐿∞- and 𝐿1-errors are considerably higher than the ones obtained with the three other optimisation procedures. The SSO approach also leads to poor convergence orders for the 𝑝= 2 and 3cases, while for 𝑝= 4,5and 6, we observe convergence orders that are in line with the (optimal) ones obtained with the cloud and fminsearch approaches, but the errors and the condition number are well above the values obtained using its SPO counterparts. Comparing the results obtained with the cloud-like and the fminsearch methods, we observe that the errors for both the numerical solution and the first-order derivatives are generally lower for the fminsearch approach. On the other hand, the condition number of 𝐺(and  𝐺) shows a weak dependency on the optimisation method used. We note that, in both cases, the use of the “equilibrate” procedure allows reducing the condition number of the system matrix up to three orders of magnitude for the finest grid. Given, on the one hand, the poor quality of the results obtained with the initial stencil approach and (to a less extent) the SSO method for the accuracy and convergence order of the numerical solution when compared with the SPO methods and, on the other hand, the fact that the two SPO algorithms lead to similar results for the numerical solution and first-order derivatives, being the fminsearch method much faster than the cloud-like algorithm (see further below), the numerical simulations for the remaining three sets of non-uniform grids are performed exclusively with the fminsearch method. GU15-1, GU30-1 and Delaunay tests. The results obtained for clouds GP15-1, GP30-1 and GD-1 (cf. Tables 22 for 𝑝= 2 and 86–89) for 𝑝= 3,…,6present an optimal convergence order for both the 𝐿∞- and the 𝐿1-norm, while we report an average value for the ratio between the 𝐿∞-error and the 𝐿1-error is in the 95% confidence interval 7.3±3.0, indicating a good distribution of the error throughout the computational domain. As far as the condition number of 𝐺and  𝐺are concerned we observe, as expected, that they increase with the number of degrees of freedom, while using the “equilibrate” formulation results in a reduction of the condition number between one and two orders of magnitude for the fminsearch method. Finally, we address, on the one hand, the computation time associated to the preprocessing computations, that is, to the choice of the stencils and the computation of the corresponding local matrices, and, on the other hand, the time needed to compute the problem solution, a process that involves the assembly of the local matrices as well as the solution of the sparse system of equations. In the former case, we consider, for the set of uniform grids GU-1, the average preprocessing CPU time per degree of freedom needed for the different optimisation approaches considered in the present work for values of 𝑝ranging from 2to 6. The results obtained, presented in Table 23 and depicted in Fig. 6, show that the logarithm of the preprocessing CPU time scales roughly in a linear manner with the approximation degree independently of the optimisation procedure involved. Additionally, and as expected, the absence of optimisation results in the fastest method, while the homemade cloud-like optimisation method is clearly the slowest. It can also be seen that the spread of preprocessing CPU times regarding the optimisation method used tends to increase with 𝑝. At last, it is worth mentioning that the optimisation based on the fminsearch procedure is, on average, 3 to 5 times slower than the (rather basic) SSO approach. Turning now our attention to the latter case, we consider the relation between the accuracy of the numerical solution and the CPU time needed to determine the problem solution. For that we consider once again the set of uniform grids GU-1 and the two extreme values of the approximation degree, 𝑝= 2 and 𝑝= 6, for illustration purposes. The results obtained are presented in Tables 24 and 25 and depicted in Fig. 7. The results show that for a given CPU time the accuracy provided by the cloud-like and fminsearch optimisation methods is rather similar and up to two orders of magnitude better than the no optimisation procedure, while up to one order of magnitude better than the accuracy provided by the SSO technique for the finer grids. Mathematics and Computers in Simulation 218 (2024) 49–78 67 S. Clain and J. Figueiredo Table 21 Results for the test case 1 considering four uniform grids and 𝑝= 2. The error and convergence order for the numerical solution, as well as the condition number of the system matrix (𝐺) and its regularised counterpart (  𝐺), are presented for different optimisation procedures. For the cloud and the fminsearch optimisation procedures the error and convergence order for the first-order derivatives are also presented. Optimisation DOF 117 437 1677 7565 ACO no opt. Solution L∞-error 7.40e−01 3.02e−01 7.78e−02 2.66e−01 Order – 1.4 2.0 −1.6 0.5 L1-error 1.38e−01 6.01e−02 1.24e−02 1.01e−02 Order – 1.3 2.3 0.3 1.3 cond. nb. G 1.3e04 4.7e05 1.7e07 2.1e10 cond. nb. ˜ G 3.3e02 2.8e03 2.6e04 6.8e06 stencil size opt. Solution L∞-error 1.57e−01 5.04e−02 2.26e−02 8.89e−03 Order – 1.7 1.2 1.2 1.4 L1-error 2.40e−02 1.04e−02 5.69e−03 2.29e−03 Order – 1.3 0.9 1.2 1.1 cond. nb. G 1.9e02 1.0e03 5.0e03 3.5e04 cond. nb. ˜ G 1.3e01 5.1e01 2.0e02 9.5e02 cloud opt. Solution L∞-error 1.16e−01 1.62e−02 3.38e−03 1.15e−03 Order – 3.0 2.3 1.4 2.2 L1-error 1.25e−02 9.85e−04 5.62e−04 1.98e−04 Order – 3.9 0.8 1.4 2.1 𝜕𝑥 L∞-error 8.85e−01 2.45e−01 8.10e−02 2.13e−02 Order – 1.9 1.6 1.8 1.8 L1-error 8.80e−02 2.17e−02 5.16e−03 1.21e−03 Order – 2.1 2.1 1.9 2.1 𝜕𝑦 L∞-error 1.33e00 5.93e−01 1.99e−01 5.15e−02 Order – 1.2 1.5 1.8 1.6 L1-error 1.80e−01 4.56e−02 1.07e−02 2.27e−03 Order – 2.1 2.2 2.1 2.1 cond. nb. G 2.2e02 1.1e03 6.2e03 4.2e04 cond. nb. ˜ G 1.9e01 7.4e01 2.9e02 1.3e03 fminsearch opt. Solution L∞-error 5.93e−02 9.34e−03 2.11e−03 4.39e−04 Order – 2.8 2.2 2.1 2.4 L1-error 1.10e−02 2.60e−03 6.26e−04 1.33e−04 Order – 2.2 2.1 2.1 2.1 𝜕𝑥 L∞-error 6.94e−01 2.52e−01 7.48e−02 1.75e−02 Order – 1.5 1.8 1.9 1.8 L1-error 5.93e−02 1.58e−02 4.02e−03 8.77e−04 Order – 2.0 2.0 2.0 2.0 𝜕𝑦 L∞-error 9.99e−01 3.62e−01 8.79e−02 2.62e−02 Order – 1.5 2.1 1.6 1.7 L1-error 1.28e−01 3.34e−02 8.52e−03 1.86e−03 Order – 2.0 2.0 2.0 2.0 cond. nb. G 2.0e02 1.1e03 6.1e03 4.1e04 cond. nb. ˜ G 1.9e01 7.5e01 3.0e02 1.4e03 6.2.2. Test case 2 The computational domain is an annulus having internal and external radius 1∕3 and 1, respectively, that is 𝛺={(𝑥, 𝑦)=(𝜌cos 𝜃, 𝜌 sin 𝜃) ∈ R2∶𝜌∈]1∕3,1[, 𝜃 ∈[0,2𝜋[}, where the inner and outer boundaries are given by 𝛤𝐷={(𝑥, 𝑦) = (cos 𝜃, sin 𝜃) ∈ R2∶𝜃∈[0,2𝜋[}, 𝛤𝑁={(𝑥, 𝑦) = (cos 𝜃, sin 𝜃)∕3 ∈ R2∶𝜃∈[0,2𝜋[}, We seek the numerical solution of BVP1 with the Dirichlet and Neumann boundary conditions corresponding to set 𝜙=𝑒𝑥+2𝑦, for all (𝑥, 𝑦) ∈ 𝛤𝐷, 𝑛𝑥𝜕𝑥𝜙+𝑛𝑦𝜕𝑦𝜙=𝑒𝑥+2𝑦(𝑛𝑥+ 2𝑛𝑦),for all (𝑥, 𝑦) ∈ 𝛤𝑁. The geometrical setting and the exact solution of the BVP, 𝜙=𝑒𝑥+2𝑦, are illustrated in Fig. 8. To assess the numerical accuracy, we use a set of four uniform clouds (U-2) with 190,740,2938, and 11 257 points. A second set of four clouds (P15-2) is obtained by performing an up to 15% random perturbation of a uniform grid, while a set of four grids (P30-2) with an up to 30% random perturbation is also used. The number of points of the four successive clouds (P15-2) (and also (P30-2)) are 168,711,2860, and 11 442 points. At last, a set of four grids (GD-2) with 228,880,3369, and 11 888 points is obtained thanks to a Delaunay-like triangulation procedure. Fig. 9 illustrates the different sets of grids used in the test case 2. All the simulations are carried out with the SPO fminsearch procedure. Mathematics and Computers in Simulation 218 (2024) 49–78 68 S. Clain and J. Figueiredo Fig. 6. Test case 1. Average preprocessing time (in CPU seconds) per degree of freedom as a function of 𝑝and optimisation procedure (see Table 23). Fig. 7. Test case 1. Relation between the numerical solution accuracy (L1-error) and the computation time (in CPU seconds) for different optimisation procedures for four uniform grids with 117,437,1677 and 7565 points. The left panel stands for 𝑝= 2 while the right panel corresponds to 𝑝= 6 (see Tables 24 and 25, respectively). Fig. 8. Test case 2. Geometrical setting and exact solution. Mathematics and Computers in Simulation 218 (2024) 49–78 69 S. Clain and J. Figueiredo Table 22 Results for the test case 1 for the fminsearch optimisation procedure considering nonuniform grids and 𝑝= 2. The error and convergence order of the numerical solution, as well as the condition number of the system matrix (𝐺) and its regularised counterpart (  𝐺), are presented for three sets of nonuniform grids: two sets of grids with perturbation levels up to 15% and 30% with respect to the uniform case, and a set of grids obtained with a Delaunay triangulation procedure. Grid type DOF 117 437 1677 7565 ACO <15% perturbation Solution L∞-error 1.15e−01 1.87e−02 2.11e−03 7.22e−04 Order – 3.2 1.5 2.1 2.4 L1-error 1.21e−02 2.63e−03 6.26e−04 1.76e−04 Order – 2.2 1.6 2.1 2.0 cond. nb. G 2.3e02 1.3e03 6.1e03 5.4e04 cond. nb. ˜ G 2.1e01 8.0e01 3.0e02 1.6e03 <30% perturbation Solution L∞-error 6.33e−02 2.15e−02 2.11e−03 9.80e−04 Order – 2.4 2.0 2.1 2.0 L1-error 1.14e−02 4.10e−03 6.26e−04 1.85e−04 Order – 2.0 2.3 2.1 2.0 cond. nb. G 2.7e02 1.9e03 6.1e03 7.7e04 cond. nb. ˜ G 2.2e01 9.0e01 3.0e02 1.9e03 DOF 209 815 3193 12378 ACO Delaunay Solution L∞-error 6.93e−02 5.48e−03 1.30e−03 3.35e−04 Order – 3.7 2.1 2.0 2.6 L1-error 5.61e−03 1.41e−03 3.89e−04 9.84e−05 Order – 2.0 1.9 2.0 2.0 cond. nb. G 7.8e02 4.9e03 2.8e04 1.8e05 cond. nb. ˜ G 4.5e01 1.9e02 7.9e02 3.1e03 Table 23 Results for the average preprocessing time (in CPU seconds) per degree of freedom as a function of 𝑝and optimisation procedure for the test case 1 considering four uniform grids with 117, 437, 1677 and 7565 points. Optimisation 𝑝 2 3 4 5 6 no opt. (2.9 ±1.0)e−04 (4.4 ±2.1)e−04 (4.9 ±1.4)e−04 (5.9 ±1.6)e−04 (9.9 ±2.5)e−04 stencil size opt. (4.1 ±1.4)e−04 (6.7 ±1.6)e−04 (1.5 ±0.2)e−03 (3.7 ±0.2)e−03 (9.3 ±0.3)e−03 fminsearch opt. (1.4 ±0.5)e−03 (2.8 ±0.5)e−03 (7.3 ±0.5)e−03 (1.9 ±0.1)e−02 (1.1 ±0.2)e−01 cloud opt. (3.0 ±0.1)e−02 (8.0 ±0.7)e−02 (5.0 ±0.2)e−01 (1.6 ±0.1)e−00 (3.8 ±0.1)e−00 Table 24 Results for the computation time (in CPU seconds) and solution L1-error as a function of the optimisation procedure for the test case 1 considering 𝑝= 2 and four uniform grids with 117, 437, 1677 and 7565 points. Optimisation DOF 117 437 1677 7565 no opt. L1-error 1.38e−01 6.01e−02 1.24e−02 1.01e−02 CPU time 1.36e−02 3.33e−02 1.83e−01 2.22e−01 stencil size opt. L1-error 2.40e−02 1.04e−02 5.69e−03 2.29e−03 CPU time 1.35e−02 3.16e−02 1.59e−01 1.69e−01 cloud opt. L1-error 1.25e−02 9.85e−04 5.62e−04 1.98e−04 CPU time 1.32e−02 3.19e−02 1.58e−01 1.70e−01 fminsearch opt. L1-error 1.10e−02 2.60e−03 6.26e−04 1.33e−04 CPU time 1.31e−02 3.18e−02 1.58e−01 1.67e−01 Errors and convergence orders, as well as for the condition number of 𝐺and  𝐺, are reported in Tables 26 (𝑝= 2) and 90– 93 (𝑝= 3,…,6). The optimal convergence is achieved for all values of 𝑝independently of the shape of the clouds. The average value for the ratio between the 𝐿∞-error and the 𝐿1-error is higher than in the previous test case, with a 95% confidence interval 13.0±3.4, reflecting the increase in the complexity of the BVP due to the curved geometry and mixed boundary conditions. The error distribution over the computational domain can still be considered good despite the moderate augmentation of this ratio. Comparing the results obtained for each set of grids, one observes that, in general, the smallest errors correspond to the uniform and Delaunay-like grids, while the highest errors are associated with the (P30-2) set of grids. As in the previous test case, we report that the condition number of 𝐺and  𝐺increases with the number of degrees of freedom by up to 3 orders of magnitude. As to the effect of the “equilibrate” procedure, we observe that it can lead to a reduction of the condition number up to 3 orders of magnitude. Mathematics and Computers in Simulation 218 (2024) 49–78 70 S. Clain and J. Figueiredo Table 25 Results for the computation time (in CPU seconds) and solution L1-error for the test case 1. Similar to Table 24, but for 𝑝= 6. Optimisation DOF 117 437 1677 7565 no opt. L1-error 1.91e−01 1.29e−03 1.06e−05 4.36e−08 CPU time 1.82e−02 6.49e−02 4.44e−01 1.84e+01 stencil size opt. L1-error 5.79e−04 2.12e−05 3.38e−07 4.25e−09 CPU time 1.58e−02 5.09e−02 3.40e−01 7.75e+00 cloud opt. L1-error 2.45e−04 3.38e−06 3.90e−08 8.38e−10 CPU time 1.67e−02 5.01e−02 3.29e−01 8.15e+00 fminsearch opt. L1-error 9.59e−05 4.00e−06 5.53e−08 5.92e−10 CPU time 1.74e−02 5.17e−02 3.30e−01 7.42e+00 Table 26 Results for the test case 2 for the fminsearch optimisation procedure and 𝑝= 2. The error and convergence order of the numerical solution, as well as the condition number of the system matrix (𝐺) and its regularised counterpart (  𝐺), are presented for four sets of grid types: a set of uniform grid types, two sets of grid types with perturbation levels up to 15% and 30% with respect to the uniform case, and a set of grid types obtained with a Delaunay triangulation procedure. Grid type DOF 190 740 2938 11 257 ACO Uniform Solution L∞-error 1.95e−01 3.36e−02 5.32e−03 3.03e−03 Order – 2.6 2.7 0.8 2.0 L1-error 1.54e−02 2.17e−03 3.77e−04 9.28e−05 Order – 2.9 2.5 2.1 2.5 cond. nb. G 1.1e03 3.6e03 5.6e04 2.5e05 cond. nb. ˜ G 7.3e01 3.3e02 2.2e03 3.5e04 DOF 168 711 2860 11442 ACO <15% perturbation Solution L∞-error 5.56e−01 9.05e−02 1.62e−02 5.20e−03 Order – 2.5 2.5 1.6 2.2 L1-error 4.02e−02 6.71e−03 1.58e−03 2.17e−04 Order – 2.5 2.1 2.9 2.5 cond. nb. G 1.6e03 7.8e04 3.2e05 6.5e05 cond. nb. ˜ G 9.9e01 4.9e03 8.4e03 1.2e04 <30% perturbation Solution L∞-error 7.56e−01 1.31e−01 2.84e−02 6.77e−03 Order – 2.4 2.2 2.1 2.2 L1-error 2.37e−02 5.26e−03 1.13e−03 2.58e−04 Order – 2.1 2.2 2.1 2.1 cond. nb. G 1.8e03 6.9e04 1.2e06 1.7e06 cond. nb. ˜ G 7.8e01 3.6e03 3.5e04 2.0e04 DOF 228 880 3369 11888 ACO Delaunay Solution L∞-error 1.12e−01 8.45e−03 2.80e−03 5.95e−04 Order – 3.8 1.6 2.5 2.6 L1-error 8.76e−03 1.13e−03 3.14e−04 8.74e−05 Order – 3.1 1.9 2.0 2.3 cond. nb. G 2.7e03 8.0e03 9.7e04 2.0e05 cond. nb. ˜ G 1.4e02 3.5e02 1.4e02 5.0e03 Table 27 Results for the test case 3 for the fminsearch optimisation procedure and 𝑝= 2. The error and convergence order of the numerical solution, as well as the condition number of the system matrix (𝐺) and its regularised counterpart (  𝐺), are presented for four grids obtained with a Delaunay triangulation technique combined with a level-set procedure. DOF 217 847 3365 13 696 ACO Solution L∞-error 2.35e−01 9.41e−02 1.74e−02 5.31e−03 Order – 1.3 2.4 1.7 1.8 L1-error 5.04e−02 1.85e−02 4.58e−03 1.07e−03 Order – 1.5 2.0 2.1 1.9 cond. nb. G 1.1e03 5.7e03 4.3e04 2.4e05 cond. nb. ˜ G 3.4e01 1.4e02 6.6e02 2.8e03 Mathematics and Computers in Simulation 218 (2024) 49–78 71 S. Clain and J. Figueiredo Fig. 9. Grids for the test case 2. For 𝛺={(𝑥, 𝑦)=(𝜌cos 𝜃, 𝜌 sin 𝜃) ∈ R2∶𝜌∈ [1∕3,1], 𝜃 ∈[0,2𝜋[}the panels illustrate, from left to right, a uniform grid, two grids with perturbation levels up to 15% and 30% with respect to the uniform case, and a grid derived from a Delaunay-like triangulation procedure. The green filled circles, the red squares and the blue triangles correspond, respectively, to interior points, and boundary points where Dirichlet and Neumann conditions are imposed. Fig. 10. Test case 3. Geometrical setting and exact solution. 6.2.3. Test case 3 We look at the numerical solution of BVP2 over a complex computational domain to show that one can address numerical approximation with a curved boundary domain given by the level-set 𝛺={(𝑥, 𝑦) ∈ R2∶√𝑥4+𝑦2 16 +(3 2𝑥+𝑦2)2 +𝑦 < 1}. Dirichlet condition is prescribed on the boundary 𝛤𝐷={(𝑥, 𝑦) ∈ R2∶√𝑥4+𝑦2 16 +(3 2𝑥+𝑦2)2 +𝑦= 1}, using the exact solution 𝜙= cos(2𝜋𝑥) sin(2𝜋𝑦). The geometrical setting and the exact solution of the BVP are illustrated in Fig. 10. We carried out the simulations with use four grids (GD-3) of 217,847,3365, and 13 656 points, respectively. In all cases, the interior points are obtained with a Delaunay-like triangulation procedure, while the points on the boundary are set using a level-set technique. Fig. 11 illustrates the coarsest grids used in the test case 3. As in the test case 2, the simulations are performed using the SPO fminsearch procedure. Errors, convergence orders, and condition numbers of 𝐺and  𝐺, are depicted in Tables 27 (𝑝= 2) and 94–97 for 𝑝= 3,…,6. The errors are definitively higher in the last case than in the two previous tests independently of 𝑝. This is due to the increased complexity of both the domain shape and the exact solution of the BVP. Nevertheless, optimal convergence is still achieved for all values of 𝑝despite the increased complexity of the BVP. We also notice that the average value for the ratio between the 𝐿∞-error and the 𝐿1-error is considerably lower than in the previous tests, with a 95% confidence interval 5.0±0.7, reflecting a very good error distribution over the computational domain, particularly for 𝑝= 2,…,5. Similarly to the previous cases, the condition number of 𝐺and  𝐺increases with the number of degrees of freedom by up to 3 orders of magnitude. We also report a reduction of the condition number ranging from 1 to 3 orders of magnitude in the framework of the “equilibrate” procedure. 7. Conclusions In the framework of Generalised Finite Difference Methods, we proposed an optimisation procedure to achieve the best stencil (stencil size optimisation) and coefficients ℎ𝑀and ℎ𝑊(scaling parameters optimisation) that minimise the condition number of Mathematics and Computers in Simulation 218 (2024) 49–78 72 S. Clain and J. Figueiredo Fig. 11. Grids for the test case 3. For 𝛺= {(𝑥, 𝑦) ∈ R2∶√𝑥4+𝑦2∕16 + (3𝑥∕2 + 𝑦2)2+𝑦≤1} the panels illustrate two grids, derived from a Delaunay-like triangulation procedure, with 217 (left) and 847 (right) points. The green filled circles and the red squares correspond, respectively, to interior points and boundary points where Dirichlet conditions are imposed. Fig. 12. Stencil choice for the random distribution case with SSO for 𝑝= 2 (top left panel), 𝑝= 3 (top middle panel), 𝑝= 4 (top right panel), 𝑝= 5 (bottom left panel), 𝑝= 6 (bottom right panel). The black filled circle corresponds to the reference node, the red filled circles stand for the remaining nodes of the final stencil, while the black circles correspond to the nodes discarded during the optimisation procedure. The figures correspond to 𝑆= 1 and the relative position of points does not depend on the scale value. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) matrix 𝑀𝑊 𝑀t. In respect to the stencil, we propose a discrete greedy optimisation procedure by eliminating, one by one, the non-necessary nodes in order to reduce the stencil size while preserving or reducing the condition number. For the continuous optimisation we recurred to two derivative-free optimisation procedures, the MATLAB®function fminsearch and a homemade cloud-like optimisation routine, which delivered very similar results for the values of ℎ𝑀and ℎ𝑊and reduction of the condition number. Since the fminsearch procedure is much faster, we adopted it as the optimisation tool used in the benchmarks presented. We performed an extensive series of numerical tests to assess the condition number reduction and the accuracy of the firstand second-order derivatives discretisation for different cloud geometries. We found that the use of the full optimisation procedure leads to an optimal convergence with respect to the cloud scaling independently of the geometric scenario considered. Mathematics and Computers in Simulation 218 (2024) 49–78 73 S. Clain and J. Figueiredo Fig. 13. Stencil choice for the uniform distribution case with SSO for 𝑝= 2 (top left panel), 𝑝= 3 (top middle panel), 𝑝= 4 (top right panel), 𝑝= 5 (bottom left panel), 𝑝= 6 (bottom right panel). The black filled circle corresponds to the reference node, the red filled circles stand for the remaining nodes of the final stencil, while the black circles correspond to the nodes discarded during the optimisation procedure. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 14. Stencil choice for the half-plane distribution case with SSO for 𝑝= 2 (top left panel), 𝑝= 3 (top middle panel), 𝑝= 4 (top right panel), 𝑝= 5 (bottom left panel), 𝑝= 6 (bottom right panel). The black filled circle corresponds to the reference node, the red filled circles stand for the remaining nodes of the final stencil, while the black circles correspond to the nodes discarded during the optimisation procedure. The figures correspond to 𝑆= 1 and the relative position of points does not depend on the scale value. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Finally, we performed a set of simulations for three test cases involving the convection diffusion reaction equation supplemented by Dirichlet/Neumann boundary conditions for different geometrical settings and grid types. Using the derivatives computed with the optimal reconstruction, proposed in the present work, we obtained effective convergence rates for the numerical approximation