III.D.2.2 Programmation de I₂.

La programmation de I₂ est beaucoup plus compliquée que celle de I₁, beaucoup de situations possibles sur les paramètres de calcul interviennent. Nous allons présenter au lecteur plusieurs algorithmes possibles de calcul et pour l'un deux un résumé de la marche suivie pour minimiser l'erreur.

Rappelons que I₂ est donné par

$$I_2(\alpha,\beta) = 2S e^{i{{\mathbf k}}\vec{G}} e^{ -\frac{... ...}{\displaystyle \beta }}{\displaystyle 2i( \alpha-\beta )}) \$$ (III.D.10)

Le domaine de définition de $I_2(\alpha,\beta)$ est ${\mathbb{C}}^2$ privé des droites $\alpha=0$, $\beta=0$ et $\alpha=\beta$. On prolonge cette fonction par continuité sur ${\mathbb{C}}^2$: la fonction prolongée est analytique sur ${\mathbb{C}}^2$. Comme nous l'avons expliqué pour le calcul de I₁ (III.D.9), il n'est pas possible de calculer I₂ directement sur le plan complexe puisque la division par zéro est interdite sur un calculateur. Nous allons suivre la même démarche que pour le calcul de I₁. On remarque que si l'on pose

$$f(\alpha) = e^{ i \alpha } \frac{\displaystyle \sin{ \alpha } }{\displaystyle \alpha }$$

on calcule I₂ par

$$I_2(\alpha,\beta) = 2S e^{i{{\mathbf k}}\vec{G}} e^{ -\frac{... ...le f(\alpha)-f(\beta) }{\displaystyle 2i( \alpha-\beta )}) \ .$$

Nous savons maintenant calculer la fonction $f(\alpha)$ sur ${\mathbb{C}}$ puisque c'est le même problème que celui du calcul de $I_1(\alpha)$. Ceci implique que nous sommes capables de calculer $I_2(\alpha,\beta)$ sur les axes $\alpha=0$ et $\beta=0$. Le problème nouveau de notre étude est donc de calculer la fonction I₂ dans un voisinage de la droite $\alpha=\beta$ à la meilleure précision possible et avec le moins d'opérations possible. Le critère supplémentaire à considérer est donc $\gamma=\alpha-\beta$ proche de zéro ou pas, et dans quel sens. En effet, on a,

$$I_2 = \tilde{I}_2(\alpha,\gamma) = 2S e^{i{{\mathbf k}}\vec{... ...f(\alpha)-f(\alpha -\gamma) }{\displaystyle 2i( \gamma )}) \ .$$

Lorsque $\gamma \to 0$, la fonction $\tilde{I}_2(\alpha,\gamma)$ tend vers

$$-iS e^{i{{\mathbf k}}\vec{G}} e^{ -\frac{4}{3}i \alpha} f'(\alpha)$$

au premier ordre, avec

$$f'(\alpha) = ie^{ i \alpha } \frac{\displaystyle \sin{ \alph... ...lpha } - \sin{ \alpha } }{\displaystyle \alpha^2 } \right) \ .$$

La fonction dérivée $f'(\alpha)$ admet une limite en zéro. Nous pouvons donc proposer l'algorithme suivant.

• Si $\vert\gamma\vert \le \epsilon_2$, un réel positif à définir que l'on suppose supérieur à $\epsilon_1$, alors
  • si $\vert\alpha\vert \le \epsilon_3$, un autre réel positif à définir, alors on sait que $f'(\alpha)=i$ et l'on a directement $I_2=S e^{i{{\mathbf k}}\vec{G}}$,
  • sinon, on calcule
    
    $$I_2 = -iS e^{i{{\mathbf k}}\vec{G}} e^{ -\frac{4}{3}i \alpha... ...lpha } - \sin{ \alpha } }{\displaystyle \alpha^2 } \right) \ .$$
    
    

• Si $\vert\gamma\vert > \epsilon_2$, alors
  • si $\vert\alpha\vert \le \epsilon_1$, (d'où $\beta>\epsilon_1$) alors
    
    $$I_2 = 2S e^{i{{\mathbf k}}\vec{G}} e^{ -\frac{2}{3}i \beta }... ...splaystyle 1-f(\beta) }{\displaystyle 2i( \alpha-\beta )}) \ ,$$
    
    

  • de même si l'on inverse les rôles de $\alpha$ et $\beta$,
  • sinon, on calcule I₂ par la formule (III.D.10).

