cuadratura gaussiana

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

Comparación entre cuadratura gaussiana y trapezoidal de 2 puntos.
La curva azul muestra la función cuya integral definida en el intervalo [−1, 1] se va a 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 está representada por 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 se aproxima mejor a 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 la cuadratura gaussiana de 2 puntos incluso arroja un resultado exacto.

En 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 expresa como

$${\displaystyle \int _{-1}^{1}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i}),}$$

lo cual es exacto para polinomios de grado 2 n − 1 o menos. Esta regla exacta se conoce como regla de cuadratura de Gauss-Legendre . La regla de la cuadratura solo será una aproximación precisa a la integral anterior si f ( x ) está bien aproximada por 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 de punto final . 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 ) está bien aproximado 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. Éstas 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 es posible que desee integrar en intervalos semiinfinitos ( 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 (ver 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 pesos y nodos de cuadratura de Gauss.

Cuadratura de Gauss-Legendre

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

Para el problema de integración más simple mencionado anteriormente, es decir, f ( x ) está bien aproximado 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, xi , es la i -ésima raíz de P _n y los pesos vienen dados por 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}}}.}$$

Algunas reglas de cuadratura de bajo orden se tabulan a continuación (sobre el intervalo [−1, 1] , consulte la sección siguiente para conocer otros intervalos).

Cambio de intervalo

Una integral sobre [ a , b ] debe convertirse en una integral sobre [−1, 1] antes de aplicar la regla de cuadratura gaussiana. 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 del punto 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 la cuadratura de Gauss de dos puntos para aproximar la distancia en metros cubierta por un cohete desde a como está dada por ${\displaystyle t=8\mathrm {s} }$ ${\displaystyle t=30\mathrm {s} ,}$

$${\displaystyle x=\int _{8}^{30}{\left(2000\ln \left[{\frac {140000}{140000-2100t}}\right]-9.8t\right){dt}}}$$

Cambie los límites para que se puedan usar los pesos y abscisas que se dan en la Tabla 1. Además, encuentre el error verdadero relativo absoluto. El valor real es 11061,34 m.

Solución

Primero, cambiar los límites de integración de a da ${\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 la función de la Tabla 1 para la regla de los dos puntos,

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

Ahora podemos usar la fórmula de la cuadratura de Gauss.

$${\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 verdadero 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 es calcular

$${\displaystyle \int _{a}^{b}\omega (x)\,f(x)\,dx}$$

abωa = −1b = 1ω ( x ) = 1Abramowitz y Stegun

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.}$$

Tenga en cuenta que esto será cierto para todos los polinomios ortogonales anteriores, porque cada p _n se construye para que sea ortogonal a los otros polinomios p _j para j < n , y x ^k está en el intervalo de ese conjunto.

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

Para probar la primera parte de esta afirmación, sea h ( x ) cualquier polinomio de grado 2 n − 1 o menos. Divídelo por el polinomio ortogonal p _n para obtener

$${\displaystyle h(x)=p_{n}(x)\,q(x)+r(x).}$$

q ( x )n − 1p _nr ( x )n − 1p _nnq ( 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.}$$

Dado que el resto r ( x ) es de grado n − 1 o menos, 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 ) (consulte a continuación otras fórmulas para los pesos). Pero todas las x _i son raíces de p _n , por lo que la fórmula de división anterior nos dice que

$${\displaystyle h(x_{i})=p_{n}(x_{i})\,q(x_{i})+r(x_{i})=r(x_{i}),}$$

yo

$${\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 menos, su integral viene dada exactamente por la suma de cuadratura gaussiana.

Para probar 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 aab no cambiará de signo durante ese intervalo . Finalmente, para los factores correspondientes a las raíces x _i dentro del intervalo de a a b que son 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 ahora son de multiplicidad par. Entonces la integral

$${\displaystyle \int _{a}^{b}p_{n}(x)\,\left(\prod _{i}(x-x_{i})\right)\,\omega (x)\,dx\neq 0,}$$

ω ( x )p _nn -1

$${\displaystyle \prod _{i}(x-x_{i})}$$

np _nnab

Fórmula general para los pesos.

Los pesos se pueden expresar como

¿ Dónde está el coeficiente de in ? Para probar esto, tenga en cuenta que usando la interpolación de Lagrange se puede expresar r ( x ) en términos de como ${\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}}}}$$

r ( x )nnω ( x )ab

$${\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 están dados por tanto

$${\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 se puede expresar en términos de polinomios ortogonales 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)}}}$$

¿ Dónde está el coeficiente de in ? Tomando el límite de x para obtener resultados 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}}}}$$

Por lo tanto, podemos 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}$$

proporcionado , porque ${\displaystyle k\leq n}$

$${\displaystyle {\frac {1-\left({\frac {x}{x_{i}}}\right)^{k}}{x-x_{i}}}}$$

k − 1q ( x ) ${\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 en el lado derecho de la siguiente manera. Como es un polinomio de grado n − 1 , 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)}$$

s ( x )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 tanto es ortogonal a . Por tanto, la integral se puede escribir como ${\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}$$

