Stable computation of high order Gauss quadrature rules using discretization for measures in radiation transfer

Abstract

The solution of the radiation transfer equation for the Earth's atmosphere needs to account for the reflectivity of the ground. When using the spherical harmonics method, the solution for this term involves an integral with a particular measure that presents numerical challenges. We are interested in computing a high-order Gauss quadrature rule for this measure. We show that the two classical algorithms to compute the desired Gauss quadrature rule, namely the Stieltjes algorithm and the method using moments are unstable in this case. In their place, we present a numerically stable method to compute Gauss quadrature rules of arbitrary high order. The key idea is to discretize the measure in the integral before computing the recurrence coefficients of the orthogonal polynomials which lead to the quadrature rule. For discrete measures, one can use a numerically stable orthogonal reduction method to compute the recurrence coefficients. Refining the discretization we arrive at the nodes and weights of the Gauss quadrature rule for the continuous case in a stable fashion. This technique is completely general and can be applied to other measures whenever high-order Gauss quadrature rules are needed.