Kriging-based finite element method: element-by-element Kriging interpolation.
An enhancement of the finite element method with Kriging shape functions (K-FEM) was recently proposed. In this method, the field variables of a boundary value problem are approximated using 'element-by-element' piecewise Kriging interpolation (el-KI). For each element, the interpolation function is constructed from a set of nodes within a prescribed domain of influence comprising the element and its several layers of neighbouring elements. This paper presents a numerical study on the accuracy and convergence of the el-KI in function fitting problems. Several examples of functions in two-dimensional space are employed in this study. The results show that very accurate function fittings and excellent convergence can be attained by the el-KI..

Keywords: finite element, Kriging, function fitting.

Article Type:
Finite element method (Research)
Interpolation (Research)
Geology (Statistical methods)
Geology (Research)
Wong, F.T.
Kanok-Nukulchai, W.
Pub Date:
Name: Civil Engineering Dimension Publisher: Civil Engineering Dimension Audience: Academic Format: Magazine/Journal Subject: Engineering and manufacturing industries; Science and technology Copyright: COPYRIGHT 2009 Civil Engineering Dimension ISSN: 1979-570X
Date: March, 2009 Source Volume: 11 Source Issue: 1
Event Code: 310 Science & research
Geographic Scope: United States Geographic Code: 1USA United States
Accession Number:
Full Text:

The finite element method (FEM) is at present very widely used as a numerical method to solve various kinds of problems in engineering and science. The power and versatility of FEM have been tested for several decades of real engineering practices. One important issue in the FEM is mesh generation. Users often prefer to use the simplest elements, namely three-node triangular elements for two dimensional problems and four-node tetrahedral elements for three-dimensional problems, as they can be easily or even automatically generated and are more amenable to adaptive procedure. Nevertheless, it is well-known that these elements often give solutions of poor accuracy, in particular for the gradients of field variables such as stresses or stress resultants in solid and structural mechanics.

Motivated by the desire to eliminate the need for a mesh in numerical analysis, in the last two decades a large variety of mesh-free (or meshless) methods have been proposed as alternatives to the FEM. These methods can be categorized, in view of their formulation bases, into strong-form based mesh-free methods such as a generalized finite difference method [1] and the finite point method [2], and into weak-form based mesh-free methods.

The latter can further be categorized into the methods using a global weak form such as the element-free Galerkin method (EFGM) [3] and point interpolation methods [4, 5] and those using a local weak form such as the meshless local Petrov-Galerkin method [6-9]. The common advantages of this new class of methods are: (1) No mesh is required for construction of the approximate function. (2) High order approximate function can be easily achieved. (3) The solutions are usually more accurate and smoother than the FEM. A detailed review of the methods is presented in [10-12].

Among countless proposed mesh-free methods, to the mesh-free methods having the same basic formulation as the FEM, i.e. those using global weak forms (the EFGM and its variants) were of interest. Even though this class of methods were claimed to be "element free" or "mesh-free", actually elements or background cells are still needed for geometric modeling and numerical integration. Another disadvantage of the EFGM and its variants is that the computational procedures are difficult to incorporate in existing general purpose FEM codes. Due to these inconveniences, these methods up to present do not find wide acceptance in real engineering practices.

In order to eliminate the aforementioned disadvantages, Plengkhom and Kanok-Nukulchai [13] proposed some modifications to the EFGM with moving Kriging interpolation (KI) [14]. The problem domain is subdivided into elements like in the conventional FEM. The KI is constructed for each element using a set of nodes in a domain of influence (DOI) composed of several layers of elements (the DOI is in the form of polygon for 2D problems). Combining the KI of all elements, the global field variable is thus approximated by piecewise KI. For evaluating the integration in the Galerkin weak form, the elements are employed as integration cells. The modified method can be viewed as an enhancement of the FEM using Kriging shape functions. Thus, this method was subsequently referred to as Kriging-based FEM (K-FEM) [15].

With K-FEM, highly accurate field variables and their gradients can be obtained even using the simplest form of elements. Solutions refinements can be achieved without any change to the element or mesh structure. The formulation and coding are very similar to the conventional FEM so that an existing general-purpose FE program can be easily extended to incorporate the enhanced method. Thus, the K-FEM has a higher chance to be accepted in real engineering practices.

