Cuadratura gaussiana

Aproximación de la integral definida de una función

Comparación entre la cuadratura gaussiana de 2 puntos y la cuadratura trapezoidal.
La curva azul muestra la función cuya integral definida en el intervalo [−1, 1] se debe calcular (el integrando). La regla trapezoidal aproxima la función con una función lineal que coincide con el integrando en los puntos finales del intervalo y se representa mediante una línea discontinua naranja. La aproximación aparentemente no es buena, por lo que el error es grande (la regla trapezoidal da una aproximación de la integral igual a y (–1) + y (1) = –10 , mientras que el valor correcto es 2 ⁄ 3 ). Para obtener un resultado más preciso, el intervalo debe dividirse en muchos subintervalos y luego debe usarse la regla trapezoidal compuesta
, que requiere muchos más cálculos. La cuadratura gaussiana elige puntos más adecuados, por lo que incluso una función lineal aproxima mejor la función (la línea discontinua negra). Como el integrando es el polinomio de grado 3 ( y ( x ) = 7 x ³ – 8 x ² – 3 x + 3 ), la regla de cuadratura gaussiana de 2 puntos incluso devuelve un resultado exacto.

En el análisis numérico , una regla de cuadratura gaussiana de n puntos , llamada así en honor a Carl Friedrich Gauss , ^[1] es una regla de cuadratura construida para producir un resultado exacto para polinomios de grado 2 n − 1 o menos mediante una elección adecuada de los nodos x _i y pesos w _i para i = 1, ..., n .

La formulación moderna que utiliza polinomios ortogonales fue desarrollada por Carl Gustav Jacobi en 1826. ^[2] El dominio de integración más común para dicha regla se toma como [−1, 1] , por lo que la regla se enuncia como $${\displaystyle \int _{-1}^{1}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i}),}$$

que es exacta para polinomios de grado 2 n − 1 o menos. Esta regla exacta se conoce como regla de cuadratura de Gauss-Legendre . La regla de cuadratura solo será una aproximación precisa a la integral anterior si f ( x ) se aproxima bien mediante un polinomio de grado 2 n − 1 o menos en [−1, 1] .

La regla de cuadratura de Gauss- Legendre no se suele utilizar para funciones integrables con singularidades en los puntos finales . En cambio, si el integrando se puede escribir como

$${\displaystyle f(x)=\left(1-x\right)^{\alpha }\left(1+x\right)^{\beta }g(x),\quad \alpha ,\beta >-1,}$$

donde g ( x ) se aproxima bien mediante un polinomio de bajo grado, entonces los nodos alternativos x _i ' y los pesos w _i ' generalmente darán reglas de cuadratura más precisas. Estas se conocen como reglas de cuadratura de Gauss-Jacobi , es decir,

$${\displaystyle \int _{-1}^{1}f(x)\,dx=\int _{-1}^{1}\left(1-x\right)^{\alpha }\left(1+x\right)^{\beta }g(x)\,dx\approx \sum _{i=1}^{n}w_{i}'g\left(x_{i}'\right).}$$

Los pesos comunes incluyen ( Chebyshev–Gauss ) y . También se puede querer integrar sobre intervalos semi-infinitos ( cuadratura de Gauss–Laguerre ) e infinitos ( cuadratura de Gauss–Hermite ). ${\textstyle {\frac {1}{\sqrt {1-x^{2}}}}}$ ${\textstyle {\sqrt {1-x^{2}}}}$

Se puede demostrar (véase Press et al., o Stoer y Bulirsch) que los nodos de cuadratura x _i son las raíces de un polinomio que pertenece a una clase de polinomios ortogonales (la clase ortogonal con respecto a un producto interno ponderado). Esta es una observación clave para calcular los nodos de cuadratura y los pesos de Gauss.

Cuadratura de Gauss-Legendre

Gráficas de polinomios de Legendre (hasta n = 5)

Para el problema de integración más simple indicado anteriormente, es decir, f ( x ) se aproxima bien mediante polinomios en , los polinomios ortogonales asociados son polinomios de Legendre , denotados por P _n ( x ) . Con el n -ésimo polinomio normalizado para dar P _n (1) = 1 , el i -ésimo nodo de Gauss, x _i , es la i -ésima raíz de P _n y los pesos se dan mediante la fórmula ^[3] ${\displaystyle [-1,1]}$ $${\displaystyle w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)\left[P'_{n}(x_{i})\right]^{2}}}.}$$

A continuación se tabulan algunas reglas de cuadratura de orden bajo (en el intervalo [−1, 1] ; consulte la sección a continuación para otros intervalos).

Cambio de intervalo

Una integral sobre [ a , b ] debe transformarse en una integral sobre [−1, 1] antes de aplicar la regla de cuadratura de Gauss. Este cambio de intervalo se puede realizar de la siguiente manera: $${\displaystyle \int _{a}^{b}f(x)\,dx=\int _{-1}^{1}f\left({\frac {b-a}{2}}\xi +{\frac {a+b}{2}}\right)\,{\frac {dx}{d\xi }}d\xi }$$

con ${\displaystyle {\frac {dx}{d\xi }}={\frac {b-a}{2}}}$

La aplicación de la regla de cuadratura gaussiana de puntos da como resultado la siguiente aproximación: ${\displaystyle n}$ ${\displaystyle (\xi ,w)}$ $${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{2}}\sum _{i=1}^{n}w_{i}f\left({\frac {b-a}{2}}\xi _{i}+{\frac {a+b}{2}}\right).}$$

