Calcul pratique de l'approximation pseudo-spectrale des dérivées

On calcule les $\frac{d^n}{d x^n} \widetilde{P_N} u_j$ en les considérant comme le résultat de produits linéaires avec les valeurs de $u_k$. Soit :

$$U = U^{(0)} = \left[ \begin{array}{c} u_0\\ \vdots\\ u_{N-... ...{c} u_0^{(r)}\\ \vdots\\ u_{N-1}^{(r)}\\ \end{array}\right]$$

où $u_{k}^{(r)}$ est l'approximation pseudo-spectrale de $\frac{d^r}{dx^r} u_k$. On souhaite écrire l'opérateur différentiel $D$ tel que $U^(r) = D U^{(r-1)} = \cdots = D^{(r)} U$. L'opérateur $D$ est nommé matrice de différentiation pseudo-spectrale. On le construit en trois étapes :

1. A partir des $u_i$ $(i=0,\cdots,N-1)$, on calcule les $\widetilde{a_k}$ $(k=-N/2+1,\cdots,N/2)$ en utilisant la transformée de Fourier. En prenant soin de convertir les indexes, cette opération est strictement équivalente à une multiplication matricielle $A=F_1U$. Le nombre d'opérations est alors en $\bigcirc(N^2)$.
2. Calcul des termes spectraux des dérivées $1^{eres}$ en multipliant chaque $a_k$ par $2 \pi i k / l$. Cette opération est équivalente à un produit scalaire $A' = V A$ avec $V=[2 \pi i (-N/2+1) / l, \cdots,2 \pi i (N/2) / l]$. Le nombre d'opérations est ici en $\bigcirc(N)$.
3. A partir du résultat obtenu à l'étape précédente, on calcule les $\frac{d}{d x} \widetilde{P_N} u_i$ en utilisant une transformée de Fourier inverse. Comme à la première étape, cette opération est strictement équivalente à une multiplication matricielle par une matrice $\frac{d}{d x} \widetilde{P_N} U = F_2 A'$. Le nombre d'opérations est encore en $\bigcirc(N^2)$.

L'opérateur $D$ vaut alors :

$$D= F_2 W F_1$$

avec $W$ la matrice dont la diagonale est égale à $V$ et les autres termes sont nuls.

L'opération est au final de complexité algorithmique $N^2$ ce qui est rédhibitoire pour de grandes valeurs de $N$. Comme cela a été évoqué dans la section précédente, en pratique, on réduit considérablement la complexité algorithmique de l'opération en utilisant la Transformée de Fourier Rapide (en Anglais Fast Fourier Transform ou FFT). Pour cela, $N$ doit être une puissance de $2$. Une transformée de Fourier passe d'une complexité $\bigcirc(N)$ à une complexitée $\bigcirc(log_2N)$. La procédure de calcul de $\frac{d^r}{dx^r} u_k$ est alors la suivante :

1. A partir des $u_i$ $(i=0,\cdots,N-1)$, on calcule les $\widetilde{a_k}$ $(k=-N/2+1,\cdots,N/2)$ en utilisant la FFT. Le nombre d'opérations est maintenant en $\bigcirc(N log_2 N)$ avec $N$ pour les $N$ $u_i$ et $\log_2 N$ pour la FFT.
2. Calcul des termes spectraux des dérivées $r^{\textrm{\\lq emes}}$ en multipliant chaque $a_k$ par $(2 \pi i k / l)^r$. L'opération est toujours ici en $\bigcirc(N)$.
3. A partir du résultat obtenu à l'étape précédente, calcule des $\frac{d^r}{d x^r} \widetilde{P_N} u_j$ en utilisant une FFT inverse. Comme à la première étape, le nombre d'opérations est en $\bigcirc(N log_2 N)$.

Au final, la complexité algorithmique est maintenant en $\bigcirc(N log_2 N)$ ce qui est beaucoup plus raisonnable. On ne peut plus par contre considérer cette opération comme un produit matriciel mais plutôt comme une boîte noire, ce qui n'est pas bien génant tant que les résultats sont les mêmes.