In the pioneering work of Plengkhom and Kanok-Nukulchai [13], the K-FEM was applied to solve 1D bar and 2D elasticity problems. Subsequently, it was improved through the use of adaptive correlation parameters and by introducing the quartic spline correlation function and developed for analyses of Reissner-Mindlin plates [15, 16]. A drawback of the present method is that the interpolation function is discontinuous at the inter-element boundaries. In spite of this discontinuity the K-FEM, with appropriate choice of shape functions, passes weak patch tests and therefore the convergence is guaranteed [17]. The basic concepts and advances of the K-FEM have been recently presented [18, 19]. The current development of the K-FEM is the extension and application to different problems in engineering. The authors have recently developed the K-FEM for analyses of general shell structures [20].

As mentioned above, the manner of approximating the field variables in the K-FEM is "element-by-element" Kriging interpolation (el-KI). The previous papers on the K-FEM presented its developments and applications to different problems of solid and structural mechanics. However, direct application of the el-KI for approximation of mathematical functions has not been presented. It is the purpose of this paper to present examinations of the accuracy and convergence of the el-KI in function fitting problems. Several examples of 2D-functions (surfaces) are considered.

Kriging Interpolation

Named after Danie G. Krige, a South African mining engineer, Kriging is a well-known geostatistical technique for spatial data interpolation in geology and mining [21, 22]. Using this interpolation, every unknown value at a point can be interpolated from known values at scattered points in its specified neighbourhood.


Consider a continuous field variable u(x) defined in a domain [OMEGA]. The domain is represented by a set of properly scattered nodes [x.sub.I], I=1, 2, ..., N, where N is the total number of nodes in the whole domain. Given N field values, u([x.sub.1]), ..., u([x.sub.N]), the problem is to obtain an estimate value of u at a point [x.sub.0] [member of] [OMEGA].

Let us consider a set of nodes [x.sub.i], i = 1, 2, ... , n, surrounding point [x.sub.0] inside a sub-domain [[OMEGA].sub.x0] [subset or equal to] [OMEGA]. Here small letter index is used in place of the capital latter to emphasize that the numbering is referred to the sub-domain [[OMEGA].sub.x0]. The Kriging estimated value [u.sup.h](x.sub.0]) is a linear combination of u([x.sub.1]), ..., u([x.sub.n]), i.e.

[u.sup.h]([x.sub.0]) = [[summation].sup.n.sub.i=1] [[lambda].sub.i]u([x.sub.i]) (1)

where [[lambda].sub.i]'s are the (Kriging) weights and n is the number of nodes inside [[OMEGA].sub.x0]. This sub-domain is referred to as DOI in this paper. Considering individual function values, u([x.sub.1]), ..., u([x.sub.n]), as the realization of random variables U([x.sub.1], ..., U([x.sub.n]), Eq. (1) can be written as

[U.sup.h] ([x.sub.0]) = [[summation].sup.n.sub.i=1] [[lambda].sub.i]U([x.sub.i]) (2)

The Kriging weights are determined by requiring that the estimator [U.sup.h]([x.sub.0]) is unbiased, i.e.

E[[U.sup.h]([x.sub.0]) - U ([x.sub.0])] = 0 (3)

and minimizing the variance of estimation error, var [[U.sup.h] ([x.sub.0]) - U ([x.sub.0])]. Using the method of Lagrange for constraint optimization problems, the requirements of minimum variance and unbiased estimator lead to the following Kriging equation system:

R[lambda] = P[mu] = r([x.sub.0]) [P.sup.T][lambda] = p ([x.sub.0]) (4a)

in which


[lambda] = [[[lambda].sub.1] ... [[lambda].sub.n]].sup.T]; [mu] = [[[mu].sub.1] ... [[mu].sub.m]].sup.T] (4c)

r([x.sub.0]) = [[C([h.sub.10]) C([h.sub.20]) ... C([h.sub.n0])]].sup.T]; p([x.sub.0]) = [[p.sub.1] ([x.sub.0]) ... [p.sub.m] ([x.sub.0])]].sup.T] (4d)

Here, R is n x n matrix of covariance between U(x) at nodes [x.sub.1], ..., [x.sub.n]; P is n x m matrix of polynomial values at the nodes; [lambda] is 1 x n vector of Kriging weights; [mu] is 1 x m vector of Lagrange multipliers; r([x.sub.0]) is n x 1 vector of covariance between the nodes and the node under consideration, [x.sub.0]; and p([x.sub.0]) is m x 1 vector of polynomial basis at [x.sub.0]. While in Eqs. (4b) and (4d), C([h.sub.ij]) = cov [U([x.sub.i]), U ([x.sub.j])]. Kriging weights [[lambda].sub.i] can be obtained by solving the Kriging equations, Eq. (4a).