Ejemplo de regla de cuadratura de Gauss de dos puntos

Utilice la regla de cuadratura de Gauss de dos puntos para aproximar la distancia en metros recorrida por un cohete desde a como se indica en ${\displaystyle t=8\mathrm {s} }$ ${\displaystyle t=30\mathrm {s} ,}$ $${\displaystyle s=\int _{8}^{30}{\left(2000\ln \left[{\frac {140000}{140000-2100t}}\right]-9.8t\right){dt}}}$$

Cambie los límites de modo que se puedan utilizar los pesos y las abscisas que se indican en la Tabla 1. Además, encuentre el error verdadero relativo absoluto. El valor verdadero se indica como 11061,34 m.

Solución

En primer lugar, al cambiar los límites de integración de a se obtiene ${\displaystyle \left[8,30\right]}$ ${\displaystyle \left[-1,1\right]}$

$${\displaystyle {\begin{aligned}\int _{8}^{30}{f(t)dt}&={\frac {30-8}{2}}\int _{-1}^{1}{f\left({\frac {30-8}{2}}x+{\frac {30+8}{2}}\right){dx}}\\&=11\int _{-1}^{1}{f\left(11x+19\right){dx}}\end{aligned}}}$$

A continuación, obtenga los factores de ponderación y los valores de los argumentos de función de la Tabla 1 para la regla de dos puntos.

• ${\displaystyle c_{1}=1.000000000}$
• ${\displaystyle x_{1}=-0.577350269}$
• ${\displaystyle c_{2}=1.000000000}$
• ${\displaystyle x_{2}=0.577350269}$

Ahora podemos utilizar la fórmula de cuadratura de Gauss ya que $${\displaystyle {\begin{aligned}11\int _{-1}^{1}{f\left(11x+19\right){dx}}&\approx 11\left[c_{1}f\left(11x_{1}+19\right)+c_{2}f\left(11x_{2}+19\right)\right]\\&=11\left[f\left(11(-0.5773503)+19\right)+f\left(11(0.5773503)+19\right)\right]\\&=11\left[f(12.64915)+f(25.35085)\right]\\&=11\left[(296.8317)+(708.4811)\right]\\&=11058.44\end{aligned}}}$$ $${\displaystyle {\begin{aligned}f(12.64915)&=2000\ln \left[{\frac {140000}{140000-2100(12.64915)}}\right]-9.8(12.64915)\\&=296.8317\end{aligned}}}$$ $${\displaystyle {\begin{aligned}f(25.35085)&=2000\ln \left[{\frac {140000}{140000-2100(25.35085)}}\right]-9.8(25.35085)\\&=708.4811\end{aligned}}}$$

Dado que el valor real es 11061,34 m, el error verdadero relativo absoluto es ${\displaystyle \left|\varepsilon _{t}\right|}$ $${\displaystyle \left|\varepsilon _{t}\right|=\left|{\frac {11061.34-11058.44}{11061.34}}\right|\times 100\%=0.0262\%}$$

Otras formas

El problema de integración se puede expresar de una manera ligeramente más general introduciendo una función de peso positiva ω en el integrando y permitiendo un intervalo distinto de [−1, 1] . Es decir, el problema consiste en calcular algunas opciones de a , b y ω . Para a = −1 , b = 1 y ω ( x ) = 1 , el problema es el mismo que el considerado anteriormente. Otras opciones conducen a otras reglas de integración. Algunas de ellas se tabulan a continuación. Se dan los números de ecuación para Abramowitz y Stegun (A y S). $${\displaystyle \int _{a}^{b}\omega (x)\,f(x)\,dx}$$

Teorema fundamental

Sea p _n un polinomio no trivial de grado n tal que $${\displaystyle \int _{a}^{b}\omega (x)\,x^{k}p_{n}(x)\,dx=0,\quad {\text{for all }}k=0,1,\ldots ,n-1.}$$

Nótese que esto será cierto para todos los polinomios ortogonales anteriores, porque cada p _n está construido para ser ortogonal a los otros polinomios p _j para j < n , y x ^k está en el espacio de ese conjunto.

Si elegimos los n nodos x _i como los ceros de p _n , entonces existen n pesos w _i que hacen que la integral calculada en cuadratura gaussiana sea exacta para todos los polinomios h ( x ) de grado 2 n − 1 o menor. Además, todos estos nodos x _i estarán en el intervalo abierto ( a , b ) . ^[4]

Para demostrar la primera parte de esta afirmación, sea h ( x ) un polinomio cualquiera de grado 2 n − 1 o menor. Dividúzcalo por el polinomio ortogonal p _n para obtener donde q ( x ) es el cociente, de grado n − 1 o menor (porque la suma de su grado y la del divisor p _n debe ser igual a la del dividendo), y r ( x ) es el resto, también de grado n − 1 o menor (porque el grado del resto es siempre menor que el del divisor). Como p _n es por supuesto ortogonal a todos los monomios de grado menor que n , debe ser ortogonal al cociente q ( x ) . Por lo tanto $${\displaystyle h(x)=p_{n}(x)\,q(x)+r(x).}$$ $${\displaystyle \int _{a}^{b}\omega (x)\,h(x)\,dx=\int _{a}^{b}\omega (x)\,{\big (}\,p_{n}(x)q(x)+r(x)\,{\big )}\,dx=\int _{a}^{b}\omega (x)\,r(x)\,dx.}$$