Nous laissons au lecteur le soin de déterminer $\epsilon_2$ et $\epsilon_3$ de façon à minimiser l'erreur globale du calcul. Pour cela, il faudra effectuer des estimations pointues des restes des fonctions approchant I₂.

Nous avons préféré utiliser un autre algorithme qui exploite les propriétés particulières de la fonction $I_2(\alpha,\beta)$ issues de considérations géométriques. En effet, rappelons que l'on a défini $\alpha = {{\mathbf k}}\, \frac{\displaystyle (\vec{x}_3-\vec{x}_2)}{\displaystyle 2}$, $\beta = {{\mathbf k}}\, \frac{\displaystyle (\vec{x}_1-\vec{x}_2)}{\displaystyle 2}$ et $\gamma = {{\mathbf k}}\, \frac{\displaystyle (\vec{x}_3-\vec{x}_1)}{\displaystyle 2}$. Si on remplace $\alpha$ par $-\beta$ et $\beta$ par $\gamma$, on sait que

$$I_2 = I_2(\beta,-\gamma) = 2S e^{i{{\mathbf k}}\vec{G}} e^{ ... ... }{\displaystyle \beta }}{\displaystyle 2i(\beta+\gamma)}) \ .$$

La fonction $I_2(\beta,-\gamma)$ est une alternative au calcul de I₂ donné par la formule (III.D.10).

Nous pouvons donc développer l'algorithme plus simple suivant.

• Si $\vert\gamma\vert \le \epsilon_2$, un réel positif à définir que l'on suppose supérieur à $\epsilon_1$, alors nous avons le choix suivant.
  • Si $\vert\alpha\vert \geq \epsilon_2$, alors,
    • si $\vert\beta\vert \geq \epsilon_1$, alors,
      
      $$I_2 = I_2(\beta,-\gamma) = 2S e^{i{{\mathbf k}}\vec{G}} e^{ ... ...\beta } }{\displaystyle \beta }}{\displaystyle 2i \alpha}) \ ,$$
      
      

    • sinon,
      
      $$I_2 = 2S e^{i{{\mathbf k}}\vec{G}} e^{ -\frac{2}{3}i \left(\... ...a } }{\displaystyle \gamma } -1}{\displaystyle 2i \alpha}) \ .$$
      
      

  • Si $\vert\alpha\vert < \epsilon_2$, alors on effectue un développement de Taylor par rapport à $\alpha$ et $\beta$, et l'on peut calculer I₂ par
    
    $$I_2 = 2S e^{i{{\mathbf k}}\vec{G}} e^{ -\frac{2}{3}i \left(\gamma \right)} T_n(\alpha,\beta)$$
    
    
    où $T_n(\alpha,\beta)$ est le polynôme approchant le développement limité de la différence des deux sinus cardinaux de $\alpha$ et $\beta$ à l'ordre n-2 par rapport aux deux variables: n est l'ordre du développement de Taylor à effectuer sur la fonction I₂. Par exemple,
    

    $$\null\,\vcenter{\openup\jot \let\\ =\@ \ialign{\strut\hfil... ... \beta ^3) + \beta ^4 )}{\displaystyle 45} \ . \cr \crcr}}\,$$ (III.D.11)

    
    
    Notons que si l'on veut éviter le calcul et la multiplication par $e^{ -\frac{2}{3}i \left(\gamma \right)}$, on peut aussi calculer I₂ par
    
    $$I_2 = 2S e^{i{{\mathbf k}}\vec{G}} D_n(\alpha,\beta)$$
    
    
    où $D_n(\alpha,\beta)$ est le polynôme approchant le développement limité de I₂ à l'ordre n par rapport aux deux variables $\alpha$ et $\beta$. Par exemple,
    

    $$\null\,\vcenter{\openup\jot \let\\ =\@ \ialign{\strut\hfil... ...alpha^3 + 3\beta^2\alpha^2}{\displaystyle 405} \cr \crcr}}\,$$ (III.D.12)

    
    

• Si $\vert\gamma\vert > \epsilon_2$, alors,
  • si $\vert\alpha\vert \le \epsilon_1$, (d'où $\beta \ge \epsilon_2-\epsilon_1$) alors
    
    $$I_2 = 2S e^{i{{\mathbf k}}\vec{G}} e^{ -\frac{2}{3}i \beta }... ...\displaystyle \beta } }{\displaystyle 2i( \alpha-\beta )}) \ ,$$
    
    

  • de même si l'on inverse les rôles de $\alpha$ et $\beta$,
  • sinon, on calcule I₂ par la formule (III.D.10).

