1.15: Cuadratura Gaussiana - el Algoritmo

La cuadratura gaussiana es un método alternativo de integración numérica que a menudo es mucho más rápido y espectacular que la regla de Simpson. Cuadratura gaussiana permite llevar a cabo la integración

\[\int_{-1}^1 f(x) dx. \label{1.15.1} \tag{1.15.1}\]

Pero, ¿qué pasa si tus límites de integración no lo son\(\pm 1\)? ¿Qué pasa si quieres integrar

\[\int_a^b F(t) dt? \label{1.15.2} \tag{1.15.2}\]

Eso no es ningún problema en absoluto — solo haces un cambio de variable. Así, vamos

\[x = \frac{2t-a-b}{b-a} , \quad t= \frac{1}{2} [(b-a)x+a+b], \label{1.15.3} \tag{1.15.3}\]

y los nuevos límites son entonces\(x = \pm 1\).

A riesgo de ser poco sólido pedagógicamente voy a describir primero, sin ningún desarrollo teórico, justo lo que haces, con un ejemplo —siempre y cuando prometes mirar la derivación después, en la Sección 1.16.

Para nuestro ejemplo, intentemos evaluar

\[I = \int_0^{\pi/2} \sin \theta d \theta. \label{1.15.4} \tag{1.15.4}\]

Hagamos el cambio de variable dado por la ecuación\(\ref{1.15.3}\) (con\(t = \theta\),\(a = 0, \ b = \pi/2\)), y ahora tenemos que evaluar

\[I = \int_{-1}^1 \frac{\pi}{4} \sin \frac{\pi}{4} (x+1) dx . \label{1.15.5} \tag{1.15.5}\]

Para una cuadratura gaussiana de 5 puntos, se evalúa el integrando a cinco valores de\(x\), donde estos cinco valores de\(x\) son las soluciones\(P_5(x) = 0\) dadas en la Sección 1.14,\(P_5\) siendo el polinomio Legendre. Es decir, evaluamos el integrando en\(x = \pm 0.906 \ 469 \ 514 \ 203, \ \pm 0.538 \ 469 \ 310 \ 106\) y\(0\).

Ahora afirmo, sin derivación (hasta más tarde), que

\[I = \sum_{i=1}^5 c_{5,i} \ f(x_{5,i}), \label{1.15.6} \tag{1.15.6}\]

donde los coeficientes\(c_{l , i}\) (todos positivos) se enumeran con las raíces de los polinomios de Legendre en la Sección 1.14.

Vamos a probarlo.

\ begin {array} {c c c}
x_ {5, i} & f (x_ {5, i}) & c_ {5, i}\
\\
+0.906\ 179\ 845\ 939 & 0.783\ 266\ 908\ 39 & 0.236\ 926\ 885\ 06\
+0.538\ 469\ 310\ 106 & 0.734\ 361\ 739\ 69 & 0.478\ 628\ 670\ 50\\
0 .000\ 000\ 000\ 000 & 0.555\ 360\ 367\ 27 & 0.568\ 888\ 888\ 89\\
-0.538\ 469\ 310\ 006 & 0.278\ 501\ 544\ 60 & 0.478\ 628\ 670\ 50\
-0.906\ 179\ 845\ 939 & 0.057\ 820\ 630\ 35 & 0.236\ 926\ 885\ 06\
\ end {array}

y la expresión\(\ref{1.15.6}\) llega a\(1.000 \ 000 \ 000 \ 04\), y presumiblemente podría haberse acercado aún más a 1 si hubiéramos dado\(x_{l,i}\), y\(c_{l,i}\), a cifras más significativas.

Ahora deberías escribir un programa de computadora para cuadratura gaussiana — tendrás que almacenar el\(x_{l,i}\) y por\(c_{l,i}\) supuesto. Presumiblemente ya has escrito un programa para la regla de Simpson.

En un texto sobre integración, el autor invitó al lector a evaluar las siguientes integrales por cuadratura gaussiana:

\ begin {array} {c c c c}
(a) &\ int_1^ {1.5} x^2\ ln x dx & & & (e) &\ int_0^ {\ pi/4} e^ {3x}\ sin 2x dx\
(b) &\ int_0^1 x^2 x^ {-x} dx && (f) &\ int_1^ {1.6} frac {2x} {x^2 - 4} dx\\
(c) &\ int_0^ {0.35}\ frac {2} {x^2-4} dx & & (g) &\ int_3^ {3.5}\ frac {x} {\ sqrt {x^2 - 4}} dx\\
(d) &\ int_0^ {\ pi/4} x^2\ sin x dx && (h) &\ int_0^ {\ pi/4}\ cos^2 x dx\
\ end {array}