Como el resto r ( x ) es de grado n − 1 o menor, podemos interpolarlo exactamente usando n puntos de interpolación con polinomios de Lagrange l _i ( x ) , donde $${\displaystyle l_{i}(x)=\prod _{j\neq i}{\frac {x-x_{j}}{x_{i}-x_{j}}}.}$$

Tenemos $${\displaystyle r(x)=\sum _{i=1}^{n}l_{i}(x)\,r(x_{i}).}$$

Entonces su integral será igual $${\displaystyle \int _{a}^{b}\omega (x)\,r(x)\,dx=\int _{a}^{b}\omega (x)\,\sum _{i=1}^{n}l_{i}(x)\,r(x_{i})\,dx=\sum _{i=1}^{n}\,r(x_{i})\,\int _{a}^{b}\omega (x)\,l_{i}(x)\,dx=\sum _{i=1}^{n}\,r(x_{i})\,w_{i},}$$

donde w _i , el peso asociado con el nodo x _i , se define como igual a la integral ponderada de l _i ( x ) (ver más abajo otras fórmulas para los pesos). Pero todos los x _i son raíces de p _n , por lo que la fórmula de división anterior nos dice que para todos los i . Por lo tanto, finalmente tenemos $${\displaystyle h(x_{i})=p_{n}(x_{i})\,q(x_{i})+r(x_{i})=r(x_{i}),}$$ $${\displaystyle \int _{a}^{b}\omega (x)\,h(x)\,dx=\int _{a}^{b}\omega (x)\,r(x)\,dx=\sum _{i=1}^{n}w_{i}\,r(x_{i})=\sum _{i=1}^{n}w_{i}\,h(x_{i}).}$$

Esto demuestra que para cualquier polinomio h ( x ) de grado 2 n − 1 o menor, su integral está dada exactamente por la suma de la cuadratura gaussiana.

Para demostrar la segunda parte de la afirmación, considere la forma factorizada del polinomio p _n . Cualquier raíz conjugada compleja producirá un factor cuadrático que es estrictamente positivo o estrictamente negativo en toda la línea real. Cualquier factor para raíces fuera del intervalo de a a b no cambiará de signo en ese intervalo. Finalmente, para los factores correspondientes a raíces x _i dentro del intervalo de a a b que sean de multiplicidad impar, multiplique p _n por un factor más para formar un nuevo polinomio. $${\displaystyle p_{n}(x)\,\prod _{i}(x-x_{i}).}$$

Este polinomio no puede cambiar de signo en el intervalo de a a b porque todas sus raíces son ahora de multiplicidad par. Por lo tanto, la integral, ya que la función de peso ω ( x ), es siempre no negativa. Pero p _n es ortogonal a todos los polinomios de grado n -1 o menor, por lo que el grado del producto debe ser al menos n . Por lo tanto, p _n tiene n raíces distintas, todas reales, en el intervalo de a a b . $${\displaystyle \int _{a}^{b}p_{n}(x)\,\left(\prod _{i}(x-x_{i})\right)\,\omega (x)\,dx\neq 0,}$$ $${\displaystyle \prod _{i}(x-x_{i})}$$

Fórmula general para los pesos

Los pesos se pueden expresar como

donde es el coeficiente de en . Para demostrarlo, observe que utilizando la interpolación de Lagrange se puede expresar r ( x ) en términos de como porque r ( x ) tiene grado menor que n y, por lo tanto, está fijado por los valores que alcanza en n puntos diferentes. Multiplicando ambos lados por ω ( x ) e integrando de a a b se obtiene ${\displaystyle a_{k}}$ ${\displaystyle x^{k}}$ ${\displaystyle p_{k}(x)}$ ${\displaystyle r(x_{i})}$ $${\displaystyle r(x)=\sum _{i=1}^{n}r(x_{i})\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {x-x_{j}}{x_{i}-x_{j}}}}$$ $${\displaystyle \int _{a}^{b}\omega (x)r(x)dx=\sum _{i=1}^{n}r(x_{i})\int _{a}^{b}\omega (x)\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {x-x_{j}}{x_{i}-x_{j}}}dx}$$

Los pesos w _i vienen dados por $${\displaystyle w_{i}=\int _{a}^{b}\omega (x)\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {x-x_{j}}{x_{i}-x_{j}}}dx}$$

Esta expresión integral para se puede expresar en términos de los polinomios ortogonales y de la siguiente manera. ${\displaystyle w_{i}}$ ${\displaystyle p_{n}(x)}$ ${\displaystyle p_{n-1}(x)}$

Podemos escribir $${\displaystyle \prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}\left(x-x_{j}\right)={\frac {\prod _{1\leq j\leq n}\left(x-x_{j}\right)}{x-x_{i}}}={\frac {p_{n}(x)}{a_{n}\left(x-x_{i}\right)}}}$$

donde es el coeficiente de en . Tomando el límite de x a obtenemos usando la regla de L'Hôpital ${\displaystyle a_{n}}$ ${\displaystyle x^{n}}$ ${\displaystyle p_{n}(x)}$ ${\displaystyle x_{i}}$ $${\displaystyle \prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}\left(x_{i}-x_{j}\right)={\frac {p'_{n}(x_{i})}{a_{n}}}}$$

