Alvise Sommariva
Alvise Sommariva
Punti per interpolazione, mesh polinomiali e formule di cubatura in Matlab e Python
Interpolation pointsets, polynomial meshes and cubature rules in Matlab and Python
Purpose of this homepage. Cubature rules in Python, orthogonal polynomials in Python, weakly admissible meshes, cubature on some domains not included in this homepage (e.g. polygons, some circular domains, planar domains with boundary defined by splines/NURBS, spherical polygons, polyhedra, etc).
Citation. In case you use software from this homepage, please provide reference in your work. Thank you.
License. These programs are open-source. You can use, redistribute and/or modify them under the terms of the GNU/General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your opinion) any later version.
In this homepage we have stored several pointsets suitable for interpolation or cubature on intervals, simplex (triangle), square, disk, sphere and tetrahedron.
A comprehensive list of cubature rules in Phyton by Nico Schlömer can be found at
https://pypi.org/project/quadpy/
https://github.com/sigma-py/quadpy
A comprehensive list of orthogonal polynomials in Python by Nico Schlömer can be found at
https://github.com/sigma-py/orthopy
A comprehensive list of weakly admissible meshes (WAM) is available at
https://www.math.unipd.it/~marcov/wam.html
For cubature on domains that are not included in this page, as
trigonometric integration on subintervals of [0,2Π]
general polygons,
domains reated to circular arcs (e.g. sectors, circular segments, lunes, lenses),
bivariate domains whose boundary is defined by splines or NURBS,
polygonal domains with a circular edge,
union of disks,
spherical polygons,
spherical geographic rectangles,
union of balls (via QMC),
polyhedra,
see the software section
https://sites.google.com/view/alvisesommarivaunipd/home-page/software
» Interpolation
General set with low Lebesgue constant: [matlab]
Extended Chebyshev: [matlab]
Gauss-Legendre-Lobatto (Fekete points in the interval): [matlab]
» Quadrature
Personal routines
Fejer I, Fejer II, Clenshaw-Curtis w.r.t.
Legendre weight function [matlab]
Gegenbauer (ultraspherical) weight function [matlab]
Zip file containing the routines above and demos: [matlab]
Routine based on J. Waldvogel work
Clenshaw-Curtis (w.r.t. Legendre weight function) [matlab]
Routines supplied by Dirk Laurie and Walter Gautschi
Jacobi weight function (includes the cases of Legendre, Chebyshev, Gegenbauer weight functions)
Recurrence coefficients for monic Jacobi polynomials [matlab]
Gauss quadrature rule from recurrence coefficients [matlab]
Gauss-Jacobi-Lobatto
Recurrence coefficients for monic Jacobi polynomials [matlab]
Gauss-Lobatto quadrature rule for Jacobi weight function from Jacobi recurrence coefficients [matlab]
Gauss-Jacobi-Radau
Gauss-Radau quadrature rule for Jacobi weight function [matlab]
Recurrence coefficients for monic Jacobi polynomials [matlab]
Hermite weight function
Recurrence coefficients for monic generalized Hermite polynomials [matlab]
Gauss quadrature rule from recurrence coefficients [matlab]
Laguerre weight function
Recurrence coefficients for monic generalized Laguerre polynomials [matlab]
Gauss quadrature rule from recurrence coefficients [matlab]
Useful routines:
Modified Chebyshev algorithm [matlab]
Clenshaw algorithm (evaluation of a polynomial from three-term recurrence relation) [matlab]
Gauss-Kronrod quadrature formula (requires "recurrence coefficients") [matlab]
More details of the impressive work of W. Gautschi can be found in original homepage (now offline).
Routines of the Chebfun project
Gaussian rules (Chebfun repository) w.r.t.:
Legendre weight function [matlab]
Jacobi weight function [matlab]
Gegenbauer (ultraspherical) weight function [matlab]
Chebyshev weight function [matlab]
Gauss-Jacobi-Lobatto [matlab]
Gauss-Jacobi-Radau [matlab]
Hermite weight function [matlab]
Laguerre weight function [matlab]
Routines of the "Cagliari Numerical Analysis" GitHub repository
Approximation of 1D and 2D integrals using anti-Gauss quadrature rules [matlab]
» Polynomial basis
→ Note: The reference interval is [-1,1].
Vandermonde matrix defined by orthogonal polynomial basis w.r.t. Jacobi weight functions (several normalizations are taken into account) [matlab]
Vandermonde matrix defined by orthonormal polynomial basis w.r.t. Legendre weight function (several normalizations are taken into account) [matlab], [matlab]
References
Numer. Math., vol. 2, (1960), 197–205.
C. W. Clenshaw and A. R. Curtis
A new method for the numerical integration of a "well-behaved" function over a finite range of argument is described. It consists essentially of expanding the integrand in a series of Chebyshev polynomials, and integrating this series term by term. Illustrative examples are given, and the method is compared with the most commonly-used alternatives, namely Simpson's rule and the method of Gauss.
J. Approx. Soft., 2 issue 1, 2025,
P. Díaz de Alba, L. Fermo and G. Rodriguez
Anti-Gauss formulae have attracted a great deal of interest in recent years. There are two main reasons for this. On the one hand, they make it possible to construct new rules, namely averaged or stratified formulae, which are characterized by a higher accuracy and a lower computational cost, compared to formulae with the same number of nodes. On the other hand, they provide numerical estimates of the error of the Gaussian rule for a fixed number of points. This allows one to determine the number of nodes required to achieve a prescribed accuracy in the integral approximation. In this paper, we provide a Matlab toolbox, including an interactive graphical user interface (GUI), to assist the interested reader in performing the computation of 1D integrals using anti-Gauss formulae and one of their generalization. When coupled to Gauss rules, they lead to averaged and generalized averaged scheme, respectively. We also consider the computation of 2D integrals by Gauss/anti-Gauss rules. With the exception of the GUI, the software also runs on Octave.
ORTHPOL: A package of routines for generating orthogonal polynomials and Gauss-type quadrature rules
ACM Transactions on Mathematical Software (TOMS), Volume 20, Issue 1 (1994), pp. 21 - 62.
W. Gautschi
A collection of subroutines and examples of their uses, as well as the underlying numerical methods, are described for generating orthogonal polynomials relative to arbitrary weight functions. The object of these routines is to produce the coefficients in the three-term recurrence relation satisfied by the orthogonal polynomials. Once these are known, additional data can be generated, such as zeros of orthogonal polynomials and Gauss-type quadrature rules, for which routines are also provided.
Acta Numerica , Volume 5 , January 1996 , pp. 45 - 119
We give examples of problem areas in interpolation, approximation, and quadrature, that call for orthogonal polynomials not of the classical kind. We then discuss numerical methods of computing the respective Gauss-type quadrature rules and orthogonal polynomials. The basic task is to compute the coefficients in the three-term recurrence relation for the orthogonal polynomials. This can be done by methods relying either on moment information or on discretization procedures. The effect on the recurrence coefficients of multiplying the weight function by a rational function is also discussed. Similar methods are applicable to computing Sobolev orthogonal polynomials, although their recurrence relations are more complicated. The paper concludes with a brief account of available software.
Science China Mathematics Volume 55 (2012), pp.1749–1760.
N. Hale, N.L. Trefethen
Chebfun is a Matlab-based software system that overloads Matlab’s discrete operations for vectors and matrices to analogous continuous operations for functions and operators. We begin by describing Chebfun’s fast capabilities for Clenshaw–Curtis and also Gauss–Legendre, –Jacobi, –Hermite, and –Laguerre quadrature, based on algorithms of Waldvogel and Glaser, Liu, and Rokhlin. Then we consider how such methods can be applied to quadrature problems including 2D integrals over rectangles, fractional derivatives and integrals, functions defined on unbounded intervals, and the fast computation of weights for barycentric interpolation.
Journal of Computational and Applied Mathematics,Volume 127, Issues 1–2, 15 January 2001, pp. 201-217.
D. P. Laurie
Gaussian formulas for a linear functional L (such as a weighted integral) are best computed from the recursion coefficients relating the monic polynomials orthogonal with respect to L. In Gauss-type formulas, one or more extraneous conditions (such as pre-assigning certain nodes) replace some of the equations expressing exactness when applied to high-order polynomials. These extraneous conditions may be applied by modifying the same number of recursion coefficients. We survey the methods of computing formulas from recursion coefficients, methods of obtaining recursion coefficients and modifying them for Gauss-type formulas, and questions of existence and numerical accuracy associated with those computations.
Mathematics of Computation, Volume 66, number 219, July 1997, pp. 1133–1145.
D. P. Laurie
The Jacobi matrix of the (2n+1)-point Gauss-Kronrod quadrature rule for a given measure is calculated efficiently by a five-term recurrence relation. The algorithm uses only rational operations and is therefore also useful for obtaining the Jacobi-Kronrod matrix analytically. The nodes and weights can then be computed directly by standard software for Gaussian quadrature formulas.
Mathematics of Vomputation, Volume 65, number 214, April 1996, pp. 739–747.
D. P. Laurie
An anti-Gaussian quadrature formula is an (n + 1)-point formula of degree 2n − 1 which integrates polynomials of degree up to 2n + 1 with an error equal in magnitude but of opposite sign to that of the n-point Gaussian formula. Its intended application is to estimate the error incurred in Gaussian integration by halving the difference between the results obtained from the two formulas. We show that an anti-Gaussian formula has positive weights, and that its nodes are in the integration interval and are interlaced by those of the corresponding Gaussian formula. Similar results for Gaussian formulas with respect to a positive weight are given, except that for some weight functions, at most two of the nodes may be outside the integration interval. The antiGaussian formula has only interior nodes in many cases when the Kronrod extension does not, and is as easy to compute as the (n + 1)-point Gaussian formula.
Electron. Trans. Numer. Anal, vol. 45, (2016), 371–404.
S. E. Notaris.
A. S. Kronrod [Sov. Phys., Dokl. 9, 17–19 (1964); translation from Dokl. Akad. Nauk SSSR 154, 283–286 (1964; Zbl 0131.15003)], trying to estimate economically the error of the n-point Gauss quadrature formula for the Legendre weight function, developed a new formula by adding to the n Gauss nodes n+1 new ones, which are determined, together with all weights, such that the new formula has maximum degree of exactness. It turns out that the new nodes are zeros of a polynomial orthogonal with respect to a variable-sign weight function, considered by T. J. Stieltjes [C. R. Acad. Sci., Paris 118, 1315–1317 (1894; JFM 25.0326.02)], without though making any reference to quadrature. We survey the considerable research work that has been emerged on this subject, during the past fifty years, after Kronrod’s original idea.
S. Olver, R.M. Slevinsky, A. Alex Townsend
This paper includes observations on the state of the art of numerical cubature.
Computers Mathematics with Applications, Volume 65, Issue 4, February 2013, Pages 682-693.(Open Access)
A. Sommariva
The main purpose of this paper is to compute the weights of Clenshaw-Curtis and Fejer type quadrature rules via DCT and DST, for general weight functions w. The approach is different from that used by Waldvogel in [25], where the author considered the computation of these sets only for the Legendre weight w ≡ 1 using DFT arguments.
BIT Numerical Mathematics Vol. 43, No. 1, pp. 001–018
J. Waldvogel
We present an elegant algorithm for stably and quickly generating the weights of Fejer’s quadrature rules and of the Clenshaw-Curtis rule. The weights for an arbitrary number of nodes are obtained as the discrete Fourier transform of an explicitly defined vector of rational or algebraic numbers. Since these rules have the capability of forming nested families, some of them have gained renewed interest in connection with quadrature over multi-dimensional regions.
» Interpolation
General set with low Lebesgue constant: [matlab] (last update: Jan 08, 2017, old version: [matlab]).
General set with low Lebesgue constant and Gauss-Legendre-Lobatto distribution on the side: [matlab]
General set with high absolute value of the Vandermonde matrix, i.e. (quasi-) Fekete points: [matlab]
Symmetric set with low Lebesgue constant: [matlab]
Approximate Fekete points (with degree from 1 to 70): [matlab]
» Polynomial mesh
Weakly Admissible Mesh: [matlab]
» Cubature
→ Note: The reference simplex has vertices (0,0), (1,0), (0,1).
Near minimal cubature sets on the triangle (up to degree 50): [matlab]
Slobodkins-Tausch cubature sets on the triangle (up to degree 50): [matlab]
Gauss-Legendre product rule: [matlab]
» Comparisons
For a comparison on several interpolation sets see also: [html].
For a comparison on several cubature sets see also: [html].
» Polynomial basis
You may find useful the M-file for the evaluation of the Vandermonde matrix w.r.t.Dubiner Legendre basis [matlab] as orthonormal basis; see also: [zip] where we test orthogonality of the basis.
Proriol-Dubiner:
a) [matlab] as orthogonal basis (more recent version, no derivatives);
b) [matlab] as orthogonal basis with its derivatives (in a form suitable for cubature, i.e. it is the transpose of the Vandermonde matrix for interpolation purposes).
» Interpolation
General set with low Lebesgue constant: [.m] (last update: Jan 08, 2017, old version: [matlab]).
General set with high absolute value of the Vandermonde matrix, i.e. (quasi-) Fekete points: [matlab]
Padua-Jacobi points with low Lebesgue constant: [matlab]
Padua-Jacobi points with high absolute value of the Vandermonde matrix: [matlab]
Padua points: [matlab]
» Cubature
→ Note: The reference square is [-1,1] x [-1,1].
Sommariva-Festa rules (Legendre weight): [matlab]
Slobodkins-Tausch rules (Legendre weight): [matlab]
Padua points rules (Legendre weight): [matlab]
Near minimal rules (Legendre weight, lowest known cardinality): [matlab]
Morrow-Patterson-Xu (minimal w.r.t. a Chebyshev measure): [matlab]
Tensorial (Jacobi tensor product): [matlab]
Anti Gauss (by CaNA group): [matlab]
» Polynomial basis
→ Note: The reference square is [-1,1] x [-1,1].
Vandermonde matrix defined by tensorial Chebyshev basis: [matlab]
Vandermonde matrix defined by tensorial orthonormal Chebyshev basis: [matlab]
Vandermonde matrix defined by tensorial orthonormal Legendre basis: [matlab]
Polynomial degree ordering (necessary for the definition of the total degree basis ordering): [matlab]
» Polynomial mesh
WAM on Lissajous curves of the cube: [matlab]
» Cubature
Morrow-Patterson like (a certain tensorial Chebyshev measure): [matlab]
Tensor product rule (tensorial Jacobi measure): [matlab]
» Polynomial basis
→ Note: The reference cube is [-1,1] x [-1,1] x [-1,1].
Vandermonde matrix relative to the tensorial Chebyshev basis: [matlab]
Vandermonde matrix relative to the tensorial orthonormal Chebyshev basis: [matlab]
Vandermonde matrix relative to the tensorial orthonormal Legendre basis: [matlab]
Polynomial degree ordering (necessary for the definition of the total degree basis ordering): [matlab]
» Cubature
→ Note: The reference tetrahedron has vertices is [0 0 0], [1 0 0], [0 1 0], [0 0 1]. In case of need we propose also the version with barycentrical coordinates, in which the sum of the weights is normalised to 1.
» Interpolation
Maximum Determinant (Fekete, Extremal) points on the sphere S2 (R. Womersley homepage).
Minimum Energy points on the sphere S2 (R. Womersley homepage).
Recursive Zonal Equal Area (EQ) Sphere Partitioning (P. Leopardi homepage).
» Cubature
Spherical designs: Efficient Spherical Designs with Good Geometric Properties (R. Womersley homepage):
1. Spherical t-designs on S2 with N = t2/2 + t + O(1) points: [matlab zip]
2. Symmetric (antipodal) spherical t-designs on S2 with N = t2/2 + t/2 + O(1) points: [matlab zip]
Spherical designs: Quadrature Rules on Manifolds: Putatively Optimal Quadrature Rules on the Sphere S2 (M. Graf homepage).
Spherical designs: Spherical Designs (R. H. Hardin and N. J. A. Sloane homepage). Each file stores a unique vector v, so that x=v(1:3:end), y=v(2:3:end), z=v(3:3:end) and weights are equal.
Albrecht-Collatz rules: [matlab].
Heo-Xu rules: [matlab].
McLaren rules: [matlab].
» Polynomial basis
Vandermonde matrix defined by Spherical Harmonics (based on a code by Paul Leopardi): [matlab]
Purpose: Computation of quadrature/cubature rules on some less standard domains and special pointsets for least-squares approximation.
» Quadrature
Trigauss 2011 (original version): [matlab]
Trigauss 2020 (faster version): [matlab]
Trigauss 2017 (via specific Gauss-Legendre rule): [matlab]
Trigauss 2013 (version via Fejer I, Fejer II, Clenshaw-Curtis rules): [matlab]
Trigauss Antigauss 2020 (useful for error estimates): [matlab]
Trigauss Gauss-Kronrod rule (useful for error estimates): [matlab]
Zip file containing all the routines above and demos: [matlab (zip)]
Notes:
This demos may require Chebfun routines. See https://www.chebfun.org/.
The software of Trigauss 2013 (version via Fejer I, Fejer II, Clenshaw-Curtis rules) requires Matlab "Signal Processing Toolbox".
» Polynomial basis
Vandermonde matrix defined by subperiodic trigonometric orthonormal basis: [matlab]
Zip file containing the routine above and demo: [matlab (zip)]
Reference papers
G. Da Fies and M. Vianello
Electron. Trans. Numer. Anal. 39 (2012), 102--112.
In this paper routine "Trigauss 2011" have been introduced. The authors construct a quadrature formula with n +1 angles and positive weights, exact in the (2n+1)-dimensional space of trigonometric polynomials of degree ≤ n on intervals with length smaller than 2π. We apply the formula to the construction of product Gaussian quadrature rules on circular sectors, zones, segments and lenses.
Numerical Algorithms, 67, Issue 3 (2014), 491-506.
G. Meurant and A. Sommariva
In this paper the authors provided routines useful to propose "Trigauss 2020" that is in general faster than "Trigauss 2011" . In particular, we investigate variants of the well-known Golub and Welsch algorithm for computing nodes and weights of Gaussian quadrature rules for symmetric weights w in intervals (−a,a) (not necessarily bounded). The purpose is to reduce the complexity of the Jacobi eigenvalue problem stemming from Wilf’s theorem and show the effectiveness of these methods for reducing the computer times. Numerical examples on three test problems show the benefits of these variants.
A. Sommariva and M. Vianello
Appl. Anal. Discrete Math. 11 (2017), 470--483.
In this paper the authors provided the routines of "Trigauss 2017" based on specific Gauss-Legendre rules. In some instances they behave as "Trigauss 2011", "Trigauss 2020", sometimes even with less nodes/weights. We show that Gauss-Legendre quadrature applied to trigonometric polynomials on subintervals of the period can be competitive with subperiodic trigonometric Gaussian quadrature. For example with intervals corresponding to few angular degrees, relevant for regional scale models on the earth surface, we see a subsampling ratio of one order of magnitude already at moderate trigonometric degrees.
A. Sommariva
Computers Mathematics with Applications, Volume 65, Issue 4, February 2013, Pages 682-693.
In this paper the author provided the routines of "Trigauss 2013" based on Fejer I, Fejer II, Clenshaw-Curtis type weighted rules. The advantage is the fast computation of nodes and weights, especially for rules with high trigonometric degree of precision. The main purpose of this paper is to compute the weights of Clenshaw-Curtis and Fejer type quadrature rules via DCT and DST, for general weight functions w. The approach is di erent from that used by Waldvogel, where the author considered the computation of these sets only for the Legendre weight w 1 using DFT arguments.
G. Da Fies, A. Sommariva and M. Vianello.
"Festschrift for the 80th Birthday of Ian Sloan" (invited paper), J. Dick, F.Y. Kuo, H. Wozniakowski Eds., Springer, 2018, pp. 283-304.
In this paper we investigate how to determine an orthonormal basis for the space of trigonometric polynomials in some intervals and show applications to hyperinterpolation. Using recent results on subperiodic trigonometric Gaussian quadrature and the construction of subperiodic trigonometric orthogonal bases, we extend Sloan’s notion of hyperinterpolation to trigonometric spaces on subintervals of the period. The result is relevant, for example, to function approximation on spherical or toroidal rectangles.
» Polynomial mesh
Weakly admissible mesh (via triangulation on a polygonal region) [matlab]
Weakly admissible mesh (via quadrangulation on a polygon without self-intersections) [matlab]
Zip file containing the routines above and demo: [matlab (zip)]
» Cubature
Polygauss (2007): [matlab]
Cubature via triangulation (2020): [matlab]
Cubature adaptive rule: [matlab]
Zip file containing all the routines above and demos: [matlab (zip)]
Reference papers
BIT Numerical Mathematics, 47 (2007), 441-453 (Open Access)
A. Sommariva and M. Vianello
In this work we propose the Polygauss code for cubature over polygons. The rule, may have internal nodes if the domain has the so called "axis-property", as in the case of convex domains. The original code has been refurbished in 2013.
We have implemented in Matlab a Gauss-like cubature formula over convex, nonconvex or even multiply connected polygons. The formula is exact for polynomials of degree at most "2n − 1" using "N ∼ m n^2" nodes, "m" being the number of sides that are not orthogonal to a given line, and not lying on it. It does not need any preprocessing like triangulation of the domain, but relies directly on univariate Gauss-Legendre quadrature via Green’s integral formula. Several numerical tests are presented.
Journal of Computational and Applied Mathematics, Volume 370, 2020, 112658.
B. Bauman, A. Sommariva and M. Vianello.
In this work we propose the code Cubature via triangulation (2020) for cubature over polygons. Differently from Polygauss, the rule, has always internal nodes and positive weights.
In this paper we propose an algorithm to determine cubature rules of algebraic degree of exactness δ on general polygons P, by means of Matlab polyshape objects and near minimal rules on triangles, obtaining by Caratheodory-Tchakaloff subsampling a PI (Positive Interior) final formula with cardinality at most N_δ = (δ + 1)(δ + 2)/2. We test our algorithm on polygons with different shape, and we also discuss an application to the computation of the RMSWE (Root Mean Square Wavefront Error) on obscured and vignetted pupils, in the framework of optical design by numerical ray tracing for the LSST (Large Synoptic Survey Telescope).
Polynomial interpolation and cubature over polygons,
Journal of Computational and Applied Mathematics 235 (2011), 5232-5239.
M. Gentile, A. Sommariva and M. Vianello.
In this paper we introduced WAM via quadrangulations. The WAMs via triangulation were defined in the thesis by F.Basaglia constructed as soon as a triangulation was at hand. We modified the code since the first release of Polyshape toolbox.
We have implemented a Matlab code to compute Discrete Extremal Sets (of Fekete and Leja type) on convex or concave polygons, together with the corresponding interpolatory cubature formulas. The method works by QR and LU factorizations of rectangular Vandermonde matrices on Weakly Admissible Meshes (WAMs) of polygons, constructed by polygon quadrangulation.
Circular sectors, circular segments, circular zones, circular lenses, circular symmetric lenses, linear blending of elliptic arcs, intersection of disks, union of disks.
» Polynomial mesh
Linear blending of elliptic arcs (see initial comments on the routine for application to a specific domain): [matlab]
Circular lune: [matlab]
Zip file containing the routines above and demo: [matlab (zip)]
» Cubature
Circular sector: [matlab]
Circular segment: [matlab]
Circular zone: [matlab]
Circular lens: [matlab]
Circular lune: [matlab]
Circular symmetric lens: [matlab]
Linear blending of elliptic arcs: [matlab]
Zip file containing all the routines above and demos: [matlab (zip)]
Linear blending of elliptic arcs (adaptive): [matlab]
Zip file containing the routine above and demo: [matlab (zip)]
Intersection of many disks [matlab]
Zip file containing the routine above and demo: [matlab (zip)]
Union of many disks [matlab] (previous slower version [matlab])
Zip file containing the routines above and demos: [matlab (zip)]
Reference papers
Applied Numerical Mathematics, 74 (2013), 49-61. (Open Access)
G. Da Fies, A. Sommariva and M. Vianello.
In this paper we introduce cubature rules in the circular domains mentioned above. The routine about "Linear blending of elliptic arcs" allows the description of all those regions and many more. The specific routines are written to simplify the application of the rules to the more common circular planar domains.
We construct a cubature formula of algebraic degree of exactness n with n^2/2 + O(n) nodes, on the bidimensional domains generated by linear blending of two arcs of ellipses corresponding to the same angular interval. The construction is based on recent results on “subperiodic” trigonometric quadrature. Our formula generalizes several recent cubature formulas on standard circular sections. Among its numerous possible applications, we quote for example integration of functions with singularities, and integration on nonstandard circular sections arising in optical design or in meshfree methods with compactly supported radial bases
FILOMAT 31:13 (2017), 4105--4115 (open access)
A. Sommariva and M. Vianello.
In this paper we introduced algorithms to compute rules over intersection and union of planar disks.
We provide an algorithm that computes algebraic quadrature formulas with cardinality not exceeding the dimension of the exactness polynomial space, on the intersection of any number of planar disks with arbitrary radius. Applications arise for example in computational optics and in wireless networks analysis. By the inclusion-exclusion principle, we can also compute algebraic formulas for the union of a small number of disks. The algorithm is implemented in Matlab, via subperiodic trigonometric Gaussian quadrature and compression of discrete measures.
Dolomites Research Notes on Approximation, Volume 15, Issue 4, 2022, pp. 73-81.
A. Sommariva and M. Vianello.
In this paper we suggest an algorithm to compute rules over union of planar disks. The main difference is that it is faster than that proposed in "Numerical quadrature on the intersection of planar disks".
In this work we present a new algorithm that computes cubature formulas with positive weights, interior nodes and fixed algebraic degree of precision, over domains Ω that are arbitrary union of disks. This novel approach first determines the boundary ∂ Ω and then defines a decomposition of Ω by means of nonoverlapping circular segments and polygons, where algebraic positive interior rules can be locally constructed. The resulting global Positive Interior (PI) formula is finally compressed by CaratheodoryTchakaloff subsampling implemented via NonNegative Least-Squares.
Product Gaussian quadrature on circular lunes
Numer. Math. Theory Methods Appl. 7 (2014), 251--264.
G. Da Fies and M. Vianello.
Resorting to recent results on subperiodic trigonometric quadrature, we provide three product Gaussian quadrature formulas exact on algebraic polynomials of degree n on circular lunes. The first works on any lune, and has n^2+O(n) cardinality. The other two have restrictions on the lune angular intervals, but their cardinality is n^2/2 + O(n).
Algebraic cubature on planar lenses and bubbles
Dolomites Res. Notes Approx. DRNA 5 (2012), 7--12.
G. Da Fies and M. Vianello
Study on circular segments, planar lenses and double bubbles, introducing the pertinent software.
By a recent result on subperiodic trigonometric Gaussian quadrature, we construct a cubature formula of algebraic degree of exactness n on planar circular lenses (intersection of two overlapping disks) and “double bubbles” (union of two overlapping disks), with n 2 /2 + O(n) nodes. An application is shown to RBF projection methods.
Journal of Inequalities and Special Functions, Volume 9, 2018, pp. 113-121.
A. Sommariva and M. Vianello.
Study on polynomial meshes on linear blending of elliptical arcs, introducing the pertinent software. Such domains include circular sectors, circular segments, circular zones, circular lenses, circular symmetric lenses.
By discrete trigonometric norming inequalities on subintervals of the period, we construct norming meshes with optimal cardinality growth for algebraic polynomials on sections of sphere, ball and torus.
General spherical rectangles, including caps, slices as special cases.
» Polynomial mesh
Spherical rectangles (cardinality O(2n^2)): [matlab]
Zip file containing the routine above and a demo: [matlab (zip)]
» Cubature
Spherical rectangles: [matlab]
Spherical rectangle (adaptive rule): [matlab]
Zip file containing all the routines above and demos: [matlab (zip)]
Reference papers
Polynomial approximation and quadrature on geographic rectangles,
Applied Mathematics and Computation 297 (2017), 159-179.
M. Gentile, A. Sommariva and M. Vianello
Using some recent results on subperiodic trigonometric interpolation and quadrature, and the theory of admissible meshes for multivariate polynomial approximation, we study product Gaussian quadrature, hyperinterpolation and interpolation on some regions of S^d , d ≥ 2. Such regions include caps, zones, slices and more generally spherical rectangles defined by longitudes and (co)latitudes (geographic rectangles). We provide the corresponding Matlab codes and discuss several numerical examples on S^2.
» Polynomial Mesh
Admissible mesh on spherical triangle [matlab]
Zip file containing the routine above and demo: [matlab (zip)]
» Cubature
Spherical triangles: [matlab]
Spherical polygons: [matlab]
Spherical polygons (adaptive rule): [matlab]
Zip file containing all the routines above and demos: [matlab (zip)]
Notes:
This software requires Chebfun routines. See https://www.chebfun.org/.
The spherical polygon is centered in the origin and contained in a cap with polar angle inferior to "pi".
Reference papers
Near-algebraic Tchakaloff-like quadrature on spherical triangles,
Applied Mathematics Letters, Volume 120, October 2021, 107282.
A. Sommariva and M. Vianello
We present a numerical code for the computation of nodes and weights of a low-cardinality positive quadrature formula on spherical triangles, nearly exact for polynomials of a given degree. The algorithm is based on subperiodic trigonometric gaussian quadrature for planar elliptical sectors and on Caratheodory-Tchakaloff quadrature compression via NNLS.
Mathematics and Computers in Simulation, Volume 190, December 2021, pp. 15-22.
A. Sommariva and M. Vianello
We present a numerical method (implemented in Matlab) for computing an orthogonal polynomial basis on spherical triangles, via a recent near-algebraic quadrature formula, and constructing the corresponding weighted orthogonal projection (hyperinterpolation) of a function sampled at the quadrature nodes.
Mediterranean Journal of Mathematics, 9, 2022, article 68.
A. Sommariva and M. Vianello
By the fundamental notion of Dubiner distance on a compact set, we construct Chebyshev polynomial norming grids in the sup-norm on spherical triangles. These grids can be used to extract Fekete-like interpolation points with slowly increasing Lebesgue constant.
Numerical cubature and hyperinterpolation over Spherical Polygons,
Applied Mathematics and Computation, Volume 495, 15 June 2025, 129335,
A. Sommariva
The purpose of this work is to introduce a strategy for determining the nodes and weights of a low-cardinality positive cubature formula nearly exact for polynomials of a given degree over spherical polygons. In the numerical section we report the results about numerical cubature over a spherical polygon P approximating Australia and reconstruction of functions over such P, also affected by perturbations, via hyperinterpolation and some of its variants.
» Polynomial mesh
Toroidal cap: [matlab]
Toroidal rectangle (longitude-latitude): [matlab]
Toroidal slice: [ma
tlab]
Toroidal zone: [matlab]
Zip file containing the routines above and demos: [matlab (zip)]
» Cubature
Toroidal rectangle (longitude-latitude): [matlab]
Toroidal rectangle (longitude-latitude) adaptive rule: [matlab]
Zip file containing the routines above and demos: [matlab (zip)]
Reference papers
Journal of Inequalities and Special Functions, Volume 9, 2018, pp. 113-121.
A. Sommariva and M. Vianello.
By discrete trigonometric norming inequalities on subintervals of the period, we construct norming meshes with optimal cardinality growth for algebraic polynomials on sections of sphere, ball and torus.
Purpose: Computation of Approximate Fekete Points, Discrete Leja points, compression of cubature rules, determination of well-conditioned Vandermonde basis.
Important: these sources are taken from Marco Vianello homepage.
1. Extraction of discrete extremal sets of Fekete and Leja type from a WAM, 2. Discrete Orthogonal Polynomials on a WAM, 3. Evaluation of DOPs,
4. Polynomial fitting on a WAM, 5. Lebesgue constant of the fitting/interpolation operator.
dexsets (extraction of discrete extremal sets of Fekete and Leja type from a WAM)
wamdop (Discrete Orthogonal Polynomials on a WAM)
wamdopeval (evaluation of DOPs)
wamfit (polynomial fitting on a WAM; note: interpolates on discrete extremal sets)
wamleb (Lebesgue constant of the fitting/interpolation operator)
1. Polynomial bases and Vandermonde-like matrices, 2. Discrete measure compression, 3. Regression design optimization, 4. Polynomial fitting.
>> Polynomial bases and Vandermonde-like matrices
mono_next_grlex.m (d-variate monomial ordering, by J. Burkardt)
dCHEBVAND (d-variate Chebyshev-Vandermonde matrix)
dORTHVAND (d-variate Vandermonde matrix in a discrete orthogonal basis)
vandermonde_jacobi 1-D Vandermonde matrix defined by orthogonal polynomial basis w.r.t. Jacobi weight functions (several normalizations are taken into account)
>> Discrete measure compression
dCATCH (d-variate CAratheodory-TCHakaloff discrete measure compression)
LHDM (Lawson-Hanson NNLS solver accelerated by Deviation Maximization)
>> Regression design optimization
dNORD (d-variate Near (G-)Optimal Regression Designs)
>> Polynomial fitting
dPOLYFIT (d-variate polynomial fitting construction)
dPOLYVAL (d-variate polynomial fitting evaluation)
dLEB (d-variate Lebesgue constant evaluation, i.e. uniform norm of the fitting operator)
demo5bubble (compression of grid points in a 5bubble: union of 5 balls, d=3)
demoChebgrids (compression of Chebyshev grids in d-cubes)
demoHalton (compression of Halton points in d-cubes)
haltonseq (generation of d-variate Halton sequences, by D. Dougherty)
Compression of multivariate discrete measures and applications
A. Sommariva and M. Vianello
Numer. Funct. Anal. Optim. 36 (2015), 1198--1223
dCATCH: a numerical package for d-variate near G-optimal Tchakaloff regression via fast NNLS
M. Dessole, F. Marcuzzi and M. Vianello
MDPI-Mathematics 8 (2020) - Special Issue "Numerical Methods"