The expression for the estimated value [u.sup.h] given by Eq. (1) can be rewritten in matrix form,

[u.sup.h]([x.sub.0]) = [[lambda].sup.T]d (5)

where d = [[u([x.sub.1]) ... u([x.sub.n])].sup.T] is n x 1 vector of nodal values. Since the point [x.sub.0] is an arbitrary point in the DOI, [x.sub.0] can be replaced by x, i.e. the position of any point in the DOI. Thus, using the usual finite element terminology, Eq. (5) can be expressed as

[u.sup.h] (x) = N (x) d = [[summation].sup.n.sub.i=1] [N.sub.i] (x) [u.sub.i] (6)

in which N(x)= [[lambda].sup.T](x) is the matrix of shape functions.

Kriging shape functions resulting from Eq. (4a) have Kronecker delta (or interpolation) and consistency properties [13, 14]. Due to the former, KI passes through all the nodes thus requiring no special treatment for boundary conditions. The consequence of the latter ensures reproduction of a linear function if the polynomial basis includes the constant and linear term.

These properties make KI very appropriate to be employed as the trial function in a Galerkin method.

Layered-Element Domain of Influence

Consider a 2D domain meshed with triangular elements, such as illustrated in Fig. 1. For each element, KI is constructed based upon a set of nodes in a polygonal DOI encompassing a predetermined number of layers of elements. The KI function over the element is given by Eq. (6). Combining the KI of all elements in the domain, the global field variable is thus approximated by el-KI. This way of approximation is very similar to the approximation in the conventional FEM.

The number of layers for each element must cover a minimum number of nodes in such a way that the Kriging equation, Eq. (4a), can be solved (not singular). If an m-order polynomial basis is employed, the DOI is required to cover a number of nodes, n, that is equal or greater than the number of terms in the basis function [13].

Within each element the el-KI is naturally continuous. However, along the element edges between two adjacent elements the function is not continuous because the KI for the edge of each neighboring element is constructed using different set of nodes. Therefore, K-FEM is a nonconforming method. The issue of non-conformity and its effects on the convergence of the solutions have been addressed by the authors in [17].


Polynomial Basis and Correlation Function

Constructing Kriging shape functions in Eq. (6) requires a polynomial basis function and a model of covariance function. For the basis function, in addition to complete polynomial bases, it is also possible to use incomplete polynomial bases such as bi-linear, bi-quadratic and bi-cubic bases.

Covariance between a pair of random variables U(x) and U(x+h) can be expressed in terms of coefficient of correlation function or shortly, correlation function, i.e. [rho](h) = C(h) / [[sigma].sup.2], where [[sigma].sup.2] = var [U(x)]. According to Gu [14], [[sigma].sup.2] has no effect on the final results and can be taken equals to 1. One of the widely used correlation model in the area of computational mechanics is the Gaussian correlation function [13-15, 17, 18], viz.

[rho](h) = [rho](h) = exp(-[[theta] h / d).sup.2]) (7)

where [theta] > 0 is the correlation parameter, h = [parallel]h[parallel], i.e. the Euclidean distance between points x and x+h, and d is a scale factor to normalize the distance. In this study, d is taken to be the largest distance between any pair of nodes in the DOI. Besides the Gaussian, the authors introduced the quartic spline (QS) correlation function [15] as follows:


The authors found that Kriging shape functions with this correlation function were not very sensitive to the change in parameter [theta]. Moreover, the convergence characteristics of the K-FEM with the QS in many cases were better than the Gaussian function [17].

Correlation Parameter

A proper choice of parameter [theta] is important as it affects the quality of KI. In order to obtain reasonnable results in the K-FEM, Plengkhom and Kanok-Nukulchai [13] suggested a rule of thumb for choosing [theta], i.e. [theta] should be selected so that it satisfies the lower bound,

[absolute value of [[summation].sup.n.sub.i=1] [N.sub.i] - 1] [less than or equal to] 1 x [10.sup.-10+a] (9)

where [N.sub.i], i=1, ...., n are Kriging shape functions and a is the order of basis function, and also satisfies the upper bound,