Podemos entonces escribir la expresión integral para los pesos como

En el integrando, escribiendo $${\displaystyle {\frac {1}{x-x_{i}}}={\frac {1-\left({\frac {x}{x_{i}}}\right)^{k}}{x-x_{i}}}+\left({\frac {x}{x_{i}}}\right)^{k}{\frac {1}{x-x_{i}}}}$$

rendimientos $${\displaystyle \int _{a}^{b}\omega (x){\frac {x^{k}p_{n}(x)}{x-x_{i}}}dx=x_{i}^{k}\int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx}$$

siempre que , porque es un polinomio de grado k − 1 que es entonces ortogonal a . Por lo tanto, si q ( x ) es un polinomio de grado n como máximo tenemos ${\displaystyle k\leq n}$ $${\displaystyle {\frac {1-\left({\frac {x}{x_{i}}}\right)^{k}}{x-x_{i}}}}$$ ${\displaystyle p_{n}(x)}$ $${\displaystyle \int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx={\frac {1}{q(x_{i})}}\int _{a}^{b}\omega (x){\frac {q(x)p_{n}(x)}{x-x_{i}}}dx}$$

Podemos evaluar la integral del lado derecho de la siguiente manera. Como es un polinomio de grado n − 1 , tenemos donde s ( x ) es un polinomio de grado . Como s ( x ) es ortogonal a tenemos ${\displaystyle q(x)=p_{n-1}(x)}$ ${\displaystyle {\frac {p_{n}(x)}{x-x_{i}}}}$ $${\displaystyle {\frac {p_{n}(x)}{x-x_{i}}}=a_{n}x^{n-1}+s(x)}$$ ${\displaystyle n-2}$ ${\displaystyle p_{n-1}(x)}$ $${\displaystyle \int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx={\frac {a_{n}}{p_{n-1}(x_{i})}}\int _{a}^{b}\omega (x)p_{n-1}(x)x^{n-1}dx}$$

Entonces podemos escribir $${\displaystyle x^{n-1}=\left(x^{n-1}-{\frac {p_{n-1}(x)}{a_{n-1}}}\right)+{\frac {p_{n-1}(x)}{a_{n-1}}}}$$

El término entre paréntesis es un polinomio de grado , que por lo tanto es ortogonal a . La integral puede escribirse así: ${\displaystyle n-2}$ ${\displaystyle p_{n-1}(x)}$ $${\displaystyle \int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx={\frac {a_{n}}{a_{n-1}p_{n-1}(x_{i})}}\int _{a}^{b}\omega (x)p_{n-1}(x)^{2}dx}$$

De acuerdo con la ecuación ( 2 ), los pesos se obtienen dividiendo esto por y que da como resultado la expresión de la ecuación ( 1 ). ${\displaystyle p'_{n}(x_{i})}$

${\displaystyle w_{i}}$también se puede expresar en términos de los polinomios ortogonales y ahora . En la relación de recurrencia de 3 términos , el término con se anula, por lo que en la ecuación (1) se puede reemplazar por . ${\displaystyle p_{n}(x)}$ ${\displaystyle p_{n+1}(x)}$ ${\displaystyle p_{n+1}(x_{i})=(a)p_{n}(x_{i})+(b)p_{n-1}(x_{i})}$ ${\displaystyle p_{n}(x_{i})}$ ${\displaystyle p_{n-1}(x_{i})}$ ${\textstyle {\frac {1}{b}}p_{n+1}\left(x_{i}\right)}$

Prueba de que los pesos son positivos

Considere el siguiente polinomio de grado donde, como se indicó anteriormente, las x _j son las raíces del polinomio . Claramente . Dado que el grado de es menor que , se aplica la fórmula de cuadratura gaussiana que involucra los pesos y nodos obtenidos de . Dado que para j no es igual a i , tenemos ${\displaystyle 2n-2}$ $${\displaystyle f(x)=\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {\left(x-x_{j}\right)^{2}}{\left(x_{i}-x_{j}\right)^{2}}}}$$ ${\displaystyle p_{n}(x)}$ ${\displaystyle f(x_{j})=\delta _{ij}}$ ${\displaystyle f(x)}$ ${\displaystyle 2n-1}$ ${\displaystyle p_{n}(x)}$ ${\displaystyle f(x_{j})=0}$ $${\displaystyle \int _{a}^{b}\omega (x)f(x)dx=\sum _{j=1}^{n}w_{j}f(x_{j})=\sum _{j=1}^{n}\delta _{ij}w_{j}=w_{i}>0.}$$

Como tanto y son funciones no negativas, se deduce que . ${\displaystyle \omega (x)}$ ${\displaystyle f(x)}$ ${\displaystyle w_{i}>0}$

Cálculo de las reglas de cuadratura gaussiana

Existen muchos algoritmos para calcular los nodos x _i y los pesos w _i de las reglas de cuadratura gaussiana. Los más populares son el algoritmo de Golub-Welsch que requiere O ( n ² ) operaciones, el método de Newton para resolver utilizando la recurrencia de tres términos para la evaluación que requiere O ( n ² ) operaciones y fórmulas asintóticas para n grandes que requieren O ( n ) operaciones. ${\displaystyle p_{n}(x)=0}$

