Meshless quadrature formulas arising from numerical differentiation
Oleg Davydov Bruno Degli Esposti




Meshless quadrature formulas 

arising from numerical differentiation
    
Oleg DavydovDepartment of Mathematics,
	University of Giessen, Arndtstrasse 2, 35392 Giessen, Germany,
	<oleg.davydov@math.uni-giessen.de>, 
Bruno Degli EspostiDepartment of Mathematics
	and Computer Science, University of Florence,
	Viale Morgagni 67, 50134 Firenze, Italy,
	<bruno.degliesposti@unifi.it>

    September 9, 2024
=======================================================================================================================================================================================================================================================================================================================



§ ABSTRACT 
We suggest a method for simultaneously generating high order
quadrature weights for integrals over 
Lipschitz domains and their boundaries that requires
neither meshing nor moment computation.
The weights are determined on pre-defined scattered nodes as a
minimum norm solution of a sparse underdetermined linear system
arising from a discretization of a suitable boundary value
problem by either collocation or meshless finite differences.
The method is easy to implement independently of the domain's
representation, since it only requires as inputs the
position of all quadrature nodes and the direction
of outward-pointing normals at each node belonging to the boundary.
Numerical experiments demonstrate the
robustness and high accuracy of the method on a number of smooth
and piecewise smooth domains in 2D and
3D, including some with reentrant corners and edges.

 

./figures/



§ INTRODUCTION



Efficient high order numerical integration over
piecewise smooth domains Ω in ^2 and ^3 and on surfaces is
needed in many numerical methods and engineering applications,
see recent extensive literature reviews in
<cit.>. In this paper we suggest
a new technique for generating quadrature formulas on scattered
nodes that only requires as input the normals to the boundary at the
boundary nodes, in addition to the positions of the nodes
themselves, and relies on the discretization of a boundary value problem
for the divergence operator by numerical differentiation 
and the least squares solution of a sparse
underdetermined linear system.

Recall that the classic way of obtaining a weight vector =(_i)_i=1^N∈^N of a quadrature formula

    ∫_Ωf(x)(x)≈∑_i=1^N_if(y_i)

for the integral with respect to some measure
μ in a domain Ω, with given knots
Y={y_1,…,y_N}⊂Ω, is to choose a finite
dimensional linear space ={ s_1,…, s_M}⊂
C(Ω) with good approximation properties, such as
(piecewise) polynomials or radial basis functions, and require
the exactness of the formula for all  f∈, which is
equivalent to the exact reproduction of the moments
∫_Ωs_j,

    ∫_Ωs_j(x)(x)=∑_i=1^N_is_j(y_i),  
    j=1,…,M,

which leads to a linear system for the weight vector ,
as soon as the moments are known. In particular, the book
<cit.> presents many quadrature formulas
for special multivariate regions such as cubes, simplices, balls
or spheres, and special node sets, where exactness is enforced
for spaces of polynomials, and the moments are computed by
analytic methods. With the help of the substitution rule these
formulas can be adapted to the smooth images of such elementary
regions. A more complicated domain may be partitioned into mapped
elementary regions leading to a compound quadrature formula
generated by summing up the formulas over this partition. This
approach is used in particular in the finite element method. For
the high order integration it however requires challenging high
order mesh generation,  often prohibitively costly and non-robust
in applications
<cit.>.
Instead of meshing into images of elementary regions,
different types of partitioning may be used followed by the dimension
reduction that replaces volume or surface integrals by repeated 
univariate integration, such as the standard iterated
integration on generalized rectangles or sectors
<cit.>.
Recent versions of the dimension reduction method
<cit.>
rely on the divergence theorem
over arbitrary subregions with known parametrizations
of their piecewise smooth boundaries.
However, a significant drawback of these 
methods is that many quadrature nodes must be placed
outside of  in its bounding box, unless
the shape of the integration domain satisfies some restrictive
assumptions (see e.g. the concept of base-line
in <cit.>).


When the quadrature nodes in Y are fixed a priori 
(for example, supplied by the user)
and the task is to generate the weight vector ,
the problem setting is related to scattered data fitting,
because the integral
∫_Ωf can be approximated by
∫_Ωs if s is an approximation of f, and
this leads to a quadrature formula as soon as s is obtained for
each fixed Y by a linear scattered data fitting operator 
_Y:^N→, with s=_Y f|_Y, where 
f|_Y:=(f(y_i))_i=1^N. The most established tool
for scattered data fitting is kernel-based interpolation, in
particular with radial basis functions <cit.>. The
meshless
quadrature formulas obtained this way provide optimal recovery 
of the integration functional <cit.>, in
particular polyharmonic and thin plate spline RBF give rise to
optimal quadrature formulas for multivariate Sobolev spaces 
<cit.>. However, numerical implementation of
these methods for general domains is problematic because the
corresponding linear system (<ref>) is not sparse, and the
moments of the radial basis functions have to be provided. The 
first problem may be tackled by either partitioning Ω or
using  spaces  with locally supported bases such as
B-splines, but the computation of the moments remains a major
obstacle unless Ω is a polygonal domain
<cit.>. A dimension reduction approach for the
moments over curved triangles has been suggested in
<cit.>. When using polynomial spaces 
on scattered nodes, a feasible approach is to choose the
polynomial degree such that the  dimension M of  is
significantly smaller than the number of nodes N, which
typically makes the system (<ref>) underdetermined and
solvable. Then a weight vector  satisfying (<ref>) 
may be found as a minimum norm solution 
<cit.> or by
subset selection <cit.>. These methods require the
moments of the polynomials, in particular <cit.>
suggests a sophisticated algorithm for computing them on
bivariate domains bounded by B-spline curves.
















The method that we suggest does not require any meshing, and bypasses
the moment computation problem. 
It relies on the discretization of 
two linear differential operators  and  satisfying 
an “integration by parts” formula of the form

    ∫_Ω u=∫_∂Ω u .

Many identities of this type may be obtained from the generalized Stokes
theorem. 
For example, a natural choice is =Δ, the Laplace
operator (or Laplace-Beltrami operator when Ω
is a Riemannian manifold), and =∂_ν,
the outward normal derivative at the boundary.
This pair is suitable for smooth domains, but unfortunately
not for piecewise smooth ones, as will be demonstrated in the paper. 
Instead, we recommend =, the divergence operator,
and =γ_ν, the normal trace operator, as a general
purpose choice, which works well in all of our numerical tests
on Lipschitz domains in 2D and 3D with either smooth
or piecewise smooth boundary.

Given two functions f:Ω→ and
g:∂Ω→, we look for a combined quadrature
to approximate the difference of two integrals

    ∫_Ωf(x)(x)-∫_∂Ωg(x)(x) 
    ≈∑_i=1^N_Y_if(y_i)- ∑_i=1^N_Z_ig(z_i),

with knots
Y={y_1,…,y_N_Y}⊂Ω and
Z={z_1,…,z_N_Z}⊂∂Ω and weight
vectors =(_i)_i=1^N_Y∈^N_Y and
=(_i)_i=1^N_Z∈^N_Z.
Clearly, (<ref>) provides in particular a quadrature
formula with knots in Y and weight vector 
for ∫_Ωf if we set g≡ 0, and a
quadrature formula with knots in Z and weight vector 
for ∫_∂Ωg if we set f≡ 0.
Assuming that the nodes in Y and Z are provided by the user,
or placed by a node generation algorithm <cit.>,
we propose a method for the computation of the weight
vectors  and . 

As in the classical approach based on the reproduction of moments
(<ref>), we may employ a finite
dimensional space of functions ={ s_1,…, s_M}. 
However, we rely on the simultaneous approximation
of f by s and g by s,
rather than on the direct approximation of f by s∈.
Therefore, we require the exactness of the combined
formula (<ref>) for all
pairs f,g obtained as f= s and g= s from
the same s∈.  According to (<ref>), the left hand
side of (<ref>) is zero in this case, so that the exactness
condition is equivalent to

    ∑_i=1^N_Y_i∑_j=1^Mc_j s_j(y_i)-
    ∑_i=1^N_Z_i∑_j=1^Mc_j s_j(z_i)=0  for all  c_j∈.

In addition, we require the exactness of (<ref>) for a single
pair
of functions f̂,ĝ with known integrals 
∫_Ωf̂ and ∫_∂Ωĝ, 
such that 

    ∫_Ωf̂(x)(x)-∫_∂Ωĝ(x)(x)0.

In particular, we may use f̂≡ 1, ĝ≡ 0 if
we know the measure of Ω, or f̂≡ 0, ĝ≡ 1 if we know the boundary measure of ∂Ω. 
If neither is known, then we may employ f̂≡ 0,
ĝ=∂_νΦ with Φ being a fundamental solution
of the Laplace(-Beltrami) equation centered in Ω,
see Section <ref>.
In any case, we get the following linear system for
the weight vectors  and ,

    ∑_i=1^N_Yℓ_ij_i-∑_i=1^N_Zb_ij_i
       = 0,   j=1,…,M,
    ∑_i=1^N_Y_if̂(y_i)- ∑_i=1^N_Z_iĝ(z_i) 
       = ∫_Ωf̂-∫_∂Ωĝ,

where ℓ_ij and b_ij are the entries of the collocation matrices

    L=( s_j(y_i))_i=1,j=1^N_Y,M,   B=( s_j(y_i))_i=1,j=1^N_Z,M

for the operators
, w.r.t. the sets Y, Z and the space . 
If Ω is a closed manifold, then ∂Ω=∅
and (<ref>) simplifies in that everything related to the
integral over ∂Ω disappears.


Moreover, we may move away from using any function space 
and replace L and B by numerical differentiation matrices on
irregular nodes, as those used in the meshless finite difference
method, see for example <cit.>, which leads to
a purely meshless approach with full flexibility of arbitrarily
distributed degrees of freedom  associated with an additional
finite set X={x_1,…,x_N_X} discretizing either
Ω or some larger set D⊃Ω.
In this case we again find  and  by solving the
linear system (<ref>), with the coefficients
ℓ_ij,b_ij of matrices L,B obtained as the weights of
suitable numerical differentiation  formulas

    u(y_i)≈∑_k=1^N_Xλ_ik^Tu(x_k),
       i=1,…, N_Y,
         u(z_i)≈∑_k=1^N_Xβ_ik^Tu(x_k),
       i=1,…,N_Z,

where

    L=(ℓ_ij)_i=1,j=1^N_Y,M
    =(λ_ik^T)_i=1,k=1^N_Y,N_X,
       B=(b_ij)_i=1,j=1^N_Z,M
    =(β_ik^T)_i=1,k=1^N_Z,N_X.

Note that λ_ij are scalars if u is scalar-valued, but
in general they are vectors, in particular d-vectors when
 is the divergence operator on a d-dimensional manifold
and thus u is a vector field. 

We choose the functional space  or the set of nodes X
such that the system (<ref>) is underdetermined,
and compute the combined quadrature weight vector (,)
as its solution with the smallest 2-norm,

    _2^2+_2^2   →min
    subject to  [     L^T    -B^T; f̂|_Y^T -ĝ|_Z^T ][ ;  ]   = [            0; ∫_Ωf̂ -∫_∂Ωĝ ],

where 

    f̂|_Y=(f̂(y_i))_i=1^N_Y,  ĝ|_Z=(ĝ(z_i))_i=1^N_Z,

which is a linearly constrained least squares problem. Efficient
algorithms for this problem are available and can handle very big
node sets Y,Z if the constraints are given by a sparse matrix
<cit.>. In particular, 
we use a direct
method based on the sparse QR decomposition provided by
the open source library SuiteSparseQR
<cit.>.
The sparsity of the matrices
L and B is achieved by using a locally supported basis of
, such as tensor-product B-splines, or by choosing small sets of
influence in the numerical differentiation formulas (<ref>),
with λ_ik 0 and β_ik 0 only for small
distances between y_i,z_i and x_k.




This leads to a scheme that uniquely combines the following features.



  * High order quadrature weights are provided for 
user-supplied nodes on 
Lipschitz domains,
with the only additional input of outward-pointing normals
at the boundary nodes.
Therefore the method is not tied to any specific representation
of the integration domain, e.g. parametric or implicit.


  * No meshing is needed, which makes this method particularly
well-suited for the meshless workflow in numerical computations
that attracts growing attention as an alternative to mesh-based
methods, because mesh generation remains challenging in more
complicated applications <cit.>.


  * Moment computation, the typical bottleneck of meshless
quadrature methods, is avoided.


  * Simple and easy to implement algorithms rely on a
discretization of a boundary value problem for the divergence
operator and a minimum-norm solution to an underdetermined sparse
linear system, for which  efficient linear algebra routines are
available.



We develop the details of both  the meshless finite difference
approach based on the numerical differentiation matrices
(<ref>), and the functional approach based on the
collocation matrices (<ref>), and formulate two
practical algorithms for domains in ^d that use polyharmonic
kernel numerical differentiation (Algorithm <ref>) and
collocation with tensor-product B-splines (Algorithm
<ref>), respectively.

Extensive numerical experiments demonstrate the effectiveness
of both algorithms, with the largest tests involving more than
100k quadrature nodes, and convergence orders up to h^7 for
both interior and boundary quadrature on several 2D and 3D domains,
including some with reentrant corners and edges.
In addition to SuiteSparseQR, our
MATLAB code relies on the numerical differentiation formulas
provided by mFDlab <cit.>, and on several MATLAB
toolboxes (optimization, curve fitting, statistics).

