La transformée de Fourier discrète 1D

On rappelle la définition de transformée de Fourier discrète.

Définition: Soit un vecteur $\textbf{u}=(u_0,\cdots,u_{N-1})^t$ de $\mathbb C^N$, on note $\hat{\textbf{u}}$ le vecteur de $\mathbb C^N$, image de $\textbf{u}$ par la transformée de Fourier discrète défini par : $$\hat{\textbf{u}}_k = \frac{1}{N}\sum_{j=0}^{N-1}u_j e^{-2i\pi j k / N}, \quad \forall \, k\in \textbf{k}=0, 1, \ldots, N-1,$$ Cette transformation admet une transformation réciproque appelée transformée de Fourier inverse donnée par : $$ u_j = \sum_{k=0}^{N-1}\hat{\textbf{u}}_k e^{2i\pi j k / N}, \quad \forall \, j\in\textbf{j}=0, 1, \ldots, N-1.$$

Dans le but d'utiliser cette transformation pour des EDPs, Il est intéressant de voir le vecteur $\textbf{u}$ comme l'évaluation sur une grille régulière d'une fonction périodique et ainsi voir les coefficients de la transformée de Fourier comme les coefficients de la fonction décomposée dans une base de fonction trigonométrique périodique. Nous allons préciser cela.