Relación de recurrencia

Los polinomios ortogonales con para un producto escalar , grado y coeficiente principal uno (es decir, polinomios ortogonales mónicos ) satisfacen la relación de recurrencia ${\displaystyle p_{r}}$ ${\displaystyle (p_{r},p_{s})=0}$ ${\displaystyle r\neq s}$ ${\displaystyle (\cdot ,\cdot )}$ ${\displaystyle (p_{r})=r}$ $${\displaystyle p_{r+1}(x)=(x-a_{r,r})p_{r}(x)-a_{r,r-1}p_{r-1}(x)\cdots -a_{r,0}p_{0}(x)}$$

y producto escalar definido $${\displaystyle (f(x),g(x))=\int _{a}^{b}\omega (x)f(x)g(x)dx}$$

donde n es el grado máximo que puede tomarse como infinito, y donde . En primer lugar, los polinomios definidos por la relación de recurrencia que comienza con tienen coeficiente principal uno y grado correcto. Dado el punto de partida por , la ortogonalidad de se puede demostrar por inducción. Para uno tiene ${\displaystyle r=0,1,\ldots ,n-1}$ ${\textstyle a_{r,s}={\frac {\left(xp_{r},p_{s}\right)}{\left(p_{s},p_{s}\right)}}}$ ${\displaystyle p_{0}(x)=1}$ ${\displaystyle p_{0}}$ ${\displaystyle p_{r}}$ ${\displaystyle r=s=0}$ $${\displaystyle (p_{1},p_{0})=(x-a_{0,0})(p_{0},p_{0})=(xp_{0},p_{0})-a_{0,0}(p_{0},p_{0})=(xp_{0},p_{0})-(xp_{0},p_{0})=0.}$$

Ahora bien, si son ortogonales, entonces también , porque en todos los productos escalares se anulan excepto en el primero y en el que cumple el mismo polinomio ortogonal. Por lo tanto, ${\displaystyle p_{0},p_{1},\ldots ,p_{r}}$ ${\displaystyle p_{r+1}}$ $${\displaystyle (p_{r+1},p_{s})=(xp_{r},p_{s})-a_{r,r}(p_{r},p_{s})-a_{r,r-1}(p_{r-1},p_{s})\cdots -a_{r,0}(p_{0},p_{s})}$$ ${\displaystyle p_{s}}$ $${\displaystyle (p_{r+1},p_{s})=(xp_{r},p_{s})-a_{r,s}(p_{s},p_{s})=(xp_{r},p_{s})-(xp_{r},p_{s})=0.}$$

Sin embargo, si el producto escalar satisface (que es el caso de la cuadratura gaussiana), la relación de recurrencia se reduce a una relación de recurrencia de tres términos: Para es un polinomio de grado menor o igual a r − 1 . Por otro lado, es ortogonal a todo polinomio de grado menor o igual a r − 1 . Por lo tanto, se tiene y para s < r − 1 . La relación de recurrencia entonces se simplifica a ${\displaystyle (xf,g)=(f,xg)}$ ${\displaystyle s<r-1,xp_{s}}$ ${\displaystyle p_{r}}$ ${\displaystyle (xp_{r},p_{s})=(p_{r},xp_{s})=0}$ ${\displaystyle a_{r,s}=0}$ $${\displaystyle p_{r+1}(x)=(x-a_{r,r})p_{r}(x)-a_{r,r-1}p_{r-1}(x)}$$

o $${\displaystyle p_{r+1}(x)=(x-a_{r})p_{r}(x)-b_{r}p_{r-1}(x)}$$

(con la convención ) donde ${\displaystyle p_{-1}(x)\equiv 0}$ $${\displaystyle a_{r}:={\frac {(xp_{r},p_{r})}{(p_{r},p_{r})}},\qquad b_{r}:={\frac {(xp_{r},p_{r-1})}{(p_{r-1},p_{r-1})}}={\frac {(p_{r},p_{r})}{(p_{r-1},p_{r-1})}}}$$

(el último debido a , ya que difiere de en un grado menor que r ). ${\displaystyle (xp_{r},p_{r-1})=(p_{r},xp_{r-1})=(p_{r},p_{r})}$ ${\displaystyle xp_{r-1}}$ ${\displaystyle p_{r}}$

El algoritmo de Golub-Welsch

La relación de recurrencia de tres términos se puede escribir en forma matricial donde , es el vector base estándar, es decir, , y J es la siguiente matriz tridiagonal , llamada matriz de Jacobi: ${\displaystyle J{\tilde {P}}=x{\tilde {P}}-p_{n}(x)\mathbf {e} _{n}}$ ${\displaystyle {\tilde {P}}={\begin{bmatrix}p_{0}(x)&p_{1}(x)&\cdots &p_{n-1}(x)\end{bmatrix}}^{\mathsf {T}}}$ ${\displaystyle \mathbf {e} _{n}}$ ${\displaystyle n}$ ${\displaystyle \mathbf {e} _{n}={\begin{bmatrix}0&\cdots &0&1\end{bmatrix}}^{\mathsf {T}}}$ $${\displaystyle \mathbf {J} ={\begin{bmatrix}a_{0}&1&0&\cdots &0\\b_{1}&a_{1}&1&\ddots &\vdots \\0&b_{2}&\ddots &\ddots &0\\\vdots &\ddots &\ddots &a_{n-2}&1\\0&\cdots &0&b_{n-1}&a_{n-1}\end{bmatrix}}.}$$

