Intégrales de cosinus de polynômes

Dans ce billet, on cherche à calculer des intégrales de la forme

\begin{equation*} \int_a^b Q(t) \cos(\alpha t^2 + \beta t + \gamma) {\rm d} t \end{equation*}

ou, de façon similaire,

\begin{equation*} \int_a^b Q(t) \sin(\alpha t^2 + \beta t + \gamma) {\rm d} t \end{equation*}

pour un polynôme \(Q\) quelconque. On résoudra dans un premier temps le problème de façon analytique avant de donner une implémentation C++ du calcul de la solution.

Réduction du problème

On donne ici la résolution dans le cas du cosinus, le cas du sinus procédant de même mutatis mutandis. Commençons par nous ramener à un problème plus simple. Par linéarité de l'intégrale, on peut se ramener au cas où \(Q\) est un monôme et \(a = 0\). De plus, en utilisant la parité du cosinus (respectivement l'imparité du sinus), on peut supposer \(\alpha > 0\).

Si \(Q\) est de degré \(n + 1\) avec \(n \geq 1\), une intégration par parties nous ramène au calcul d'une constante et du cas où le monôme est de degré \(n - 1\). Si \(Q\) est de degré 1, on écrit l'intégrale comme partie réelle de son pendant complexe :

\begin{equation*} \Re\left(\int_0^b t e^{it}{\rm d}t\right) \end{equation*}

Un calcul de la partie complexe permet de réécrire l'expression ainsi :

\begin{equation*} \frac{1}{2 \alpha} \left[\sin(\alpha b^2 + \beta b + \gamma) - \sin \gamma - \beta \int_a^b \sin(\alpha t^2 + \beta t + \gamma) {\rm d}t\right] \end{equation*}

On s'est ainsi ramené au cas où Q est la multiplication par une constante, ce qui résout le calcul du polynôme « externe » de l'expression. Reste celui du polynôme interne, c'est-à-dire à l'intérieur des cosinus/sinus de l'intégrande, qui est de degré 2. Éliminons son monôme de degré 1 par un changement de variable :

\begin{equation*} \int_0^b \cos(\alpha t^2 + \beta t + \gamma) {\rm d}t = \frac{1}{\sqrt{\alpha}} \int_{\frac{\beta}{2\sqrt{\alpha}}}^{\sqrt{\alpha} b + \frac{\beta}{2\sqrt{\alpha}}} \cos\left[u^2 + \left(\gamma - \frac{\beta^2}{4 \alpha}\right)\right] {\rm d}u \end{equation*}

En appliquant les formules trigonométriques des cosinus et sinus de sommes (\(\cos(a + b) = \cos(a) \cos(b) - \sin(a) \sin(b)\), etc.), on élimine alors le facteur constant (\(\gamma - \frac14 \beta / \alpha\)), modulo quelques multiplications par des constantes. Il ne nous reste plus alors qu'à intégrer un \(\cos(u^2)\).

Résolution du problème simplifié

Considérons le développement en série entière de l'intégrande :

\begin{equation*} \int_0^{b} \cos(u^2) {\rm d} u = \int_0^{b} \sum_{n=0}^{+\infty} (-1)^n \frac{u^{4n}}{(2n)!} {\rm d}u \end{equation*}

Par un critère de convergence normale, on peut intervertir somme et intégrale dans cette expression, ce qui nous donne :

\begin{equation*} \int_0^{b} \cos(u^2) {\rm d} u = \sum_{n=0}^{+\infty} \frac{(-1)^n}{(2n)!} \frac{b^{4n+1}}{4n+1} \end{equation*}

Le terme général de cette série est facile à calculer et son reste facile à estimer la série est alternée. Pour peu que \(b\) ne soit pas « trop grand », on peut approximer l'intégrale en calculant la série jusqu'à un certain terme \(N\), le reste d'ordre \(N\) nous donnant une indication sur l'erreur commise.

Toutefois, cette méthode ne fonctionne pas quand \(b\) est trop grand : les termes de la série ont alors des ordres de grandeur très différents, et les termes d'ordre élevé occultent ceux d'ordre inférieur lorsqu'on calcule la somme avec la représentation habituelle des réels sur un nombre fini de bits (par exemple un long double en C). Une solution (non-triviale) à ce problème consiste à regrouper les termes de même ordre de grandeur : la série étant alternée, la somme des termes d'ordre k est au-moins d'ordre k-1. Nous ne rentrerons pas dans les détails d'implémentation de cette solution ici.

Implémentation

Voici une implémentation itérative en C++ de la méthode que nous venons de présenter. Commençons par le problème réduit :

long double integrateCosX2(long double b) {
    const long double precision = 1e-20;
    const long double stopThres = precision / b;

    if (b > 1)
       warn("Precision issue due to high bound.");

    long double curFact = 1., curTerm = 1., sum = curTerm;
    for (int i = 1; curTerm < 0 || curTerm > stopThres; i++) {
       curFact *= (-1.) * (b*b*b*b) / (2.*i) / (2.*i - 1.);
       curTerm = curFact / (4.*i + 1.);
       sum += curTerm;
    }

    return b * sum;
}

Real integrateCosX2(Real from, Real to) {
    return integrateCosX2(to) - integrateCosX2(from);
}

La fonction integrateSinX2(long double) s'implémente de façon similaire. Maintenant, en supposant qu'on a défini par ailleurs un type polynôme Poly où l'opérateur [] permet la lecture des coefficients, on peut coder la fonction de calcul générale ainsi :

long double integrateTrigPoly(TrigFun trig, Poly P, long double from, long double to) {
    if (P[2] < 0)
        switch (trig) {
           case COS:
              return integrateTrigPoly(trig, -P, from, to);
           case SIN:
           default:
              return -integrateTrigPoly(trig, -P, from, to);
        }

    long double I1, I2, dcSqrt = sqrtl(P[2]);
    long double newFrom = dcSqrt * from + P[1] / 2. / dcSqrt;
    long double newTo = dcSqrt * to + P[1] / 2. / dcSqrt;
    long double delta = P[0] - P[1]*P[1] / 4. / dcSqrt;

    switch (trig) {
       case COS:
          I1 = cosl(delta) * integrateCosX2(newFrom, newTo);
          I2 = sinl(delta) * integrateSinX2(newFrom, newTo);
          return I1 - I2;
       case SIN:
       default:
          I1 = cosl(delta) * integrateSinX2(newFrom, newTo);
          I2 = sinl(delta) * integrateCosX2(newFrom, newTo);
          return I1 + I2;
    }
}

Bilan de ce billet : par des simplifications analytiques et du calcul de séries entières, nous sommes parvenus à calculer en peu de code des intégrales de cosinus et sinus de polynômes. Toutefois, notre algorithme devient instable à mesure que le domaine d'intégration s'agrandit (i.e., pour \(b\) croissant) et nous n'avons pas résolu cette instabilité de manière satisfaisante.

Edit: page mise-à-jour le 18 juin 2014.

Edit: page mise-à-jour le 29 janvier 2016. Les formules sont maintenant générées par MathJax.

Pages of this website are under the CC-BY 4.0 license.