Todos estos pueden integrarse analíticamente, así que voy a invitar al lector a evaluarlos primero analíticamente, y luego numéricamente por la regla de Simpson y nuevamente por cuadratura gaussiana, y para ver en cuántos puntos el integrando tiene que ser evaluado por cada método para lograr una precisión de nueve o diez cifras. Lo intenté, y los resultados son los siguientes. La primera columna es la respuesta, la segunda columna es el número de puntos requeridos por la regla de Simpson, y la tercera columna es el número de puntos requeridos por la cuadratura gaussiana.

\ begin {array} {c c c c c}
(a) & 0.192\ 259\ 358 && 33 & 4\\
(b) & 0.160\ 602\ 794 && 99 & 5\\
(c) & −0.176\ 820\ 020 && 19 & 4\\
(d) & 0.088\ 755\ 284\ 4 && 111 & 5\\
(e) & 2.588\ 628\ 633 && 453 & 7\\
(f) & -0.733\ 969\ 175 && 143 & 8\\
(g) & 0.636\ 213\ 346 && 31 & 5\\
(h) & 0.642\ 699\ 082 && 59 & 5\\
\ end {array}

Echemos ahora un vistazo a cuatro de las integrales que conocimos en la Sección 1.2.

1. \(\int_0^1 \frac{x^4 dx}{\sqrt{2(1+x^2)}}\). Esto fue sencillo. Cuenta con una solución analítica de\(\frac{\sqrt{18} \ln (1 + \sqrt{2}) -2 }{16} = 0.108 \ 709 \ 465 \). Necesitaba evaluar la integral en 89 puntos para poder obtener esta respuesta a nueve figuras significativas usando la regla de Simpson. Para usar cuadratura gaussiana, observamos que integrand contiene solo potencias pares de\(x\) y por lo tanto es simétrico alrededor\(x = 0\), y por lo tanto la integral es igual a\(\frac{1}{2} \int_{-1}^1 \frac{x^4 dx}{\sqrt{2(1+x^2)}}\), ¡lo que la hace inmediatamente conveniente para cuadratura gaussiana! Doy a continuación las respuestas que obtuve para cuadratura gaussiana de 3 a 7 puntos.

\ begin {array} {c c c}
& 3 & 0.108\ 667\ 036\\
& 4 & 0.108\ 711\ 215\
& 5 & 0.108\ 709\ 441\
& 6 & 0.108\ 709\ 463\
& 7 & 0.108\ 709\ 465\
\ text {Respuesta correcta} & & 0.108\ 709\ 465\\
\ end {array}

2. \(\int_0^2 \frac{y^2 dy}{\sqrt{2-y}}\). Esto tenía la dificultad de que el integrando es infinito en el límite superior. Nos volvimos a esto por medio de la sustitución\(y = 2 \sin^2 \theta\), y la integral se vuelve\(\sqrt{128} \int_0^{\pi/2} \sin^5 \theta d \theta\). Esto tiene una solución analítica de\(\sqrt{8192}/15 = .6 033977866\). Necesitaba 59 puntos para obtener esta respuesta a diez figuras significativas usando la regla de Simpson. Para usar cuadratura gaussiana podemos dejar\(y = 1 + x\), para que la integral se convierta,\(\int_{-1}^1 \frac{(1+x)^2 dy}{\sqrt{1-x}},\) lo que parece ser inmediatamente adecuado para cuadratura gaussiana. Antes de continuar, recordamos que el integrando se vuelve infinito en el límite superior, y todavía lo hace después de nuestro cambio de variable. Observamos, sin embargo, que con la cuadratura gaussiana, no evaluamos el integrando en el límite superior, por lo que esto parecería ser una gran ventaja del método sobre el método de Simpson. ¡Ay! — esto resulta no ser el caso. Si, por ejemplo, utilizamos una cuadratura de 17 puntos, el mayor valor\(x\) para el cual evaluamos el integrando es igual a la solución más grande de\(P_{17} (x) = 0\), que es\(0.9906\). Simplemente no podemos pasar por alto el hecho de que el integrando dispara hasta el infinito más allá de esto, así que hemos dejado atrás una gran parte de la integral. En efecto, con una cuadratura gaussiana de 17 puntos, obtuve una respuesta de\(5.75\), que está muy lejos de la respuesta correcta de\(6.03\).

Por lo tanto tenemos que hacer un cambio de variable, como hicimos para el método de Simpson, para que el límite superior sea finito. Elegimos\(y = 2 \sin^2 \theta\) cuál cambió la integral a\(\sqrt{128} \int_0^{\pi/2} \sin^5 \theta d \theta\). Para que esto sea adecuado para cuadratura gaussiana, ahora debemos hacer la sustitución adicional (ver Ecuación\(\ref{1.15.3}\))\(x = 4\theta/\pi-1\),\(\theta = \frac{\pi}{4} (x+1)\). Si queremos impresionar, podemos hacer las dos sustituciones en un solo paso, así: Vamos\(y = 2 \sin^2 \frac{\pi}{4}(1+x)\),\(x = \frac{4}{\pi} \sin^{-1} \sqrt{\frac{y}{2}} - 1\). La integral se vuelve\(\sqrt{8} \pi \int_{-1}^1 \sin^5 \frac{\pi}{4}(1+x)dx\), y no hay más dificultades. Con una integración de 9 puntos, obtuve la respuesta, correcta a diez cifras significativas,\(6.033 \ 977 \ 866\). La regla de Simpson requirió 59 puntos.