Los ceros de los polinomios de grado n como nodos de la cuadratura gaussiana se pueden hallar calculando los valores propios de esta matriz. Este procedimiento se conoce como algoritmo de Golub-Welsch . ${\displaystyle x_{j}}$

Para calcular los pesos y nodos, es preferible considerar la matriz tridiagonal simétrica con elementos ${\displaystyle {\mathcal {J}}}$ $${\displaystyle {\begin{aligned}{\mathcal {J}}_{k,i}=J_{k,i}&=a_{k-1}&k&=1,2,\ldots ,n\\[2.1ex]{\mathcal {J}}_{k-1,i}={\mathcal {J}}_{k,k-1}={\sqrt {J_{k,k-1}J_{k-1,k}}}&={\sqrt {b_{k-1}}}&k&={\hphantom {1,\,}}2,\ldots ,n.\end{aligned}}}$$

Eso es,

$${\displaystyle {\mathcal {J}}={\begin{bmatrix}a_{0}&{\sqrt {b_{1}}}&0&\cdots &0\\{\sqrt {b_{1}}}&a_{1}&{\sqrt {b_{2}}}&\ddots &\vdots \\0&{\sqrt {b_{2}}}&\ddots &\ddots &0\\\vdots &\ddots &\ddots &a_{n-2}&{\sqrt {b_{n-1}}}\\0&\cdots &0&{\sqrt {b_{n-1}}}&a_{n-1}\end{bmatrix}}.}$$

J yson matrices similares y por lo tanto tienen los mismos valores propios (los nodos). Los pesos se pueden calcular a partir de los vectores propios correspondientes: Sies un vector propio normalizado (es decir, un vector propio con norma euclidiana igual a uno) asociado con el valor propio x _j , el peso correspondiente se puede calcular a partir del primer componente de este vector propio, es decir: ${\displaystyle {\mathcal {J}}}$ ${\displaystyle \phi ^{(j)}}$ $${\displaystyle w_{j}=\mu _{0}\left(\phi _{1}^{(j)}\right)^{2}}$$

¿Dónde está la integral de la función peso? ${\displaystyle \mu _{0}}$ $${\displaystyle \mu _{0}=\int _{a}^{b}\omega (x)dx.}$$

Véase, por ejemplo, (Gil, Segura y Temme 2007) para más detalles.

Estimaciones de errores

El error de una regla de cuadratura gaussiana se puede expresar de la siguiente manera. ^[5] Para un integrando que tiene 2 n derivadas continuas, para algún ξ en ( a , b ) , donde p _n es el polinomio ortogonal mónico (es decir, el coeficiente principal es 1 ) de grado n y donde $${\displaystyle \int _{a}^{b}\omega (x)\,f(x)\,dx-\sum _{i=1}^{n}w_{i}\,f(x_{i})={\frac {f^{(2n)}(\xi )}{(2n)!}}\,(p_{n},p_{n})}$$ $${\displaystyle (f,g)=\int _{a}^{b}\omega (x)f(x)g(x)\,dx.}$$

En el caso especial importante de ω ( x ) = 1 , tenemos la estimación del error ^[6] $${\displaystyle {\frac {\left(b-a\right)^{2n+1}\left(n!\right)^{4}}{(2n+1)\left[\left(2n\right)!\right]^{3}}}f^{(2n)}(\xi ),\qquad a<\xi <b.}$$

Stoer y Bulirsch señalan que esta estimación del error es incómoda en la práctica, ya que puede resultar difícil estimar la derivada de orden 2 n y, además, el error real puede ser mucho menor que un límite establecido por la derivada. Otro enfoque consiste en utilizar dos reglas de cuadratura gaussiana de órdenes diferentes y estimar el error como la diferencia entre los dos resultados. Para este propósito, pueden resultar útiles las reglas de cuadratura de Gauss-Kronrod.

Reglas de Gauss-Kronrod

Si se subdivide el intervalo [ a , b ] , los puntos de evaluación de Gauss de los nuevos subintervalos nunca coinciden con los puntos de evaluación anteriores (excepto en cero para números impares), y por lo tanto el integrando debe evaluarse en cada punto. Las reglas de Gauss-Kronrod son extensiones de las reglas de cuadratura de Gauss generadas al agregar n + 1 puntos a una regla de n puntos de tal manera que la regla resultante sea de orden 2 n + 1 . Esto permite calcular estimaciones de orden superior mientras se reutilizan los valores de función de una estimación de orden inferior. La diferencia entre una regla de cuadratura de Gauss y su extensión de Kronrod se usa a menudo como una estimación del error de aproximación.

Reglas de Gauss-Lobatto

También conocida como cuadratura de Lobatto , ^[7] llamada así por el matemático holandés Rehuel Lobatto . Es similar a la cuadratura gaussiana con las siguientes diferencias:

1. Los puntos de integración incluyen los puntos finales del intervalo de integración.
2. Es preciso para polinomios hasta grado 2 n – 3 , donde n es el número de puntos de integración. ^[8]

Cuadratura de Lobatto de la función f ( x ) en el intervalo [−1, 1] : $${\displaystyle \int _{-1}^{1}{f(x)\,dx}={\frac {2}{n(n-1)}}[f(1)+f(-1)]+\sum _{i=2}^{n-1}{w_{i}f(x_{i})}+R_{n}.}$$