Soit $u$ une fonction définie de $\mathbb R$ dans $\mathbb C$ supposée $L-$périodique. On note $\textbf{x}=(0, h, \ldots, (N-1)h)$ avec $h=\frac L N$. On note $\phi_k$ ($0\le k<N$) les fonctions : $$ \phi_k(x)=e^{2i\pi k\frac x L},~x\in[0,L[.$$ On reécrit la transformée de Fourier discrète pour $\textbf{u}=u(\textbf{x})$ (au sens de Hadamard) et on a $$ \textbf{u}=u(\textbf{x})=\sum_{k=0}^{N-1}\hat{\textbf{u}}_k \phi_k(\textbf{x}). $$ Si de plus $u\in Vect(\phi_k, ~0\le k<N)$ alors $u=\sum_{k=0}^{N-1}\hat{\textbf{u}}_k \phi_k$.

Le choix de la base $Vect(\phi_k, ~0\le k<N)$ est un choix pertinent de base pour une méthode de Galerkine pour des EDPs à coefficients constants et des conditions limites périodiques. Ce point sera détaillé plus tard. Il faut d'abord évaluer le coût d'un changement de base pour passer des coefficients de l'évaluation de la fonction sur la grille aux coefficients de Fourier.

Coût d'une transformée de Fourier

Une application directe des formules de transformation montre une complexité en $\mathcal O(N^2)$ puisque chacun des $N$ coefficients nécessite de l'ordre de $N$ sommes et $N$ produits. Il existe cependant un algorithme dit "FFT" pour "Fast Fourier Transform" qui repose sur un découpage récursif de la transformation en exploitant les propriétés de l'exponentielle complexe. Il s'agit de l'algorithme de Cooley-Tukey dont la construction est triviale si on se limite à $N$ une puissance de $2$.

$$\hat{\textbf{u}}_k = \frac{1}{N}\sum_{j=0}^{N-1}u_j e^{-2i\pi j k / N}=\frac{1}{N}\sum_{j=0}^{N-1}u_j W_N^{kj},$$

où $W_{N}=e^{-2i\pi / N}$. On note ici $\mathcal F_N(\textbf{u})=\hat{\textbf{u}}$ et $\textbf{u}^{odd}\in \mathbb C^{N/2}$ le vecteur des indices impairs et $\textbf{u}^{even}\in \mathbb C^{N/2}$ le vecteur des indices pairs. On supposera ici $N$ une puissance de $2$. Ainsi, comme $W_N=W_{N/2}^2$, \begin{split} \mathcal F_N(\textbf{u})_k &=\frac{1}{N}\sum_{j=0}^{N-1}u_j W_N^{kj}=\frac{1}{N}\sum_{j=0}^{N/2-1}u_{2j} W_{N}^{2kj}+\frac{1}{N}\sum_{j=0}^{N/2-1}u_{2j+1} W_N^{k(2j+1)}\\ &=\frac 1 2\frac 1 {N/2}\sum_{j=0}^{N/2-1}u^{even}_{j} W_{N/2}^{kj}+\frac 1 2\frac 1 {N/2}\sum_{j=0}^{N/2-1}u^{odd}_{j} W_{N/2}^{kj}W_N^k \end{split} Finalement, pour $k<N/2$, on a directement : $$ \mathcal F_N(\textbf{u})_k=\frac 1 2\mathcal F_{N/2}(\textbf{u}^{even})_{k}+\frac 1 2 W_N^k\mathcal F_{N/2}(\textbf{u}^{odd})_{k}. $$ Si $k\ge N/2$, on retranche $N/2$ à $k$ et on l'ajoute pour écrire la transformée discrète d'indice supérieure à $N/2$ en fonction des indices admissibles de la transformée discrète des vecteurs d'indices paires et impaires de taille $N/2$. Pour $k\ge N/2$, \begin{split} \mathcal F_N(\textbf{u})_k &=\frac 1 2\frac 1 {N/2}\sum_{j=0}^{N/2-1}u^{even}_{j} W_{N/2}^{(k-N/2)j}W_{N/2}^{jN/2}+\frac 1 2\frac 1 {N/2}\sum_{j=0}^{N/2-1}u^{odd}_{j} W_{N/2}^{(k-N/2)j}W_N^kW_{N/2}^{jN/2} \\ &=\frac 1 2\mathcal F_{N/2}(\textbf{u}^{even})_{k-N/2}+\frac 1 2W_N^k\mathcal F_{N/2}(\textbf{u}^{odd})_{k-N/2},\\ &=\frac 1 2\mathcal F_{N/2}(\textbf{u}^{even})_{k-N/2}-\frac 1 2W_N^{k-N/2}\mathcal F_{N/2}(\textbf{u}^{odd})_{k-N/2} \end{split} puisque $W_{N/2}^{jN/2}=1$ et $W_{N}^{N/2}=-1$.

Le calcul de complexité s'obtient par récurrence sur $l$ avec $N=2^l$. En notant $c_l$ la complexité en fonction de $l$, on a $c_{l+1}=2c_l+2^{l+2}$. Ainsi $c_l=(c_0+4l)2^l$, soit $c_l=(1+4\log_2(N))N=\mathcal O(N\log_2(N))$. Ce coût est drastiquement inférieur au coût d'un calcul direct en $\mathcal O(N^2)$.

Méthode de Galerkine

Pour comprendre l'intérêt de la FFT pour la résolution de certaines EDPs, il suffit de réaliser que la base $Vect(\phi_k, ~0\le k<N)$ est une base de vecteurs propres pour les opérateurs de dérivation avec condition limite périodique. En dimension supérieure, c'est encore vrai; pour d'autres conditions limites, Dirichlet et Neumann, c'est encore vrai pour les opérateurs de dérivation paire avec des bases trigonométriques modifiées construites avec une complexité encore en $\mathcal O(N\log_2(N))$.

Exemple: On considère la solution $1$ périodique de $-u''(x)=f(x)$ pour une fonction $f$ $1-$périodique donnée à moyenne nulle et par exemple dans $L^2(0,1)$. On note $P_N$ la projection sur $V_N=Vect(\phi_k, ~0\le k<N)$. On cherche alors $u_N\in V_N$ tel que $-u_N''=P_N f$. Ainsi, on cherche $N$ scalaires $u_i$, $0\le i<N$ tel que $$ u_N=\sum_{k=0}^{N-1} u_k \phi_k,~P_N f=\sum_{k=0}^{N-1} f_k \phi_k $$ et l'équation différentielle donne, comme $\phi_k''=-(2\pi k/L)^2\phi_k$ et après identification des coefficients dans la base: $$ (2\pi k/L)^2 u_k=f_k. $$ Ceci revient à résoudre un système diagonal. La clé est donc de choisir la bonne base alors même que le changement de base a un coût très réduit. On remarque que la condition nécessaire et suffisante d'existence d'une solution est $f_0=0$ ce qui est équivalent à dire que $f$ est à moyenne nulle. On comprend de suite, avec $u_0$ quelconque, que la solution est définie à une constante près.

Une autre façon d'écrire la même caractérisation est d'écrire une formulation faible de l'EDP: on cherche $u_N\in V_N$ tel que $$ a(u_N,v)=\int_0^1 u_N'(x) v'(x)\,dx=l(v)=\int_0^1 f(x)v(x)\,dx,~\forall v \in V_N. $$ De façon générale pour les méthodes de Galerkine, le système linéaire s'écrit alors: $$ \sum_{k=0}^{N-1} a(\phi_k,\phi_i)u_k=l(\phi_i),~\forall 0\le i<N. $$ La matrice de ce système linéaire (coefficient $a(\phi_k,\phi_i)$ en ligne $i$ colonne $k$) est diagonale puisque les éléments de la base $V_N$ sont orthogonaux dans $L^2(0,1)$ (et aussi $H^1(0,1)$ puisque la dérivée des fonctions de base est proportionnelle à elle-même).

Remarque : On se rappelle que les méthodes "éléments finis conformes" rentrent dans ce formalisme des méthodes de Galerkine. On voit aussi que dans ce cas, les fonctions de base sont choisies à support le plus réduit possible pour rendre la matrice la plus creuse possible. Mais il n'y a ancun espoir d'obtenir une matrice diagonale. Avec les méthodes spectrales (décomposition sur les fonctions propres de l'opérateur), les supports des fonctions de base couvrent tout le domaine, mais c'est une propriété d'orthogonalité de ces fonctions dans les Sobolev qui trivialise le système avec une matrice diagonale au lieu d'une matrice pleine. En revanche, si l'EDP n'est pas à coefficients constants, la propriété d'othogonalité est perdue et la matrice est potentiellement pleine.

Vérification des bases par transformation

Transformée rapide réelle

La FFT opère en complexe. C'est le bonne outil pour résoudre l'équation de Schrödinger sur domaine périodique. Pour une EDP 1D réelle, c'est la fonction "rfft" qu'on utilisera en python depuis la librairie scipy. La fonction "rfft" existe aussi dans la librairie numpy, elle renvoie un tableau de taille $\frac N 2+1$ de complexes associée à une base de fonction complexes. Toutes les fonctions qui suivront peuvent être chargées pour fortran ou C à l'aide du module proposé sur le site http://www.fftw.org/.

L'exemple suivant permet de comprendre la normalisation choisie des fonctions de base réelles $L-$périodique $\phi_0=\frac 1 {2N}$, $\phi_1(x)=\frac 1 {N} \cos(2\pi\frac x L)$, $\phi_2(x)=-\frac 1 {N} \sin(2\pi\frac x L)$, $\cdots$, $\phi_{N-2}(x)=-\frac 1 {N} \sin(2\pi\frac{N-2}2\frac x L)$, $\phi_{N-1}(x)=\frac 1 {2N} \cos(2\pi\frac N 2\frac x L)$, associée à la fonction "rfft" de scipy.

Transformée rapide en cosinus

Pour choisir une base de fonction réelle satisfaisant une dérivée nulle aux extrémités, une base en cosinus est adaptée. On remarque que la grille de pas $dx$ est placée au milieu des segments $N$ segments de longueur $dx$ qui couvrent l'intervalle $[0,L]$. Ansi on évalue la fonction en $\frac {dx} 2$ puis par pas de $dx$ avant d'en calculer la transformée de Fourier "dct". La base construite est $\phi_0=\frac 1 {2N}$, $\phi_1(x)=\frac 1 {N} \cos(\pi\frac x L)$, $\phi_2(x)=\frac 1 {N} \cos(\pi 2\frac x L)$, $\cdots$, $\phi_{N-2}(x)=\frac 1 {N} \cos(\pi(N-2)\frac x L)$, $\phi_{N-1}(x)=\frac 1 {N} \cos(\pi(N-1)\frac x L)$.

Contrairement à la fonction "irfft" qui est bien la transformée inverse de la fonction "rfft", la fonction "idct" doit être renormalisée par $\frac 1 {2N}$.

Transformée rapide en sinus

Pour choisir une base de fonction réelle satisfaisant une condition de Dirichlet homogène aux extrémités, une base en sinus est adaptée. On remarque que la grille de pas $dx$ est placée au milieu des segments $N$ segments de longueur $dx$ qui couvrent l'intervalle $[0,L]$. Ansi on évalue la fonction en $\frac {dx} 2$ puis par pas de $dx$ avant d'en calculer la transformée de Fourier "dst". La base construite est $\phi_0=\frac 1 {N} \sin(\pi\frac x L)$, $\phi_1(x)=\frac 1 {N} \sin(\pi 2\frac x L)$, $\phi_2(x)=\frac 1 {N} \sin(\pi 3\frac x L)$, $\cdots$, $\phi_{N-2}(x)=\frac 1 {N} \sin(\pi(N-1)\frac x L)$, $\phi_{N-1}(x)=\frac 1 {2N} \sin(\pi N\frac x L)$.

Résolution du problème de Neumann homogène

L'exemple qui suit compare la solution analytique à la solution obtenue par méthode spectrale en exploitant les bases orthogonales identifiées. On rappelle que la fonction "idct" n'est pas exactement la transformée de Fourier inverse. Elle a besoin d'être normalisée avec le facteur $\frac 1{2N}$.

Calcul de l'erreur en fonction de $N$

On propose ici de calculer l'erreur évaluée dans la cellule précédente pour différentes valeurs de $N$ ainsi que la vérification de la complexité du calcul en fonction de $N$.

Le temps de simulation se comporte grossièrement comme une constante fois $N$ pour les valeurs de $N$ testées. Pour les petites valeurs de $N$ un surcôut est probablement dû à des temps systèmes. De ce fait, le comportement en $\mathcal O(N\log(N))$ plutôt que $\mathcal O(N)$ n'est pas observable, sauf à exclure les 8 premières données.

Problème de Dirichlet 1D

On propose ici une vérification du problème de Dirichlet 1D.

Extension en dimension 2

La transformée de Fourier discrète se définit en toute dimension. On présentera les calculs similaires au 1D en dimension 2 en étendant le calcul discret aux matrices de taille $N_0,N_1$ avec $\omega_m=j_m k_m /N_m$ ($m=0,1$): \begin{split}\hat{u}_{k_0,k_1} &= \frac{1}{N_0}\sum_{j_0 \in \textbf{j}_0}\Big( e^{-2\pi i \omega_0} \frac{1}{N_1} \sum_{j_1\in \textbf{j}_1} \Big( e^{-2\pi i \omega_1} u_{j_0,j_1}\Big) \Big), \quad \forall \, (k_0, k_1) \in \textbf{I}_0 \times \textbf{I}_1, \\ u_{j_0, j_1} &= \sum_{k_1\in \textbf{k}_1} \Big( e^{2\pi i \omega_1} \sum_{k_0\in\textbf{k}_0} \Big( e^{2\pi i \omega_0} \hat{u}_{k_0, k_1} \Big) \Big), \quad \forall \, (j_0, j_1) \in \textbf{I}_0 \times \textbf{I}_1,\end{split} avec, comme un 1D, $\textbf{I}_m=\{0, 1, \ldots, N_m-1\}$.

On va comparer le calcul d'un problème de Neumann homogène calculé par transformée de Fourier rapide à l'aide de la généralisation de la fonction "dct" à "dctn" quelle que soit le nombre d'indices du tableau en argument (ici 2). On verra que le temps de calcul est drastiquement réduit par cette méthode comparée à un système creux résolu pa factorisation $LU$ ou par méthode itérative. Le ratio de temps est de plusieurs centaine pour $10^6$ inconnues. Seul les méthodes de résolution multigrilles peuvent tenter de rivaliser.

Résolution d'un problème de Poisson 3D

On observe ici l'extension 3D du problème de Neumann homogène sur un cube pour s'apercevoir que ce problème 3D est accessible au calcul séquentiel sur grille fine (par exemple $10^8$ inconnues). Avec quelques millions d'inconnues, la résolution est largement sous la seconde.

Equation de la chaleur 1D

On propose ici de résoudre l'équation de la chaleur par transformée de Fourier discrète.

Remarque: On a observé la superposition de solutions au cours du temps. On remarque que tant qu'il n'y a pas de second membre, le calcul de la solution peut se faire sans itération en temps et le résulat au temps final sera le même sans aucune erreur lié au temps.

Résolution de Ginzburg-Landau 1D

La résolution se fait par splitting de Strang en résolvant la partie EDO non linéaire dans l'espace des phases avec un schéma inconditionnellement A-stable et la partie équation de la chaleur linéaire en Fourier.

Résolution de Ginzburg-Landau 2D

La résolution se fait par splitting de Strang en résolvant la partie EDO non linéaire dans l'espace des phases avec un schéma inconditionnellement A-stable et la partie équation de la chaleur linéaire en Fourier. Les conditons limites sont ici de type Neumann et non pas des conditions limites périodiques.