La cuadratura de Clenshaw-Curtis y la cuadratura de Fejér son métodos para la integración numérica , o "cuadratura", que se basan en una expansión del integrando en términos de polinomios de Chebyshev . De manera equivalente, emplean un cambio de variables y use una aproximación de transformada de coseno discreta (DCT) para la serie de coseno . Además de tener una precisión de convergencia rápida comparable a las reglas de cuadratura de Gauss, la cuadratura de Clenshaw-Curtis conduce naturalmente a reglas de cuadratura anidadas (donde diferentes órdenes de precisión comparten puntos), lo cual es importante tanto para la cuadratura adaptativa como para la cuadratura multidimensional ( cubatura ).
Brevemente, la función a integrarse se evalúa en el extremos o raíces de un polinomio de Chebyshev y estos valores se utilizan para construir una aproximación polinomial para la función. Este polinomio se integra entonces exactamente. En la práctica, las ponderaciones de integración para el valor de la función en cada nodo se calculan previamente, y este cálculo se puede realizar entiempo mediante algoritmos rápidos relacionados con la transformada de Fourier para el DCT. [1] [2]
Método general
Una forma sencilla de entender el algoritmo es darse cuenta de que la cuadratura de Clenshaw-Curtis (propuesta por esos autores en 1960) [3] equivale a integrar mediante un cambio de variable x = cos (θ). El algoritmo se expresa normalmente para la integración de una función f ( x ) sobre el intervalo [-1,1] (cualquier otro intervalo puede obtenerse mediante el reajuste apropiado). Para esta integral, podemos escribir:
Es decir, hemos transformado el problema de integrar a uno de integrar . Esto se puede realizar si conocemos la serie del coseno para:
en cuyo caso la integral se convierte en:
Por supuesto, para calcular los coeficientes de la serie del coseno
hay que volver a realizar una integración numérica, por lo que al principio puede parecer que esto no ha simplificado el problema. Sin embargo, a diferencia del cálculo de integrales arbitrarias, las integraciones en serie de Fourier para funciones periódicas (como, por construcción), hasta la frecuencia de Nyquist , son calculados con precisión por el puntos igualmente espaciados y ponderados por igual por (excepto que los puntos finales están ponderados por 1/2, para evitar el doble conteo, equivalente a la regla trapezoidal o la fórmula de Euler-Maclaurin ). [4] [5] Es decir, aproximamos la integral de la serie del coseno por la transformada de coseno discreta de tipo I (DCT):
por y luego use la fórmula anterior para la integral en términos de estos . Porque soloes necesario, la fórmula se simplifica aún más en un DCT de tipo I de orden N / 2, asumiendo que N es un número par :
A partir de esta fórmula, queda claro que la regla de cuadratura de Clenshaw-Curtis es simétrica, ya que pondera f ( x ) yf (- x ) por igual.
Debido al aliasing , solo se calculan los coeficienteshasta k = N / 2, ya que el muestreo discreto de la función hace que la frecuencia de 2 k sea indistinguible de la de N –2 k . De manera equivalente, elson las amplitudes del polinomio de interpolación trigonométrica limitada de banda única que pasa por los puntos N +1 donde se evalúa f (cos θ ), y aproximamos la integral por la integral de este polinomio de interpolación. Hay cierta sutileza en la forma en que se trata a los coeficiente en la integral, sin embargo, para evitar la doble contabilización con su alias, se incluye con peso 1/2 en la integral aproximada final (como también se puede ver al examinar el polinomio de interpolación):
Conexión a polinomios de Chebyshev
La razón por la que esto está conectado a los polinomios de Chebyshev es que, por definición, , por lo que la serie de coseno anterior es realmente una aproximación de por polinomios de Chebyshev:
y así estamos "realmente" integrando integrando su expansión aproximada en términos de polinomios de Chebyshev. Los puntos de evaluacióncorresponden a los extremos del polinomio de Chebyshev.
El hecho de que tal aproximación de Chebyshev sea solo una serie de coseno bajo un cambio de variables es responsable de la rápida convergencia de la aproximación a medida que más términosestán incluidos. Una serie de coseno converge muy rápidamente para funciones que son uniformes , periódicas y suficientemente suaves. Esto es cierto aquí, ya que es uniforme y periódico en por construcción, y es k veces diferenciable en todas partes sies k -veces diferenciable en. (Por el contrario, la aplicación directa de una expansión en serie de coseno a en vez de generalmente no convergerá rápidamente porque la pendiente de la extensión periódica par generalmente sería discontinua).
Cuadratura de Fejér
Fejér propuso dos reglas de cuadratura muy similares a la cuadratura de Clenshaw-Curtis, pero mucho antes (en 1933). [6]
De estos dos, la "segunda" regla de cuadratura de Fejér es casi idéntica a la de Clenshaw-Curtis. La única diferencia es que los puntos finales y se ponen a cero. Es decir, Fejér solo usó los extremos interiores de los polinomios de Chebyshev, es decir, los verdaderos puntos estacionarios.
La "primera" regla de cuadratura de Fejér evalúa la evaluando en un conjunto diferente de puntos igualmente espaciados, a medio camino entre los extremos: por . Estas son las raices de, y se conocen como los nodos de Chebyshev . (Estos puntos medios igualmente espaciados son la única otra opción de puntos de cuadratura que preservan tanto la simetría par de la transformada del coseno como la simetría de traslación de la serie periódica de Fourier). Esto conduce a una fórmula:
que es precisamente el DCT tipo II. Sin embargo, la primera regla de cuadratura de Fejér no está anidada: los puntos de evaluación para 2 N no coinciden con ninguno de los puntos de evaluación para N , a diferencia de la cuadratura Clenshaw-Curtis o la segunda regla de Fejér.
A pesar de que Fejér descubrió estas técnicas antes que Clenshaw y Curtis, el nombre "cuadratura de Clenshaw-Curtis" se ha convertido en estándar.
Comparación con la cuadratura gaussiana
El método clásico de cuadratura gaussiana evalúa el integrando enpuntos y está construido para integrar exactamente polinomios hasta el grado . En contraste, la cuadratura de Clenshaw-Curtis, arriba, evalúa el integrando en puntos e integra exactamente polinomios solo hasta el grado . Por tanto, puede parecer que Clenshaw-Curtis es intrínsecamente peor que la cuadratura gaussiana, pero en realidad no parece ser así.
En la práctica, varios autores han observado que Clenshaw-Curtis puede tener una precisión comparable a la de la cuadratura gaussiana para el mismo número de puntos. Esto es posible porque la mayoría de los integrandos numéricos no son polinomios (especialmente porque los polinomios pueden integrarse analíticamente) y la aproximación de muchas funciones en términos de polinomios de Chebyshev converge rápidamente (ver Aproximación de Chebyshev ). De hecho, resultados teóricos recientes [7] argumentan que tanto la cuadratura de Gauss como la de Clenshaw-Curtis tienen un error delimitado porpara un integrando diferenciable k veces.
Una ventaja que se cita a menudo de la cuadratura de Clenshaw-Curtis es que los pesos de cuadratura se pueden evaluar en tiempo por algoritmos de transformada rápida de Fourier (o sus análogos para el DCT), mientras que la mayoría de los algoritmos para pesos en cuadratura gaussiana requeríantiempo para calcular. Sin embargo, los algoritmos recientes han logradocomplejidad para la cuadratura de Gauss-Legendre. [8] Como cuestión práctica, la integración numérica de alto orden rara vez se realiza simplemente evaluando una fórmula de cuadratura para muy grandes. En su lugar, generalmente se emplea un esquema de cuadratura adaptativa que primero evalúa la integral a un orden bajo y luego refina sucesivamente la precisión aumentando el número de puntos de muestra, posiblemente solo en regiones donde la integral es inexacta. Para evaluar la precisión de la cuadratura, se compara la respuesta con la de una regla de cuadratura de orden aún menor. Idealmente, esta regla de cuadratura de orden inferior evalúa el integrando en un subconjunto de los N puntos originales , para minimizar las evaluaciones de integrando. Esto se llama una regla de cuadratura anidada , y aquí Clenshaw-Curtis tiene la ventaja de que la regla de orden N utiliza un subconjunto de los puntos de orden 2 N . En contraste, las reglas de cuadratura de Gauss no están anidadas de forma natural, por lo que se deben emplear fórmulas de cuadratura de Gauss-Kronrod o métodos similares. Las reglas anidadas también son importantes para las cuadrículas dispersas en cuadratura multidimensional, y la cuadratura de Clenshaw-Curtis es un método popular en este contexto. [9]
Integración con funciones de peso
De manera más general, se puede plantear el problema de integrar un contra una función de peso fijo que se conoce de antemano:
El caso más común es , como anteriormente, pero en ciertas aplicaciones es deseable una función de peso diferente. La razón básica es que, dado quepuede tenerse en cuenta a priori , se puede hacer que el error de integración dependa sólo de la precisión en la aproximación, independientemente de lo mal que se comporte la función de ponderación.
La cuadratura de Clenshaw-Curtis se puede generalizar a este caso de la siguiente manera. Como antes, funciona encontrando la expansión de la serie coseno dea través de un DCT, y luego integrando cada término en la serie del coseno. Ahora, sin embargo, estas integrales tienen la forma
Para la mayoría , esta integral no se puede calcular analíticamente, a diferencia de antes. Dado que la misma función de ponderación se usa generalmente para muchos integrandos, sin embargo, uno puede darse el lujo de calcular estos numéricamente a alta precisión de antemano. Además, dado que generalmente se especifica analíticamente, a veces se pueden emplear métodos especializados para calcular .
Por ejemplo, se han desarrollado métodos especiales para aplicar la cuadratura de Clenshaw-Curtis a integrandos de la forma con función de peso que es muy oscilatoria, por ejemplo, una función sinusoide o de Bessel (véase, por ejemplo, Evans y Webster, 1999 [10] ). Esto es útil para el cálculo de la serie de Fourier de alta precisión y la serie de Fourier-Bessel , dondeLos métodos de cuadratura son problemáticos debido a la alta precisión requerida para resolver la contribución de las oscilaciones rápidas. Aquí, la parte de oscilación rápida del integrando se tiene en cuenta mediante métodos especializados para, mientras que la función desconocida suele comportarse mejor.
Otro caso en el que las funciones de ponderación son especialmente útiles es si el integrando es desconocido pero tiene una singularidad conocida de alguna forma, por ejemplo, una discontinuidad conocida o divergencia integrable (como 1 / √ x ) en algún punto. En este caso, la singularidad se puede incorporar a la función de peso. y sus propiedades analíticas se pueden utilizar para calcular con precisión de antemano.
Tenga en cuenta que la cuadratura gaussiana también se puede adaptar para varias funciones de peso, pero la técnica es algo diferente. En la cuadratura de Clenshaw-Curtis, el integrando siempre se evalúa en el mismo conjunto de puntos independientemente de, correspondiente a los extremos o raíces de un polinomio de Chebyshev. En la cuadratura gaussiana, diferentes funciones de peso conducen a diferentes polinomios ortogonales y, por lo tanto, a diferentes raíces donde se evalúa el integrando.
Integración en intervalos infinitos y semi-infinitos
También es posible utilizar la cuadratura de Clenshaw-Curtis para calcular integrales de la forma y , utilizando una técnica de reasignación de coordenadas. [11] La alta precisión, incluso la convergencia exponencial para integrandos suaves, se puede conservar siempre quedecae con la suficiente rapidez como | x | se acerca al infinito.
Una posibilidad es utilizar una transformación de coordenadas genérica como x = t / (1− t 2 )
para transformar un intervalo infinito o semi-infinito en uno finito, como se describe en Integración numérica . También existen técnicas adicionales que se han desarrollado específicamente para la cuadratura de Clenshaw-Curtis.
Por ejemplo, se puede utilizar la reasignación de coordenadas , donde L es una constante especificada por el usuario (se podría simplemente usar L = 1; una elección óptima de L puede acelerar la convergencia, pero depende del problema [11] ), para transformar la integral semi-infinita en:
El factor que multiplica sen (θ), f (...) / (...) 2 , puede luego expandirse en una serie de coseno (aproximadamente, usando la transformada discreta del coseno) e integrarse término por término, exactamente como era hecho para f (cos θ) arriba. Para eliminar la singularidad en θ = 0 en este integrando, simplemente se requiere que f ( x ) vaya a cero lo suficientemente rápido cuando x se aproxima al infinito, y en particular f ( x ) debe decaer al menos tan rápido como 1 / x 3/2 . [11]
Para un intervalo de integración doblemente infinito, se puede utilizar la reasignación de coordenadas (donde L es una constante especificada por el usuario como arriba) para transformar la integral en: [11]
En este caso, hemos utilizado el hecho de que el integrando reasignado f ( L cotθ) / sen 2 (θ) ya es periódico y, por lo tanto, puede integrarse directamente con una precisión alta (incluso exponencial) utilizando la regla trapezoidal (asumiendo que f es suficientemente suave y decayendo rápidamente); no es necesario calcular la serie del coseno como paso intermedio. Tenga en cuenta que la regla de la cuadratura no incluye los puntos finales, donde hemos asumido que el integrando va a cero. La fórmula anterior requiere que f ( x ) decaiga más rápido que 1 / x 2 cuando x llega a ± ∞. (Si f decae exactamente como 1 / x 2 , entonces el integrando pasa a un valor finito en los puntos finales y estos límites deben incluirse como términos de punto final en la regla trapezoidal. [11] ). Sin embargo, si f decae solo polinomialmente rápidamente, entonces puede ser necesario usar un paso adicional de la cuadratura de Clenshaw-Curtis para obtener una precisión exponencial de la integral reasignada en lugar de la regla trapezoidal, dependiendo de más detalles de las propiedades limitantes de f : El problema es que, aunque f ( L cotθ) / sen 2 (θ) es de hecho periódica con período π, no es necesariamente suave en los puntos finales si todas las derivadas no desaparecen allí [por ejemplo, la función f ( x ) = tanh ( x 3 ) / x 3 decae como 1 / x 3 pero tiene una discontinuidad de salto en la pendiente de la función reasignada en θ = 0 y π].
Se sugirió otro enfoque de reasignación de coordenadas para integrales de la forma , en cuyo caso se puede utilizar la transformación transformar la integral en la forma dónde , en cuyo punto se puede proceder de manera idéntica a la cuadratura de Clenshaw-Curtis para f como se indicó anteriormente. [12] Sin embargo, debido a las singularidades de los extremos en esta reasignación de coordenadas, se usa la primera regla de cuadratura de Fejér [que no evalúa f (−1)] a menos que g (∞) sea finito.
Calcular previamente los pesos en cuadratura
En la práctica, es inconveniente realizar una DCT de los valores de función muestreados f (cos θ) para cada nuevo integrando. En su lugar, normalmente se calculan previamente los pesos en cuadratura(para n de 0 a N / 2, asumiendo que N es par) de modo que
Estos pesos también se calculan mediante un DCT, como se ve fácilmente al expresar el cálculo en términos de álgebra matricial . En particular, calculamos los coeficientes de la serie de coseno a través de una expresión de la forma:
donde D es la forma matricial de la DCT de tipo I de ( N / 2 + 1) punto de arriba, con entradas (para índices de base cero ):
y es
Como se discutió anteriormente, debido al aliasing , no tiene sentido calcular coeficientes más allá de, entonces D es unmatriz. En términos de estos coeficientes c , la integral es aproximadamente:
desde arriba, donde c es el vector de coeficientesarriba yd es el vector de integrales para cada coeficiente de Fourier:
(Tenga en cuenta, sin embargo, que estos factores de ponderación se alteran si se cambia la matriz D de DCT para usar una convención de normalización diferente. Por ejemplo, es común definir la DCT de tipo I con factores adicionales de 2 o √ 2 factores en la primera y últimas filas o columnas, lo que conduce a las alteraciones correspondientes en las entradas d ). la suma se puede reorganizar para:
donde w es el vector de los pesos deseados arriba, con:
Dado que la matriz transpuestatambién es una DCT (por ejemplo, la transposición de una DCT de tipo I es una DCT de tipo I, posiblemente con una normalización ligeramente diferente según las convenciones que se empleen), los pesos de cuadratura w se pueden calcular previamente en O ( N log N ) tiempo para un N dado usando algoritmos DCT rápidos.
Los pesos son positivos y su suma es igual a uno. [13]
Ver también
- Fórmula de Euler-Maclaurin
- Fórmula de cuadratura de Gauss-Kronrod
Referencias
- ^ W. Morven Gentleman, "Implementación de la cuadratura I de Clenshaw-Curtis: Metodología y experiencia", Comunicaciones del ACM 15 (5), p. 337-342 (1972).
- ^ Jörg Waldvogel, " Construcción rápida de las reglas de cuadratura de Fejér y Clenshaw-Curtis ", BIT Numerical Mathematics 46 (1), p. 195-202 (2006).
- ^ CW Clenshaw y AR Curtis " Un método para la integración numérica en una computadora automática Numerische Mathematik 2 , 197 (1960).
- ^ JP Boyd, Chebychev y Fourier Spectral Methods , 2ª ed. (Dover, Nueva York, 2001).
- ^ Véase, por ejemplo, SG Johnson, " Notas sobre la convergencia de la cuadratura de la regla trapezoidal ", notas del curso en línea del MIT (2008).
- ^ Leopold Fejér, " Sobre las secuencias infinitas que surgen en las teorías del análisis armónico, de la interpolación y de las cuadraturas mecánicas ", Boletín de la American Mathematical Society 39 (1933), págs. 521–534. Leopold Fejér, " Mechanische Quadraturen mit positiven Cotesschen Zahlen , Mathematische Zeitschrift 37 , 287 (1933).
- ^ Trefethen, Lloyd N. (2008). "¿Es la cuadratura de Gauss mejor que la de Clenshaw-Curtis?". Revisión SIAM . 50 (1): 67–87. CiteSeerX 10.1.1.157.4174 . doi : 10.1137 / 060659831 .
- ^ Ignace Bogaert, Computación libre de iteraciones de Gauss - Nodos y pesos en cuadratura de Legendre, SIAM Journal on Scientific Computing vol. 36 , págs. A1008 – A1026 (2014)
- ^ Erich Novak y Klaus Ritter, "Alta integración dimensional de funciones suaves sobre cubos", Numerische Mathematik vol. 75 , págs. 79-97 (1996).
- ^ GA Evans y JR Webster, "Una comparación de algunos métodos para la evaluación de integrales altamente oscilatorias", Journal of Computational and Applied Mathematics , vol. 112 , pág. 55 - 69 (1999).
- ^ a b c d e John P. Boyd, " Esquemas de cuadratura exponencialmente convergentes de Fourier-Chebshev [ sic ] en intervalos acotados e infinitos", J. Scientific Computing 2 (2), p. 99-109 (1987).
- ^ Nirmal Kumar Basu y Madhav Chandra Kundu, "Algunos métodos de integración numérica en un intervalo semi-infinito", Aplicaciones de las matemáticas 22 (4), p. 237 - 243 (1977).
- ^ JP Imhof, "Sobre el método de integración numérica de Clenshaw y Curtis", Numerische Mathematik 5 , p. 138-141 (1963).