Méthode de quadrature de Clenshaw-Curtis
En analyse, les méthodes de quadrature de Clenshaw–Curtis et de quadrature de Fejér sont des méthodes d'intégration numérique s'appuyant sur le développement de la fonction à intégrer en polynômes de Tchebychev. De façon équivalente, ils emploient un changement de variable  et utilisent une approximation de la transformée en cosinus discrète pour un développement en cosinus. En plus d'avoir des résultats de convergence rapide comparables à la quadrature de Gauss, la quadrature de Clenshaw–Curtis mène naturellement à des méthodes de quadrature imbriquée (où des points se retrouvent dans plusieurs ordres de précision), ce qui devient intéressant pour la quadrature adaptative et les méthodes de quadrature multidimensionnelles (cubature).

En résumé, la fonction  à intégrer est évaluée aux  extrema ou racines d'un polynôme de Tchebychev et ces valeurs sont utilisées pour construire une approximation polynomiale de la fonction. Ce polynôme est ensuite intégré de façon exacte. En pratique, les poids d'intégration en chaque nœud sont pré-calculés, en un temps en  par des algorithmes de transformée de Fourier rapide adaptés à la TCD.

Méthode générale 

Une façon de comprendre l'algorithme est de voir que la quadrature de Clenshaw–Curtis (telle que proposée par ses auteurs en 1960) revient à intégrer par un changement de variable . L'algorithme est normalement exprimé pour l'intégration d'une fonction  sur l'intervalle [-1 ; 1] (tout autre intervalle peut être obtenu par un changement d'échelle approprié). On utilise alors l'égalité :

En utilisant le développement en cosinus de  :

ce qui permet de réécrire l'intégrale par inversion série-intégrale :

Afin de calculer les coefficients de la série en cosinus

on utilise souvent une intégration numérique, ce qui ne semble pas indiquer une simplification du problème. Cependant, l'intégration des séries de Fourier de fonctions périodiques (comme , par construction), au-delà de la fréquence de Nyquist , sont calculés précisément pour les  points  pour , uniformément pondérés (sauf les extrémités 1/2, afin d'éviter le double comptage, comme pour la méthode des trapèzes ou la formule d'Euler-Maclaurin). Ainsi, on approche l'intégrale de la série en cosinus par la transformée en cosinus discrète de type I :

pour  et on utilise cette formule précédente pour l'intégrale en termes de ces . Puisque seul  est nécessaire, la formule se simplifie d'autant en une TCD de type I d'ordre N/2, si  est un nombre pair :

De cette formule, il apparait clairement que la quadrature de Clenshaw–Curtis est une formule symétrique, dans le sens où il pondère  et  également.

À cause de l'aliasing, on ne calcule les coefficients  que jusqu'à k=N/2, car l'échantillonnage discret de la fonction rend la fréquence de 2k indistinguable de celle de N–2k. De façon équivalente, les  sont les amplitudes de l'unique polynôme trigonométrique d'interpolation à bande limitée passant par les N+1 points où  est évalué, et on approche l'intégrale par l'intégrale de ce polynôme d'interpolation. Il y a une certaine subtilité dans la façon de considérer le coefficient  dans l'intégrale, cependant — pour éviter un double comptage dans le signal — il est inclus avec un poids moitié dans l'intégrale approchée finale (ce qui peut être vu aussi en examinant le polynôme d'interpolation) :

Connexion aux polynômes de Tchebychev 

En utilisant l'égalité caractérisant les polynômes de Tchebychev , soit , on peut voir la série de cosinus comme une approximation de  en polynômes de Tchebychev :

et l'intégration de  devient "exacte" en intégrant le développement approchant en termes de polynômes de Tchebychev. Les points d'évaluation  correspondent aux extrema du polynôme .

Le fait que l'approximation de Tchebychev est une simple série de cosinus par un changement de variable permet la convergence rapide de l'approximation car plus de termes de  sont inclus. Une série de cosinus converge très rapidement pour des fonctions paires, périodiques et assez lisses. C'est le cas ici, dans le sens où  est paire et périodique en  par construction, et de classe C partout si  est de classe C sur [-1 ; 1]. (A l'inverse, appliquer directement un développement en série de cosinus à  au lieu de  ne converge pas dans la plupart des cas rapidement car la partie paire du développement est en général discontinu).

Quadrature de Fejér 

Fejér a proposé deux méthodes de quadrature très similaires à la quadrature de Clenshaw–Curtis, mais plus tôt (vers 1933).

Sur les deux, la deuxième méthode de quadrature de Fejér est très proche de celle de Clenshaw–Curtis, à la seule différence que les extrémités  et  sont fixées à zéro et qu'elle n'utilise que les extrema intérieurs des polynômes de Tchebyschev, soit les vrais points stationnaires.

La première méthode de quadrature de Fejér évalue les coefficients  en calculant  sur un ensemble de points également espacés, à mi-chemin entre les extrema : . Ce sont les racines de , parfois appelés "nœuds de Tchebychev". Ces points sont le seul autre choix permettant de conserver la symétrie paire de la transformée en cosinus et la symétrie de translation de la série de Fourier périodique. On obtient ainsi la formule :