det(R) [less than or equal to] 1 x [10.sup.-b] (10)

where R is the covariance matrix and b is the dimension of the problem. For 2D problem with cubic basis function, for example, b=2 and a=3.

Numerical investigations on the upper and lower bound values of [theta] [15] revealed that the parameter bounds vary with respect to the number of nodes in the DOI. The authors proposed explicit parameter functions for practical implementation of the K-FEM as follows [15-18]:

For the Gaussian correlation function, the parameter function is

[theta] = (1 - f)[[theta].sup.low] + f [[theta].sup.up], 0 [less than or equal to] f [less than or equal to] 0.8 (11a)

where f is a scale factor, [[theta].sup.low] and [[theta].sup.up] are the lower and upper bound functions as follows:



For the QS correlation function, the parameter function is


With these functions, adaptive values of [theta] can be used now in place of a uniform value of [theta]. Here, "adaptive" means that the correlation parameters used in an analysis are adjusted to the number of nodes in the DOI of each element. The advantage of the use of adaptive [theta] from practical viewpoint is that a user of K-FEM program is not required to input a value of [theta] in an analysis since the parameter functions can be embedded in the program.

Numerical Tests

To study the accuracy and convergence of the el-KI, the following functions are considered:

z = - [x.sup.2] - [y.sup.2] + 1,

[OMEGA] = {(x,y)| [x.sup.2] + [y.sup.2] = 1,0 [less than or equal to] x [less than or equal to] 1, 0 [less than or equal to] y [less than or equal to] 1} (13)

z = [square root of 625 - [x.sup.2]],

[OMEGA] = {(x,y)| 0 [less than or equal to] x [less than or equal to] 25 sin 40[degrees], 0 [less than or equal to] y [less than or equal to] 25} (14)

z = [square root of 100 - [x.sup.2] - [y.sup.2]],

[OMEGA] = {(x,y)| [x.sup.2] + [y.sup.2] = 100, 0 [less than or equal to] x [less than or equal to] 10, 0 [less than or equal to] y [less than or equal to] 10} (15)

Here [OMEGA] is the domain of the function. The first function is a 2D quadratic function, the graph of which is a paraboloid above a quarter of the unit circle on the x-y plane (Fig 2). The second and third functions are adopted from the widely-used benchmark problems for shell structures [23], namely a quarter of the cylindrical surface of Scordelis-Lo roof (Fig. 3) and a quarter of the hemispherical dome (Fig. 4).




To measure the accuracy of the approximation, the following relative L2 norm of error is employed,


in which [[OMEGA].sup.h], denotes a domain for a mesh with the element characteristic size, h and [z.sup.h] denotes the approximation function. The domain [[OMEGA].sup.h] for the mesh with the number of elements on each side M=4 are shown in Fig. 5 for the paraboloid and hemispherical functions and in Fig. 6 for the Scordelis-Lo roof function. The meshes are generated automatically using Delaunay triangulation in Matlab version 6.5 [24]. The integration is performed numerically on each element using the six-point quadrature rule for triangles.



In the following, abbreviations in the form of P*-*-G or P*-*-QS, in which the asterisk represents a number, are used to designate various options of the el-KI. The first syllable denotes the polynomial basis with the order indicated by a number next to letter P. The middle asterisk denotes the number of layers. The last syllable denotes the Gaussian or QS correlation function. The adaptive correlation parameters given by Eq. (11a) for the Gaussian function and by Eq. (12) for the QS function are employed. The scale factor f in Eq. (11a) is taken to be 0.5.

Paraboloid Function

The domain is approximated using meshes with different number of elements on each side, i.e. M=2, 4, 8, 16, and 32. The options P1-2-G, P1-2-QS, P2-3G, and P2-3-QS are considered. The relative L2 norm errors are presented in Table 1 and plotted in Fig. 7. The average of convergence rates are calculated for each option and shown in the legend of Fig. 7.

It can be observed that the accuracy and convergence rate for the el-KI with Gaussian and QS correlation functions are comparable. The Kriging with quadratic basis, as expected, can reproduce the quadratic function exactly (for P2-3-QS) or almost exactly (for P2-3-G). The el-KI with P2-3-G for M=8, 16, and 32 (Table 1) does not yield the relative error of order -16 probably because the Kriging equation system, Eq. 4a, is close to singular for a number of elements.