3. \(\int_0^{\pi/2} \sqrt{\sec \theta} d \theta\). Esta integral ocurre en la teoría de un simple péndulo que se balancea a través\(90^\circ\). Por lo que puedo decir no tiene una solución analítica simple a menos que tengamos que recurrir a integrales elípticas desconocidas, que tendríamos que evaluar numéricamente en cualquier caso. La integral tiene la dificultad de que el integrando es infinito en el límite superior. Obtenemos alrededor de esto por medio de una sustitución. Así vamos\(\sin \phi = \sqrt{2} \sin \frac{1}{2} \theta\). (¿No pensaste en esto?) La integral se vuelve\(\sqrt{2} \int_0^{\pi/2} \frac{d\phi}{\sqrt{1 - \frac{1}{2} \sin^2 \phi}}.\) necesaria 13 puntos por regla de Simpson para obtener la respuesta a diez cifras significativas,\(2.622 \ 057 \ 554\). Para hacer los límites\(\pm 1\), adecuados para cuadratura gaussiana, podemos hacer la segunda sustitución (como en el ejemplo 2),\(\phi = \frac{\pi}{4}(x+1)\). Si realmente deseamos impresionar a nuestros amigos, podemos hacer las dos sustituciones en un solo paso, así: Vamos\(\sin \frac{\pi}{4} (1+x) = \sqrt{2} \sin \frac{1}{2} \theta\). (¡Nadie va a adivinar cómo pensamos de eso!) La integral se convierte en la\(\frac{\pi}{2} \int_{-1}^1 \frac{dx}{\sqrt{2 - \sin^2 \frac{\pi}{4} (x+1)}},\) que ahora está lista para la cuadratura gaussiana. Obtuve la respuesta\(2.622 \ 057 \ 554\) en una cuadratura gaussiana de 10 puntos, que sólo es un poco más rápida que los 13 puntos requeridos por la regla de Simpson.

4. \(\int_0^\infty \frac{dy}{y^5\left(e^{1/y} - 1 \right)}.\)Esta integral ocurre en la teoría de la radiación de cuerpo negro. Tiene la dificultad de un límite superior infinito. Obtenemos alrededor de esto por medio de una sustitución. Así vamos\(y = \tan \theta\). La integral se convierte\(\int_0^{\pi/2} \frac{c^3(c^2+1)}{e^c - 1} d\theta\), donde\(c = \cot \theta\). Cuenta con una solución analítica de\(\pi^4 /15 = 6.493 \ 939 \ 402\). Necesitaba 261 puntos por regla de Simpson para obtener la respuesta a diez cifras significativas. Para prepararlo para cuadratura gaussiana, podemos dejar\(\theta = \frac{\pi}{4}(x+1)\), como hicimos en el ejemplo 2, para que la integral se convierta, donde\(\frac{\pi}{4} \int_{-1}^1 \frac{c^3 (c^2+1)}{e^c -1}dx,\) donde\(c = \cot \frac{\pi}{4} (x+1)\). Usando cuadratura gaussiana de 16 puntos, obtuve 6.48. Por lo tanto, necesitaríamos extender nuestra tabla de constantes para el método gaussiano a un orden mucho más alto para poder utilizar el método con éxito. Sin duda, el método gaussiano sería entonces más rápido que el método Simpson, pero no necesitamos un conjunto extenso (y difícil de calcular) de constantes para este último. Otro pequeño punto: Es posible que hayas notado que no es inmediatamente obvio que el integrando sea cero en los puntos finales, y que se necesita algún trabajo para demostrarlo. Pero con el método gaussiano no se evalúa el integrando en los puntos finales, ¡así que eso es una cosa menos de lo que preocuparse!

Así hemos encontrado que en la mayoría de los casos el método gaussiano es mucho más rápido que el método Simpson. En algunos casos sólo es marginalmente más rápido. En otros más probablemente sería más rápido que la regla de Simpson, pero se necesitan constantes de orden superior para aplicarla. Ya sea que usemos la regla de Simpson o la cuadratura gaussiana, tenemos que llevar a cabo la integración con órdenes sucesivamente superiores hasta que ir a órdenes superiores no dé como resultado ningún cambio adicional en el número de cifras significativas deseadas.