Il nous reste à choisir $\epsilon_2$ de façon à optimiser la précision du calcul de I₂ puis à minimiser le nombre d'opérations à effectuer. Il s'agit d'étudier les précisions des calculs de I₂ par les différentes formules présentées, puis d'obtenir les mêmes précisions sur tous les cas limites afin que la majoration globale de l'erreur sur le calcul de I₂ soit aussi faible que possible. Sous l'hypothèse $\epsilon_1 \le \epsilon_2/2$, les différentes formules de calcul utilisent les trois fonctions suivantes:

1. la fonction h(x,y)
   

   $$h(x,y) = \frac{\displaystyle e^{ i x } \frac{\displaystyle \... ...ystyle \sin{ y } }{\displaystyle y }}{\displaystyle 2i( x-y )}$$ (III.D.13)

   
   
   définie sur $\vert x-y\vert > \epsilon_2$ et $\vert x\vert >\epsilon_1$ et $\vert y\vert >\epsilon_1$,
2. la fonction g(x,y)
   

   $$g(x,y) = \frac{\displaystyle e^{ i x } \frac{\displaystyle \sin{ x } }{\displaystyle x } -1}{\displaystyle 2i( x-y )}$$ (III.D.14)

   
   
   définie sur $\vert x-y\vert > \epsilon_2$ et $\vert y\vert \le \epsilon_1$,
3. la fonction d(x,y)
   

   $$\null\,\vcenter{\openup\jot \let\\ =\@ \ialign{\strut\hfil$... ...2 y ^2) + ( x y ^3) + y ^4 )}{\displaystyle 45} \cr \crcr}}\,$$ (III.D.15)

   
   
   définie sur $\vert x-y\vert \le \epsilon_2$ et $\vert y\vert \le \epsilon_2$.

La fonction d(x,y) permet d'approcher I₂ en effectuant une erreur de l'ordre du premier terme oublié dans la série de Taylor, terme de la forme

$$\frac{\displaystyle 2i}{\displaystyle 315} \left( x^5 +x^4y +x^3y^2 +x^2y^3 +xy^4 +y^5 \right)$$

dont le module est maximal pour $y=\epsilon_2$ et $x=2\epsilon_2$, et de valeur $\frac{2}{5} \epsilon_2^5$. Pour $\epsilon_2$ assez petit, nous ferons l'hypothèse que l'erreur est donc $\frac{2}{5} \epsilon_2^5$. Par exemple, pour $\epsilon_2=1/10$, le terme suivant est 50 fois inférieur à ce terme.

Etudions comment se comporte la fonction g(x,y) (III.D.14) autour du point limite $x=\epsilon_2$ et y=0. On a, puisque $\epsilon_2$ est petit,

$$g(\epsilon_2,0) \approx \frac{\displaystyle 1 + i \epsilon_2 - 2/3 \epsilon_2^2 -1}{\displaystyle 2i( \epsilon_2 )} \ .$$

On constate que cette fonction a bien une limite informatique pour $\epsilon_2$ petit non nul et que cette limite est 1/2. En revanche, autour du point $x=-i\epsilon_2$ et y=0, on a

$$g(-i\epsilon_2,0) \approx \frac{\displaystyle 1 + \epsilon_2 + 2/3 \epsilon_2^2 -1}{\displaystyle 2( \epsilon_2 )} \ .$$

Pour $\epsilon_2$ assez petit, l'ordinateur calcule 1-1 qui est nul, soit une erreur de l'ordre de 1/2, ce qui est inacceptable. L'erreur sur cette fonction de $\epsilon_2$ est

$$p_m \times \frac{1}{\epsilon_2} \ .$$

La même analyse peut être effectuée sur la fonction h(x,y) différence de deux fonctions évaluées à la précision machine et dont l'erreur sur la différence revient à étudier la fonction g(x,y). Nous choisissons donc $\epsilon_2$ de façon à avoir

$$\frac{2}{5} \epsilon_2^5 = p_m \times \frac{1}{\epsilon_2}$$

soit,

$$\epsilon_2 = \left( \frac{5}{2} p_m \right)^{1/6}$$

où l'on vérifie bien l'hypothèse $\epsilon_1 < \epsilon_2/2$ à condition d'avoir $p_m \le 1/10$. L'erreur relative sur le calcul de I₂ est alors majorée par

p_m^5/6

qui est proche de la précision machine.

L'erreur effectuée sur le calcul de I₂ est donc faible, largement inférieure à toute méthode de calcul par intégration numérique (sur des points de Gauss par exemple).