cuadratura gaussiana

En análisis numérico , una regla de cuadratura gaussiana $de n$ puntos , llamada así en honor a Carl Friedrich Gauss , ^[¹^] es una regla de cuadratura construida para producir un resultado exacto para polinomios de grado $2$ $n$ $- 1$ o menor 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 $\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 menor. Esta regla exacta se conoce como la regla de cuadratura de Gauss-Legendre . La regla de cuadratura solo será una aproximación precisa de la integral anterior si $f (x)$ se aproxima bien mediante un polinomio de grado $2 n - 1$ o menor en $[-1, 1]$ .

La regla de cuadratura de Gauss- Legendre no se usa normalmente para funciones integrables con singularidades en los extremos . En cambio, si el integrando se puede escribir como

$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,

$\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 ${\textstyle {\frac {1}{\sqrt {1-x^{2}}}}}$ ( Chebyshev–Gauss ) y ${\estilo de texto {\sqrt {1-x^{2}}}}$ También se puede querer integrar sobre intervalos semiinfinitos ( cuadratura de Gauss-Laguerre ) e infinitos ( cuadratura de Gauss-Hermite ).

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 perteneciente a una clase de polinomios ortogonales (la clase ortogonal respecto a un producto interno ponderado). Esta es una observación clave para el cálculo de los nodos y pesos de la 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)$ se aproxima bien mediante polinomios en $[-1,1]$ , 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 vienen dados por la fórmula ^{[ 3 ]} $w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)\left[P'_{n}(x_{i})\right]^{2}}}.$

A continuación se presentan en una tabla algunas reglas de cuadratura de bajo orden (en el intervalo $[-1, 1]$ ; consulte la sección siguiente 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 gaussiana. Este cambio de intervalo puede realizarse de la siguiente manera: $\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 ${\frac {dx}{d\xi }}={\frac {b-a}{2}}$

Aplicando el $n$ cuadratura gaussiana puntual $(\xi ,w)$ La regla da como resultado la siguiente aproximación: $\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 la 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 $t=8\mathrm {s}$ a $t=30\mathrm {s} ,$ según lo indicado por $s=\int _{8}^{30}{\left(2000\ln \left[{\frac {140000}{140000-2100t}}\right]-9.8t\right){dt}}$

Modifique los límites para poder utilizar los pesos y las abscisas que se muestran en la Tabla 1. Además, calcule el error relativo absoluto verdadero. El valor verdadero es 11061,34 m.

Solución

Primero, cambiar los límites de integración de $\left[8,30\right]$ a $\left[-1,1\right]$ da

${\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 dos puntos,

$c_{1}=1.000000000$
$x_{1}=-0.577350269$
$c_{2}=1.000000000$
$x_{2}=0.577350269$

Ahora podemos usar la fórmula de cuadratura de Gauss. ${\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}}$ desde ${\begin{aligned}f(12.64915)&=2000\ln \left[{\frac {140000}{140000-2100(12.64915)}}\right]-9.8(12.64915)\\&=296.8317\end{aligned}}$ ${\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, $\left|\varepsilon _{t}\right|$ es $\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 un poco 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 $\int _{a}^{b}\omega (x)\,f(x)\,dx$ para algunas elecciones de $a$ , $b$ y $ω$ . Para $a = -1$ , $b = 1$ y $ω (x) = 1$ , el problema es el mismo que el considerado anteriormente. Otras elecciones conducen a otras reglas de integración. Algunas de ellas se tabulan a continuación. Los números de ecuación se proporcionan para Abramowitz y Stegun (A & S).

Teorema fundamental

Sea $p n$ un polinomio no trivial de grado $n$ tal que $\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 ser ortogonal a los otros polinomios $p j$ para $j < n$ , y $x k$ está en el espacio generado por ese conjunto.

Si elegimos los $n$ nodos $x i$ como los ceros de $p n$ , $p_{n}\propto \prod _{i=1}^{n}(x-x_{i})$ Entonces existen $n$ pesos $w i$ que hacen que la integral calculada mediante 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)$ cualquier polinomio de grado $2 n - 1$ o menor. Dividiéndolo por el polinomio ortogonal $p n$ se obtiene $h(x)=p_{n}(x)\,q(x)+r(x).$ donde $q (x)$ es el cociente, de grado $n - 1$ o menor (porque la suma de su grado y el del divisor $p n$ debe ser igual al del dividendo), y $r (x)$ es el resto, también de grado $n - 1$ o menor (porque el grado del resto siempre es menor que el del divisor). Dado que $p n$ es, por hipótesis, ortogonal a todos los monomios de grado menor que $n$ , debe ser ortogonal al cociente $q (x)$ . Por lo tanto $\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 menor, podemos interpolarlo exactamente usando $n$ puntos de interpolación con polinomios de Lagrange $l i (x)$ , donde $l_{i}(x)=\prod _{j\neq i}{\frac {x-x_{j}}{x_{i}-x_{j}}}.$

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

Entonces su integral será igual a $\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 $h(x_{i})=p_{n}(x_{i})\,q(x_{i})+r(x_{i})=r(x_{i}),$ para todo $i$ . Por lo tanto, finalmente tenemos $\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 viene dada exactamente por la suma de cuadratura gaussiana.

Para demostrar la segunda parte de la afirmación, consideremos la forma factorizada del polinomio $p n$ . Cualquier raíz compleja conjugada producirá un factor cuadrático que es estrictamente positivo o estrictamente negativo en toda la recta real . Cualquier factor para raíces fuera del intervalo de $a$ a $b$ no cambiará de signo en ese intervalo. Finalmente, para 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 obtener un nuevo polinomio. $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 allí ahora tienen multiplicidad par. Por lo tanto, la integral $\int _{a}^{b}p_{n}(x)\,\left(\prod _{i}(x-x_{i})\right)\,\omega (x)\,dx\neq 0,$ ya que la función de peso $ω (x)$ siempre es no negativa. Pero $p n$ es ortogonal a todos los polinomios de grado $n - 1$ o menor, por lo que el grado del producto $\prod _{i}(x-x_{i})$ debe ser al menos $n$ . Por lo tanto, $p n$ tiene $n$ raíces distintas, todas reales, en el intervalo de $a$ a $b$ .

Fórmula general para los pesos

Los pesos se pueden expresar como

dónde $a_{k}$ es el coeficiente de $x^{k}$ en $p_{k}(x)$ Para demostrar esto, observe que utilizando la interpolación de Lagrange se puede expresar $r (x)$ en términos de $r(x_{i})$ como $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}}}$ porque $r (x)$ tiene grado menor que $n$ y, por lo tanto, está determinado por los valores que alcanza en $n$ puntos diferentes. Multiplicando ambos lados por $ω (x)$ e integrando de $a$ a $b$ se obtiene $\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 lo tanto $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 $w_{i}$ puede expresarse en términos de polinomios ortogonales $p_{n}(x)$ y $p_{n-1}(x)$ como sigue.

