Computational aspects of simultaneous Gaussian quadrature

In this paper, we derive a new method to compute the nodes and weights of simultaneous n-point Gaussian quadrature rules. The method is based on the eigendecomposition of the banded lower Hessenberg matrix that contains the coefficients of the recurrence relations for the corresponding multiple orthogonal polynomials. The novelty of the approach is that it uses the property of total nonnegativity of this matrix associated with the particular considered multiple orthogonal polynomials, in order to compute its eigenvalues and eigenvectors in a numerically stable manner. The overall complexity of the computation of all the nodes and weights is O(n^2).