Abscisas: x _i es el 1er cero de , aquí denota el polinomio de Legendre estándar de m -ésimo grado y el guión denota la derivada. ${\displaystyle (i-1)}$ ${\displaystyle P'_{n-1}(x)}$ ${\displaystyle P_{m}(x)}$

Pesos: $${\displaystyle w_{i}={\frac {2}{n(n-1)\left[P_{n-1}\left(x_{i}\right)\right]^{2}}},\qquad x_{i}\neq \pm 1.}$$

Resto: $${\displaystyle R_{n}={\frac {-n\left(n-1\right)^{3}2^{2n-1}\left[\left(n-2\right)!\right]^{4}}{(2n-1)\left[\left(2n-2\right)!\right]^{3}}}f^{(2n-2)}(\xi ),\qquad -1<\xi <1.}$$

Algunos de los pesos son:

Una variante adaptativa de este algoritmo con 2 nodos interiores ^[9] se encuentra en GNU Octave y MATLAB como quadly integrate. ^{[10] [11]}

Referencias

Citas

1. Gauss 1815
2. Jacobi 1826
3. Abramowitz y Stegun 1983, pág. 887
4. Stoer y Bulirsch 2002, págs. 172-175
5. Stoer y Bulirsch 2002, Thm 3.6.24
6. Kahaner, Moler y Nash 1989, §5.2
7. Abramowitz y Stegun 1983, pág. 888
8. Quarteroni, Sacco y Saleri 2000
9. Gander y Gautschi 2000
10. MathWorks 2012
11. Eaton y otros, 2018

Bibliografía

• Abramowitz, Milton ; Stegun, Irene Ann , eds. (1983) [junio de 1964]. "Capítulo 25.4, Integración". Manual de funciones matemáticas con fórmulas, gráficos y tablas matemáticas . Serie de Matemáticas Aplicadas. Vol. 55 (Novena reimpresión con correcciones adicionales de la décima impresión original con correcciones (diciembre de 1972); primera ed.). Washington DC; Nueva York: Departamento de Comercio de los Estados Unidos, Oficina Nacional de Normas; Publicaciones de Dover. ISBN 978-0-486-61272-0. LCCN  64-60036. MR  0167642. LCCN  65-12253.
• Anderson, Donald G. (1965). "Fórmulas de cuadratura gaussiana para ∫ 0 1 − ln ⁡ ( x ) f ( x ) d x {\displaystyle \int _{0}^{1}-\ln(x)f(x)dx}". Matemáticas. Comp . 19 (91): 477–481. doi : 10.1090/s0025-5718-1965-0178569-1 .
• Danloy, Bernard (1973). "Construcción numérica de fórmulas de cuadratura gaussiana para y ". Math. Comp . Vol. 27, núm. 124. págs. 861–869. doi :10.1090/S0025-5718-1973-0331730-X. MR  0331730. ${\displaystyle \int _{0}^{1}(-\log x)x^{\alpha }f(x)dx}$ ${\displaystyle \int _{0}^{\infty }E_{m}(x)f(x)dx}$
• Eaton, John W.; Bateman, David; Hauberg, Søren; Wehbring, Rik (2018). «Funciones de una variable (GNU Octave)» . Consultado el 28 de septiembre de 2018 .
• Gander, Walter; Gautschi, Walter (2000). "Cuadratura adaptativa - Revisada". BIT Matemáticas numéricas . 40 (1): 84–101. doi :10.1023/A:1022318402393.
• Gauss, Carl Friedrich (1815). Methodus nova integralium valores per approximationem inveniendi. Com. Soc. Ciencia. Matemáticas de Gotinga. vol. 3. págs. 29–76.fecha de 1814, también en Werke, Band 3, 1876, págs. 163–196. Traducción al inglés por Wikisource.
• Gautschi, Walter (1968). "Construcción de fórmulas de cuadratura de Gauss-Christoffel". Math. Comp . Vol. 22, núm. 102. págs. 251–270. doi :10.1090/S0025-5718-1968-0228171-0. MR  0228171.
• Gautschi, Walter (1970). "Sobre la construcción de reglas de cuadratura gaussianas a partir de momentos modificados". Math. Comp . Vol. 24. págs. 245–260. doi :10.1090/S0025-5718-1970-0285117-6. MR  0285177.
• Gautschi, Walter (2020). Un repositorio de software para cuadraturas gaussianas y funciones de Christoffel . SIAM. ISBN 978-1-611976-34-2.
• Gil, Amparo; Segura, Javier; Temme, Nico M. (2007), "§5.3: Cuadratura de Gauss", Métodos numéricos para funciones especiales , SIAM, ISBN 978-0-89871-634-4
• Golub, Gene H. ; Welsch, John H. (1969). "Cálculo de las reglas de cuadratura de Gauss". Matemáticas de la computación . 23 (106): 221–230. doi : 10.1090/S0025-5718-69-99647-1 . JSTOR  2004418.
• Jacobi, CGJ (1826). "Ueber Gauß' neue Methode, die Werthe der Integrale näherungsweise zu finden". Journal für die Reine und Angewandte Mathematik . 1 . S. 301–308und Werke, banda 6.{{cite journal}}: CS1 maint: postscript (link)
• Kabir, Hossein; Matikolaei, Sayed Amir Hossein Hassanpour (2017). "Implementación de una solución de cuadratura gaussiana generalizada precisa para encontrar el campo elástico en un medio anisotrópico homogéneo". Revista de la Sociedad Serbia de Mecánica Computacional . 11 (1): 11–19. doi :10.24874/jsscm.2017.11.01.02.
• Kahaner, David; Moler, Cleve ; Nash, Stephen (1989). Métodos numéricos y software . Prentice-Hall . ISBN. 978-0-13-627258-8.
• Laudadio, Teresa; Mastronardi, Nicola; Van Dooren, Paul (2023). "Cálculo de reglas de cuadratura gaussiana con alta precisión relativa". Algoritmos numéricos . 92 : 767–793. doi : 10.1007/s11075-022-01297-9 .
• Laurie, Dirk P. (1999), "Recuperación precisa de coeficientes de recursión a partir de fórmulas de cuadratura gaussiana", J. Comput. Appl. Math. , 112 (1–2): 165–180, doi : 10.1016/S0377-0427(99)00228-9
• Laurie, Dirk P. (2001). "Cálculo de fórmulas de cuadratura de tipo Gauss". J. Comput. Appl. Math . 127 (1–2): 201–217. Bibcode :2001JCoAM.127..201L. doi :10.1016/S0377-0427(00)00506-9.
• MathWorks (2012). "Integración numérica - Integral de MATLAB".
• Piessens, R. (1971). "Fórmulas de cuadratura gaussiana para la integración numérica de la integral de Bromwich y la inversión de la transformada de Laplace". J. Eng. Math . Vol. 5, no. 1. págs. 1–9. Bibcode :1971JEnMa...5....1P. doi :10.1007/BF01535429.
• Press, WH ; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Sección 4.6. Cuadraturas gaussianas y polinomios ortogonales", Recetas numéricas: el arte de la computación científica (3.ª ed.), Nueva York: Cambridge University Press, ISBN 978-0-521-88068-8
• Quarteroni, Alfio ; Sacco, Ricardo; Saleri, Fausto (2000). Matemáticas Numéricas. Nueva York: Springer-Verlag . págs. 425–478. doi :10.1007/978-3-540-49809-4_10. ISBN 0-387-98959-5.
• Riener, Cordian; Schweighofer, Markus (2018). "Enfoques de optimización para la cuadratura: Nuevas caracterizaciones de la cuadratura gaussiana en la línea y la cuadratura con pocos nodos en curvas algebraicas planas, en el plano y en dimensiones superiores". Journal of Complexity . 45 : 22–54. arXiv : 1607.08404 . doi :10.1016/j.jco.2017.10.002.
• Sagar, Robin P. (1991). "Una cuadratura gaussiana para el cálculo de integrales generalizadas de Fermi-Dirac". Comput. Phys. Commun . 66 (2–3): 271–275. Bibcode :1991CoPhC..66..271S. doi :10.1016/0010-4655(91)90076-W.
• Stoer, Josef; Bulirsch, Roland (2002), Introducción al análisis numérico (3.ª ed.), Springer , ISBN 978-0-387-95452-3
• Temme, Nico M. (2010), "§3.5(v): Cuadratura de Gauss", en Olver, Frank WJ ; Lozier, Daniel M.; Boisvert, Ronald F.; Clark, Charles W. (eds.), Manual del NIST de funciones matemáticas , Cambridge University Press, ISBN 978-0-521-19225-5, Sr.  2723248.
• Yakimiw, E. (1996). "Cálculo preciso de pesos en las reglas clásicas de cuadratura de Gauss-Christoffel". J. Comput. Phys . 129 (2): 406–430. Bibcode :1996JCoPh.129..406Y. doi :10.1006/jcph.1996.0258.

Enlaces externos

• "Fórmula de cuadratura de Gauss", Enciclopedia de Matemáticas , EMS Press , 2001 [1994]
• ALGLIB contiene una colección de algoritmos para integración numérica (en C# / C++ / Delphi / Visual Basic / etc.)
• Biblioteca científica GNU: incluye la versión C de los algoritmos QUADPACK (consulte también Biblioteca científica GNU )
• De la cuadratura de Lobatto a la constante de Euler e
• Regla de integración de la cuadratura gaussiana: notas, PPT, Matlab, Mathematica, Maple, Mathcad en el Instituto de métodos numéricos holísticos
• Weisstein, Eric W. "Cuadratura Legendre-Gauss". MundoMatemático .
• Cuadratura Gaussiana por Chris Maes y Anton Antonov, Proyecto de Demostraciones Wolfram .
• Pesos tabulados y abscisas con código fuente de Mathematica, pesos de cuadratura Legendre-Gaussiana de alta precisión (16 y 256 decimales) y abscisas, para n = 2 a n = 64, con código fuente de Mathematica.
• Código fuente de Mathematica distribuido bajo la GNU LGPL para la generación de abscisas y pesos para funciones de ponderación arbitrarias W(x), dominios de integración y precisiones.
• Cuadratura gaussiana en Boost.Math, para precisión arbitraria y orden de aproximación
• Cuadratura de Gauss-Kronrod en Boost.Math
• Nodos y pesos de la cuadratura gaussiana Archivado el 14 de abril de 2021 en Wayback Machine