The paper is organized as follows. Section <ref> presents
the algorithms for the computation of the quadrature weights and
their theoretical background, while
Section <ref>  is devoted to numerical
experiments. We start Section <ref> by an abstract error
analysis that motivated our approach 
(Section <ref>), and
then describe in Section <ref> a general framework
for the development of quadrature formulas of the new type,
followed by the details of the options available for each of
its aspect, namely the choice of the functions  and
 satisfying (<ref>) in 
Section <ref>, differential
operators  and  in 
Section <ref>, methods for
numerical differentiation by either collocation or meshless finite
differences in  Section <ref>, and
the methods for solving the underdetermined linear system 
(<ref>) in 
Section <ref>. Two
practical, ready to implement algorithms are presented in
Section <ref>. In the numerical testing part of the
paper, Section <ref> is devoted to
the validation of the settings suggested in
Algorithms <ref> and <ref>, and 
Section <ref> to convergence order tests.



./figures/



§ QUADRATURE FORMULAS


Let  be a C^0 manifold of dimension d,
and Ω⊂ a domain (that is, a connected
and open subset) with compact closure ,
and let μ,σ
be finite, strictly positive Borel measures defined
on  and ∂Ω, respectively.
We suppose that μ(∂Ω) = 0, so that
the space L^1(Ω) of absolutely integrable functions 
can be identified with L^1().

For any given set Y of N_Y scattered nodes in ,
and, in the case of nonempty boundary ∂Ω∅, 
any given set Z of N_Z scattered nodes in 
∂Ω, we seek quadrature weights
∈^N_Y and ∈^N_Z such that

    ∫_Ω f(x) ≈∑_i=1^N_Y_i f(y_i)
      and  ∫_∂Ω g(x) ≈∑_i=1^N_Z_i g(z_i)

for all sufficiently regular integrands 
f → and g ∂Ω→.
To make the pointwise evaluation f(y_i) and g(z_i)
a meaningful operation we assume that f∈ L^1_Y(Ω) and
g∈ L^1_Z(∂Ω), where L^1_Y(Ω) and
L^1_Z(∂Ω) are the subsets of
L^1(Ω) and L^1(∂Ω) whose elements
(equivalence classes of functions) have at least one
representative that is continuous at each point of Y and Z,
respectively.
We will denote the two quadrature
formulas (Y, ) and (Z, ), respectively.

We will compute the quadrature formulas (Y, ) and (Z, )
simultaneously. Moreover, they arise as a combined
quadrature

    Q(f,g) ∑_i=1^N_Y_i f(y_i)-∑_i=1^N_Z_i g(z_i)

for the difference[Note that we could use the sum of the
two integrals instead of their difference, and respectively the
sum of the two quadrature formulas  in the definition of the
combined quadrature, which would look more natural in a sense.
However we preferred the difference in order to match the
usual form of the compatibility condition
(<ref>) and to emphasize the
"cancellation" nature of the equation I(f,g)=0 that plays a
crucial role in the method.] 
of the two integrals

    I(f,g) ∫_Ω f(x) -∫_∂Ω g(x).


Let → L^1(Ω) and
→ L^1(∂Ω)
be operators such that

    ∫_Ω u 
    = ∫_∂Ω u   for all u ∈,

with  being a common domain of definition for the two.
In particular, we have in mind the
situation where  is a linear space of functions on
, both  and   are linear differential
operators, measures μ,σ are induced by a Riemannian
metric in , and identity
(<ref>) holds by the divergence
theorem or another corollary of the generalized Stokes theorem.
When Ω= is a closed manifold with ∂Ω=∅, 
the measure σ and the operator  are
clearly irrelevant, and the right-hand
side in  (<ref>)
is understood to be identically zero.

The underlying idea of the method may be roughly explained as follows.
The identity (<ref>) implies that 

    I(f,g)=I( u, u)=0
 
as soon as f= u and g= u for the same function
u∈. If u is discretized by a vector c∈^M,
and matrices L ∈^N_Y × M
and B ∈^N_Z × M are such that

    Lc≈ u Y=fY,   
    Bc≈ u Z=gZ,

then 

    Q(f,g)
    ≈^T L c - ^T B c
    = (^T L - ^T B)   c,

and hence requiring that the weight vectors , satisfy

    ^TL-^TB=0

implies 

    Q(f,g)≈ 0=I(f,g)  for all pairs (f,g) such that  
    f= u and  g= u  for some  u ∈.

In addition, we require from , that 

    Q(,)=I(,)  for a single pair (,) 
    with I(,)0.

As we will demonstrate, this combination of
conditions produces accurate quadrature
formulas (Y, ) and (Z, ) for appropriate choices of
the operators , and numerical
differentiation matrices L,B.

We first make these considerations more precise in 
Section <ref>, then
formulate in <ref> a general framework 
for obtaining quadrature formulas of this type, 
and explore in Sections <ref>–<ref>
the options available for each of its steps.
The final subsection <ref> presents two
specific algorithms that we recommend for practical
use on piecewise smooth Lipschitz domains in ^d,
and which are extensively tested in Section <ref>.



 §.§ Error analysis in an abstract setting




For any given operators  and  and finite sets 
Y⊂ and Z⊂∂Ω, we say that a
triple (Ψ,L,B), where Ψ→^M,
L ∈^N_Y × M and B ∈^N_Z × M,
is a numerical differentiation scheme
with discretization operator Ψ and
differentiation matrices L,B.
We denote by _L(u) and _B(u) the recovery error
with which the matrices L and B approximate the
pointwise values of u and u
using the information about u∈ contained in Ψ(u):

    ε_L(u)   := u Y - LΨ(u)_∞=
    max_i=1,…,N_Y u (y_i) - (LΨ(u))_i,
    ε_B(u)   := u Z - BΨ(u)_∞=
    max_i=1,…,N_Z u (z_i) - ( BΨ(u))_i.

The vector Ψ(u) ∈^M 
may correspond to pointwise values of u over a scattered
set of nodes, coefficients of splines or radial basis functions 
that approximate u,
Fourier coefficients, or other kinds of discrete information about u.
Tools of approximation theory may then be used to define
matrices L and B that accurately reconstruct
the pointwise values of u and u
from a given vector Ψ(u).


Our approach is motivated by the following proposition that
estimates the combined quadrature error

    δ(f,g)  I(f,g)-Q(f,g).




Let (,) ∈ L^1_Y(Ω) × L^1_Z(∂Ω)
be a fixed pair of functions such that

    I(,)
    ≠ 0.

Given any pair (f,g) ∈ L^1_Y(Ω) × L^1_Z(∂Ω), 
let _f,g⊂ be the set of solutions u
to the boundary value problem

    u = f - α   in Ω
     u = g - α   on ∂Ω,

where the coefficient α∈ is chosen so that 
(<ref>) satisfies
the compatibility condition (<ref>):

    α=I(f,g)/I(,).

Then for any pair of quadrature formulas (<ref>), 
the following estimate of the combined error holds:


    δ(f,g)≤αε̂
    + inf_u ∈_f,g{_1 ε_L(u)
    	+ _1 ε_B(u)
    	+ Ψ(u)_∞L^T -B^T_1
    },

where :=|δ(,)|.



For all u ∈_f,g, by the linearity of δ,

    δ(f,g)
    = δ(f - α, g - α)
    	+ δ(α,α)
    = δ( u,  u) + αδ(,)≤δ( u,  u) + αε̂.

Moreover, by the compatibility condition
(<ref>) and by the definitions
of ε_L(u) and ε_B(u),

    δ( u,  u)   = ∑_i=1^N_Y_i ( u)(y_i)
    	- ∑_i=1^N_Z_i ( u)(z_i)
       = ∑_i=1^N_Y_i ( u)(y_i)
    	±∑_i=1^N_Y_i ( LΨ(u) )_i
    	±∑_i=1^N_Z_i ( BΨ(u) )_i
    	- ∑_i=1^N_Z_i ( u)(z_i)
       ≤_1 ε_L(u)
    	+  ^T L Ψ(u) -  ^T B Ψ(u)
    	+ q_1 ε_B(u) 
       ≤_1 ε_L(u)
    	+ _1 ε_B(u)
    	+ Ψ(u)_∞L^T -B^T_1.

Taking the infimum over u ∈_f,g completes the proof.


By choosing either f≡0 or g≡0, we 
obtain separate error bounds for both quadrature formulas
(<ref>).




By the previous proposition, the special cases

    δ(f,0)
    = ∫_Ω f  - ∑_i=1^N_Y_i f(y_i),
        δ(0,g)
    = -∫_∂Ω g  + ∑_i=1^N_Z_i g(z_i),

immediately imply the following error bounds for
quadrature formulas (Y, ) and (Z, ),

    ∫_Ω f  - ∑_i=1^N_Y_i f(y_i)   ≤αε̂
    + inf_u ∈_f{_1 ε_L(u)
    	+ _1 ε_B(u)
    	+ Ψ(u)_∞L^T -B^T_1
    }, 
    ∫_∂Ω g  - ∑_i=1^N_Z_i g(z_i)   ≤βε̂
    + inf_u ∈_g{_1 ε_L(u)
    	+ _1 ε_B(u)
    	+ Ψ(u)_∞L^T -B^T_1
    },

with _f and _g being the sets of solutions to the
boundary value problems

    u = f - α   in Ω
     u =   - α   on ∂Ω,
        and     u =   - β   in Ω
     u = g - β   on ∂Ω,

respectively, where  α=I(f,0)/I(,) and
β=I(0,g)/I(,)
are again determined
by the compatibility condition
(<ref>).


The case of a closed manifold simplifies to just one quadrature
formula.


If ∂Ω = ∅, then there is only one
quadrature formula (Y, ), and the argumentation of
Proposition <ref>
leads to the error estimate

    ∫_Ω f  - ∑_i=1^N_Y_i f(y_i)≤αε̂
    + inf_u ∈_f{_1 ε_L(u)
    	+ Ψ(u)_∞L^T_1
    },

where _f is the set of solutions to the
equation u = f - α in Ω, and
α=∫_Ω f /∫_Ω, with
∈ L_Y^1(Ω) such that ∫_Ω0.




 §.§ A general framework


Motivated by Proposition <ref> we
intend to generate quadrature weights , by trying to
minimize the error |δ(f,g)| of the combined quadrature.
It is clear from (<ref>) that 
|δ(f,g)| is small if 

    α, _L(u), _B(u) and L^T -B^T_1 are small

and 

    _1, _1  and Ψ(u)_∞ are bounded.

Moreover, if we choose the weight vectors  and 
such that L^T -B^T =0, then the size of
Ψ(u)_∞ is irrelevant. We may
also require that =|δ(,)|=0,
which is just another linear equation for  and 
to satisfy. Assuming that _L(u) and _B(u)
are small thanks to (a) the smoothness of u enforced by
the appropriate choice of the operators ,, and (b)
the choice of the numerical differentiation scheme, we may use any
remaining degrees of freedom in   and 
in order to minimize the stability constants
_1 and _1. 
We may minimize the joint 1-norm _1 + _1,
or some other norm ( , )_♯
that may be preferable for computational or other reasons.
All these ideas lead us to the following algorithm.

 
Given Ω, , , Y, Z
as in Section <ref>, 
compute quadrature weights ^* ∈^N_Y
and ^* ∈^N_Z as follows:




  * Choose auxiliary functions ∈ L^1_Y(Ω)
and ∈ L^1_Z(∂Ω) such that

    ∫_Ω
    - ∫_∂Ω≠ 0.



  * Choose operators  and  such that

    ∫_Ω u 
    = ∫_∂Ω u   for all u ∈,

where the set  contains at least one solution 
of the boundary value problem
(<ref>) for each pair of functions f ∈ L^1_Y(Ω)
and g ∈ L^1_Z(∂Ω)
whose integrals we wish to approximate using
formulas (Y, ^*) and (Z, ^*).


  * Choose a numerical differentiation scheme (Ψ,L,B)
for the linear reconstruction of the pointwise values
u Y and u Z
from the vector Ψ(u) in the form of matrix-vector
products LΨ(u) and BΨ(u).


  * Assemble the non-homogeneous linear system
with unknowns ( , )

    L^T   - B^T   = 0,     ∑_i=1^N_Y_i (y_i)
    - ∑_i=1^N_Z_i (z_i)
    = ∫_Ω
    - ∫_∂Ω,

which can be written more compactly as 

    A x = b,   with  A = [    L^T   -B^T;  |_Y^T -|_Z^T ],   
    x = [ ;  ],   
    b = [          0; ∫_Ω - ∫_∂Ω ].



  * Compute a solution ( ^*, ^*) of
(<ref>) that minimizes
some norm ( , )_♯^N_Y×^N_Z→^+ defined on the space of
all possible weights ( , ).



The algorithm terminates successfully as long as the linear system
(<ref>) is consistent,
and Corollary <ref> implies that the
quadrature formulas (Y, ^*) and (Z, ^*)
satisfy the error bounds

    ∫_Ω f  - ∑_i=1^N_Y_i^* f(y_i)   ≤inf_u ∈_f{ ^*_1 ε_L(u)
    	+  ^*_1 ε_B(u)
    }, 
    ∫_∂Ω g  - ∑_i=1^N_Z_i^* g(z_i)   ≤inf_u ∈_g{ ^*_1 ε_L(u)
    	+  ^*_1 ε_B(u)
    }.


The consistency of (<ref>) may be
characterized in terms of a certain
discrete incompatibility condition as follows.
Recall that  and  satisfy I(,) ≠ 0
and hence in view of (<ref>) 
are incompatible as the right hand sides of the boundary value
problem for the operators  and , that is, the problem

    u = λ f̂, 
     u = λ ĝ,

is not solvable for u∈ whenever
λ∈∖{0}.


The linear system Ax=b given by (<ref>) is
consistent if and only if the
following discrete incompatibility condition holds: every solution
(c,λ) ∈^M × to the linear system

    L c = λ f̂|_Y, 
    
    B c = λ ĝ|_Z,