Scordelis-Lo Roof and Hemispherical Functions

In order to observe the effects of using different number of layers and different basis functions, the functions are firstly approximated using the el-KI of options P1-1-G, P1-2-G, and P1-3-G, and secondly with options P1-3-G, P2-3-G, and P3-3-G. The domains are represented by meshes with M=2, 4, 8, 16, and 32. The relative [L.sub.2] norm errors for the Scordelis-Lo roof function are plotted in Figs. 8a and 8b, while for the hemisphere function the errors are plotted in Figs. 9a and 9b.


The figures show that the accuracy improves quite significantly when the number of layers increases from one to two, but it practically does not improve when the number of layers increases from two to three. For the Scordelis-Lo roof function, the convergence rates are a little bit faster as the number of layers increases. The accuracy and convergence rates improve significantly as the degree of polynomial basis raises. However, this is not the case for the hemispherical function; the accuracy and convergence rates are nearly equal although the degree of polynomial basis rises. The possible reason for this is that at the circular boundary of the hemispherical function, the gradient of the function is singular. Better approximation can be achieved by using finer mesh at the region close to the circular boundary.





The manner of approximation in the K-FEM is element-by-element piecewise Kriging interpolation (el-KI). A series of 2D function fittings using the el-KI have been carried out. The results show that the el-KI can give highly accurate function fittings and very good convergence characteristics. The accuracy is generally higher as the order of polynomial basis raises and as the number of layers increases. Raising the order of polynomial basis may result in significant increase in the convergence rate. For functions containing points of singular gradient, however, the accuracy and convergence rate may not be improved by raising the polynomial order or increasing the number of layers.


[1.] Liszka, T. and Orkisz, J., The Finite Difference Method at Arbitrary Irregular Grids and Its Application in Applied Mechanics, Computers and Structures, Elsevier, 11 (1-2), 1980, pp. 83-95.

[2.] Onate, E., Idelsohn, S., Zienkiewicz, O. C., and Taylor, R. L., A Finite Point Method in Computational Mechanics. Applications to Convective Transport and Fluid Flow, International Journal for Numerical Methods in Engineering, John Wiley and Sons, 39 (22), 1996, pp. 3839-3866.

[3.] Belytschko, T., Lu, Y. Y., and Gu, L., Element-Free Galerkin Methods, International Journal for Numerical Methods in Engineering, John Wiley and Sons, 37 (2), 1994, pp. 229-256.

[4.] Liu, G. R. and Gu, Y. T., A Point Interpolation Method for Two-dimensional Solid, International Journal for Numerical Methods in Engineering, John Wiley and Sons, 50 (4), 2001, pp. 937-951.

[5.] Wang, J. G. and Liu, G.,R., A Point Interpolation Meshless Method Based on Radial Basis Functions, International Journal for Numerical Methods in Engineering, John Wiley and Sons, 54 (11), 2002, pp. 1623-1648.

[6.] Atluri, S. N. and Zhu, T., A New Meshless Local Petrov-Galerkin Method (MLPG) Approach in Computational Mechanics, Computational Mechanics, Springer, 22 (2), 1998, pp. 117-127.

[7.] Atluri, S.N., Kim, H. G., and Cho, J. Y., A Critical Assessment of the Truly Meshless Local Petrov-Galerkin (MLPG) and Local Boundary Integral Equation (LBIE) Methods, Computational Mechanics, Springer, 24 (5), 1999, pp. 348-372.

[8.] Pudjisuryadi, P., Introduction to Meshless Local Petrov-Galerkin Method, Civil Engineering Dimension (Dimensi Teknik Sipil), Petra Christian University, 4 (2), 2002, pp.112-116.

[9.] Pudjisuryadi, P., Adaptive Meshless Local Petrov-Galerkin Method with Variable Domain of Influence in 2D Elastostatic Problems, Civil Engineering Dimension, Petra Christian University, 10 (2), 2008, pp. 99-108.

[10.] Belytschko, T., Krongauz, Y., Organ, D., Fleming, M., and Krysl, P., Meshless Methods: An Overview and Recent Developments, Computer Methods in Applied Mechanics and Engineering, Elsevier, 139 (1-4), 1996, pp. 3-47.

[11.] Fries, T. P. and Matthies, H. G., Classification and Overview of Meshfree Methods, Institute of Scientific Computing, Technical University Braunschweig, Brunswick, Germany, 2004.