soit la TCD de type II. Cependant, cette quadrature n'est pas localisable : les points d'évaluation pour 2N ne coïncident avec aucun des points d'évaluation pour N, contrairement à la quadrature de Clenshaw–Curtis et la deuxième méthode de Fejér.

Bien que Fejér ait publié ces techniques avant Clenshaw et Curtis, le nom de "quadrature de Clenshaw–Curtis" est plus courant.

Comparaison à la quadrature gaussienne 

La méthode classique de la quadrature de Gauss évalue l'intégrande en  points et les utilise pour obtenir une valeur exacte pour des polynômes de degré au plus . En revanche, la quadrature de Clenshaw–Curtis n'est construite que pour atteindre l'égalité pour des polynômes de degré au plus . On pourrait alors considérer que la version de Clenshaw–Curtis est intrinsèquement moins bonne que la Gaussienne, mais la pratique montre que cet a priori semble faux.

Plusieurs auteurs ont pu observer par l'application que la méthode de Clenshaw–Curtis a une précision comparable à celle de Gauss pour un même nombre de points. La cause est que plusieurs intégrandes numériques ne sont pas des polynômes (d'autant que les polynômes peuvent être intégrés analytiquement), et l'approximation de nombreuses fonctions en termes de polynômes de Tchebychev converge rapidement (par l'approximation de Tchebychev). En fait, des résultats théoriques récents avancent que les quadratures de Gauss et de Clenshaw–Curtis ont toutes deux des erreurs en  pour une intégrande k-fois différentiable.

Un avantage souvent cité de la quadrature de Clenshaw–Curtis est que les pondérations peuvent être évaluées en un temps en  par un algorithme de transformée de Fourier rapide (ou les analogues pour la DCT), alors que la plupart des algorithmes pour les poids d'une quadrature gaussienne quadrature sont en un temps en  pour les calculer. Cependant, il existe aujourd'hui des algorithmes atteignent une complexité en  pour une quadrature de Gauss–Legendre. Dans les faits, une intégration numérique d'ordre élevé est rarement effectuée par une simple évaluation d'une formule de quadrature pour un  très grand. On préférera plutôt employer une quadrature adaptative qui évalue d'abord l'intégrale à un ordre petit avant d'augmenter le nombre de points de calcul, éventuellement dans certaines zones où l'erreur sur l'intégrale est plus grande. Afin d'évaluer la précision sur la quadrature, on compare la réponse avec celle obtenue par une règle de quadrature d'ordre plus faible. Idéalement, cette quadrature d'ordre plus faible évalue l'intégrande sur un sous-ensemble des  points originaux, afin de minimiser les évaluations des intégrandes. On parle ainsi de quadrature imbriquée, et la règle de Clenshaw–Curtis a l'avantage que les points utilisés pour l'ordre  sont un sous-ensemble des points de l'ordre . Ce n'est pas le cas en général pour la quadrature gaussienne, qu'on peut modifier par la méthode de quadrature de Gauss–Kronrod, par exemple. Les règles d'imbrication sont aussi importantes pour les grilles creuses en quadrature multidimensionnelles, et la quadrature de Clenshaw–Curtis est populaire dans ce contexte.

Intégration avec des fonctions poids 

On peut également considérer le cas de l'intégration d'une fonction avec une fonction poids donnée  fixe et connue :

Le cas commun est la fonction constante , mais dans certaines applications, une fonction poids différente peut être utile. Une raison basique est que, puisque  peut être prise en compte a priori, l'erreur d'intégration peut être générée pour ne dépendre que de la précision sur l'approximation de , peu importe le comportement de la fonction poids.

La quadrature de Clenshaw–Curtis peut donc être généralisée, sur la même méthode, mais on en arrive à calculer les intégrales

Dans le cas général, une intégrale de ce genre ne peut pas être calculée analytiquement, ce qui était le cas pour la version simple. Comme la même fonction poids est utilisée pour plusieurs intégrandes, on peut envisager de calculer les  numériquement pour les cas de haute précision. De plus que, puisque  est généralement spécifiée analytiquement, il est possible dans certains cas d'utiliser des méthodes spécialisées pour calculer les coefficients.

Par exemple, des méthodes spéciales ont été développées pour la quadrature de Clenshaw–Curtis aux fonctions de la forme  avec une fonction poids  hautement oscillatoire, e.g. une sinusoïde ou une fonction de Bessel. On obtient alors des calculs de séries de Fourier et de Fourier-Bessel de haute précision, alors que la version classique aurait nécessité un ordre élevé pour prendre en compte l'impact des oscillations rapides.

Un autre cas où des fonctions poids ont un intérêt si l'intégrande est inconnue mais a une singularité, e.g. une discontinuité ponctuelle connue, ou pour le cas d'une intégrale impropre (comme 1/). Dans ce cas, la singularité peut être intégrée à la fonction poids et ses propriétés analytiques seront utilisés pour le calcul précis des .

