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 Modèle:Math 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 Modèle:Math à intégrer est évaluée aux Modèle:Mvar 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 Modèle:Math par des algorithmes de transformée de Fourier rapide adaptés à la TCDModèle:,.

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 Modèle:Math. L'algorithme est normalement exprimé pour l'intégration d'une fonction Modèle:Math sur l'intervalle [-1 ; 1] (tout autre intervalle peut être obtenu par un changement d'échelle approprié). On utilise alors l'égalité :

$\underset{-1}{\overset{1}{\int }}f(x)\mathrm{d}x=\underset{0}{\overset{\pi }{\int }}f(\cos \theta )\sin (\theta )\mathrm{d}\theta .$

En utilisant le développement en cosinus de Modèle:Math :

$f(\cos \theta )=\frac{a_0}{2}+\sum _{k=1}^{\mathrm{\infty }}{a}_k\cos (k\theta )$

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

$\underset{0}{\overset{\pi }{\int }}f(\cos \theta )\sin (\theta )\mathrm{d}\theta =a_0+\sum _{k=1}^{\mathrm{\infty }}\frac{2a_{2k}}{1-(2k{)}^2}.$

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

$a_k=\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}f(\cos \theta )\cos (k\theta )\mathrm{d}\theta ,k=0,1,2,\dots ,$

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 Modèle:Math, par construction), au-delà de la fréquence de Nyquist Modèle:Mvar, sont calculés précisément pour les Modèle:Math points Modèle:Math pour Modèle:Math, 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)Modèle:,. Ainsi, on approche l'intégrale de la série en cosinus par la transformée en cosinus discrète de type I :

$a_k\approx \frac{2}{N}[\frac{f(1)}{2}+\frac{f(-1)}{2}(-1{)}^k+\sum _{n=1}^{N-1}f(\cos [n\pi /N])\cos (nk\pi /N)]$

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

$a_{2k}\approx \frac{2}{N}[\frac{f(1)+f(-1)}{2}+f(0)(-1{)}^k+\sum _{n=1}^{N/2-1}\{f(\cos [n\pi /N])+f(-\cos [n\pi /N])\}\cos \left(\frac{nk\pi }{N/2}\right)]$

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

À cause de l'aliasing, on ne calcule les coefficients Modèle:Math 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 Modèle:Math sont les amplitudes de l'unique polynôme trigonométrique d'interpolation à bande limitée passant par les N+1 points où Modèle:Math 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 Modèle:Mvar 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) :

$\underset{0}{\overset{\pi }{\int }}f(\cos \theta )\sin (\theta )\mathrm{d}\theta \approx a_0+\sum _{k=1}^{N/2-1}\frac{2a_{2k}}{1-(2k{)}^2}+\frac{a_N}{1-N^2}.$

Connexion aux polynômes de Tchebychev

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

$f(x)=\frac{a_0}{2}{T}_0(x)+\sum _{k=1}^{+\mathrm{\infty }}{a}_k{T}_k(x),$

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

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 Modèle:Math 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ù Modèle:Math est paire et périodique en Modèle:Mvar par construction, et de classe CModèle:Exp partout si Modèle:Math est de classe CModèle:Exp sur [-1 ; 1]. (A l'inverse, appliquer directement un développement en série de cosinus à Modèle:Math au lieu de Modèle:Math 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 Modèle:Math et Modèle:Math 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 Modèle:Mvar en calculant Modèle:Math sur un ensemble de points également espacés, à mi-chemin entre les extrema : ${\theta }_n=(n+0,5)\pi /N,0\le n<N$. Ce sont les racines de Modèle:Math, 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 :

$a_k\approx \frac{2}{N}\sum _{n=0}^{N-1}f(\cos \left[\frac{(2n+1)\pi }{2N}\right])\cos \left(\frac{(2n+1)k\pi }{2N}\right)$

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 Modèle:Math points et les utilise pour obtenir une valeur exacte pour des polynômes de degré au plus Modèle:Math. En revanche, la quadrature de Clenshaw–Curtis n'est construite que pour atteindre l'égalité pour des polynômes de degré au plus Modèle:Mvar. 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 Modèle:Math 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 Modèle:Math 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 Modèle:Math pour les calculer. Cependant, il existe aujourd'hui des algorithmes atteignent une complexité en Modèle:Math 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 Modèle:Mvar 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 Modèle:Mvar 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 Modèle:Mvar sont un sous-ensemble des points de l'ordre Modèle:Math. 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 Modèle:Math fixe et connue :

$\underset{-1}{\overset{1}{\int }}f(x)w(x)\mathrm{d}x=\underset{0}{\overset{\pi }{\int }}f(\cos \theta )w(\cos \theta )\sin (\theta )\mathrm{d}\theta .$

Le cas commun est la fonction constante Modèle:Math, mais dans certaines applications, une fonction poids différente peut être utile. Une raison basique est que, puisque Modèle:Math 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 Modèle:Math, 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

$W_k=\underset{0}{\overset{\pi }{\int }}w(\cos \theta )\cos (k\theta )\sin (\theta )\mathrm{d}\theta .$

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 Modèle:Mvar numériquement pour les cas de haute précision. De plus que, puisque Modèle:Math 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 Modèle:Math avec une fonction poids Modèle:Math 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/Modèle:Racine). 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 Modèle:Mvar.

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 Modèle:Math, à 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 $\underset{0}{\overset{\mathrm{\infty }}{\int }}f(x)\mathrm{d}x$ et $\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}f(x)\mathrm{d}x$, 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 Modèle:Math ait une décroissance suffisante forte pour $\left|x\right|$ tendant vers l'infini.

