Cuadratura de Gauss-Legendre, cálculo de las abscisas y pesos

Question

Me gustaría escribir un programa para calcular las abscisas y los pesos de la cuadratura de Gauss-Legendre. He encontrado la siguiente fuente

http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature .

No entiendo por qué la primera suposición $x_0$ para el $i$ -raíz de a $n$ -polinomio de orden $P_n$ se define como

$x_0=\cos\left(\pi\frac{i-\frac{1}{4}}{n+\frac{1}{2}}\right)$

¿Podría alguien explicarme esto?

Answer 1

H. Bruns demostró en 1881 que la arcocosina del $i$ raíz de $P_n$ se encuentra en el intervalo $$ \left[ \frac{2i-1}{2n+1}\pi, \ \frac{2i}{2n+1}\pi\right] $$ (Cuando se trabaja con polinomios de Legendre, a menudo es conveniente componerlos con el coseno, es decir, mirar $P_n(\cos\theta)$ en lugar de $P_n(x)$ .)

Así que es natural tomar el coseno del punto medio de este intervalo como la primera conjetura para la raíz. Una fuente de lo anterior es El documento de Szegő mencionado por uranix; otro es el libro de Szegő Polinomios ortogonales Teorema 6.21.2.

Dicho esto, creo que utilizar un procedimiento de búsqueda de raíces para cada raíz individualmente no es particularmente eficiente, ya que esto no aprovecha la estructura especial de los polinomios $P_n$ . La propiedad de recurrencia de $P_n$ conduce a una familia de matrices tridiagonales para las que $P_n$ es (un múltiplo de) el polinomio característico. Es más fácil encontrar los valores propios de una matriz tridiagonal que las raíces de un polinomio genérico: véase el artículo Cálculo de las reglas de cuadratura de Gauss por G. H. Golub y J. H. Welsch.

Spoiler: el papel ¿Es mejor la cuadratura de Gauss que la de Clenshaw-Curtis? de Lloyd N. Trefethen tiene la siguiente implementación de 7 líneas de cuadratura gaussiana en $[-1,1]$ en Matlab.

function I = gauss(f,n)               % (n+1)-pt Gauss quadrature of f
beta = .5./sqrt(1-(2*(1:n)).ˆ(-2));   % 3-term recurrence coeffs
T = diag(beta,1) + diag(beta,-1);     % Jacobi matrix
[V,D] = eig(T);                       % eigenvalue decomposition
x = diag(D); [x,i] = sort(x);         % nodes (= Legendre points)
w = 2*V(1,i).ˆ2;                      % weights
I = w*feval(f,x);                     % the integral