Méthode de quadrature de Clenshaw-Curtis

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)[3] revient à intégrer par un changement de variable x = cos(θ). L'algorithme est normalement exprimé pour l'intégration d'une fonction f(x) sur l'intervalle [-1 ; 1] (tout autre intervalle peut être obtenu par un changement d'échelle approprié). On utilise alors l'égalité :

${\displaystyle \int _{-1}^{1}f(x)\,\mathrm {d} x=\int _{0}^{\pi }f(\cos \theta )\sin(\theta )\,\mathrm {d} \theta .}$

En utilisant le développement en cosinus de f(cos θ) :

${\displaystyle f(\cos \theta )={\frac {a_{0}}{2}}+\sum _{k=1}^{\infty }a_{k}\cos(k\theta )}$

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

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

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

${\displaystyle a_{k}={\frac {2}{\pi }}\int _{0}^{\pi }f(\cos \theta )\cos(k\theta )\,\mathrm {d} \theta ,\quad 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 f(cos θ), par construction), au-delà de la fréquence de Nyquist k = N, sont calculés précisément pour les N+1 points θ_n = nπ/N pour n = 0 , ... , N, 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)[4] - [5]. Ainsi, on approche l'intégrale de la série en cosinus par la transformée en cosinus discrète de type I :

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

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

${\displaystyle a_{2k}\approx {\frac {2}{N}}\left[{\frac {f(1)+f(-1)}{2}}+f(0)(-1)^{k}+\sum _{n=1}^{N/2-1}\left\{f(\cos[n\pi /N])+f(-\cos[n\pi /N])\right\}\cos \left({\frac {nk\pi }{N/2}}\right)\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 f(x) et f(−x) également.

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

${\displaystyle \int _{0}^{\pi }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}}}.}$

En utilisant l'égalité caractérisant les polynômes de Tchebychev T_k(x), soit T_k(cos θ) = cos(k θ), on peut voir la série de cosinus comme une approximation de f(x) en polynômes de Tchebychev :

${\displaystyle f(x)={\frac {a_{0}}{2}}T_{0}(x)+\sum _{k=1}^{+\infty }a_{k}T_{k}(x),}$