satisfies λ = 0.



The characterization essentially follows
from the identity (A)^⊥ = (A^T)
and the usual assumptions on (,).
We can reason as follows:

    b ∈(A)
         ⇔  [              0; ∫_Ωf̂
    	- ∫_∂Ωĝ ]∈(A^T)^⊥
         ⇔  [              0; ∫_Ωf̂
    	- ∫_∂Ωĝ ]·[ c; λ ]
    	= 0   for all
    		(c,λ) ∈(A^T)
         ⇔  λ	= 0   for all
    		(c,λ) ∈(A^T)
         ⇔  λ = 0   for all
    		(c,λ) that satisfy
    	(<ref>).



Since (<ref>) is a
discretization of (<ref>), we expect that the
discrete incompatibility condition holds whenever the
discretization obtained via the choice of Y,Z and the numerical 
differentiation scheme is sufficiently accurate for the linear
system to inherit this important feature of the continuous
problem.



Even though we cannot guarantee the consistency of 
(<ref>) in general, in particular for
arbitrary user-supplied quadrature nodes Y,Z, a simple
strategy of choosing the size M of the discretization
vector Ψ(u) such that the linear system
(<ref>)
is underdetermined with significantly fewer rows than columns
makes the system solvable in all of our numerical experiments,
leaving enough room for enforcing reasonable size of
the stability constants _1 and _1
through minimization, see Section <ref>.





In the following subsections we illustrate and compare
several options for the choices to be made at each step of
Algorithm <ref>.



 §.§ Auxiliary functions f̂
and ĝ



The non-homogeneous constraint

    ∑_i=1^N_Y_i (y_i)
    - ∑_i=1^N_Z_i (z_i)
    = ∫_Ω
    - ∫_∂Ω≠ 0

plays a fundamental role in the definition of linear system
(<ref>),
because it rules out the trivial solution ( ^*, ^*) = (0,0).
In practice, choices for  and  need not
be complicated: 

    (,) ≡ (1,0)  or  (,) ≡ (0,1)

are usually effective,
as long as at least one of the moments

    Ω∫_Ω 1 ,     ∂Ω∫_∂Ω 1

is either known in advance, or can be numerically approximated in
a straightforward way. In particular, this can be easily done
using standard quadrature formulas when the boundary
∂Ω is defined by explicit parametric patches over
simple regions in the parameter domain. 

As an alternative, we suggest the following approach that only
requires knowing sufficiently many points on the boundary
∂Ω to generate the set Z of quadrature nodes,
as well as the outward-pointing normals
on ∂Ω at Z.
This applies for instance to implicit surfaces as used in the
level-set method, and to trimmed multipatch surfaces as used in 
Computer Aided Design. 

Assume that Ω is a domain in ^d and
μ and σ are the standard
Lebesgue and hypersurface measures.
Let Φ(x,x_0) be the fundamental solution of the
Laplace equation centered at a point x_0 ∈Ω:

    Φ(x,x_0)
    = 1/2πlog(x-x_0)
      if d = 2,    Φ(x,x_0)
    = -x-x_0^2-d/d(d-2)ω_d  if d = 1 or d > 2,

with ω_d being the measure of
the d-dimensional unit ball. Let ν be
the outward-pointing unit normal field along ∂Ω,
and let ∂_ν be the normal derivative operator
on ∂Ω. By Green's third identity,

    ∫_∂Ω∂_νΦ(x,x_0)  = 1,

and so the choice 

    (x)=0  and  (x)= ∂_νΦ(x,x_0)=ν(x)^T(x-x_0)/dω_dx-x_0^d

satisfies

    ∫_Ω
    - ∫_∂Ω = -1

independently of the shape of the domain and the position
of x_0 ∈Ω.
Observe that the singularity of Φ(x,x_0) in
Ω plays no role in the solution of boundary
value problems (<ref>),
because ≡ 0.
In practice, care must be taken to avoid
points x_0 too close to the boundary,
for otherwise the smoothness of the solution u of
(<ref>) may be affected by the
near-singularity of , and as a consequence
the quantities ε_L(u), ε_B(u)
in the quadrature errors estimates may blow up.






 §.§ Operators  and 



The error bounds in
Section <ref>
depend on the accurate reconstruction of pointwise
values of u and u from a given vector
Ψ(u), with u being a solution to a
boundary value problem of the form

    u = f    in Ω, 
     u = g    on ∂Ω,
        with  ∫_Ω f  = ∫_∂Ω g .

No matter which numerical differentiation method is used
to define matrices L and B, the accuracy of the
reconstruction, and thus the quantities _L(u) and
_B(u) in (<ref>), depend crucially
on the regularity of u as measured by e.g. Sobolev norms.
On the contrary, the accuracy of the quadrature formulas should
rather depend on the regularity of the integrands f and g,
as is natural in numerical integration.
Therefore, we are led to consider operators
 and  for which the boundary value problem
(<ref>) admits solutions of sufficiently high
regularity, depending on the smoothness of f and g.



  §.§.§ Neumann problem for Laplace-Beltrami operator



The  classical results in the theory of elliptic partial
differential equations, see e.g. <cit.>,
suggest the use of the Laplace-Beltrami
operator 
and the normal derivative trace operator.

Let Ω be a domain in a C^∞ Riemannian
manifold , with compact closure  and (possibly empty)
smooth boundary ∂Ω, and let μ and
σ be the measures on Ω and ∂Ω
induced by the Riemannian metric of . Denote by
H^s(Ω) and H^s(∂Ω) the L^2-norm Sobolev
spaces on Ω and ∂Ω of (possibly fractional)
order s ≥ 0, and by ν be the outward-pointing unit
normal field along ∂Ω. Then it follows from the 
divergence theorem that 

    ∫_ΩΔ u  = ∫_∂Ω∂_ν u   for all u ∈ H^2(Ω),

where Δ is the Laplace-Beltrami operator
on the manifold, and ∂_ν is the normal derivative
trace operator. In other words, identity
(<ref>) is satisfied for

    =Δ, 
       = ∂_ν,

with = H^2(Ω). 



Under the above assumptions, for all f ∈ H^s(Ω) and 
g ∈ H^s+1/2(∂Ω), the Neumann boundary value problem

    Δ u = f    in Ω
    ∂_ν u = g    on ∂Ω  with constraint  ∫_Ω u  = 0

has a unique solution u ∈ H^s+2(Ω)
that depends continuously on f and g, if and only if

    ∫_Ω f  = ∫_∂Ω g ,

a constraint known as compatibility condition
for the Neumann boundary value problem.


As a direct consequence of this result, the choice
(<ref>)
is fully satisfactory for a bounded domain Ω in a smooth
Riemannian manifold  with (possibly empty) smooth boundary
∂Ω, because solutions to the boundary value
problem (<ref>)
are guaranteed to exist, and to be more regular
than the corresponding data in the right-hand side.

However, if ∂Ω is not smooth, as is
often the case in practical applications, 
elliptic regularity results as in Proposition
<ref> do not apply anymore, and so
the solution u to the Neumann boundary value problem
(<ref>) may not be regular enough
for the numerical differentiation scheme

    Δ u Y≈ LΨ(u) ,   ∂_ν u Z≈ BΨ(u).

to be sufficiently accurate even for highly smooth f and g,
which indicates that the choice
(<ref>) may be problematic in this
case because the numerical
differentiation errors ε_L(u) and ε_B(u)
featuring in the estimate (<ref>)
are presumably affected by the loss of smoothness of u.

Indeed, let Ω be a bounded Lipschitz domain in ^d,
and μ,σ the standard Lebesgue and
hypersurface measures. Then Δ is
the usual Laplace operator. The Lipschitz regularity of Ω
is sufficient to define an outward-pointing
unit normal field ν almost everywhere on ∂Ω,
and (<ref>) still holds.
Moreover, for all f ∈ L^2(Ω) and
all g ∈ L^2(∂Ω), the Neumann boundary problem
(<ref>) has a unique
solution u ∈ H^3/2(Ω) if and only if
the compatibility condition (<ref>)
is satisfied, as implied by the results in <cit.>.
However, standard examples in <cit.>
for the corresponding Dirichlet problem on piecewise
smooth domains Ω⊂^2 with reentrant corners can be
modified to apply to the Neumann problem, and generate for any
s>3/2 a solution u ∉ H^s(Ω), even when f is
infinitely differentiable and g=0.
This means that the regularity of u is not only limited
by the regularity of f and g, but also by the regularity
of ∂Ω. 




The following example, adapted from
<cit.>, shows that this may happen
even on elementary geometries, such as a square.

 
Let Ω = (0,1)^2, and let Γ_1,…,Γ_4 be the
left, bottom, right, and top sides of ∂Ω,
respectively. For all i = 1,…,4 and s ≥ 0,
let γ_i be the trace operator from H^s+1/2(Ω)
to H^s(Γ_i).
Functions in H^s(Γ_i) are parametrized
by the restrictions of Cartesian coordinates (x,y)
to ∂Ω, that is, x ∈ [0,1] for i = 2,4
and y ∈ [0,1] for i = 1,3. The following function,
defined using polar coordinates (r,θ)
centered at the origin, is harmonic in Ω:

    v(r,θ) = r^2 (
    log(r) cos(2θ) - θsin(2θ)
    ).

Let φ [0,1] → be a smooth function
such that φ(r) ≡ 1 for all r ∈ [0,1/3],
and φ(r) ≡ 0 for all r ∈ [2/3,1].
Then, the function u(r,θ) φ(r)   v(r,θ)
is the unique solution (up to a constant) of

    Δ u = Δ (φ v)
    	   in Ω
    ∂_ν u = 0
    	   on ∂Ω∖Γ_1
    ∂_ν u = γ_1(-∂_x u) = π y φ(y)
    	   on Γ_1.

The right hand side functions f,g of this Neumann problem 
are smooth in the sense
that they are restrictions to Ω and ∂Ω
of suitable functions that belong to C^∞(^2).
Nevertheless, it is clear that u ∉ H^3(Ω),
because u and v coincide for r < 1/3 and

    ∂^3 v(x,y)/∂ x^3
    = 2x/x^2+y^2
    = 2cos(θ)/r∉ L^2(Ω).



In principle, the estimate (<ref>) that
relates the error of our quadrature to the numerical
differentiation errors ε_L(u) and ε_B(u),
is just an upper bound that is not guaranteed to be tight.
Moreover, on Lipschitz domains in ^2 with piecewise
smooth boundary, it can be proved that
u still has regularity H^s+2 away from corners,
so one could hope that a few inaccurate approximations
of Δ u(y_i) and ∂_ν u(z_i)
around corners would not affect the overall
accuracy of the quadrature formulas. However, our numerical
experiments in Section <ref>
demonstrate that a catastrophic loss of accuracy
in the quadrature formulas may occur on a Lipschitz domain
with piecewise smooth boundary, whereas we do not register 
this behavior for  the alternative divergence approach suggested
below.



  §.§.§ The divergence approach




The most basic instance of the identity
(<ref>) is clearly the divergence
theorem, 
which suggests the choice 

    = ,
       = γ_ν,

where  is the divergence operator on the manifold
and γ_ν is the normal trace operator.
This leads us to the boundary value problem

    F  = f    in Ω, 
    γ_ν(F) = g    on ∂Ω    with  ∫_Ω f  = ∫_∂Ω g ,

where the measures μ and σ are determined by the
Riemannian metric, the compatibility condition on the right must
be satisfied due to the divergence theorem, and the solution F 
is sought in an  appropriate space  of vector fields on
Ω. In the case of the manifold =^d, vector fields
may be identified with vector valued functions
F:Ω→^d, and the normal trace is just the
projection γ_ν(F)=ν^TF of F to the
outward-pointing unit normal ν on the boundary ∂Ω
for sufficiently smooth F.

Any solution u of the Neumann problem (<ref>)
gives rise to a solution F=∇ u of
(<ref>). Hence, no smooth solutions are ever
lost by using operators  and γ_ν
instead of Δ and ∂_ν. However,
unlike (<ref>), the boundary value problem
(<ref>) is underdetermined,
with highly non-unique solutions F.
Recall that the error bound (<ref>)
involves the infimum over all solutions of
(<ref>), so the non-uniqueness is an
advantage as (<ref>) may rely on the
solution F with the smallest numerical differentiation errors
_L(F) and _B(F).
In particular, one may hope that on Lipschitz domains
with piecewise smooth boundary there still exists a solution F
of (<ref>) that inherits a high order of
smoothness from smooth functions f,g in the right hand side,
even when the same is not true about (<ref>).

Although we have not found results of this type in the
literature, considerable attention has been paid to the corresponding
Dirichlet problem, where the full trace of F on the boundary is
prescribed instead of the normal trace ν^TF in
(<ref>). The Dirichlet problem, however, requires 
additional compatibility conditions between f  and the boundary
data when Ω is not sufficiently smooth.  In particular, if
F vanishes on the boundary of a polygon Ω⊂^2,
then its divergence f must vanish at the vertices of Ω,
which is a quite unnatural condition to impose on integrands f.
Nevertheless, several results on the Dirichlet problem for domains in
=^d imply that a solution
F of (<ref>) belongs to H^s+1(Ω)^d
when f∈ H^s(Ω) and g satisfies
certain smoothness and compatibility assumptions. In particular, this follows from
<cit.> for Lipschitz domains in
^d, with s∈, when g=0 and f vanishes on the
boundary   ∂Ω to order s-1. Furthermore,
<cit.> gives the same regularity result for
polygons in ^2 when g=0 and f vanishes at all vertices of
the polygon, and also for non-homogeneous g, with smoothness
and compatibility conditions that are difficult to recast to our
setting, in which g is the normal projection of the Dirichlet
boundary data studied in <cit.>. These results
are generalized in <cit.> to 
arbitrary bounded Lipschitz domains in smooth
manifolds, see in particular <cit.>
for the case =^d.


