La cuadratura de Clenshaw-Curtis y la cuadratura de Fejér son métodos de 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 usan una aproximación de transformada discreta del coseno (DCT) para la serie de cosenos . Además de tener una precisión de convergencia rápida comparable a las reglas de cuadratura gaussiana , la cuadratura de Clenshaw-Curtis conduce naturalmente a reglas de cuadratura anidadas (donde diferentes órdenes de precisión comparten puntos), lo que es importante tanto para la cuadratura adaptativa como para la cuadratura multidimensional ( cubatura ).
Brevemente, la función a integrar se evalúa en los extremos o raíces de un polinomio de Chebyshev y estos valores se utilizan para construir una aproximación polinómica para la función. Luego, este polinomio se integra exactamente. En la práctica, los pesos de integración para el valor de la función en cada nodo se calculan previamente y este cálculo se puede realizar a tiempo mediante algoritmos relacionados con la transformada rápida de Fourier para la DCT. [1] [2]
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 ) en el intervalo [−1,1] (cualquier otro intervalo se puede obtener mediante un reescalado apropiado). Para esta integral, podemos escribir:
Es decir, hemos transformado el problema de integración a uno de integración . Esto se puede realizar si conocemos la serie de cosenos para :
En cuyo caso la integral se convierte en:
Por supuesto, para calcular los coeficientes de la serie del coseno
es necesario realizar nuevamente 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 de series de Fourier para funciones periódicas (como , por construcción), hasta la frecuencia de Nyquist , se calculan con precisión mediante los puntos igualmente espaciados e igualmente ponderados para (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 mediante la transformada discreta del coseno (DCT) de tipo I :
para y luego use la fórmula anterior para la integral en términos de estos . Como solo se necesita , la fórmula se simplifica aún más en una DCT de tipo I de orden N /2 , suponiendo 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 ) y f (− x ) por igual.
Debido al aliasing , solo se calculan los coeficientes hasta 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, son las amplitudes del único polinomio de interpolación trigonométrica de banda limitada que pasa por los N +1 puntos donde se evalúa f (cos θ ) , y aproximamos la integral mediante la integral de este polinomio de interpolación. Sin embargo, hay cierta sutileza en cómo se trata el coeficiente en la integral: para evitar el doble conteo con su alias, se incluye con el peso 1/2 en la integral aproximada final (como también se puede ver al examinar el polinomio de interpolación):
La razón por la que esto está conectado con los polinomios de Chebyshev es que, por definición, , y por lo tanto la serie de cosenos anterior es realmente una aproximación de por los polinomios de Chebyshev:
y por lo tanto estamos "realmente" integrando al integrar su expansión aproximada en términos de polinomios de Chebyshev. Los puntos de evaluación corresponden a los extremos del polinomio de Chebyshev .
El hecho de que dicha aproximación de Chebyshev sea simplemente una serie de cosenos bajo un cambio de variables es responsable de la rápida convergencia de la aproximación a medida que se incluyen más términos. Una serie de cosenos converge muy rápidamente para funciones que son pares , periódicas y suficientemente suaves. Esto es cierto aquí, ya que es par y periódica en por construcción, y es k -veces diferenciable en todas partes si es k -veces diferenciable en . (Por el contrario, aplicar directamente una expansión de serie de cosenos a en lugar de por lo general no convergerá rápidamente porque la pendiente de la extensión par-periódica generalmente sería discontinua).
Fejér propuso dos reglas de cuadratura muy similares a la cuadratura de Clenshaw-Curtis, pero mucho antes (en 1933). [6]
De estas 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 extremos y se establecen en cero. Es decir, Fejér solo utilizó 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 mediante la evaluación en un conjunto diferente de puntos igualmente espaciados, a medio camino entre los extremos: para . Estas son las raíces 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 traslacional de la serie periódica de Fourier). Esto conduce a una fórmula:
que es precisamente la DCT de 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 de 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 Clenshaw-Curtis" se ha convertido en estándar.
El método clásico de cuadratura gaussiana evalúa el integrando en los puntos y está construido para integrar exactamente polinomios hasta el grado . En cambio, la cuadratura de Clenshaw-Curtis, antes mencionada, evalúa el integrando en los puntos e integra exactamente polinomios solo hasta el grado . Por lo 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 la cuadratura de 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 (véase aproximación de Chebyshev ). De hecho, resultados teóricos recientes [7] sostienen que tanto la cuadratura gaussiana como la de Clenshaw-Curtis tienen un error acotado por para 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 el tiempo mediante algoritmos de transformada rápida de Fourier (o sus análogos para la DCT), mientras que la mayoría de los algoritmos para los pesos de cuadratura gaussiana requerían tiempo para calcularlos. Sin embargo, los algoritmos recientes han alcanzado complejidad 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 . En cambio, generalmente se emplea un esquema de cuadratura adaptativa que primero evalúa la integral a 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 del integrando. Esto se denomina regla de cuadratura anidada, y aquí Clenshaw–Curtis tiene la ventaja de que la regla para el orden N utiliza un subconjunto de los puntos de orden 2 N . Por el contrario, las reglas de cuadratura gaussiana 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]
De manera más general, se puede plantear el problema de integrar una función de peso arbitraria frente a una función de peso fija que se conoce de antemano:
El caso más común es , como el anterior, pero en ciertas aplicaciones es deseable una función de peso diferente. La razón básica es que, dado que se puede tener en cuenta a priori , se puede hacer que el error de integración dependa solo de la precisión en la aproximación de , independientemente de lo mal que se comporte la función de peso.
La cuadratura de Clenshaw-Curtis se puede generalizar a este caso de la siguiente manera. Como antes, funciona hallando la expansión en serie de cosenos de mediante una DCT y luego integrando cada término en la serie de cosenos. Ahora, sin embargo, estas integrales tienen la forma
Para la mayoría de las , esta integral no se puede calcular analíticamente, a diferencia de lo que sucedía antes. Sin embargo, dado que generalmente se utiliza la misma función de ponderación para muchos integrandos , se puede calcular de antemano de forma numérica con gran precisión. 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 una función de peso que es altamente oscilatoria, por ejemplo, una función sinusoide o de Bessel (véase, por ejemplo, Evans & Webster, 1999 [10] ). Esto es útil para el cálculo de series de Fourier y series de Fourier-Bessel de alta precisión , donde los métodos de cuadratura simples 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 a través de 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 cuando el integrando es desconocido pero tiene una singularidad conocida de algún tipo, por ejemplo, una discontinuidad conocida o una divergencia integrable (como 1/ √ x ) en algún punto. En este caso, la singularidad se puede incorporar a la función de ponderación y sus propiedades analíticas se pueden utilizar para realizar cálculos precisos 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.
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] Se puede conservar una alta precisión, incluso la convergencia exponencial para integrandos suaves, siempre que decaiga lo suficientemente rápido a medida que | 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 semiinfinito 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 (uno 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 entonces expandirse en una serie de cosenos (aproximadamente, usando la transformada discreta del coseno) e integrarse término por término, exactamente como se hizo para f (cos θ ) anteriormente. Para eliminar la singularidad en θ = 0 en este integrando, uno simplemente requiere que f ( x ) vaya a cero lo suficientemente rápido a medida que x se acerca 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 se indicó anteriormente) para transformar la integral en: [11]
En este caso, hemos utilizado el hecho de que el integrando reasignado f ( L cot θ )/sin 2 ( θ ) ya es periódico y, por lo tanto, se puede integrar directamente con alta precisión (incluso exponencial) utilizando la regla trapezoidal (asumiendo que f es suficientemente suave y decae rápidamente); no hay necesidad de calcular la serie de cosenos como un paso intermedio. Tenga en cuenta que la regla de cuadratura no incluye los puntos finales, donde hemos asumido que el integrando tiende a cero. La fórmula anterior requiere que f ( x ) decaiga más rápido que 1/ x 2 a medida que x tiende a ±∞. (Si f decae exactamente como 1/ x 2 , entonces el integrando tiende 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 sólo polinomialmente rápido, entonces puede ser necesario usar un paso adicional de cuadratura de Clenshaw-Curtis para obtener precisión exponencial de la integral remapeada 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 θ )/sin 2 ( θ ) es de hecho periódica con período π, no es necesariamente suave en los puntos finales si todas las derivadas no se desvanecen 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 remapeada en θ=0 y π].
Se sugirió otro enfoque de reasignación de coordenadas para integrales de la forma , en cuyo caso se puede usar la transformación para transformar la integral en la forma donde , 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 puntos finales 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.
En la práctica, no resulta práctico realizar una DCT de los valores de la función muestreada f (cos θ) para cada nuevo integrando. En cambio, normalmente se calculan previamente los pesos de cuadratura (para n desde 0 hasta N /2, suponiendo que N es par) de modo que
Estos pesos también se calculan mediante una DCT, como se puede ver fácilmente al expresar el cálculo en términos de álgebra matricial . En particular, calculamos los coeficientes de la serie de cosenos mediante una expresión de la forma:
donde D es la forma matricial de la DCT de tipo I de ( N /2+1) puntos de arriba, con entradas (para índices basados en cero ):
y es
Como se explicó anteriormente, debido al aliasing , no tiene sentido calcular coeficientes más allá de , por lo que D es una matriz. En términos de estos coeficientes c , la integral es aproximadamente:
desde arriba, donde c es el vector de coeficientes anteriores y d es el vector de integrales para cada coeficiente de Fourier:
(Sin embargo, tenga en cuenta que estos factores de peso se modifican si se cambia la matriz DCT D para utilizar 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 las primeras y últimas filas o columnas, lo que conduce a modificaciones correspondientes en las entradas d ). La suma se puede reorganizar a:
donde w es el vector de los pesos deseados arriba, con:
Dado que la matriz transpuesta también es una DCT (por ejemplo, la transpuesta de una DCT de tipo I es una DCT de tipo I, posiblemente con una normalización ligeramente diferente dependiendo de las convenciones que se empleen), los pesos de cuadratura w se pueden precalcular en tiempo O ( N log N ) para un N dado utilizando algoritmos DCT rápidos.
Los pesos son positivos y su suma es igual a uno. [13]