Podemos escribir $\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 $a_{n}$ es el coeficiente de $x^{n}$ en $p_{n}(x)$ . Tomando el límite de $x$ a $x_{i}$ rendimientos utilizando la regla de L'Hôpital $\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 escribir así la expresión integral para los pesos como

En el integrando, escribiendo ${\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 $\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$

proporcionó $k\leq n$ , porque ${\frac {1-\left({\frac {x}{x_{i}}}\right)^{k}}{x-x_{i}}}$ es un polinomio de grado $k - 1$ que es entonces ortogonal a $p_{n}(x)$ . Entonces, si $q (x)$ es un polinomio de grado como máximo n, tenemos $\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 para $q(x)=p_{n-1}(x)$ de la siguiente manera. Porque ${\frac {p_{n}(x)}{x-x_{i}}}$ es un polinomio de grado $n - 1$ , tenemos ${\frac {p_{n}(x)}{x-x_{i}}}=a_{n}x^{n-1}+s(x)$ donde $s (x)$ es un polinomio de grado $n-2$ . Dado que $s (x)$ es ortogonal a $p_{n-1}(x)$ tenemos $\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 $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 $n-2$ , que por lo tanto es ortogonal a $p_{n-1}(x)$ La integral se puede escribir, por lo tanto, como $\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 $p'_{n}(x_{i})$ y eso produce la expresión en la ecuación ( 1 ).

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

Prueba de que los pesos son positivos

Consideremos el siguiente polinomio de grado $2n-2$ $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}}}$ donde, como se indicó anteriormente, las $x j$ son las raíces del polinomio $p_{n}(x)$ . Claramente $f(x_{j})=\delta _{ij}$ . Dado que el grado de $f(x)$ es menor que $2n-1$ , la fórmula de cuadratura gaussiana que involucra los pesos y nodos obtenidos de $p_{n}(x)$ Se aplica. Dado que $f(x_{j})=0$ para $j$ distinto de $i$ , tenemos $\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.$