et l'intégration de f(x) devient "exacte" en intégrant le développement approchant en termes de polynômes de Tchebychev. Les points d'évaluation x_n = cos(nπ/N) correspondent aux extrema du polynôme T_N(x).

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 T_k(x) 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ù f(cos θ) est paire et périodique en θ par construction, et de classe C^k partout si f(x) est de classe C^k sur [-1 ; 1]. (A l'inverse, appliquer directement un développement en série de cosinus à f(x) au lieu de f(cos θ) ne converge pas dans la plupart des cas rapidement car la partie paire du développement est en général discontinu).

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

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 f(-1) et f(1) 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 a_k en calculant f(cos θ) sur un ensemble de points également espacés, à mi-chemin entre les extrema : ${\displaystyle \theta _{n}=(n+0,5)\pi /N,0\leq n<N}$. Ce sont les racines de T_N(cos θ), 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 :

${\displaystyle a_{k}\approx {\frac {2}{N}}\sum _{n=0}^{N-1}f\left(\cos \left[{\frac {(2n+1)\pi }{2N}}\right]\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.

La méthode classique de la quadrature de Gauss évalue l'intégrande en N + 1 points et les utilise pour obtenir une valeur exacte pour des polynômes de degré au plus 2N + 1. En revanche, la quadrature de Clenshaw–Curtis n'est construite que pour atteindre l'égalité pour des polynômes de degré au plus N. 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[7] avancent que les quadratures de Gauss et de Clenshaw–Curtis ont toutes deux des erreurs en O([2N]^-k/k) 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 O(N log N) 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 O(N²) pour les calculer. Cependant, il existe aujourd'hui des algorithmes atteignent une complexité en O(N) pour une quadrature de Gauss–Legendre[8]. 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 N 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 N 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 N sont un sous-ensemble des points de l'ordre 2N. 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[9].

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

${\displaystyle \int _{-1}^{1}f(x)w(x)\,\mathrm {d} x=\int _{0}^{\pi }f(\cos \theta )w(\cos \theta )\sin(\theta )\,\mathrm {d} \theta .}$

Le cas commun est la fonction constante w(x) = 1, mais dans certaines applications, une fonction poids différente peut être utile. Une raison basique est que, puisque w(x) 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 f(x), 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

${\displaystyle W_{k}=\int _{0}^{\pi }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 W_k numériquement pour les cas de haute précision. De plus que, puisque w(x) 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 f(x) w(x) avec une fonction poids w(x) hautement oscillatoire, e.g. une sinusoïde ou une fonction de Bessel[10]. 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/√x). 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 W_k.

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 w(x), à 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.

Il est également possible d'adapter la quadrature de Clenshaw–Curtis pour le calcul d'intégrales de la forme ${\displaystyle \int _{0}^{\infty }f(x)\,\mathrm {d} x}$ et ${\displaystyle \int _{-\infty }^{\infty }f(x)\,\mathrm {d} x}$, par une technique de remappage[11]. Les propriétés de précision et de convergence exponentielle pour des intégrandes lisses peuvent être conservés pourvu que f(x) ait une décroissance suffisante forte pour ${\displaystyle \left|x\right|}$ tendant vers l'infini.

Une possibilité est d'utiliser une transformation de coordonnées telle que x = t/(1−t²)

${\displaystyle \int _{-\infty }^{\infty }f(x)\,\mathrm {d} x=\int _{-1}^{1}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 x = L cot²(θ/2), où L est une constante à spécifier ou optimiser pour accélérer la convergence selon le problème[11]), ce qui donne alors :

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

Le facteur multipliant sin(θ), f(...)/(...)², 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 f(cos θ). Afin d'éliminer la singularité en θ = 0 de cet intégrande, il suffit que f(x) tende vers 0 suffisamment vers quand x tend vers l'infini, au moins aussi vite que 1/x^3/2[11].

Pour un intervalle infini d'intégration, on peut utiliser le changement de variable x = L cot(θ), ce qui mène à[11]:

${\displaystyle \int _{-\infty }^{\infty }f(x)\,\mathrm {d} x=L\int _{0}^{\pi }{\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 f(L cot θ)/sin²(θ) 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ù f 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 f(x) tende vers 0 plus vite que 1/x² en ±∞ (si f décroit en 1/x², 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[11]). Cependant, si f 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 f : le problème est que, même si f(L cot θ)/sin²(θ) 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 f(x) = tanh(x³)/x³ décroit en 1/x³ mais a une discontinuité de saut dans la pente de la fonction remappée en θ = 0 et π].

Un autre changement de variable a été suggéré pour les intégrales de la forme ${\displaystyle \int _{0}^{\infty }\mathrm {e} ^{-x}g(x)\,\mathrm {d} x}$, en utilisant la transformation x = -ln[(1 + cosθ)/2] pour transformer l'intégrale vers la forme ${\displaystyle \int _{0}^{\pi }f(\cos \theta )\sin \theta \,\mathrm {d} \theta }$ avec f(u) = g(-ln[(1+u)/2])/2, à partir de là on peut réutiliser le schéma de la quadrature de Clenshaw–Curtis pour f comme vue plus haut[12]. À 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 f(−1)] sauf si g(∞) est fini.

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

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

Ces poids w_n 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 a_2k par une expression de la forme :

${\displaystyle c={\begin{pmatrix}a_{0}\\a_{2}\\a_{4}\\\vdots \\a_{N}\end{pmatrix}}=D{\begin{pmatrix}y_{0}\\y_{1}\\y_{2}\\\vdots \\y_{N/2}\end{pmatrix}}=Dy,}$

où D est la matrice de la TCD de type I en (N/2+1) points, avec :

${\displaystyle D_{k,n}={\frac {2}{N}}\cos \left({\frac {nk\pi }{N/2}}\right)\times {\begin{cases}1/2&n=0,N/2\\1&\mathrm {sinon} \end{cases}}}$

et y_n tel que :

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

Comme vu auparavant, à cause de l'aliasing, il n'y a aucun intérêt à calculer les coefficients au-delà de a_N, donc D est une matrice (N/2+1) × (N/2+1). En termes de ces coefficients c, l'intégrale est approchée par :

${\displaystyle \int _{-1}^{1}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ù c est le vecteur des coefficients a_2k et d est le vecteur des intégrales de chaque coefficient de Fourier :

${\displaystyle d={\begin{pmatrix}1\\2/(1-4)\\2/(1-16)\\\vdots \\2/(1-[N-2]^{2})\\1/(1-N^{2})\end{pmatrix}}.}$

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

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

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

${\displaystyle w=D^{T}d.\!}$

Puisque la matrice transposée D^T 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 w peuvent être pré-calculés en un temps en O(N log N) pour un N par des algorithmes de transformées de Fourier rapide.

Les poids w_n sont positifs et leur somme vaut 1[13].

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

• (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Clenshaw–Curtis quadrature » (voir la liste des auteurs).