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 x = cos θ 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 f(x) à intégrer est évaluée aux N 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 O(N log N) par des algorithmes de transformée de Fourier rapide adaptés à la TCD[1],[2].

Méthode générale modifier

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é :

 

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

 

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 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 :

 

pour k = 0,...,N et on utilise cette formule précédente pour l'intégrale en termes de ces ak. Puisque seul a2k 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 :

 

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 a2k 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 a2k 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 aN 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 modifier

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

 

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 xn = cos(nπ/N) correspondent aux extrema du polynôme TN(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 Tk(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 Ck partout si f(x) est de classe Ck 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).

Quadrature de Fejér modifier

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 ak en calculant f(cos θ) sur un ensemble de points également espacés, à mi-chemin entre les extrema :  . Ce sont les racines de TN(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 :

 

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 modifier

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(N2) 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].

Intégration avec des fonctions poids modifier

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 :

 

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

 

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 Wk 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 Wk.

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.

Intégration sur des intervalles semi-infinis ou infinis modifier

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[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   tendant vers l'infini.

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

 

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 cot2(θ/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 :

 

Le facteur multipliant sin(θ), f(...)/(...)2, 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/x3/2[11].

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

 

Dans ce cas, on utilise le fait que le terme f(L cot θ)/sin2(θ) 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/x2 en ±∞ (si f décroit en 1/x2, 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 θ)/sin2(θ) 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(x3)/x3 décroit en 1/x3 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  , en utilisant la transformation x = -ln[(1 + cosθ)/2] pour transformer l'intégrale vers la forme   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.

Pré-calcul des poids de quadrature modifier

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 wn (pour n de 0 à N/2, en supposant N pair) de sorte que

 

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

 

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

 

et yn tel que :

 

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

 

c est le vecteur des coefficients a2k et d 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 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 dT c peut être réarrangée en :

 

w est le vecteur des poids désirés wn, avec:

 

Puisque la matrice transposée DT 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 wn sont positifs et leur somme vaut 1[13].

Voir aussi modifier

Références modifier

  1. (en) W. Morven Gentleman, « Implementing Clenshaw-Curtis quadrature I: Methodology and experience », Communications of the ACM, vol. 1, no 5,‎ , p. 337-342.
  2. (en) Jörg Waldvogel, « Fast construction of the Fejér and Clenshaw-Curtis quadrature rules », BIT Numerical Mathematics, vol. 46, no 1,‎ , p. 195-202 (lire en ligne).
  3. (en) C. W. Clenshaw et A. R. Curtis, « A method for numerical integration on an automatic computer », Numerische Mathematik, vol. 2, no 197,‎ (lire en ligne).
  4. J. P. Boyd, Chebychev and Fourier Spectral Methods, Dover, New York, , 2e éd..
  5. Voir S. G. Johnson, "Notes on the convergence of trapezoidal-rule quadrature," online MIT course notes (2008).
  6. Leopold Fejér, « On the infinite sequences arising in the theories of harmonic analysis, of interpolation, and of mechanical quadratures », Bulletin of the American Mathematical Society, no 39,‎ , p. 521–534 (lire en ligne). Leopold Fejér, « Mechanische Quadraturen mit positiven Cotesschen Zahlen », Mathematische Zeitschrift, vol. 37, no 287,‎ (lire en ligne).
  7. (en) Lloyd N. Trefethen, « Is Gauss quadrature better than Clenshaw-Curtis? », SIAM Review, vol. 50, no 1,‎ , p. 67–87 (DOI 10.1137/060659831, CiteSeerx 10.1.1.157.4174)
  8. (en) Ignace Bogaert, « Iteration-Free Computation of Gauss--Legendre Quadrature Nodes and Weights », SIAM Journal on Scientific Computing, vol. 36,‎ , A1008–A1026
  9. (en) Erich Novak et Klaus Ritter, « High dimensional integration of smooth functions over cubes », Numerische Mathematik, vol. 75,‎ , p. 79–97.
  10. (en) G. A. Evans et J. R. Webster, « A comparison of some methods for the evaluation of highly oscillatory integrals », Journal of Computational and Applied Mathematics, vol. 112,‎ , p. 55-69.
  11. a b c d et e (en) John P. Boyd, « Exponentially convergent Fourier–Chebshev [sic] quadrature schemes on bounded and infinite intervals », J. Scientific Computing, vol. 2, no 2,‎ , p. 99-109.
  12. (en) Nirmal Kumar Basu et Madhav Chandra Kundu, « Some methods of numerical integration over a semi-infinite interval », Applications of Mathematics, vol. 22, no 4,‎ , p. 237-243.
  13. (en) J. P. Imhof, « On the Method for Numerical Integration of Clenshaw and Curtis », Numerische Mathematik, vol. 5,‎ , p. 138-141.