[12.] Gu, Y. T., Meshfree Methods and Their Comparisons, International Journal of Computational Methods, World Scientific, 2 (4), 2005, pp. 477-515.

[13.] Plengkhom, K. and Kanok-Nukulchai, W., An Enhancement of Finite Element Method with Moving Kriging Shape Functions, International Journal of Computational Methods, World Scientific, 2 (4), 2005, pp. 451-475.

[14.] Gu, L., Moving Kriging Interpolation and Element-Free Galerkin Method, International Journal for Numerical Methods in Engineering, John Wiley and Sons, 56 (1), 2003, pp. 1-11.

[15.] Wong, F. T. and Kanok-Nukulchai, W., Kriging-Based Finite Element Method for Analyses of Reissner-Mindlin Plates, Proceedings of the Tenth East-Asia Pacific Conference on Structural Engineering and Construction, Emerging Trends: Keynote Lectures and Symposia, 3-5 August 2006, Bangkok, Thailand, Asian Institute of Technology, pp. 509-514.

[16.] Wong, F. T. and Kanok-Nukulchai, W., On Alleviation of Shear Locking in the Kriging-Based Finite Element Method, Proceedings of International Civil Engineering Conference "Towards Sustainable Engineering Practice", 25-26 August 2006, Surabaya, Indonesia, Petra Christian University, pp. 39-47.

[17.] Wong, F. T. and Kanok-Nukulchai, W., On the Convergence of the Kriging-Based Finite Element Method, International Journal of Computational Methods, World Scientific (accepted on 25 October 2008).

[18.] Kanok-Nukulchai, W. and Wong, F. T., A Finite Element Method Using Node-Based Interpolation, Proceedings of the Third Asia-Pacific Congress on Computational Mechanics and the Eleventh International Conference on Enhancement and Promotion of Computational Methods in Engineering and Sciences, 3-6 December 2007, Kyoto, Japan, APACM and EPMESC, Paper No. PL3-3.

[19.] Kanok-Nukulchai, W. and Wong, F. T., A Break-Through Enhancement of FEM Using Node-Based Kriging Interpolation, IACM Expressions, International Association for Computational Mechanics, 23, June 2008, pp. 24-29.

[20.] Wong, F. T. and Kanok-Nukulchai, W., A Kriging-Based Finite Element Method for Analyses of Shell Structures, Proceedings of the Eighth World Congress on Computational Mechanics and the Fifth European Congress on Computational Methods in Applied Sciences and Engineering, 30 June-4 July 2008, Venice, Italy, IACM and ECCOMAS, Paper No. a1247.

[21.] Olea, R. A., Geostatistics for Engineers and Earth Scientists, Kluwer Academic Publishers, Boston, USA, 1999.

[22.] Wackernagel, H., Multivariate Geostatistics, 2nd edition, Springer, Berlin, Germany, 1998.

[23.] MacNeal, R. H. and Harder, R. L., A Proposed Standard Set of Problems to Test Finite Element Accuracy, Finite Elements in Analysis and Design, Elsevier, 1 (1), 1985, pp. 3-20.

[24.] Chapman, S. J., MATLAB Programming for Engineers, 2nd edition, Brooks/Cole, Pacific Grove, California, USA, 2002.

Wong, F. T. (1) and Kanok-Nukulchai, W. (2)

(1) Department of Civil Engineering, Petra Christian University, Surabaya, Indonesia


(2) School of Engineering and Technology, Asian Institute of Technology, Pathumthani, Thailand


Note: Discussion is expected before June, 1st 2009, and will be published in the "Civil Engineering Dimension" volume 11, number 2, September 2009.

Received 13 October 2008; revised 5 December 2008; accepted 14 January 2009.
Table 1. Relative [L.sub.2] norm approximation error for the
paraboloid function

   M        P1-2-G    P1-2-QS     P2-3-G    P2-3-QS

   2       1.92E-02   4.15E-02     N.A.       N.A.
   4       1.00E-02   1.06E-02   2.87E-16   2.26E-16
   8       3.84E-03   3.48E-03   3.20E-15   2.96E-16
   16      1.14E-03   9.60E-04   2.16E-11   3.10E-16
   32      3.09E-04   2.52E-04   2.55E-10   3.06E-16

M: the number of line segments on each side
Gale Copyright:
Copyright 2009 Gale, Cengage Learning. All rights reserved.