Dado que ambos $\omega (x)$ y $f(x)$ son funciones no negativas, por lo tanto, se deduce que $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 2)$ operaciones, el método de Newton para resolver $p_{n}(x)=0$ utilizando la recurrencia de tres términos para la evaluación que requiere $O (n 2)$ operaciones, y fórmulas asintóticas para n grande que requieren $O (n)$ operaciones.

Relación de recurrencia

Polinomios ortogonales $p_{r}$ con $(p_{r},p_{s})=0$ para $r\neq s$ para un producto escalar $(\cdot ,\cdot )$ , grado $(p_{r})=r$ y el coeficiente principal uno (es decir, los polinomios ortogonales mónicos ) satisfacen la relación de recurrencia $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 $(f(x),g(x))=\int _{a}^{b}\omega (x)f(x)g(x)dx$

para $r=0,1,\ldots ,n-1$ donde $n$ es el grado máximo que puede tomarse como infinito, y donde ${\textstyle a_{r,s}={\frac {\left(xp_{r},p_{s}\right)}{\left(p_{s},p_{s}\right)}}}$ . En primer lugar, los polinomios definidos por la relación de recurrencia que comienza con $p_{0}(x)=1$ tienen coeficiente principal uno y grado correcto. Dado el punto de partida por $p_{0}$ , la ortogonalidad de $p_{r}$ puede demostrarse por inducción. $r=s=0$ uno tiene $(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 $p_{0},p_{1},\ldots ,p_{r}$ son ortogonales, entonces también $p_{r+1}$ , porque en $(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})$ Todos los productos escalares se anulan excepto el primero y el que $p_{s}$ encuentra el mismo polinomio ortogonal. Por lo tanto, $(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 $(xf,g)=(f,xg)$ (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 $s<r-1,xp_{s}$ es un polinomio de grado menor o igual a $r - 1$ . Por otro lado, $p_{r}$ es ortogonal a todo polinomio de grado menor o igual a $r - 1$ . Por lo tanto, se tiene $(xp_{r},p_{s})=(p_{r},xp_{s})=0$ y $a_{r,s}=0$ para $s < r - 1$ . La relación de recurrencia se simplifica entonces a $p_{r+1}(x)=(x-a_{r,r})p_{r}(x)-a_{r,r-1}p_{r-1}(x)$

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

(con la convención $p_{-1}(x)\equiv 0$ ) dónde $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 $(xp_{r},p_{r-1})=(p_{r},xp_{r-1})=(p_{r},p_{r})$ , desde $xp_{r-1}$ difiere de $p_{r}$ por un grado menor que $r$ ).

El algoritmo de Golub-Welsch

La relación de recurrencia de tres términos se puede escribir en forma matricial. $J{\tilde {P}}=x{\tilde {P}}-p_{n}(x)\mathbf {e} _{n}$ dónde ${\tilde {P}}={\begin{bmatrix}p_{0}(x)&p_{1}(x)&\cdots &p_{n-1}(x)\end{bmatrix}}^{\mathsf {T}}$ , $\mathbf {e} _{n}$ es el $n$ vector base estándar , es decir, $\mathbf {e} _{n}={\begin{bmatrix}0&\cdots &0&1\end{bmatrix}}^{\mathsf {T}}$ y $J$ es la siguiente matriz tridiagonal , llamada matriz de Jacobi: $\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 $x_{j}$ Los polinomios de 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 .

Para calcular los pesos y los nodos, es preferible considerar la matriz tridiagonal simétrica . ${\mathcal {J}}$ con elementos ${\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,

${\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$ y ${\mathcal {J}}$ son 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 : Si $\phi ^{(j)}$ es 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: $w_{j}=\mu _{0}\left(\phi _{1}^{(j)}\right)^{2}$

dónde $\mu _{0}$ es la integral de la función de peso $\mu _{0}=\int _{a}^{b}\omega (x)dx.$

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

Estimaciones de error

El error de una regla de cuadratura gaussiana se puede enunciar de la siguiente manera. ^{[ 5 ]} Para un integrando que tiene $2 n$ derivadas continuas, $\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})$ 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 $(f,g)=\int _{a}^{b}\omega (x)f(x)g(x)\,dx.$

En el caso especial importante de $ω (x) = 1$ , tenemos la estimación de error ^{[ 6 ]} ${\frac {\left(ba\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 de error resulta inconveniente en la práctica, ya que puede ser difícil estimar la derivada de orden $2n$ $,$ 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 ambos resultados. Para ello, pueden resultar útiles las reglas de cuadratura de Gauss-Kronrod.

Reglas de Gauss-Kronrod

Si el intervalo $[a, b]$ se subdivide, 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), por lo que 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 $2n + 1. Esto permite calcular estimaciones de orden superior mientras se reutilizan los valores de$ la 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

En algunas aplicaciones, es deseable contar con reglas de cuadratura que tengan la alta precisión de las fórmulas de Gauss, pero que también incluyan los extremos del intervalo entre los puntos de evaluación. Dichas reglas se conocen como Gauss-Lobatto , o simplemente cuadratura de Lobatto , ^{[ 7 ]} nombradas en honor al matemático neerlandés Rehuel Lobatto . Dado que para una regla de n puntos ya no se pueden elegir libremente las ubicaciones de todos los puntos de cuadratura (dos de los puntos están fijos en los extremos), cabe esperar que la regla sea menos precisa que la cuadratura gaussiana regular. De hecho, una regla de Gauss-Lobatto $de$ n puntos solo es precisa para polinomios de hasta grado $2n -3$ . ^{[ 8 ]}

Cuadratura de Lobatto de la función $f (x)$ en el intervalo $[-1, 1]$ : $\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 la $(i-1)$ primer cero de $P'_{n-1}(x)$ , aquí $P_{m}(x)$ denota el polinomio de Legendre estándar de grado $m$ y el guion denota la derivada.

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

Resto: $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

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 edición original con correcciones (diciembre de 1972); primera ed.). Washington D. C.; Nueva York: Departamento de Comercio de los Estados Unidos, Oficina Nacional de Normas; Dover Publications. ISBN 978-0-486-61272-0. LCCN 64-60036 . MR 0167642 . LCCN 65-12253 .
Anderson, Donald G. (1965). "Fórmulas de cuadratura gaussiana para $\int _{0}^{1}-\ln(x)f(x)dx$ " . Math. 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 $\int _{0}^{1}(-\log x)x^{\alpha }f(x)dx$ y $\int _{0}^{\infty }E_{m}(x)f(x)dx$ ". Math. Comp . 27 (124): 861– 869. doi : 10.1090/S0025-5718-1973-0331730-X . MR 0331730 .
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). "Cuadradura adaptativa: una revisión" . Matemáticas numéricas BIT . 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". Math. Comp . 22 (102): 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 . 24 (110): 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 mantenimiento: postscript ( enlace )
Kabir, Hossein; Matikolaei, Sayed Amir Hossein Hassanpour (2017). "Implementación de una solución precisa de cuadratura gaussiana generalizada 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" . Numerical Algorithms . 92 : 767–793 . doi : 10.1007/s11075-022-01297-9 . hdl : 2078.1/272678 .
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 . 5 (1): 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" , Numerical Recipes: The Art of Scientific Computing (3.ª ed.), Nueva York: Cambridge University Press, ISBN 978-0-521-88068-8Archivado del original el 19 de marzo de 2012 , consultado el 8 de agosto de 2011.
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.), NIST Handbook of Mathematical Functions , Cambridge University Press, ISBN 978-0-521-19225-5, MR 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 la integración numérica (en C# / C++ / Delphi / Visual Basic / etc.).
Biblioteca Científica GNU : incluye la versión en C de los algoritmos QUADPACK (véase también Biblioteca Científica GNU ).
De la cuadratura de Lobatto a la constante de Euler e
Regla de cuadratura gaussiana para la integración: 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 de Wolfram .
Pesos y abscisas tabulados con código fuente de Mathematica , pesos y abscisas de cuadratura Legendre-Gaussiana de alta precisión (16 y 256 decimales), para n =2 a n =64, con código fuente de Mathematica.
Código fuente de Mathematica distribuido bajo la licencia 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 y orden de aproximación arbitrarios.
Cuadratura de Gauss-Kronrod en Boost.Math
Nodos y pesos de la cuadratura gaussiana. Archivado el 14 de abril de 2021 en Wayback Machine.

[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 et al. 2018

[

[ 2 ]

[ 3 ]

[ 4 ]

[ 5 ]

[ 6 ]

[ 7 ]

[ 8 ]

[ 9 ]

[ 10 ]

[ 11 ]