Según la ecuación ( 2 ), los pesos se obtienen dividiendo esto por y eso produce la expresión en la ecuación ( 1 ). ${\displaystyle p'_{n}(x_{i})}$

${\displaystyle w_{i}}$También se puede expresar en términos de polinomios ortogonales y ahora . En la relación de recurrencia de 3 términos, el término con desaparece, por lo que en la ecuación. (1) puede sustituirse 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 ${\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}}}}$$

x _jji ${\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 ambas y son funciones no negativas, se deduce que . ${\displaystyle \omega (x)}$ ${\displaystyle f(x)}$ ${\displaystyle w_{i}>0}$

Cálculo de reglas de cuadratura gaussianas.

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 operaciones O ( n ² ) , el método de Newton para resolver utilizando la recurrencia de tres términos para la evaluación que requiere operaciones O ( n ² ) y fórmulas asintóticas para n grandes que requieren operaciones O ( n ) . ${\displaystyle p_{n}(x)=0}$

Relación de recurrencia

Los polinomios ortogonales con for 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}$$

porque donde n es el grado máximo que puede considerarse infinito, y donde . En primer lugar, los polinomios definidos por la relación de recurrencia que comienza con tienen el coeficiente principal uno y el grado correcto. Dado el punto de partida de , la ortogonalidad de puede demostrarse por inducción. porque 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 ${\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: For 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 enésimo vector de 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 hasta el grado n , que se utilizan como nodos para la cuadratura gaussiana, se pueden encontrar 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: sihay 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, a saber: ${\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 de 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 error

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,

$${\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})}$$

ξ( a , b )p _n1n

$${\displaystyle (f,g)=\int _{a}^{b}\omega (x)f(x)g(x)\,dx.}$$

En el importante caso especial 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 inconveniente en la práctica, ya que puede ser 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 gaussianas de diferentes órdenes 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 sumando 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 utiliza a menudo como estimación del error de aproximación.

Reglas de Gauss-Lobatto

También conocida como cuadratura de Lobatto , ^[7] lleva el nombre del 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 exacto para polinomios hasta el 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 st cero de , aquí denota el polinomio estándar de Legendre 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. Trabajos de matemáticas 2012
11. Eaton y col. 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áficas 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. SEÑOR  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 ". Matemáticas. comp . vol. 27, núm. 124, págs. 861–869. doi :10.1090/S0025-5718-1973-0331730-X. SEÑOR  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 (Octava GNU)" . 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 por aproximación de inventos. 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". Matemáticas. comp . vol. 22, núm. 102, págs. 251–270. doi :10.1090/S0025-5718-1968-0228171-0. SEÑOR  0228171.
• Gautschi, Walter (1970). "Sobre la construcción de reglas de cuadratura gaussianas a partir de momentos modificados". Matemáticas. comp . vol. 24. págs. 245–260. doi :10.1090/S0025-5718-1970-0285117-6. SEÑOR  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". Revista para 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" (PDF) . 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). "Calcular reglas de cuadratura gaussianas con alta precisión relativa". Algoritmos Numéricos . 92 : 767–793.
• Laurie, Dirk P. (1999), "Recuperación precisa de coeficientes de recursividad a partir de fórmulas de cuadratura gaussiana", J. Comput. Aplica. Matemáticas. , 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. Computación. Aplica. Matemáticas . 127 (1–2): 201–217. Código Bib : 2001JCoAM.127..201L. doi :10.1016/S0377-0427(00)00506-9.
• MathWorks (2012). "Integración numérica - integral 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. Ing. Matemáticas . vol. 5, núm. 1. págs. 1–9. Código bibliográfico : 1971JEnMa...5....1P. doi :10.1007/BF01535429.
• Presione, WH ; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Sección 4.6. Cuadraturas gaussianas y polinomios ortogonales", Recetas numéricas: el arte de la informática 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.422, 425. ISBN 0-387-98959-5.
• Riener, Cordian; Schweighofer, Markus (2018). "Enfoques de optimización de la cuadratura: Nuevas caracterizaciones de la cuadratura gaussiana en la recta y la cuadratura con pocos nodos en curvas algebraicas planas, en el plano y en dimensiones superiores". Revista de Complejidad . 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". Computadora. Física. Comunitario . 66 (2–3): 271–275. Código bibliográfico : 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 de funciones matemáticas del NIST , Cambridge University Press, ISBN 978-0-521-19225-5, señor  2723248.
• Yakimiw, E. (1996). "Cálculo preciso de pesos en las reglas clásicas de cuadratura de Gauss-Christoffel". J. Computación. Física . 129 (2): 406–430. Código Bib : 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 cuadratura gaussiana: notas, PPT, Matlab, Mathematica, Maple, Mathcad en el Holistic Numerical Methods Institute
• Weisstein, Eric W. "Cuadratura Legendre-Gauss". MundoMatemático .
• Cuadratura gaussiana de Chris Maes y Anton Antonov, Proyecto de demostraciones Wolfram .
• Pesos tabulados y abscisas con código fuente de Mathematica, alta precisión (16 y 256 decimales) Pesos y abscisas en cuadratura Legendre-Gaussiana, para n =2 a n =64, con código fuente de Mathematica.
• Código fuente de Mathematica distribuido bajo GNU LGPL para abscisas y generación de 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