La quadrature de Gauss peut également être adaptée à des intégrales avec poids, mais la technique a ses différences. Dans la quadrature de Clenshaw–Curtis, l'intégrande est toujours évalué aux mêmes points peu importe le choix de , à savoir les extrema d'un polynôme de Tchebychev, alors qu'adapter la quadrature de Gauss amènent à l'utilisation des racines des polynômes orthogonaux associés au poids choisi.

Intégration sur des intervalles semi-infinis ou infinis 

Il est également possible d'adapter la quadrature de Clenshaw–Curtis pour le calcul d'intégrales de la forme  et , par une technique de remappage. Les propriétés de précision et de convergence exponentielle pour des intégrandes lisses peuvent être conservés pourvu que  ait une décroissance suffisante forte pour  tendant vers l'infini.

Une possibilité est d'utiliser une transformation de coordonnées telle que 

afin de ramener une intégrale sur un intervalle semi-fini ou infini vers un fini. Des techniques ont spécialement été développées pour la quadrature de Clenshaw–Curtis.

Par exemple, on peut utiliser le remappage , où  est une constante à spécifier ou optimiser pour accélérer la convergence selon le problème), ce qui donne alors :

Le facteur multipliant , , peut être développé en série de cosinus (de façon approchée, par une TCD) et intégré terme à terme, exactement comme pour le terme . Afin d'éliminer la singularité en  de cet intégrande, il suffit que  tende vers 0 suffisamment vers quand x tend vers l'infini, au moins aussi vite que .

Pour un intervalle infini d'intégration, on peut utiliser le changement de variable , ce qui mène à:

Dans ce cas, on utilise le fait que le terme  est déjà périodique et que la méthode des trapèzes permet d'obtenir une estimation très précise (même exponentielle) de l'intégration (dans le cas où  est suffisamment lisse et décroit rapidement) ; on évite ainsi le calcul intermédiaire d'un développement en série de cosinus. On remarquera que les extrémités ne sont pas dans les points de calcul, où on a supposé que l'intégrande tend vers 0. La formule précédente nécessite que  tende vers 0 plus vite que  en  (si  décroit en , alors l'intégrande tend vers une valeur finie aux extrémités et ces limites doivent être incluses dans les points de calcul dans la méthode des trapèzes). Cependant, si  a seulement une décroissance polynomiale, alors il peut être nécessaire d'utiliser une étape supplémentaire dans la quadrature de Clenshaw–Curtis pour retrouver la précision exponentielle de l'intégrale remappée au lieu de la règle des trapèzes, selon les propriétés de  : le problème est que, même si  est effectivement -périodique, il n'est pas nécessairement lisse aux extrémités si toutes les dérivées ne s'y annulent pas [e.g. la fonction  décroit en  mais a une discontinuité de saut dans la pente de la fonction remappée en  et ].

Un autre changement de variable a été suggéré pour les intégrales de la forme , en utilisant la transformation  pour transformer l'intégrale vers la forme  avec , à partir de là on peut réutiliser le schéma de la quadrature de Clenshaw–Curtis pour  comme vue plus haut. À cause des singularités aux extrémités dans le changement de variable, on peut préférer la première méthode de quadrature de Fejér [sans l'évaluation de ] sauf si  est fini.

Pré-calcul des poids de quadrature 

Dans la pratique, il est peu commode de réaliser une TCD dans la fonction échantillonnée  pour chaque nouvel intégrande. A la place, on pré-calcule les poids de quadrature  (pour n de 0 à N/2, en supposant N pair) de sorte que

Ces poids  sont aussi calculés via une TCD, comme il peut être vu facilement en exprimant le calcul avec les outils de l'algèbre linéaire. En particulier, on calcule les coefficients de séries de cosinus  par une expression de la forme :

où  est la matrice de la TCD de type I en  points, avec :

et  tel que :

Comme vu auparavant, à cause de l'aliasing, il n'y a aucun intérêt à calculer les coefficients au-delà de , donc  est une matrice . En termes de ces coefficients , l'intégrale est approchée par :

où  est le vecteur des coefficients  et  est le vecteur des intégrales de chaque coefficient de Fourier :

On note toutefois que ces facteurs de poids sont altérés si on change la matrice de la TCD  pour utiliser une convention de normalisation différente. Par exemple, il est commun de définir la TCD de type  avec des facteurs additionnels de 2 ou  dans les première et dernière lignes ou colonnes, d'où les valeurs indiquées dans la définition de ). La sommation  peut être réarrangée en :

où  est le vecteur des poids désirés , avec:

Puisque la matrice transposée  est aussi une TCD (e.g., la transposée d'une TCD de type I est aussi une TCD de type I, possiblement avec une normalisation différente selon les conventions choisies), les poids de quadrature  peuvent être pré-calculés en un temps en  pour un  par des algorithmes de transformées de Fourier rapide.

Les poids  sont positifs et leur somme vaut 1.

Voir aussi 
 Formule d'Euler-Maclaurin
 Méthode de quadrature de Gauss-Kronrod
 Transformée en cosinus discrète

Références 

Quadrature de Clenshaw-Curtis