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 3 – 8 x 2 – 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
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
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,
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]
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:
con
La aplicación de la regla de cuadratura gaussiana del punto da como resultado la siguiente aproximación:
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
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
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,
Ahora podemos usar la fórmula de la cuadratura de Gauss.
Dado que el valor verdadero es 11061,34 m, el error verdadero relativo absoluto es
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
Sea p n un polinomio no trivial de grado n tal que
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
q ( x )n − 1p nr ( x )n − 1p nnq ( x )
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
Tenemos
Entonces su integral será igual
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
yo
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.
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
ω ( x )p nn -1
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
r ( x )nnω ( x )ab
Los pesos w i están dados por tanto
Esta expresión integral se puede expresar en términos de polinomios ortogonales de la siguiente manera.
Podemos escribir
¿ Dónde está el coeficiente de in ? Tomando el límite de x para obtener resultados usando la regla de L'Hôpital
Por lo tanto, podemos escribir la expresión integral para los pesos como
En el integrando, escribiendo
rendimientos
proporcionado , porque
k − 1q ( x )
Podemos evaluar la integral en el lado derecho de la siguiente manera. Como es un polinomio de grado n − 1 , tenemos
s ( x )s ( x )
Entonces podemos escribir
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
Según la ecuación ( 2 ), los pesos se obtienen dividiendo esto por y eso produce la expresión en la ecuación ( 1 ).
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 .
Prueba de que los pesos son positivos.
Considere el siguiente polinomio de grado
x jji
Como ambas y son funciones no negativas, se deduce que .
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 2 ) , el método de Newton para resolver utilizando la recurrencia de tres términos para la evaluación que requiere operaciones O ( n 2 ) y fórmulas asintóticas para n grandes que requieren operaciones O ( n ) .
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
Ahora bien, si son ortogonales, entonces también , porque en
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
o
(con la convención ) donde
(el último debido a , ya que difiere de en 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 donde , es el enésimo vector de base estándar, es decir, y J es la siguiente matriz tridiagonal , llamada matriz de Jacobi:
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 .
Para calcular los pesos y nodos, es preferible considerar la matriz tridiagonal simétrica con elementos
Eso es,
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:
¿Dónde está la integral de la función de peso?
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,
ξ( a , b )p n1n
En el importante caso especial de ω ( x ) = 1 , tenemos la estimación del error [6]
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:
Los puntos de integración incluyen los puntos finales del intervalo de integración.
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] :
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.
Pesos:
Resto:
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
^ Gauss 1815
^ Jacobi 1826
^ Abramowitz y Stegun 1983, pág. 887
^ Stoer y Bulirsch 2002, págs. 172-175
^ Stoer y Bulirsch 2002, Thm 3.6.24
^ Kahaner, Moler y Nash 1989, §5.2
^ Abramowitz y Stegun 1983, pág. 888
^ Quarteroni, Sacco y Saleri 2000
^ Gander y Gautschi 2000
^ Trabajos de matemáticas 2012
^ 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.
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.
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
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
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.
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