Our numerical results in Section <ref>
for the divergence-based quadrature do not indicate any
loss of convergence order  on non-smooth domains, which
suggests that the smoothness of F is not lost.
Moreover, the following elementary construction produces a smooth
solution F of (<ref>) on Ω=(0,1)^2,
whenever f is smooth in Ω and g is smooth on each side
of the square, in contrast to
Example <ref> for the
elliptic Neumann problem.





  §.§.§ Smooth solution of divergence boundary
problem on the square
 

As before, let Γ_1,…,Γ_4
denote the left, bottom, right, and top sides of ∂Ω,
and let g_i, i=1,…,4, be the univariate
functions representing g|_Γ_i in the natural Cartesian
parametrizations. 
We set g̃_i  g_i-α_i,
with α_i=∫_0^1g_i(t) dt, and
consider the vector fields

    G_1(x,y)=
    [     (x-1)g̃_1(y); -∫_0^yg̃_1(t) dt ],  
    G_2(x,y)=
    [ -∫_0^xg̃_2(t) dt;     (y-1)g̃_2(x) ],


    G_3(x,y)=
    [         xg̃_3(y); -∫_0^yg̃_3(t) dt ],  
    G_4(x,y)=
    [ -∫_0^xg̃_4(t) dt;         yg̃_4(x) ]

that satisfy

    G_i=0,    ν^TG_i|_∂Ω=
    g|_Γ_i-α_i    in Γ_i,
    
    0    otherwise.

By construction, it follows that

    G̃(x,y)  G_1(x,y)+G_2(x,y)+G_3(x,y)+G_4(x,y)+
    [ α_1(x-1)+α_3x; α_2(y-1)+α_4y ]

is a solution of the problem

    G̃=α_1+α_2+α_3+α_4,  ν^TG̃|_∂Ω=g.

Moreover, let 
α=∫_Ω f(x,y) dx dy,
 h(y)=∫_0^1f(ξ,y) dξ-α and f̃(x,y)=f(x,y)-α-h(y). 
Then the vector field

    F̃(x,y)=[ ∫_0^xf̃(ξ,y) dξ;    ∫_0^yh(ζ) dζ ]

satisfies

    F̃=f-α,    ν^TF̃|_∂Ω=0.

Since α=α_1+α_2+α_3+α_4 by the
compatibility condition, we conclude that F=F̃+G̃
solves (<ref>). It is easy to see that for any
r≥0 we have F∈ C^r(Ω)^2 as soon as f∈
C^r(Ω) and g_i∈ C^r[0,1], i=1,…,4,
which even allows for g to be discontinuous at the corners of
∂Ω.



 §.§ Numerical differentiation schemes


The choice of the numerical differentiation scheme (Ψ,L,B)
is crucial for the performance of the method, because
the factors _L(u) and  _B(u) in the error bounds
(<ref>) and (<ref>), as well as the solvability of
the linear system (<ref>) with reasonable stability
constants _1 and _1
depend on this choice. In particular, the discrete 
incompatibility condition (<ref>)
must be satisfied. Moreover, for an efficient numerical
implementation it is important that the differentiation matrices 

    L = ( ℓ_ij)_i=1,j=1^N_Y,M  and  
    B = ( b_ij)_i=1,j=1^N_Z,M

are sparse.

We describe two approaches that work well in our numerical tests.
The first is based on meshless numerical differentiation formulas as used
in the generalized finite difference methods, and the second
employs tensor-product splines. 

In both cases we assume that  and  are 
linear differential operators of
integer orders on a space   of functions on a domain D 
such that Ω⊂ D. The domains Ω and D
need not have the same dimension; indeed, when 
is embedded into an ambient space  such as
a surface  in =^3, one may
prefer to choose a domain D in  rather than in ,
leading to an immersed approach. 

For the sake of notational simplicity,
but without loss of generality, we assume in this section that
the functions in  are scalar-valued.



  §.§.§ Meshless finite difference formulas


Assuming that ⊂ C(D), we choose a finite set 
X = {x_1,…,x_M}⊂ D, and define the
discretization operator Ψ→^M
by pointwise evaluation over X,

    Ψ(u) = ( u(x_1), …, u(x_M) )^T,

and approximate u (y_i) and u (z_i) at the
quadrature nodes y_i ∈ Y and z_i ∈ Z
by appropriate numerical differentiation formulas

    u (y_i) ≈∑_j=1^M ℓ_i j  u(x_j),      u (z_i) ≈∑_j=1^M b_i j  u(x_j).

Since differentiation is a local operation, sufficiently accurate
formulas can be found with the sums in
(<ref>)
restricted to small sets of indices
S_L,i and S_B,i corresponding to the nodes
in X close to y_i and z_i, known as
sets of influence. This in turn implies that
ℓ_ij=0 for j∉ S_L,i, b_ij=0 for 
j∉ S_B,i, and so the matrices L and B are sparse.

When X is a gridded set and Y∪ Z⊂ X, formulas of
this type can be obtained by  classical finite differences as
those used in the finite difference method for partial
differential equations. For irregular sets X,Y,Z, as typically
needed on more complicated domains Ω, numerical
differentiation formulas are in the core of the  meshless
finite difference methods, such as RBF-FD or GFDM, see for
example <cit.>  
and references therein.
Therefore, we generally refer to
(<ref>)
as meshless finite difference formulas.

A basic approach to obtaining accurate numerical differentiation
weights ℓ_ij, j∈ S_L,i, and b_ij,  j∈
S_B,i, is to impose exactness of formulas 
(<ref>) over a local basis of
functions that can provide a good approximation of
u∈ in a neighborhood of y_i and z_i, respectively. 
In the case when Ω is a domain in ^d, a natural local
approximation tool is given by the spaces Π^d_q of multivariate
polynomials of total degree <q. The sets of influence
are usually chosen somewhat larger than the minimum needed to admit a 
numerical differentiation formula exact for all u∈Π^d_q,
and the extra degrees of freedom are used either to minimize a
(semi-)norm of the weight vectors (ℓ_ij)_j∈ S_L,i
and (b_ij)_j∈ S_B,i, or to enhance the local
approximation space by adding radial basis functions. We refer to
<cit.> for the theory and
error bounds for these methods. Polynomial type methods are also
available on certain manifolds, such as the d-dimensional
sphere, where spherical harmonics may be employed. Otherwise,
numerical differentiation formulas on arbitrary manifolds may be
generated with the help of positive definite kernels. Suitable local error
bounds for the kernel-based numerical differentiation may be found 
in <cit.> for domains in ^d and in
<cit.> for reproducing kernels of Sobolev
spaces on manifolds.


We discuss in more detail the numerical differentiation formulas
employed in our numerical tests. They are generated by the
polyharmonic radial basis kernel K(x,y)=x-y^2m-1,
x,y∈^d, with a polynomial term in Π^d_q, q≥ m, which 
is a standard approach in RBF-FD since <cit.>. We
use m=q. Thus, the weights ℓ_ij, j∈ S_L,i, are
obtained by requiring the exactness  of the first formula in
(<ref>) for all functions
of the form

    u(x)=∑_j∈ S_L,ic_jx-x_j^2q-1+p̃(x),  
    c_j∈, p̃∈Π^d_q,  with ∑_j∈ S_L,ic_jp(x_j)=0 for all p∈Π^d_q,

