Fast and reliable high-accuracy computation of Gauss–Jacobi quadrature

Abstract

Iterative methods with certified convergence for the computation of Gauss–Jacobi quadratures are described. The methods do not require a priori estimations of the nodes to guarantee its fourth-order convergence. They are shown to be generally faster than previous methods and without practical restrictions on the range of the parameters. The evaluation of the nodes and weights of the quadrature is exclusively based on convergent processes which, together with the fourth-order convergence of the fixed point method for computing the nodes, makes this an ideal approach for high-accuracy computations, so much so that computations of quadrature rules with even millions of nodes and thousands of digits are possible on a typical laptop.

Appendix. On the conditioning of the recurrence for computing derivatives

We briefly discuss the conditioning of the computation of the derivatives with the recurrence relation (26). This is a five term recurrence relation, with a space of solutions of dimension four. For studying the conditioning as $j\rightarrow \infty$ we can divide all terms of the recurrence by j² and then all the coefficients have finite limit as $j\rightarrow +\infty$. We have

$$ \sum\limits_{n=0}^{4} C_{n} (j) a_{j+2-n}=0, $$

with

$$ \lim_{j\rightarrow \infty}C_{n}(j)=c_{n}=\frac{Q^{(n)}(x)}{n!}. $$

Then it is known (see, for instance, [8], Theorem 8.11) that the solutions of the recurrence satisfy

$$ \limsup_{n\rightarrow +\infty}\left( |a_{n}|\right)^{1/n}=W, $$

(45)

where W is the modulus of one of the solutions of the characteristic polynomial

$$ \sum\limits_{j=0}^{4} c_{n} \delta^{4-n}=0, $$

which, upon dividing by δ⁴ and denoting μ = 1/δ we can write

$$ Q(x)+Q^{\prime}(x)\mu+ Q^{\prime\prime}(x)\frac{\mu^{2}}{2}+Q^{\prime\prime\prime}(x)\frac{\mu^{3}}{3!} +Q^{(4)}(x)\frac{\mu^{4}}{4!}=0. $$

And because Q is a polynomial of degree four (Q(x) = 4(1 − x²)²) this equation is the same as Q(x + μ) = 0, which, solving for μ gives two double roots μ = x ± 1. This means that the possible values of W in (45) are W₁ = 1/|1 − x| and W₂ = 1/|1 + x| and there is a subspace of dimension 2 satisfying (45) with W = W₁ and a second subspace with W = W₂; the first space will be dominant over the second when W₁ > W₂ and the opposite when W₁ < W₂ (the case W₁ = W₂ is the degenerate case, in which no solution is exponentially dominant over the rest).

In our case, we have $a_{j}=\tilde {Y}^{(j)}(x)$ with $\tilde {Y}(x)$ given by (10). We notice that for α and β odd, $\tilde {Y}(x)$ is a polynomial of degree N = n + (α + β + 2)/2 = (L + 1)/2, and then a_j = 0 if j > N. The Taylor series have a finite number of terms in this case and the analysis of stability for a_j as $j\rightarrow \infty$ is not needed. Let us now consider that neither α nor β are odd, and we leave for later the case in which only one of the parameters (α or β) is odd.

When neither α nor β are odd, then Taylor series at x ∈ (− 1, 1) has an infinite number of terms. Because $\tilde {Y}(x)$ is a polynomial times (1 − x)^(α+ 1)/2(1 + x)^(β+ 1)/2 the radius of convergence of the series for $\tilde {Y}(x+t)$ centered at x,

$$ \tilde{Y}(x+h)=\sum\limits_{j=0}^{\infty} \frac{\tilde{Y}^{(j)}(x)}{j!} h^{j}=\sum\limits_{j=0}^{\infty} a_{j} h^{j} $$

is $R=\min \limits \{1-x,1+x\}$ (which we could expect because the ODE satisfied by $\tilde {Y}$ has singularities at x = ± 1) and then

$$ \frac{1}{R}=\limsup_{n\rightarrow\infty}\sqrt[n]{|a_{n}|}=\max\{W_{1},W_{2}\}. $$

Therefore in this case the sequence {a_j}, $a_{j}=\tilde {Y}^{(j)}(x)/j!$ is in the dominant subspace of solutions of the recurrence.

The case when either α or β is odd but not both is different. Let us for instance consider that α is odd, but not β. In this case, the convergence of the series is limited by the singularity at − 1 (but not at + 1); the radius of convergence in this case is therefore R = 1 + x and then

$$ \limsup_{n\rightarrow\infty}\sqrt[n]{|a_{n}|}=\frac{1}{|1+x|}. $$

Therefore, {a_j} is in the dominant subspace only if x ∈ (− 1, 0). In this case for large enough j the forward computation of a_j would be unstable for positive x. However, even for this case we have not observed inaccuracies in the computation of the series. For a given accuracy claim, the number of terms in the series needed are not high enough to produce stability issues.