Une possibilité est d'utiliser une transformation de coordonnées telle que Modèle:Math

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}f(x)\mathrm{d}x=\underset{-1}{\overset{1}{\int }}f\left(\frac{t}{1-t^2}\right)\frac{1+t^2}{(1-t^2{)}^2}\mathrm{d}t,$

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 Modèle:Math, où Modèle:Mvar est une constante à spécifier ou optimiser pour accélérer la convergence selon le problème), ce qui donne alors :

$\underset{0}{\overset{\mathrm{\infty }}{\int }}f(x)\mathrm{d}x=2L\underset{0}{\overset{\pi }{\int }}\frac{f[L{\cot }^2(\theta /2)]}{[1-\cos (\theta ){]}^2}\sin (\theta )\mathrm{d}\theta .$

Le facteur multipliant Modèle:Math, Modèle:Math, 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 Modèle:Math. Afin d'éliminer la singularité en Modèle:Math de cet intégrande, il suffit que Modèle:Math tende vers 0 suffisamment vers quand x tend vers l'infini, au moins aussi vite que Modèle:Math.

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

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}f(x)\mathrm{d}x=L\underset{0}{\overset{\pi }{\int }}\frac{f[L\cot (\theta )]}{{\sin }^2(\theta )}\mathrm{d}\theta \approx \frac{L\pi }{N}\sum _{n=1}^{N-1}\frac{f[L\cot (n\pi /N)]}{{\sin }^2(n\pi /N)}.$

Dans ce cas, on utilise le fait que le terme Modèle:Math 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ù Modèle:Mvar 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 Modèle:Math tende vers 0 plus vite que Modèle:Math en Modèle:Math (si Modèle:Mvar décroit en Modèle:Math, 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 Modèle:Mvar 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 Modèle:Mvar : le problème est que, même si Modèle:Math est effectivement Modèle:Math-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 Modèle:Math décroit en Modèle:Math mais a une discontinuité de saut dans la pente de la fonction remappée en Modèle:Math et Modèle:Math].

Un autre changement de variable a été suggéré pour les intégrales de la forme $\underset{0}{\overset{\mathrm{\infty }}{\int }}{\mathrm{e}}^{-x}g(x)\mathrm{d}x$, en utilisant la transformation Modèle:Math pour transformer l'intégrale vers la forme $\underset{0}{\overset{\pi }{\int }}f(\cos \theta )\sin \theta \mathrm{d}\theta$ avec Modèle:Math, à partir de là on peut réutiliser le schéma de la quadrature de Clenshaw–Curtis pour Modèle:Mvar 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 Modèle:Math] sauf si Modèle:Math 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 Modèle:Math pour chaque nouvel intégrande. A la place, on pré-calcule les poids de quadrature Modèle:Mvar (pour n de 0 à N/2, en supposant N pair) de sorte que

$\underset{-1}{\overset{1}{\int }}f(x)\mathrm{d}x\approx \sum _{n=0}^{N/2}{w}_n[f(\cos \left(\frac{n\pi }{N}\right))+f(-\cos \left(\frac{n\pi }{N}\right))].$

Ces poids Modèle:Mvar 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 Modèle:Math par une expression de la forme :

$c=\left(\begin{array}{c}{a}_0\\ a_2\\ a_4\\ ⋮\\ a_N\end{array}\right)=D\left(\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_{N/2}\end{array}\right)=Dy,$

où Modèle:Mvar est la matrice de la TCD de type I en Modèle:Math points, avec :

$D_{k,n}=\frac{2}{N}\cos \left(\frac{nk\pi }{N/2}\right)\times \{\begin{array}{cc}1/2& n=0,N/2\\ 1& \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{o}\mathrm{n}\end{array}$

et Modèle:Mvar tel que :

$y_n=f(\cos \left(\frac{n\pi }{N}\right))+f(-\cos \left(\frac{n\pi }{N}\right)).$

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

$\underset{-1}{\overset{1}{\int }}f(x)\mathrm{d}x\approx a_0+\sum _{k=1}^{N/2-1}\frac{2a_{2k}}{1-(2k{)}^2}+\frac{a_N}{1-N^2}=d^{T}c,$

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

$d=\left(\begin{array}{c}1\\ 2/(1-4)\\ 2/(1-16)\\ ⋮\\ 2/(1-[N-2{]}^2)\\ 1/(1-N^2)\end{array}\right).$

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

$\underset{-1}{\overset{1}{\int }}f(x)\mathrm{d}x\approx d^{T}c=d^{T}Dy=(D^{T}d{)}^{T}y=w^{T}y$

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

$w=D^{T}d.$

Puisque la matrice transposée Modèle:Mvar 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 Modèle:Mvar peuvent être pré-calculés en un temps en Modèle:Math pour un Modèle:Mvar par des algorithmes de transformées de Fourier rapide.

Les poids Modèle:Mvar 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