and similarly for the weights b_ij, j∈ S_B,i. In most
cases the square linear system  that arises from
these conditions is regular. However, as shown in <cit.>,
numerical differentiation weights are uniquely determined and can
be computed by a null space method also for certain `deficient'
sets of influence, such as the five point stencil for the
Laplacian on gridded nodes.


Let Ω be a domain in ^d, and Ω⊂ D⊂^d. 
Denote by k_L and k_B
the orders of the operators  and , respectively.
Suppose we use polyharmonic numerical differentiation formulas
with q=q_L>k_L for the operator , and q=q_B>k_B for .
Under appropriate assumptions on the sets of influence
X_L,i{x_j:j∈ S_L,i}, and
X_B,i{x_j:j∈ S_B,i}, such
as quasi-uniformity and boundedness of the
polynomial Lebesgue constants, the errors (<ref>) and
(<ref>) of the polyharmonic numerical differentiation
 can be estimated as 

    _L(u)= u Y - LΨ(u)_∞
    =max_1≤ i≤ N_Y| u (y_i) -
    ∑_j∈ S_L,iℓ_i j  u(x_j)|≤ 
    C_Lh_L^q_L-k_L|u|_W^q_L_∞(D),
    _B(u)= u Z - BΨ(u)_∞
    =max_1≤ i≤ N_Z| u (z_i) -
    ∑_j∈ S_B,i b_i j  u(x_j)|≤ 
    C_Bh_B^q_B-k_B|u|_W^q_B_∞(D),

where h_L is the maximum diameter of the sets {y_i}∪
X_L,i, i=1,…,N_Y, h_B is the maximum diameter of the
sets {z_i}∪ X_B,i, i=1,…,N_Z, and C_L,C_B
are some positive constants independent of h_L,h_B and u.
Note that a simple method for the selection of the sets of
influence that ensures in practice that the errors behave in
accordance with these estimates is to compose X_L,i of
2Π^d_q_L closest nodes of y_i in X, and X_B,i
of 2Π^d_q_B closest nodes of z_i, as recommended in
<cit.> for RBF-FD.
Further methods that try to optimize the selection can be
found in <cit.> and references therein.
In any case, the number of nodes in X_L,i and X_B,i
remains bounded if q_L and q_B are fixed, which in turn
implies that the matrices L and B are sparse and that
the diameters h_L and h_B shrink as the density
of the nodes increases.

Assuming that k_B=k_L-1, as in both the elliptic and the divergence
settings of Section <ref>,
we get the same approximation order for  _L(u) and
_B(u) by choosing q_B=q_L-1. 
Hence, by (<ref>) and (<ref>), the quadrature
formulas computed by Algorithm <ref>
satisfy the estimates 

    ∫_Ω f  - ∑_i=1^N_Y_i^* f(y_i)   ≤  C ( ^*_1+  ^*_1 )h^q-kinf_u ∈_fmax{|u|_W^q-1_∞(D),|u|_W^q_∞(D)}, 
    ∫_∂Ω g  - ∑_i=1^N_Z_i^* g(z_i)   ≤  C ( ^*_1+  ^*_1 )h^q-kinf_u ∈_gmax{|u|_W^q-1_∞(D),|u|_W^q_∞(D)},

where 

    C=max{C_L,C_B},   h=max{h_L,h_B},   q=q_L=q_B+1,  
    k=k_L=k_B+1.



These estimates show that, as soon as the stability constants
^*_1 and ^*_1 remain bounded and the sets
_f and _g contain functions u∈ W^q_∞(D),
the error of both quadrature formulas behaves as
(h^q-k).
In order to relate the convergence order to the smoothness of the
integrands f and g, we need regularity results for the
boundary value problems (<ref>),
such as those discussed in
Section <ref>. 

For example, in the setting of a smooth domain
Ω⊂^d and operators =Δ,
=∂_ν (with k=2), assuming that 
and  are infinitely differentiable, we may use
Proposition <ref> to infer that 
there exists u_f∈ U_f that belongs to H^s+2(Ω)
as soon as f∈ H^s(Ω).
By the Sobolev embedding theorem,  H^s+2(Ω)⊂
W^q_∞(Ω) if s+2>q+d/2, and by the Stein extension
theorem u_f may be extended to a function in W^q_∞(D).
Therefore, the estimate 

    ∫_Ω f  - ∑_i=1^N_Y_i^* f(y_i)
    =(h^q-2)

holds for all functions f∈ H^s(Ω) with s>q-2+d/2
as long as ^*_1+ ^*_1 is bounded.
Similarly, we derive from
Proposition <ref>  that  there exists
u_g∈ U_g that belongs to H^s+3/2(Ω) as soon as g∈
H^s(∂Ω). By the same arguments with the Sobolev
and Stein theorems, u_g may be extended to a  function in
W^q_∞(D) if s+3/2>q+d/2, and we obtain the estimate

    ∫_∂Ω g  - ∑_i=1^N_Z_i^* g(z_i)
    =(h^q-2)

for all functions g∈ H^s(∂Ω) with s>q-2+(d+1)/2
as long as ^*_1+ ^*_1 is bounded.





  §.§.§ Choosing discretization nodes

In practical applications, quadrature formulas often
need to be defined on fixed, user-supplied nodes,
so in this work we never assume to be in control of the placement
of quadrature nodes Y and Z. 
The set X⊂ D of discretization nodes that determines the
operator Ψ, however, can be chosen freely, and so it must be
computed by Algorithm <ref>.

First of all, we have to choose D, which may coincide with
, but it may also be a larger domain in  or in
the ambient space , for example a bounding box around
. While exploring these options numerically, we have
seen that using D= and irregular X was consistently
a better choice than a bounding box D with a Cartesian grid
X. 

Even though a discretization set X fitting  may be
obtained via mesh generation, for example by taking X
as the vertices of a triangulation of Ω,
we only need meshless nodes for our purposes.
The nodes need not be connected into grids or networks,
which simplifies node generation algorithms,
especially for complicated 3D domains or surfaces.
For a survey on the topic of meshless
node generation, see <cit.>.

The generation of X should take into account multiple
criteria that in part contradict one another.
On the one hand, for any fixed quadrature nodes Y and Z,
we want the set X to be as large as possible, in order
to improve the errors _L(u) and _B(u)
of the numerical differentiation thanks to
a higher density of X in .
On the other hand, a smaller set X produces
more stable quadrature formulas. Indeed, let us consider
two nested sets of nodes X ⊂X̃, 
and corresponding matrices L,B,L̃,B̃ so that
L,B are submatrices of L̃,B̃. Then,
the solution (^*,^*) to
(<ref>) with L̃,B̃
that minimizes the 1-norm _1 + _1
satisfies not just

    L̃^T ^* - B̃^T ^* = 0,
      but also  
    L^T ^* - B^T ^* = 0,

and so ^*_1 + ^*_1
≤^*_1 + ^*_1.
Moreover, as soon as X is large enough so that M+1>N_Y +
N_Z, the linear system (<ref>)
is overdetermined, i.e. it has more equations than unknowns,
and so it may not admit solutions anymore.
In general, X should not be unnecessarily irregular,
although a higher density of the nodes may be
advantageous near the boundary of Ω, especially in the
vicinity of its corners or fine features. In the numerical tests we 
only consider quasi-uniform nodes.
For more details on the distribution of discretization nodes X
and how the size N_X should be chosen, see
Section <ref>.









  §.§.§ Numerical differentiation by tensor-product splines



Let  be an M-dimensional linear subspace of 
spanned by a basis {s_1,…,s_M}⊂.
Any map 

    Ψ→^M,   Ψ(u) = ( c_1(u), …, c_M(u) )^T,

can be seen as a discretization operator that gives rise to 
the numerical differentiation  scheme 

    u (y_i)   ≈(
    	∑_j=1^M c_j(u) s_j ) (y_i)
    	=∑_j=1^M c_j(u) s_j(y_i)
    	=(LΨ(u))_i,   i=1,…,N_Y
     u (z_i)   ≈(
    	∑_j=1^M c_j(u) s_j ) (z_i)
    	=∑_j=1^M c_j(u) s_j(z_i)
    	=(BΨ(u))_i,   i=1,…,N_Z,

given by the matrices L ∈^N_Y × M and
B ∈^N_Z × M with entries

    ℓ_i j =  s_j(y_i)
      and  
    b_ij  =  s_j(z_i).


Since these matrices do not depend on Ψ, the weights
^*, ^* computed by Algorithm <ref> do
not depend on it either, and hence the choice of Ψ is 
irrelevant for the implementation of our quadrature in this
setting. However, numerical differentiation errors

    _L(u)=max_i=1,…,N_Y u (y_i) 
    - (∑_j=1^M c_j(u) s_j )(y_i),
      _B(u)=max_i=1,…,N_Z u (z_i) 
    - (∑_j=1^M c_j(u) s_j )(z_i)
  
of (<ref>)–(<ref>) do depend on Ψ.
Therefore we may take infimum of  (<ref>) and
(<ref>) over all possible Ψ, which is equivalent to
infimum  over all s=∑_j=1^M c_j(u) s_j∈,
and implies the estimates

    ∫_Ω f  - ∑_i=1^N_Y_i^* f(y_i)   ≤inf_u ∈_finf_s∈{ ^*_1  (u-s)|_Y_∞
    	+  ^*_1  (u-s)|_Z_∞}, 
    ∫_∂Ω g  - ∑_i=1^N_Z_i^* g(z_i)   ≤inf_u ∈_ginf_s∈{ ^*_1   (u-s)|_Y_∞
    	+  ^*_1 (u-s)|_Z_∞}.


Note that the weights ^* and ^* do not depend on the 
choice of the basis  {s_1,…,s_M} of a given space
. Indeed, if we choose another basis and let V
be the change-of-basis matrix, then the new differentiation
matrices are given by L̃ = L V
and B̃ = B V, hence the condition
L̃^T -B̃^T
= V^T (L^T-B^T )=0
is equivalent to L^T -B^T =0 because V is invertible.
Nevertheless, the matrices L and B still depend on the
choice of the basis, and their properties strongly influence the
efficiency and stability of the computation of the weights. In
particular, in order to obtain sparse  matrices L and B, we need
locally supported basis functions s_j, such that only a few of
them do not vanish in the vicinity of each point y_i or z_i.

Although the space  can be chosen freely,
and in principle it would be interesting to compare
different approaches, in this work we only consider the case
of the unfitted tensor-product spline approximation.


Suppose that the manifold  of dimension d is embedded in
^n, with n≥ d,  and let H be the n-dimensional
bounding box around  such that the lengths of its
sides are integer multiples of a parameter h_ > 0:

    ⊂ H = [a_1, b_1] ×…× [a_n, b_n],
         b_i - a_i = N_i   h_  for all i = 1,…,n.

For each dimension, let T_i be the uniform knot vector

    T_i = { a_i, …, a_i,
    a_i + h_, …, b_i - h_,
    b_i, …, b_i },

where the first and last knots are repeated q times,
and let _i be the univariate spline space
of order q (degree q-1) defined by the knot vector T_i.
We define  to be the restriction of the
tensor-product spline space
_1 ⊗…⊗_n
to  , and a natural basis
for  is given by the restriction to 
of those tensor-product B-splines s_1,…,s_M,
whose supports have non-empty intersection with Ω.
In the case when d<n this is known as the ambient,
or immersed approach <cit.>.
Clearly, an upper bound for M = is given by

    ∏_i=1^n _i
    = ∏_i=1^n ( N_i + q - 1 ),

although M will often only be a 
fraction of this upper bound, for example when
n>d, or when n=d but Ω≪H.
We denote by D the interior of the union of the supports of all
B-splines in . Then D may also be considered
as the domain of definition of the splines in .

B-splines s_j have local supports contained in cubes  with side
length qh_, and so the  matrices L and B are sparse as
each of their rows has at most q^n nonzeros.  Note that the
assembly of  L and B does not require any expensive tests
whether the  intersection of s_j with Ω is non-empty,
because  we can simply discard any basis function such that
s_j|_Y∪ Z = 0, see Algorithm <ref>.
 
Assume that n=d, such that  Ω⊂ D are domains in ^d.
As before, we denote by k_L and k_B
the orders of the operators  and , respectively,
with max{k_L,k_B}< q.
Then it follows by <cit.> that
for any u∈ W^q_∞(D) there exists s∈
such that simultaneously

    (u-s)|_Y_∞≤ C_Lh_^q-k_L|u|_W^q_∞(D)  and   (u-s)|_Z_∞≤
    C_Bh_^q-k_B|u|_W^q_∞(D),

where C_L,C_B are some positive constants independent of h_L,h_B and u.
Therefore, we obtain from (<ref>) and (<ref>)

    ∫_Ω f  - ∑_i=1^N_Y_i^* f(y_i)   ≤inf_u ∈_f{
    	C_L ^*_1 h_^q-k_L
    	+ C_B ^*_1h_^q-k_B}|u|_W^q_∞(D), 
    ∫_∂Ω g  - ∑_i=1^N_Z_i^* g(z_i)   ≤inf_u ∈_g{
    	C_L ^*_1 h_^q-k_L
    	+ C_B ^*_1h_^q-k_B}|u|_W^q_∞(D),

as long as _f and _g contain functions extensible to D
with a finite seminorm |u|_W^q_∞(D).

Similar to Section <ref>, for a smooth domain
Ω⊂^d and operators =Δ,
=∂_ν with k_L=2, k_B=1, assuming that 
and  are infinitely differentiable, we may use
Proposition <ref>, the Sobolev
embedding theorem, and the Stein extension
theorem, to obtain the same (h^q-2) estimates (<ref>)
and (<ref>), as long as the combined stability constant
^*_1+ ^*_1 is bounded. 

We refer to <cit.> for approximation 
results by ambient tensor-product splines when Ω is a
compact closed hypersurface. 




 §.§ On the choice of minimization norm



As we have seen in Section <ref>, 
assuming stability of the quadrature formulas, that is,
_1 = (1) and _1 = (1),
the regularity of the boundary value problem
(<ref>), and appropriate numerical
differentiation methods, Algorithm
<ref> computes quadrature formulas of
high convergence order. This seems to suggest to use the
combined 1-norm (,)_♯
= (,)_1 = _1+_1
in Step 5 of Algorithm <ref>.

However, the main reason for the 1-norm to appear in the
estimates is that in the error analysis it was convenient
to use the 1-norms _1, _1
of the quadrature weight vectors , ,
and the ∞-norm of the error of numerical differentiation.
In fact, the Hölder inequality with any pair of conjugated 
exponents p,p' ∈ [1,+∞] such that 1/p + 1/p' = 1,
leads to the estimates 

    ∫_Ω f 
    	- ∑_i=1^N_Y_i f(y_i)   ≤inf_u ∈_f{_p  u Y - LΨ(u)_p'
    	+ _p  u Z - BΨ(u)_p'}, 
    ∫_∂Ω g 
    	- ∑_i=1^N_Z_i g(z_i)   ≤inf_u ∈_g{_p  u Y - LΨ(u)_p'
    	+ _p  u Z - BΨ(u)_p'},

where p=1 of (<ref>) and (<ref>)
is not necessarily optimal.


Nevertheless, there are still good reasons why one may prefer
the minimization of the 1-norm of the quadrature weights.
First of all, the 1-norm is the standard way
to assess stability of a quadrature formula: any
perturbation of size ε > 0 in the values
f(y_i) and g(z_i) is potentially amplified by
_1 and _1 when computing the sums
∑_i=1^N_Y_i f(y_i) and 
∑_i=1^N_Z_i g(z_i).

Second, 1-norm minimization is useful for obtaining
positive formulas with non-negative weights
_i≥0 or _i≥0, which is a highly
desirable property in many applications, although
far from necessary for achieving stability of numerical
integration in general. If we know the measures

    Ω∫_Ω 1 ,     ∂Ω∫_∂Ω 1 ,

of the domain and its boundary,
then we can choose ≡ 1 and ≡-1
as auxiliary functions, so that

    ∑_i=1^N_Y_i +∑_i=1^N_Z_i = Ω 
    + ∂Ω

holds for all quadrature formulas (Y,) and (Z,)
satisfying (<ref>). Hence,

    Ω + ∂Ω≤_1 + _1,

and equality is attained if and only if the weights
 and  are non-negative. If a pair of non-negative
weight vectors (,) satisfying system
(<ref>) exists, then it will be found
as (^*,^*) by 1-norm minimization, and the combined
1-norm (^*,^*)_1 will be equal to
Ω + ∂Ω,
implying excellent stability.
Note that in this case we may ensure that both
formulas (Y,) and (Z,) are exact for constants by
requiring two conditions 

    ∑_i=1^N_Y_i = Ω  and  ∑_i=1^N_Z_i = ∂Ω

instead of (<ref>).  
It is however possible that positive formulas with
(<ref>) exist but none of them satisfies (<ref>).
If we enforce instead of (<ref>) just one of the two conditions in
(<ref>), that is, we use one of the pairs of auxiliary
functions (,)≡(1,0) or
(,)≡(0,1), then there is no guarantee that
1-norm minimization delivers a positive formula for either
(Y,) or (Z,), even if they exist.





From a computational point of view, any 1-norm minimization
problem can be reformulated as a linear program
with twice as many unknowns, so any general purpose
linear programming solver can be used for 1-norm
minimization. When the simplex algorithm is used,
sparse interior quadrature formulas 
with no more than m nonzero weights can be obtained,
with m being the total number of rows in the
linear system (<ref>).
For a numerical comparison of 1-norm
and 2-norm minimization, we refer to the experiments
in Section <ref>.


The choice p = 2 is very attractive from the computational
point of view: efficient linear algebra routines
based on e.g. QR decomposition are available
to find a solution of system
(<ref>)
with minimal combined 2-norm

    (,)_2
    √(_2^2 + _2^2).

Moreover, the strict convexity
of the 2-norm guarantees uniqueness of the optimal
solution (^*,^*), a property that may not
hold for the 1-norm. 




 §.§ Practical choices for effective algorithms



After discussing in
Sections <ref>–<ref> 
various options available for the realization of the framework
outlined in Algorithm <ref>, we now
describe two particular settings 
tested in the numerical experiments of
Section <ref> and recommended for practical
applications. 

We only consider 
domains Ω in ^d, even if the framework is applicable
to surfaces and other manifolds.
The evidence gathered in our experiments suggests
that these choices are essentially optimal among
the ones that will be compared in
Section <ref>, and therefore
represent a good starting point for applications
and future research.
The first algorithm is based on meshless finite difference
formulas, whereas the second one is based on tensor-product
spline spaces. This way we demonstrate that both functional and
meshless finite difference approaches have  effective
realizations. In either case, we assume that |∂Ω|
is known, and choose

    (,) ≡ (0,1),
         = ,
         = γ_ν,
        (,)_♯ = (,)_2
    = √(_2^2 + _2^2).

For simplicity, we assume that the quadrature nodes in Y and Z, 
although irregular in general, are not intentionally generated
with density varying over the domain or its boundary,
and we characterize their density by a single
spacing parameter h > 0 defined as the step size 
of a uniform Cartesian grid with N_Y nodes in a
d-dimensional cube of measure |Ω|, and similarly for Z:

    h ≈ h_Y ≈ h_Z,   with  
    h_Y ( Ω/N_Y)^1/d,   
    h_Z ( ∂Ω/N_Z)^1/(d-1).

We assume that, in addition to supplying the quadrature nodes
in Y and Z, the user provides sufficiently accurate
outward-pointing unit normals at the boundary
nodes of Z ⊂∂Ω,
and, in the case of Algorithm <ref>, is in position to
generate an additional node set X in  targeting
a prescribed spacing parameter h_X > 0.


 
Approach based on meshless finite difference formulas (MFD).



Note that in order to increase the performance of Algorithm <ref>,
one should assemble the sparse matrices L^T and B^T directly,
instead of precomputing matrices L_k and D_k for each
spatial dimension separately.
Moreover, the index sets S_L,i and S_B,i can be
determined on the fly during the assembly procedure.
Nevertheless, we have chosen to present Algorithm <ref>
in the form above to enhance readability, and to clearly
show how L^T and B^T can be put together in the divergence
case, which has not been detailed out in
Section <ref>.
Indeed, we have assumed in that section that
functions in  are scalar valued, but in the case
of = and = γ_ν
the functions in  are actually vector valued.

In the scalar case, the discretization operator
Ψ→^N_X is defined by pointwise
evaluation over X, see (<ref>), and so the elements
in the vector Ψ(u) follow the order of the nodes in X.
In the vector valued case, the discretization operator
Ψ(F) ∈^dN_X is defined by pointwise
evaluation over X of all d components of the vector field
F ∈, and so one must choose whether to order
the elements of Ψ(F) so that all d components for a
fixed node are contiguous, as in

    Ψ(F) = (
    F_1(x_1), …, F_d(x_1),
      …,   
    F_1(x_N_X), …, F_d(x_N_X) )^T,

or so that all N_X pointwise values for a fixed component
are contiguous, as in

    Ψ(F) = (
    F_1(x_1), …, F_1(x_N_X),
      …,   
    F_d(x_1), …, F_d(x_N_X) )^T.

In Algorithm <ref> we have chosen the latter ordering,
and although the discretization operator Ψ is not
used in the algorithm, the chosen ordering for Ψ(F) clearly
determines the order of the columns of L and B,
and hence their assembly.
In a high performance implementation, the order
that provides the fastest memory access should be preferred.
In any case, the final values of  and 
do not depend on the chosen ordering for Ψ(F),
because changing the order amounts to permuting the rows
of the linear system Ax=b.

 
Approach based on a tensor-product spline space (BSP).



Once again, we have written Algorithm <ref>
in a way to enhance readability rather than performance.
We conclude this section with three remarks.

First, we note that the set X of Algorithm <ref> may
also be generated as a Cartesian grid over the bounding box H
as defined in Algorithm <ref>. This version of MFD
delivered satisfactory results in our preliminary numerical 
experiments, but we decided to present it in the current form
because the theoretical analysis in Section <ref> in
this case avoids the need for the existence of the extensions of
the functions u_f and u_g  to a larger domain
D⊃Ω without loss of smoothness, which may lead to a
better performance in some settings. Moreover, being able to 
generate the node set X⊂ targeting a prescribed 
spacing parameter h_X is not a restrictive hypothesis in practice, 
because, whenever node generation is deemed expensive due to the
complexity of , one may obtain X by thinning 
the set Y ∪ Z. 

Second, Algorithms <ref> and <ref> can be 
adapted to the setting where the quadrature nodes Y and Z 
are generated with locally varying density according to a
spacing function h(x) →^+. In this
case, the set X of Algorithm <ref> may be generated by
using a scaled spacing function h_X(x) = c  h(x) for some
c>1, either from scratch, or by a suitable
subsampling algorithm, see e.g. <cit.>. 
In the case of Algorithm <ref>, the spacing function
h(x) may be used  to guide local refinement of the
tensor-product spline space.

Third, the choice of the coefficients in h_X = 1.6   h and 
h_ = 4   h is justified by the numerical experiments of
Section <ref>.
For now, we point out that these choices lead
by design to underdetermined systems Ax=b.
To see why, let us consider the case of Algorithm <ref>.
Assuming that N_Y ≫ N_Z, the matrix A
has approximately N_Y columns and dN_X rows.
By the definitions of h and h_X, 

    d N_X/N_Y≈d Ω h_X^-d/Ω h^-d
    = d (1.6)^-d,

and so the ratio of rows to columns is always smaller
than 1 for all d ∈, with a peak of about 0.8
for d=2. A similar argument shows that the linear system
Ax=b assembled in Algorithm <ref> is
also underdetermined for small enough h.






./figures/


	search path = ./tables,
	col sep = comma,
	string replace = NaN,
	empty cells with = -,
	set thousands separator = ,
	every column/.append style = 
		zerofill, sci e, precision = 2



	columns/Nsamples/.style=
		column name=N_seeds, int detect
	,
	columns/hX/.style=
		column name=h_X, sci
	,
	columns/hY/.style=
		column name=h, sci
	,
	columns/hZ/.style=
		column name=h_Z, sci
	,
	columns/NX/.style=
		column name=N_X, int detect
	,
	columns/NY/.style=
		column name=N_Y, int detect
	,
	columns/NZ/.style=
		column name=N_Z, int detect
	,
	columns/NrowsA/.style=
		column name=m, int detect
	,
	columns/NrowsA_mfd/.style=
		column name=m_MFD, int detect
	,
	columns/NrowsA_bsp/.style=
		column name=m_BSP, int detect
	,
	columns/sizeratioA/.style=
		column name=m/N, fixed
	,
	columns/sizeratioA_mfd/.style=
		column name=m_MFD/N, fixed
	,
	columns/sizeratioA_bsp/.style=
		column name=m_BSP/N, fixed
	,
	columns/nnzA/.style=
		column name=nnz(A), int detect
	,
	columns/residual/.style=
		column name=Ax-b, sci
	,
	columns/tnodegen/.style=
		column name=t_nodes, sci
	,
	columns/tmatassembly/.style=
		column name=t_mat, sci
	,
	columns/tsolve/.style=
		column name=t_solve, sci
	,
	columns/sreintrungeavg/.style=
		column name=e_avg(f_1), sci
	,
	columns/sreintrungestd/.style=
		column name=e_std(f_1), sci
	,
	columns/sreintrungerms/.style=
		column name=e_rms(f_1), sci
	,
	columns/sreintrungemax/.style=
		column name=e_max(f_1), sci
	,
	columns/srebndrungeavg/.style=
		column name=e_avg(g_1), sci
	,
	columns/srebndrungestd/.style=
		column name=e_std(g_1), sci
	,
	columns/srebndrungerms/.style=
		column name=e_rms(g_1), sci
	,
	columns/srebndrungemax/.style=
		column name=e_max(g_1), sci
	,
	columns/sreintfrankeavg/.style=
		column name=e_avg(f_2), sci
	,
	columns/sreintfrankestd/.style=
		column name=e_std(f_2), sci
	,
	columns/sreintfrankerms/.style=
		column name=e_rms(f_2), sci
	,
	columns/sreintfrankemax/.style=
		column name=e_max(f_2), sci
	,
	columns/srebndfrankeavg/.style=
		column name=e_avg(g_2), sci
	,
	columns/srebndfrankestd/.style=
		column name=e_std(g_2), sci
	,
	columns/srebndfrankerms/.style=
		column name=e_rms(g_2), sci
	,	
	columns/srebndfrankemax/.style=
		column name=e_max(g_2), sci
	,	
	columns/sreintexpsumavg/.style=
		column name=e_avg(f_3), sci
	,
	columns/sreintexpsumstd/.style=
		column name=e_std(f_3), sci
	,
	columns/sreintexpsumrms/.style=
		column name=e_rms(f_3), sci
	,
	columns/sreintexpsummax/.style=
		column name=e_max(f_3), sci
	,
	columns/srebndexpsumavg/.style=
		column name=e_avg(g_3), sci
	,
	columns/srebndexpsumstd/.style=
		column name=e_std(g_3), sci
	,
	columns/srebndexpsumrms/.style=
		column name=e_rms(g_3), sci
	,
	columns/srebndexpsummax/.style=
		column name=e_max(g_3), sci
	,
	columns/stabw/.style=
		column name=K_, fixed, dec sep align =c, precision = 2
	,
	columns/stabq/.style=
		column name=K_, fixed
	,
	columns/stabwmax/.style=
		column name=max K_, fixed,
		dec sep align =c, precision = 2
	,
	columns/stabqmax/.style=
		column name=max K_, fixed
	,
	columns/minw/.style=
		column name=min, sci
	,
	columns/minq/.style=
		column name=min, sci
	,
	columns/nnzw/.style=
		column name=(), int detect
	,
	columns/nnzq/.style=
		column name=(), int detect
	,
	columns/nnzwq/.style=
		column name=()+(), int detect
	




§ NUMERICAL TESTS



In this section we numerically evaluate the accuracy and 
stability of quadrature formulas (Y,) and (Z,)
produced by Algorithms <ref> and <ref>.
Let Ω be a bounded Lipschitz domain in ^2 or ^3
with piecewise smooth boundary.
Quadrature errors of (Y,μ) and (Z,) are computed
for two test functions f_1, f_2 ∈ C^∞(Ω), and,
respectively, their restrictions 
g_1 = f_1∂Ω and 
g_2 = f_2∂Ω.
The first test function f_1 is a multidimensional
generalization of Runge's function centered at
a domain-dependent point x_R ∈^d:

    f_1(x,x_R) = 1/1+25x-x_R_2^2.

The second test function f_2 as scaled and translated 
Franke's function <cit.> in 2D, and its 3D
generalization <cit.>. As our domains do
not fit into the square [0,1]^2 or cube [0,1]^3, for which
these test functions have been designed, we compose them with the
affine mapping 

    (x_1,…,x_d) ↦( x_1+1/2, …, x_d+1/2),  
    d=2,3,

between [0,1]^d and [-1,1]^d.
In two dimensions, Franke's function is evaluated
with the  command in MATLAB.


The accuracy of the quadrature formulas is assessed
by computing the relative errors

    e(f_1) δ(f_1,0) / I(f_1,0) ,
     … ,  
    e(g_2) δ(0,g_2) / I(0,g_2).

In all tests, the denominators are large enough
to make relative errors meaningful.
Relative errors were preferred to absolute errors to make
results comparable across different test functions and domains.
As will be clear 
from the convergence
plots of Section <ref> that go below 10^-15 in some
instances, we need very accurate reference values  of the integrals

    I(f_1,0),    I(f_2,0),    I(0,g_1),    I(0,g_2),

in order to compute at least one reliable significant digit
of the true relative errors. For every test domain,
parametrizations of all smooth pieces of ∂Ω
are known in closed form, with the parametric domain given by
either an interval for 2D domains, or a rectangle for 3D domains.
The integrals of f_1 and f_2 over Ω were
computed by finding vector fields F_1 and F_2
whose divergence is f_1 and f_2, then using
the divergence theorem to turn the integrals over Ω
into integrals over ∂Ω, and finally integrating
over the parametric domain of each smooth boundary
piece with MATLAB's adaptive quadrature routines
 and .
To meet the required accuracy target, absolute and
relative tolerances were set to 8 · 10^-16.

Note that even though the test functions f_1 and f_2
are both infinitely differentiable and of a simple shape,
we can distinguish their smoothness because the partial
derivatives of f_1 grow as the factorial of their order,
whereas those of f_2, as en entire function, grow much 
slower as an exponent of the order.
This makes f_1, that actually stems from the famous Runge
example for polynomial interpolation, a more difficult test
function for high order methods than f_2.

The stability of the quadrature formulas is assessed
by computing the normalized stability constants

    K_1/Ω∑_i=1^N_Y_i,
        
    K_1/∂Ω∑_i=1^N_Z_i,

that measure the sensitivity of the quadrature to the maximum
absolute error in the function values, such that formulas with 
smaller K_ and K_ are more stable. Again, the
normalization helps to compare the stability across all test
functions and domains. For any quadrature formula exact on
constants we have K_≥ 1 or K_≥ 1, and the
equality holds if and only if the formula is positive. Without
exactness for constants, but with the relative quadrature error
for the constant function less than some ∈(0,1), we get 
K_> 1- and K_ > 1-. Therefore, the best
stability is attained when K_ and K_ are close
to one. 
Note that Algorithms <ref> and <ref> may produce
quadrature weights such that K_ < 1, because exactness
for constants is only enforced for (Z,) by the choice
(,) ≡ (0,1). If |Ω| is known, in addition
to |∂Ω|, then one may also enforce exactness of
(Y,) for constants by adding one more row to the linear
system (<ref>), but this does not
provide any tangible benefit, as will be shown
in Section <ref>.

In what follows we first demonstrate in
Section <ref> 
the effectiveness of the settings suggested in 
Algorithms <ref> and <ref>,
which serves as a kind of a tuning step for the free
parameters in the algorithms.
We do the testing mostly for Algorithm <ref>,
because the results for Algorithm <ref>
are very similar. A second batch of numerical experiments in
Section <ref> shows that quadrature errors converge to
zero as h → 0 at the rate determined by the order of the
numerical differentiation scheme.


In all tests the weights of the  meshless finite difference formulas are
computed using the open source library mFDlab
<cit.>, and those for the tensor-product spline
differentiation by MATLAB's Curve Fitting Toolbox. 




 §.§ Validation of recommended parameters





  §.§.§ Choice of quadrature nodes



In the first test we demonstrate the robust performance of the
MFD method described in Algorithm <ref> with respect to
the choice of nodes in the sets Y and Z. These experiments are performed on 
an elliptical domain Ω_1 in ^2
with semi-axes of length 1 and 3/4 centered at the origin:

    Ω_1 = {
    (x_1,x_2) ∈^2 |
    x_1^2 + x_2^2/(3/4)^2 < 1
    }.


A domain in ^d with piecewise smooth boundary may be
discretized by scattered nodes using a wide range of techniques,
see the survey <cit.>.
We consider three important cases: a meshless advancing front
point cloud generation algorithm, and two rejection sampling
algorithms based either on the quasi-random Halton sequence
or on pseudo-random uniformly distributed samples,
in both cases mapped to a bounding box containing Ω.

The advancing front method works by first generating
a set of boundary nodes Z from
a parametric description of ∂Ω,
possibly consisting of multiple patches,
and then placing interior nodes Y_int in Ω
by advancing the front inside the domain starting from Z.
Node sets Z and Y_int are disjoint, and their nodes
are spaced according to locally varying spacing functions
h_Z ∂Ω→^+ and
h_Y Ω→^+,
although we restrict ourselves to constant spacing

    h_Z ≡ h_Y ≡ h ∈^+

for all numerical experiments, which implies that the nodes
are quasi-uniform. Outward-pointing unit normals
νZ are also computed and stored,
because they are required to assemble matrix B
in Algorithms <ref> and <ref>.
We use a custom implementation of the advancing front method,
with the set Z generated according to
<cit.>, and Y_int according to
<cit.>.

In the case of rejection sampling, a set Y_H of
N_H = (  h^-d) nodes is initially generated in a bounding
box  around Ω, 
and then only nodes inside Ω are kept by using a
level set function such that

    Ω = { x ∈|φ(x) < 0 },
        
    Y_int = { y ∈ Y_H |φ(y) < 0 }.

Nodes in a first-order approximation of a tubular
neighborhood of ∂Ω of width 2h are
projected to the zero set of φ
along the direction parallel to the gradient of φ,
and then a thinning operation is applied to the projected
points to ensure that no pair of nodes in Z
has distance smaller than h.

Regardless of how the sets Y_int and Z are generated,
we produce both open quadrature formulas, for which
Y ∩∂Ω = ∅, and closed
quadrature formulas, for which
Y ∩∂Ω≠∅.
In all of our numerical experiments, the former are obtained by
taking Y = Y_int, whereas the latter are obtained by
taking Y = Y_int∪ Z.
Note that the accuracy and stability of 
are also affected by the choice of Y, because
the weight vectors  and  are computed
simultaneously.

In all tests, unless stated otherwise, the set X of
discretization nodes used to define the numerical
differentiation scheme is generated
in the same way as the closed version of Y,
but with a larger spacing h_X = 1.6   h.


Node generation can be a stochastic process,
which applies in particular to the advancing front method
and rejection-sampling of random nodes,
or a deterministic process whose starting conditions
are arbitrary and can therefore be randomized,
like rejection-sampling of Halton nodes.
In either case, node sets X, Y, Z can be
understood as functions of a seed, an integer
variable i_s that initializes the pseudo-random number
generators used in all cases.
A naive comparison of quadrature errors based on a
single choice of seed can lead to misleading results,
because all errors are subject to random noise, and 
so any particular node distribution may overperform
or underperform compared to the others, sometimes
by even an order of magnitude. Numerical quadrature is much more
affected by this than for example the data fitting problem, where
the maximum or the root mean square error over an entire domain
is computed, which significantly reduces the influence of the
noise in the pointwise errors.

Therefore we compute average
errors as follows.
Let e(·,i_s) be the relative quadrature
error for a given seed i_s ∈.
We denote by e_rms(·) 
the root mean square values (RMS) of
e(·,i_s) for i_s = 1, …,n_s:

    e_rms(·)
    = (1/n_s∑_i_s = 1^n_s e(·,i_s)^2
    )^1/2,

where we take n_s=64.
The stability constants K_
and K_ are obtained by averaging the values

    K_(i_s)
    = 1/Ω∑_j=1^N_Y(i_s)_j(i_s)  and  
    K_(i_s)
    = 1/∂Ω∑_j=1^N_Z(i_s)_j(i_s)

over all seeds i_s = 1, …, n_s.

Table <ref> reports the errors and stability
constants for open and closed quadrature formulas
computed by Algorithm <ref>
using the three node generation methods described above
on Ω_1 and ∂Ω_1.
Numerical differentiation is performed by meshless finite
difference formulas of order q=5, and
the Runge function f_1 is centered at x_R = (0,0).
In all rows of Table <ref>, care was taken
to tweak the value of h so that N_Y ≈ 2500.
We observe that nodes generated by the advancing
front algorithm lead to the most accurate and stable
quadrature formulas, although Halton nodes and even
random nodes remain competitive. This suggests that
our quadrature formulas are robust with
respect to the exact placement of the input
nodes in Y and Z. On average, closed quadrature
formulas are more accurate and stable than
their open counterparts, and so they will be used in all
subsequent experiments. Note that the errors for the worst of 64
seeds are less than three times higher than the RMS errors in the
table in all cases except of the open quadrature on random nodes,
where they may get to almost 8 times higher. The worst stability
constant K_w is 1.71 for the closed quadrature on advancing
front nodes and is significantly higher than the average only for
random nodes, with 8.29 for the closed and 233.21 for the open
formulas. The largest constant K_v for the open formulas on
random nodes is 12.72 and never exceeds 1.07 in other cases.








  §.§.§ Choice of spacing parameters h_X and h_



In order to keep the linear system (<ref>) underdetermined
and solvable, the spacing parameters h_X and h_ must
be sufficiently large, because they determine the size of
the system matrix. Nevertheless, as seen by the 
considerations in Section <ref>,
they should remain comparable to h 
to ensure that the numerical differentiation scheme
is asymptotically accurate in the sense that

    ε_L(u) = (h^q-k)
      and  ε_B(u) = (h^q-k).

Therefore, we suggest to choose h_X and h_ as
multiples of h with some fixed coefficients, and investigate
how the errors of the quadrature formulas depend on the ratios
h_X/h and h_/h. 



In Algorithms <ref> and <ref>
we recommend to choose the spacing parameters as
h_X = 1.6   h and h_ = 4   h
based on numerical evidence that we have gathered in
our numerical experiments on several 2D and 3D domains.
The case of the ellipse Ω_1 is presented in
Figure <ref>, where quadrature errors
are plotted as a function of the ratios h_X/h and h_/h
for h = 0.025 and q = 5. We use node sets X,Y,Z generated
by the advancing front algorithm, with Z ⊂ Y.
As in the previous test, data points are the RMS of 64 quadrature
errors given by different seeds. The Runge function f_1 is
centered at x_R = (0,0).
To aid visualization, instead of the errors themselves
we plot their relative size compared
to the choices h_X = 1.6   h and h_ = 4   h.
The vertical axis is in logarithmic scale. In addition we plot 
the stability constant K_w. 



The two plots clearly show that quadrature errors
decrease as the ratios h_X/h and h_/h get smaller,
but only up to a point: when the ratios get too small,
the errors and the stability constant K_
quickly rise up again.
Reducing the ratios even further leads to overdetermined systems
for which no exact solution exists, and so their errors are not
included in the plots. The ratios
h_X = 1.6   h and h_ = 4   h
were found to be a good compromise between accuracy and stability,
as they are close the minimum of the curves in
Figure <ref>, but leaning on the side
of stability, so that a safety margin is left
in Algorithms <ref> and <ref>.
In our experiments, the same ratios were found
to also be appropriate for 3D domains. 
The stability constant for  is not included
in Figure <ref> because its
dependence on h_X/h and h_S/h is very weak.

Even though increasing the ratios h_X/h or h_/h
reduces the number of rows of the system matrix A, it does not
change the total number of nonzeros in A. Indeed, the number of
nonzeros in every column of the matrix (L^T -B^T)
depends in the MFD case only on the size n_L or n_B of the
set of influence of the corresponding node in Y or Z, as chosen
in Step 4 of Algorithm <ref>. Likewise, in the BSP case
the number of nonzeros in every column in general equals the
number q^d of B-splines of degree q-1 that include the
corresponding quadrature node in the interior of their support.
Taking into account the last row of A corresponding to the
non-homogeneous constraint, we get the following upper bound for
the number of nonzeros:

    (A)≤ N_Z+
    
    n_LN_Y+n_BN_Z,    for Algorithm <ref>,
    
    q^d(N_Y+N_Z),    for Algorithm <ref>.

Note that we do not test how the performance of the MFD
quadrature depends on the parameters n_L and n_B, leaving
this to a future work, where also more sophisticated methods for
the selection of the sets of influence could be investigated. As
mentioned in Section <ref>, we rely on the safe and
simple but possibly not optimal selection of k nearest
neighbors in X to a node in Y or Z, with k equal to twice
the dimension of the corresponding polynomial space. As a result,
(A) is somewhat larger for MFD than for BSP in the 2D case,
for the same q, and about half of it in 3D. 



  §.§.§ Choice of operators  and 



In Section <ref> we
discussed two possible ways to choose operators
 and  for our scheme, the elliptic approach
with = Δ and = ∂_ν, and the
divergence approach, with =
and = γ_ν. Theoretical considerations indicated
that the elliptic approach should be appropriate for
domains with smooth boundary.
However, the theory of the elliptic approach breaks down when
the boundary is not smooth. In the same time,
the divergence approach looks promising for piecewise smooth
domains, even though the theory to justify it is far from complete 
because of the lack of results on the regularity of the
corresponding boundary  value problem (<ref>).

In this section we compare the numerical performance of both
approaches on two 2D domains, the ellipse Ω_1 with a
smooth boundary, and the disk sector Ω_2 defined in the
polar coordinates (r,θ) centered
at the origin by

    Ω_2 = {
    (r,θ) ∈^2 |
    0 < r < 1  and  0 < θ < 3π/2
    }.

The domain Ω_2 has non-smooth boundary, with a reentrant
corner at (0,0), and serves as a standard benchmark for
testing numerical methods for the elliptic boundary value problem
(<ref>) that take into account the reduced smoothness
of the solution at the reentrant corners.


We use numerical differentiation schemes (Ψ,L,B)
based on meshless finite difference formulas for both approaches,
the same sets of quadrature nodes Y and Z generated
by the advancing front algorithm with Z ⊂ Y,
and the same non-homogeneous constraints
with (,) ≡ (0,1).
In the divergence approach we generate X with
the usual advancing front method and spacing parameter
h_X = 1.6   h, whereas in the elliptic approach
we pick X = Y. We have found this simple choice to
be more accurate and reliable compared to the generation
of a coarser set of nodes X (Figure <ref>
is specific to the divergence approach).
The Runge function f_1 is centered at x_R = (0,0)
in the case of the ellipse Ω_1, and at
x_R = (cos(3π/4)/2,sin(3π/4)/2)
in the case of the disk sector Ω_2.
We have chosen the polynomial order as q = 4+k, where
k is the order of the differential operator , so that
the errors are expected to decay as h^4 as h→0, unless the
convergence order is saturated by the low regularity of the
solution u of the boundary value problem
(<ref>).


Figure <ref> presents relative quadrature errors
as a function of spacing parameter h for both the elliptic  and
the divergence approach for the ellipse Ω_1 (top row) and
for the disk sector Ω_2 (bottom row). Relative errors are
the RMS of 8 outcomes given by distinct seeds.



Comparing the plots in the top row of
Figure <ref> for the smooth domain Ω_1,
we observe that although the errors are smaller for the
divergence approach, the decay is h^4 in both cases, as
predicted by the theory presented in Sections <ref>
and <ref>. In contrast to this, the plots in the bottom
row for the non-smooth disk sector Ω_2 show that in the case 
of the elliptic approach the quadrature errors are barely decreasing
as h → 0, whereas the divergence approach still delivers
the expected h^4 convergence order.
This is the very compelling reason why the divergence
approach is recommended in Algorithms <ref>
and <ref>, despite the additional complexity
introduced by the numerical differentiation of a  
vector field instead of a scalar function.

Note that the stability constants K_ and K_ of the
quadrature formulas generated by the elliptic approach
on the disk sector never exceed 1.1
in the tests of Figure <ref>(c).
This provides further evidence that the
large quadrature errors are not caused by the instability
of the computed formulas, but rather by the ineffectiveness
of the numerical differentiation scheme on solutions u
to the boundary value problem (<ref>)
on the disk sector due to their low regularity near the reentrant
corner, which leads to large recovery errors
ε_L(u) and ε_B(u) and influences all
estimates of Section <ref> starting
from (<ref>).



In addition to the experiments in this section, those
in Section <ref> for 4≤ q≤ 8 on the
disk sector Ω_2 and an L-shaped 3D domain (see
Figures <ref> and <ref>), also
confirm that  the divergence approach does not
suffer from saturation of the convergence order of the
quadrature formulas.
This can be interpreted as numerical evidence in support
of our conjecture that the boundary value problem
(<ref>) for the divergence operator
has solutions of high regularity order on Lipschitz
domains with piecewise smooth boundary, whenever the right hand
side is highly regular, see
Section <ref> for a more detailed discussion.




  §.§.§ Choice of non-homogeneous constraints



At least one non-homogeneous constraint
(<ref>) has to be included
in the linear system (<ref>), or else its minimum-norm
solution will be trivial.
In Table <ref> we have compared
quadrature errors and stability constants for
different combinations of the non-homogeneous
constraints introduced
in Section <ref>:

    (,) ≡ (1,0),
        
    (,) ≡ (0,1),
        
    (,) ≡ (1,-1),
        
    (,) ≡ (0,∂_νΦ(·,x_0))
     with x_0 ∈Ω.

The integration domains are the ellipse Ω_1
and the disk sector Ω_2, as defined previously.
The spacing parameter is h = 0.025 and the polynomial
order is q = 5. Nodes in X,Y,Z are generated
by the advancing front algorithm with Z ⊂ Y,
and relative errors are the RMS of 64 quadrature errors
given by different seeds. The Runge function f_1 is
centered at x_R = (0,0) in the case of the
ellipse Ω_1, and at x_R = (cos(3π/4)/2,sin(3π/4)/2)
in the case of the disk sector Ω_2.
The fundamental solution is centered at x_0 = (0.1, 0.05),
a point that does not belong to
any axis of symmetry of the two domains.



The results in the table show that quadrature errors
depend quite weakly on the choice of non-homogeneous
constraints. In particular, on the ellipse, whose boundary
is smooth, there is essentially no difference in the errors,
whereas on the disk sector, with a non-smooth boundary,
a difference of one order of magnitude or more can be
observed between enforcing exactness for
constants and enforcing exactness for
the normal derivative of the fundamental solution
centered at x_0 ∈Ω.
This is the reason to recommend the choice (,) ≡ (0,1)
in Algorithms <ref> and <ref>
whenever the measure of ∂Ω is known or can be
sufficiently accurately approximated at low cost, as in the case of domains
whose boundary is described by parametric patches with standard
parameter domains for which accurate quadrature formulas are
available.
Whenever this is challenging, the approach based on
the fundamental solution remains effective, although
a smaller spacing parameter h may be needed to achieve
the same errors on domains with non-smooth boundary.
As far as stability is concerned, all the non-homogeneous
constraints perform equally well. Boundary quadrature formulas
are close to being positive, but no constraint can consistently
deliver non-negative weights across all 64 seeds.



  §.§.§ Choice of minimization norm



In Section <ref>,
we have discussed relative merits of the norms

    (,)_1 = _1 + _1,
        (,)_2 = √(_2^2 + _2^2)

as candidates for the optimization norm
(,)_♯ to be used in the final step
of the general framework described in
Algorithm <ref>, where the 
quadrature weights ^* and ^* are computed
by minimizing (,)_♯
over all solutions to the non-homogeneous
linear system Ax=b of (<ref>).
We now present numerical experiments to compare their
performance and support our choice of the 2-norm for
the practical algorithms of Section <ref>.

The 1-norm minimization problem can be reformulated
as the linear program

    Minimize  1^T x^+ + 1^T x^-
      subject to   A x^+ - A x^-  = b
      and   x^+, x^- ≥ 0,

with x=x^+ -x^-, and solved by e.g. the dual simplex algorithm
or an interior-point method. For our numerical tests,
we have used their implementations
provided by the 2024a release of MATLAB's Optimization Toolbox.
For the solution of the 2-norm minimization problem,
we have used a direct method based on the sparse
QR factorization provided by the SuiteSparseQR library
<cit.>, version 4.3.3 available from <cit.>.

The integration domain for these tests is the ellipse Ω_1,
the polynomial order is q = 5, and 
the sets of nodes X,Y,Z are generated as usual by the
advancing front algorithm with Z ⊂ Y.
The Runge function f_1 is centered at x_R = (0,0).

Table <ref> reports quadrature errors
and stability constants of the formulas produced by
Algorithm <ref> when minimizing the norms above,
denoted in the table by L^1 and L^2.
Separate rows are devoted to the 1-norm minimization computed
by the dual simplex algorithm and the interior-point method,
because their results are significantly different.
For the same tests, Table <ref> reports
further information such as the runtime t_solve
of the solver, the number of the rows m of the system matrix A,
and the number of nonzero weights in the quadrature
formulas  and .
For each value of the spacing parameter h, results
are averaged over 8 different seeds. Values which are meant
to be integers, such as m, have been rounded to the nearest
integer, to enhance readability.


It turned out that the interior-point method for 1-norm minimization may
stall or lose feasibility when solving the linear
program (<ref>).
For this reason, we have included the column "conv" in the tables
that reports how many runs have converged out of 8,
and all statistics are computed only for the seeds for
which convergence has been achieved. 





As we see in Table <ref>,
minimization of the 1-norm leads on average to smaller
quadrature errors compared to the 2-norm,
although the advantage becomes less significant
as h decreases.
The L^1 methods also deliver more stable formulas,
especially on the boundary: the appearance of the exact
value 1 for K_ indicates that a positive formula
(Z,) has been found for all seeds, because
Algorithm <ref> relies on the non-homogeneous
constraint (,)≡(0,1).
The L^2 method only delivers positive formulas
on the boundary for a few seeds.

The L^1 methods, however, struggle with computational
efficiency and robustness. The interior-point method
does not converge for any seed when h = 0.025,
even if parameters such as ,
, and 
are increased significantly from their default values.
The dual simplex algorithm always terminates successfully
in these tests, but its running time is orders of magnitude
longer compared to SuiteSparseQR. This is the reason why
the 2-norm is recommended in Algorithms <ref>
and <ref> for general use.


The two algorithms for the minimization of the 1-norm
that we have considered so far are not the only iterative
algorithms that struggle with the system Ax=b: the ones
for the minimization of the 2-norm, such as LSQR <cit.>,
also require too many iterations to be practical,
unless an effective preconditioning strategy is employed
(the condition number of the minimum-norm
least squares problem naturally grows as h → 0).
Since out-of-the-box performance of standard left and right
preconditioners such as symmetric successive
over-relaxation and incomplete Cholesky factorization
was not found to be satisfactory, we opted for a direct solver
based on a sparse QR factorization of A^T, as implemented
in the command  provided by the
MATLAB interface to the open source library SuiteSparseQR
<cit.>. Fill-in during sparse QR factorization
is mitigated by the reordering algorithms used by SuiteSparseQR,
which were found to be quite effective, especially in the
case of 2D domains.
The QR factorization provided by SuiteSparseQR is rank-revealing,
and the numerical rank of A is taken into account when
computing the solution to Ax=b with minimal 2-norm.
To maximize accuracy, the tolerance 10^-15
was passed to  for numerical rank
computation in all tests in this paper, unless stated otherwise.

By the well-known properties of the linear program
(<ref>), its solution is not necessarily unique, and
there is at least one solution x with the number of nonzero
components not exceeding the number m of rows in the matrix A.
The dual simplex algorithm is known to terminate at a solution
with this property, and indeed we see that
()+() is less than m in all cases
when this method is used, meaning that the combined
quadrature weights x = (,)^T have a significant
amount of zeros, which may be seen as an advantage.
Moreover, the quadrature formulas (Y,) and (Z,)
are also sparse when considered individually, because
() < N_Y and () < N_Z.
In contrast to this, all weights coming from the interior-point
method and the 2-norm minimization are always nonzero.
Even sparser weights can be obtained by the dual simplex
method with a larger ratio h_X/h, although the accuracy
of the resulting quadrature formulas will be reduced.



 §.§ Convergence tests in 2D and 3D



Now that numerical evidence for the choices made
in Algorithms <ref> and <ref>
has been given, we move on to the evaluation
of relative quadrature errors
for h → 0, i.e. for an increasing number of
quasi-uniform quadrature nodes, on the six
test domains depicted in Figure <ref>.



The first two domains are the ellipse Ω_1
and the disk sector Ω_2 already defined
in Sections <ref> and <ref>.
The third and last 2D domain  is
the region enclosed by a Cassini oval,
a quartic plane curve defined as the locus
of points such that the product of the distances
to two fixed points (-a,0) and (a,0),
called foci, is a constant b^2 ∈^+:

    Ω_3 = { (x_1,x_2) ∈^2 |((x_1+a)^2+x_2)^2
    ((x_1-a)^2+x_2)^2 - b^4 < 0 }.

For a < b < a√(2), the domain Ω_3
has smooth boundary and is simply connected,
although it is not convex: its shape resembles
the number eight rotated sideways. 
For our numerical experiments, we have taken
a = 0.95 and b = 1 to stress its non-convex shape
and to ensure that its area is close to that
of Ω_1 and Ω_2.
As far as 3D domains are concerned, Ω_4 is an
ellipsoid with semiaxes (1,0.7,0.7), Ω_5 is
an L-shaped domain obtained by cutting one quadrant
from a rectangular cuboid with sides of length (2,2,2/3),
and Ω_6 is the solid torus with major radius R = 1
and minor radius r = 0.32.

For all tests we use node sets X,Y,Z generated
by the advancing front algorithm
(see Section <ref>), with Z ⊂ Y
and spacing parameter h ranging from 0.005 to 0.16
on 2D domains and from 0.025 to 0.1 on 3D domains.
Note that the size of 2D domains has been chosen so that
their areas are approximately equal, and the same has been
done for 3D domains and their volumes. This means
that the advancing front algorithm generates a similar
number N_Y of internal quadrature nodes for a given
spacing parameter h > 0, and so the plots of interior
quadrature errors can be directly compared across different
domains of the same dimension.
The exact number of quadrature nodes in Y and Z and the size
of the system matrix A are reported in Table <ref>
for the domains Ω_1 and Ω_4, exemplifying 2D
and 3D domains, respectively. The data for the other four domains is
very similar to Ω_1 or Ω_4, except that the
size of Z varies somewhat as it depends on the measure of
∂Ω.
Variables m_MFD and m_BSP denote the number
of rows of A in the case of Algorithms <ref>
and <ref>, respectively.
We report m_BSP for the order q=5, but the numbers for other
orders used in the tests are similar.
The ratios m_MFD/N and m_BSP/N are
below one, and this confirms that the linear system
Ax=b is underdetermined.




Relative quadrature errors for each test domain
and for the polynomial orders q between 4 and 8
are shown as a function of the spacing parameter h 
in Figures <ref> to <ref>.
Results obtained by Algorithms <ref>
and <ref> are shown
in logarithmic scale side by side, with the same
range for the horizontal and vertical axis to aid
visualization and comparisons.
Reference slopes corresponding to different
convergence orders q-1 are shown using dotted lines.
The center x_R of the Runge function f_1 for each test
domain is specified in the caption of its corresponding figure,
and is chosen so that x_R ∈Ω.
Boundary quadrature errors for the test function g_2 were omitted
for the sake of brevity, as these values are very similar to the
ones for g_1 in all cases.
Unlike the numerical tests in
Section <ref>,
relative quadrature errors are now reported
for a single seed i_s = 1. Because of this, the error curves
in Figures <ref>–<ref>
look somewhat noisy compared to the convergence
plots for e.g. the numerical solution of boundary
value problems, 
where typically the RMS error over the entire domain is reported, 
which averages out the noise present in individual pointwise
errors.
We have chosen to plot quadrature errors for a single seed
to stress that convergence is still evident despite the noise
induced by the node generation process,
and that no averaging with respect to the seed is needed
in practice to get accurate and reliable results.

We indeed observe in all the figures a high order of convergence
for both methods MFD and BSP, including the non-smooth domains
Ω_2 and Ω_5, with the slopes of the error plots
increasing with q and roughly corresponding to the expected
h^q-1 convergence order, in accordance with error estimates
(<ref>), (<ref>), (<ref>) and
(<ref>), if the noise induced by the node generation
process is factored out. In the 2D setting, errors close to
the accuracy limits of double-precision floating-point
numbers are reached in all tests, with values
around 10^-15 being typical for the highest order
methods on the largest number of nodes considered for the tests.
To reach such small errors in 2D tests, the rank detection
tolerance used by SuiteSparseQR (see
Section <ref>)
had to be manually set to 10^-15, or else errors would plateau
around 10^-12.
In the 3D setting, where such high precision is not reached,
the rank detection tolerance used by SuiteSparseQR
was kept to its automatically assigned value (about 10^-10
for most tests), a choice motivated by better errors
and stability constants for coarser sets of nodes.

Comparison of the results for MFD versus BSP reveals
similar behavior of the two approaches, with comparable
quadrature errors across all tests for methods with the same
polynomial order q. The main difference, although small,
is that the MFD method tends to perform better than BSP
on 3D domains with non-smooth boundary such as Ω_5,
whereas the BSP method tends to perform better than MFD
for the integration of Franke's function (the smoother
one among f_1 and f_2). Note that the numbers in the 
columns m_MFD and m_BSP of
Table <ref> indicate  that the system matrices for MFD
are much larger than  those for BSP, which increases the fill-in
for the direct  method of SuiteSparseQR. Nevertheless, as
discussed in  Section <ref>, the number of nonzeros
in the system matrix for MFD is comparable in 2D to that for BSP,
and about a half of it in 3D.

Concerning stability constants, the values of K_
and K_ are not reported in the figures,
but in all cases do not grow as h → 0, and
typically remain bounded below 2.
For the same choices of h and q, Algorithm <ref>
may produce larger stability constants, with values
up to 100 on the coarsest node distributions and
largest polynomial order q=8 (assuming that
system (<ref>) is underdetermined),
whereas Algorithm <ref> always leads to stability
constants smaller than 4 in our experiments.
For finer node distributions or smaller polynomial orders,
both methods produce stability constants smaller than 2.

In all tests, quadrature errors are dropped from the plots
whenever the linear system Ax=b is overdetermined,
i.e. the matrix A has more rows than columns.
Moreover, in all MFD tests, quadrature errors
are dropped from the plots whenever n_L > N_X, that is,
when the size of the sets of influence for the discretization
of  exceeds the size of the set X.
These two measures can be understood as rough safeguards
against ill-chosen parameters h and q.

In the case of Algorithm <ref>, the support
of some B-splines may have small intersection with Ω,
and this is known to be a source of instability in
immersed methods, see e.g.<cit.>. 
In our numerical experiments, however, we have not
observed any instability as h → 0,
and a likely explanation is that the minimization of 2-norm has
a regularizing effect on the computation of quadrature weights.
Nevertheless, stabilization techniques from the literature on immersed
methods may further improve the accuracy and stability
of the BSP approach.


























§ ACKNOWLEDGEMENTS

Bruno Degli Esposti is member
of the Gruppo Nazionale Calcolo Scientifico - Istituto
Nazionale di Alta Matematica (GNCS-INdAM).
The INdAM support through GNCS is gratefully acknowledged.

abbrv 
 

