A MATLAB package computing simultaneous Gaussian quadrature rules for Multiple Orthogonal Polynomials

TL;DR

This work introduces GaussMOP.m, a Matlab package for computing simultaneous Gaussian quadrature rules associated with multiple orthogonal polynomials (MOPs) of type II in the case $r=2$. The key challenge is the ill-conditioning of the banded lower Hessenberg matrix $H_n$ whose eigen-decomposition yields the nodes and weights; the authors implement a diagonal-balancing scheme that produces a similar matrix $ ilde{H}_n = S_n^{-1} H_n S_n$ with drastically improved conditioning and use an Ehrlich–Aberth-based solver to compute eigenvalues and left/right eigenvectors efficiently. The method delivers nodes and two sets of weights $(oldsymbol{x}, oldsymbol{ ho}^{(1)}, oldsymbol{ ho}^{(2)})$ with $ ext{O}(n^2)$ time and $ ext{O}(n)$ memory, applicable to nine MOP classes via ClassMOP.m and GaussMOP.m. Numerical tests compare GaussMOP.m against high-precision integration, demonstrating reliable accuracy and favorable performance, and the appendix provides complete Matlab code for reproduction. The approach enables reliable floating-point construction of simultaneous quadrature rules for MOPs, with potential impact in numerical analysis and related applications such as random matrix theory and approximation theory.

Abstract

The aim of this paper is to describe a Matlab package for computing the simultaneous Gaussian quadrature rules associated with a variety of multiple orthogonal polynomials. Multiple orthogonal polynomials can be considered as a generalization of classical orthogonal polynomials, satisfying orthogonality constraints with respect to $r$ different measures, with $r \ge 1$. Moreover, they satisfy $(r+2)$--term recurrence relations. In this manuscript, without loss of generality, $r$ is considered equal to $2$. The so-called simultaneous Gaussian quadrature rules associated with multiple orthogonal polynomials can be computed by solving a banded lower Hessenberg eigenvalue problem. Unfortunately, computing the eigendecomposition of such a matrix turns out to be strongly ill-conditioned and the \texttt{Matlab} function \texttt{balance.m} does not improve the condition of the eigenvalue problem. Therefore, most procedures for computing simultaneous Gaussian quadrature rules are implemented with variable precision arithmetic. Here, we propose a \texttt{Matlab} package that allows to reliably compute the simultaneous Gaussian quadrature rules in floating point arithmetic. It makes use of a variant of a new balancing procedure, recently developed by the authors of the present manuscript, that drastically reduces the condition of the Hessenberg eigenvalue problem.

Figures (1)

• Figure 1: Condition number of the eigenvalues of $H_{20}$ (left) and of the eigenvalues of $\hat{H}_{20}$ (right) for all the considered MOPs.

Theorems & Definitions (5)

• Theorem 1
• Lemma 1
• proof
• Remark 1
• Example 1