from IPython.display import display, Latex
from IPython.core.display import HTML
%reset -f
%matplotlib inline
%autosave 300
from matplotlib.pylab import *
#from time import clock
#importation des bibliothèques utiles
from scipy.sparse.linalg import spsolve
from scipy.sparse import spdiags
Autosaving every 300 seconds
Outils numériques pour des modèles d'EDP issus de problèmes physiques
Soit $\Omega\subset \mathbb R^n$ ouvert, soit $f:\Omega\to \mathbb R$, on peut choisir $f$ plus ou moins régulière, $f\in \mathcal C^\infty(\Omega)$, ou $f\notin \mathcal C^0(\Omega)$, $f\in L^1(\Omega)$...
L'appartenance à l'espace découle de la norme finie définie pour l'espace. Les espaces de fonctions sont en général de dimension infinie (sauf en numérique!), les normes ne sont pas nécessairement équivalentes ($N_2(u)\le c N_1(u)\le c'N_2(u)$). Exemples:
Remarque: Si on approche les fonctions de ces espaces par des fonctions affines par morceaux sur un maillage donné, l'espace d'approximation ainsi construit est de dimension finie et toutes les normes sont équivalentes. Mais attention, les constantes intervenants dans l'équivalence de norme dépendent du maillage et explosent par raffinement de maillage.
Dans la mesure du possible, on note en minuscule $f$ une fonction à valeur scalaire ($f:\Omega\to \mathbb R$) et $F$ ou $\bf f$ ou $\vec f$ une fonction à valeur vectorielle ($ F:\Omega\to \mathbb R^m$).
On note indifféremment $\partial_i f$, $\partial_{x_i} f$, $\frac {\partial}{\partial x_i} f$ ou encore $f_{x_i}$ la dérivée dans la direction $x_i$: $$ \partial_{x_i} f(x)=\lim_{h\to 0}\frac {f(x_1,\cdots,x_i+h,x_{i+1},\cdots)-f(x)}{h}. $$
La formule de Green s'écrit (pour des fonctions suffisament régulières pour donner un sens aux intégrales): $$ \int_{\Omega} \partial_{x_i} f(x) g(x) \,dx=-\int_{\Omega} \partial_{x_i} g(x) f(x) \,dx+\int_{\partial \Omega} f(x) g(x) {\bf n}_i \,d\gamma, $$ où $\partial \Omega$ désigne le bord du domaine $\Omega$, ${\bf n}_i$ est la $i$-ème composante de la normale unitaire sortante au domaine et $d\gamma$ désigne la mésure du bord $\partial \Omega$.
En particulier, en dimension $n=1$, on retrouve la formule de l'intégration par partie puisque ${\bf n}_1=+1$ sur le bord droit de $\Omega =]a,b[$, ${\bf n}_1=-1$ sur le bord gauche de $\Omega =]a,b[$ et la mesure du bord est l'évaluation aux points du bord: $$ \int_{\partial \Omega} f(x) g(x) {\bf n}_i \,d\gamma=f(b)g(b)-f(a)g(a). $$
\partial_3 F_1-\partial_1 F_3\\
\partial_1 F_2-\partial_2 F_1
\end {matrix}\right ]$On peut étendre les opérateurs opérant sur des scalaires à des opérateurs opérant sur des vecteurs en les définissant composante par composante. De même, les opérateurs opérant sur des vecteurs à des opérateurs opérant sur des tenseurs d'ordre 2 en les définissant composante par composante.
L'opérateur Laplacien opérant sur $\Omega\subset \mathbb R^3$ à valeur $\mathbb R^3$ peut se décomposer ainsi: $$ \Delta F =-\nabla\times \nabla\times F+\nabla \nabla \cdot F $$ </font>
#EXEMPLE D'ECOULEMENT POTENTIEL
def u1(x,y):
#phi(x,y)=x^2sin(y)
#partial_x phi=2xsin(y)
return 2*x*sin(y)
def u2(x,y):
#phi(x,y)=x^2sin(y)
#partial_y phi=x^2cos(y)
return x**2*cos(y)
X,Y=meshgrid(linspace(0.25,1.5,200),linspace(0,4,200))
U1=u1(X,Y)
U2=u2(X,Y)
figure(figsize=(10,10))
streamplot(X, Y, U1, U2,color=sqrt(U2**2+U1**2), linewidth=1, cmap='autumn');
Il faut se méfier de l'extension des résultats de fonction d'une seule variable. Si $f'=0$ sur un intervalle alors $f$ est constant sur cet intervalle. En dimension supérieur, il faut être vigilent quant à l'opérateur de dérivation considéré.
Propriété: Pour tout $f:\Omega \to \mathbb R$ régulière, alors $\nabla \times \nabla f=0$
Cette propriété est immédiate par définition de l'opérateur rotationnel. La réciproque est fort utile:
Propriété: Si $F:\Omega \to \mathbb R^3$ avec $\Omega$ simplement connexe et si $\nabla \times F=0$ sur $\Omega$ alors il existe une fonction scalaire $f:\Omega \to \mathbb R$ tel que $$ F=\nabla f. $$ Autrement dit, $F$ dérive d'un potentiel.
Remarque: Une fonction vectorielle prise au hasard n'a aucune raison de dériver d'un potentiel. C'est seulement le cas lorsque son rotationel est nul.
Propriété: Si $F:\Omega \to \mathbb R^3$ alors $\nabla \cdot \nabla \times F=0$
Cette propriété est immédiate par définition des opérateurs divergence et rotationnel. La réciproque est encore fort utile:
Propriété: Si $F:\Omega \to \mathbb R^3$ avec $\Omega$ simplement connexe et si $\nabla \cdot F=0$ sur $\Omega$ alors il existe une fonction vectorielle $G:\Omega \to \mathbb R^3$ tel que $$ F=\nabla \times G. $$
Remarque: $G$ n'est fixé qu'à un gradient près puisque la divergence d'un rotationel est nulle.
Remarque: Lorsqu'un gradient est nul sur une composante connexe, alors la fonction est constante sur cette composante. Il suffit de raisonner composante par composante.
Propriété: Il découle des propriétés précédentes qu'un champ de vitesse ${\bf u}:\Omega \to \mathbb R^3$ de $(L^{2}(\Omega))^3$ se décompose en un champ incompressible (donc un pur rotationnel) et un gradient: $$ {\bf u}=\nabla \times G +\nabla \varphi={\bf u}_{\bot}+\nabla \varphi, $$ avec $\varphi \in \dot H^1(\Omega)$. Lorsque le domaine est infini et la décroissance de ${\bf u}$ à l'infini suffisament rapide, on peut choisir une décomposition unique et orthogonale dans $(L^{2}(\Omega))^3$. Cette décomposition est très utile, elle se nomme décomposition (ou projection) de Hodge ou de Helmotz en électromagnétisme et décomposition de Leray en mécanique des fluides. On appelle ${\bf u}_{\bot}$ le projeté de ${\bf u}$ sur les champ à divergence nulle. Le calcul de cette projection est largement utilisé dans les algorithmes de résolution de Navier-Stokes en écoulement incompressible. On détaillera ce point qui repose sur le calcul de $\varphi$ comme solution d'un problème de Poisson obtenu en prenant la divergence de l'expression ci-dessus et en y ajoutant des conditions limites bien choisies.
On se sert de la seule formule de Green donnée précédemment pour les généraliser aux opérateurs divergence, gradient, rotationnel, laplacien. $$ \int_{\Omega}\nabla f\cdot G\,dx=-\int_{\Omega}f\nabla \cdot G\,dx +\int_{\partial \Omega}f G\cdot {\bf n}\,d\gamma. $$ En effet, $$ \int_{\Omega}\partial_i f G_i\,dx=-\int_{\Omega}f \partial_i G_i\,dx +\int_{\partial \Omega}f G_i {\bf n}_i\,d\gamma $$ et par sommation sur $i$, on obtient le résultat. En écrivant le laplacien comme la divergence du gradient, on a également, $$ \int_{\Omega}\Delta f \, g\,dx=\int_{\Omega}\nabla\cdot \nabla f \, g\,dx=-\int_{\Omega}\nabla f\cdot\nabla g\,dx +\int_{\partial \Omega}g \nabla f\cdot {\bf n}\,d\gamma. $$ On a encore $$ \int_{\Omega}\nabla \times F \cdot G\,dx=\int_{\Omega}\nabla \times G \cdot F\,dx +\int_{\partial \Omega}{\bf n} \times F \cdot {G}\,d\gamma. $$
Le but est de revisiter les problèmes de type diffusion comme l'opérateur laplacien sur domaine borné. On se limitera ici aux opérateurs d'ordre 2 éventuellement à coefficients non constants et pour différentes conditions limites. On rangera dans cette famile d'opérateur tous ceux qui admettent une formulation faible que satisfait le théorème de Lax-Milgram.
Théorème Lax-Milgram : On se donne
Trouver $u\in V$ tel que pour tout $v\in V$, $ a(u,v)=l(v)$.
Sous les hypothèses précédetentes, il existe une unique solution à la formuation faible.
La capacité à écrire une formulation faible sans se tromper sur le choix de l'espace des fonctions tests est supposée acquise.
Quelques clés pour se faire:
Pour s'entraîner voici quelques exercices à savoir traiter:
Trouver $u\in V$ tel que pour tout $v\in V$ $$ \int_{\Omega} \nabla u \cdot \nabla v \,dx= \int_{\Omega} f v\,dx. $$
Trouver $u\in V$ tel que pour tout $v\in V$ $$ \int_{\Omega} \nabla u \cdot \nabla v \,dx= \int_{\Omega} f v\,dx. $$
Trouver $u\in V$ tel que pour tout $v\in V$ $$ \int_{\Omega} \nabla u \cdot \nabla v \,dx+ \int_{\partial \Omega} a u v \,d\gamma= \int_{\Omega} f v\,dx+ \int_{\partial \Omega} g v \,d\gamma. $$
Quelques clés techniques pour établir les hypothèses du théorèmes de Lax-Milgram
Ce problème peut être issu d'un problème de thermique en régime stationnaire avec des parois parfaitement isolante (cas Neumann homogène). Intuitivement, on comprend de suite qu'une solution stationnaire n'est possible, alors qu'il n'y a aucun échange de chaleur avec l'extérieur, que si le terme source est globalement neutre en terme d'apport de chaleur (terme source à moyenne nulle). On verra à quel moment cette hypothèse est mathématiquement nécessaire.
Une autre motivation de ce problème est l'application en mécanique des fluides ou électromagnétisme lorsqu'il sera nécessaire de projeter un champ de vitesse sur les champ à divergence nulle.
Soit $\Omega\subset \mathbb R^d$ ($d=2$ ou $3$ en pratique). On suppose donné un champ de vitesse $\bf u\in (L^2(\Omega))^d$ . La décomposition de Hodge ou Leray stipule qu'il existe un unique $\varphi\in (\dot H^1(\Omega))^d$ tel que $$ {\bf u}={\bf u}_\bot +\nabla \varphi, \text{ dans } \Omega, $$ avec $\nabla \cdot {\bf u}_\bot=0$ dans $\Omega$ avec ${\bf u}_\bot$ le projeté orthogonal de ${\bf u}$ pour la norme de $L^2(\Omega)$ sur les champs à divergence nulle.
Ce projeté est précisément défini par : $$ \|{\bf u}-{\bf u}_\bot\|_{L^2}=\min_{v\in H}\|{\bf u}-{\bf v}\| $$ où $H=\{{\bf v} \in (L^2(\Omega))^d, \text{ tel que } \nabla \cdot {\bf v} =0 \}$.
Si on traite un écoulement ${\bf u}$ où $\partial \Omega$ sont des parois imperméables alors ${\bf u}\cdot {\bf n}=0$ et ainsi ${\bf u}$ est à divergence moyenne nulle. En faisant l'hypothèse que ${\bf u}_\bot\cdot {\bf n}=0$, on peut définir un projeté (quitte à perdre l'orthogonalité) aisément puisqu'alors sur le bord: $$ \nabla \varphi\cdot {\bf n}={\bf u}\cdot {\bf n}-{\bf u}_\bot\cdot {\bf n}=0 $$ si la décomposition est vrai au sens de la trace. En prenant la divergence de la décomposition et en ajoutant la condition limite obtenue, on a le problème de Poisson : $$ \begin {cases} &\Delta \varphi=\nabla \cdot {\bf u} \text{ sur } \Omega,\\ &\partial_n \varphi =0\text{ sur } \partial \Omega. \end {cases} $$ On remarque que la condition de compatibilité est satisfaite : $$ \int_{\Omega}\nabla \cdot {\bf u}\, dx=0=\int_{\partial\Omega}\partial_n \varphi\, d\gamma $$ Cette condition est nécessaire puisqu'en intégrant l'équation : $$ \int_{\Omega}\Delta{ \varphi}\, dx=\int_{\Omega}\nabla \cdot {\bf u}\, dx=0 $$ et que $$ \int_{\Omega}\Delta{ \varphi}\, dx= \int_{\partial\Omega}\partial_n \varphi\, d\gamma. $$
Cherchons la formulation faible du problème, on multiplie par $\psi$ et on intègre : $$ \int_{\Omega} \nabla \varphi \cdot \nabla \psi \,dx+ \int_{\partial \Omega} 0 \,d\gamma= \int_{\Omega} f \psi\,dx, $$ avec $f=-\nabla \cdot {\bf u}$ à moyenne nulle.
Le choix des fonctions tests pour assurer l'intégrabilité du premier terme est $\psi \in H^1(\Omega)$. Mais on voit que cet espace est trop large pour assurer la coercivité (ellipticité) de la forme bilinéaire (pas d'inégalité de Poincaré). On veut alors se limiter à $\dot H^1(\Omega)$ même si les fonctions régulières à support compact ne sont pas dans cet ensemble. Le théorème de Lax-Milgram va alors s'appliquer aisément, la difficulté est reporté dans l'interprétation de la solution faible pour retrouver la solution faible solution de l'EDP de départ.
Pour cela, on remarque qu'une fonction $\psi \in H^1(\Omega)$ se décompose comme $$ \psi =\bar \psi +\frac 1 {|\Omega|}{\int_{\Omega} \psi \, dx} $$ avec $\bar \psi \in \dot H^1(\Omega)$. Ainsi la solution dans $\dot H^1(\Omega)$ de, $$ \int_{\Omega} \nabla \varphi \cdot \nabla \bar\psi \,dx= \int_{\Omega} f \bar\psi\,dx, \text { pour tout }\bar\psi \in \dot H^1(\Omega), $$ est aussi solution de $$ \int_{\Omega} \nabla \varphi \cdot \nabla \psi \,dx= \int_{\Omega} f \psi\,dx, \text { pour tout }\psi \in H^1(\Omega), $$ puisque le gradient d'une constante est nul et que $f$ multiplié par une constante est à moyenne nulle. C'est ici que l'hypothèse de compatibilité intervient. Si on ne l'avait pas eu, la solution de la formulation faible résoudrait un autre problème.
On retrouve donc l'EDP $-\Delta ~\varphi=f$ au sens des distribution en restreignant encore les fonctions tests aux fonctions $\mathcal C^\infty$ à support compact puis on retrouve, $$ \int_{\Omega} -\Delta \,\varphi ~ \psi \,dx+\int_{\partial\Omega} \partial_n \,\varphi ~ \psi \,d\gamma= \int_{\Omega} f \,\psi\,dx, \text { pour tout }\psi \in H^1(\Omega). $$ Soit finalement : $$ \int_{\partial\Omega} \partial_n \,\varphi ~\psi \,d\gamma=0, \text { pour tout }\psi \in H^1(\Omega). $$ Ce qui signifie que $\partial_n \,\varphi =0$ au sens du dual de $H^{\frac 1 2}(\partial \Omega)$. On a ainsi identifié l'unique solution dans $H^1(\Omega)$ du problème de Poisson $$ \begin {cases} &-\Delta \,\varphi=f=-\nabla \cdot {\bf u} \text{ sur } \Omega,\\ &\partial_n \,\varphi =0\text{ sur } \partial \Omega. \end {cases} $$ </font>
On se donne un champ de vitesse $\bf u$ régulier ($H^1(\Omega)^2$ par exemple) sur un domaine rectangulaire $\Omega$ tel que $\bf u\cdot \bf n=0$. Montrer que la décomposition de Leray $$ {\bf u}={\bf u}_\bot +\nabla \varphi, \text{ dans } \Omega, $$ avec $\varphi$ solution du problème de Neumann homogène.
Montrer que ${\bf u}_\bot$ est bien le projeté orthogonal dans $L^2(\Omega)^2$ sur $V=\{{\bf v} \in H^1(\Omega)^2, \text{ tel que } \nabla \cdot {\bf v} =0 \text{ et } {\bf v}\cdot {\bf n}=0\}$.
Proposer une résolution numérique pour le champ ${\bf u}$ de votre choix.
</font>
On cherche à trouver l'état stationnaire d'un fluide incompressible dans un domaine $\Omega\subset \mathbb R^d$ qui satisfait un minimum de puissance associée aux contraintes visqueuses ajoutées à la puissance d'une force agissant sur le fluide. On néglige ici les effets d'inertie. On cherche ${\bf u}$ qui satisfait $$ \begin {cases} & \displaystyle \min_{{\bf v}\in V}J(v)\\ & \text {avec }J(v)=\int_{\Omega}\frac \mu 2 |\nabla {\bf v}|^2\, -{\bf f}\cdot {\bf v}\,dx\\ & \text { et }V=\{{\bf v}\in (H_0^1(\Omega)^d, \text { tel que } \nabla \cdot {\bf v}=0\}. \end{cases} $$
Si on caractérise ce minimum ${\bf u}$ par dérivation de la fonctionnelle, on obtient que ${\bf u}$ est solution de la formulation faible: $$ \mu \int_{\Omega} \nabla {\bf u}\cdot \nabla {\bf v}\,dx =\int_{\Omega} {\bf f}\cdot {\bf v}\,dx, \text{ pour tout } {\bf v}\in V. $$ Voici une formulation faible pour laquelle une méthode élément fini peut être mise en place. Il est cependant non trivial de fabriquer une base de $V$ compte tenu de la contrainte d'incompressibilité.
On va relaxer la contrainte par une méthode de pénalisation et comprendre ainsi ce qu'est la pression dans le problème de Stokes en écoulement incompressible.
On considère le problème de minimisation approché suivant : $$ \min_{{\bf v}\in (H_0^1(\Omega)^d}J_\varepsilon(v)=\int_{\Omega}\frac \mu 2 |\nabla {\bf v}|^2\,+\frac 1 {2\varepsilon}|\nabla \cdot {\bf v}|^2\, -{\bf f}\cdot {\bf v}\,dx, $$ dont la solution sera notée ${\bf u}_\varepsilon$. On parle de pénalisation de la contrainte car, dans la fonctionnelle de coût, on rajoute à $J(v)$ le coût, très cher avec $\varepsilon>0$ petit, de l'écart à la divergence nulle. Autrement dit, on a rajouté artificiellement dans la fonctionelle de puissance, une puissance qui mesure chèrement la dilatation ou contraction du fluide. L'optimal sera ainsi atteint pour une divergence du fluide petite avec $\varepsilon>0$ petit dans tout le domaine $\Omega$.
La caractérisation du minimum donne la formulation faible $$ \mu \int_{\Omega} \nabla {\bf u}_\varepsilon\cdot \nabla {\bf v}\,dx+\frac 1 {\varepsilon}\int_{\Omega} \nabla \cdot{\bf u}_\varepsilon\,\nabla \cdot{\bf v}\,dx =\int_{\Omega} {\bf f}\cdot {\bf v}\,dx, \text{ pour tout }{\bf v}\in (H_0^1(\Omega)^d. $$ On interprète alors ${\bf u}_\varepsilon$ solution de l'EDP : $$ \begin{cases} & -\mu \Delta {\bf u}_\varepsilon -\frac 1 {\varepsilon} \nabla \nabla \cdot {\bf u}_\varepsilon =f \text{ dans } \Omega,\\ & {\bf u}_\varepsilon = 0 \text{ sur } \partial \Omega. \end {cases} $$ Remarque : Ici, l'interptétation a été facilitée par l'espace de minimisation qui n'est pas contraint par la divergence.
Si on pose $p_\varepsilon=-\frac 1 {\varepsilon}\nabla \cdot {\bf u}_\varepsilon$ et si ce scalaire a le bon goût de converger lorsque $\varepsilon\to 0$, on pourra caractériser, à la limite, l'EDP de Stokes en écoulement incompressible faisant intervenir la pression $p$ limite de $p_\varepsilon$.
Remarque : On a immédiatement que $J_{\varepsilon}({\bf u}_\varepsilon)\le J_{\varepsilon}({\bf u})=J({\bf u})$ ne dépend pas de $\varepsilon$. On en déduit qu'il existe une constant $K$ indépendante de $\varepsilon$ tel que:
(preuve: controler le terme du travail de la force à l'aide du terme visqueux et de l'inégalité de Poincaré.)
On en déduit alors que :
Par passage à la limite quand $\varepsilon$ tend vers zéro, on obtient $$ \begin{cases} & -\mu \Delta {\bf u}_0 + \nabla {p} =f \text{ dans } \Omega,\\ &\nabla \cdot {\bf u}_0 =0\text{ dans } \Omega,\\ & {\bf u}_0= 0 \text{ sur } \partial \Omega. \end {cases} $$ La limite obtenue en vitesse réalise le minimum de $J$ sur $V$ car $J({\bf u }_\varepsilon)\to J({\bf u })$ et $J({\bf u }_\varepsilon)\to J({\bf u }_0)$ d'où $ \bf u =\bf u _0$.
Le couple $( {\bf u},p)$ solution de l'EDP de Stokes incompressible peut être interprété comme le point selle d'un Lagrangien.
Défintion: On appelle $( {\bf u},p)\in V\times H$ le point selle du Lagrangien $\mathcal L$ si $$ \mathcal L({\bf u},q)\le \mathcal L({\bf u},p)\le \mathcal L({\bf v},p), \text { pour tout } ({\bf v},q)\in V\times H. $$ avec $\mathcal L$ convexe en la première variable et concave en la deuxième. Le couple $( {\bf u},p)\in V\times H$ réalise le minimum en la première variable du maximum en la deuxième variable de $\mathcal L$.
Propriété: Le couple $( {\bf u},p)$ solution de l'EDP de Stokes incompressible est l'unique pour selle du Lagrangien $\mathcal L$ sur $V\times H=(H_0^1(\Omega)^d\times \dot L^2(\Omega)$ avec $$ \mathcal L( {\bf u},p)=J({\bf u})+\int_\Omega \nabla p \cdot {\bf u}\,dx. $$ Remarque: Ce formalisme de Lagrangien permet de définir un problème de dualité min-max sur des espaces non contraints et ouvre la porte aux méthodes numériques d'optimisation par dualité. On alterne une recherche de min (méthode de descente ou résolution de l'EDP qui caractérise le min) puis une recherche de max (partielle). Il s'agit de l'algoritme d'Uzawa (recherche du min en ${\bf u}$ à $p$ fixé puis une itération de méthode de gradient à pas fixe dans la deuxième variable) et d'Harrow-Hurvicz (une itération de descente à pas fixe dans la première variable puis dans la deuxième variable).
Remarque: Pour garantir une plus grande vitesse de convergence dans l'algorithme d'Uzawa pour le problème de Stokes, on augmente le Lagrangien par ($r>0$) : $$ \mathcal L_r( {\bf u},p)=\mathcal L( {\bf u},p)+\frac r 2 \int_\Omega |\nabla \cdot {\bf u} |^2\, dx. $$ L'alogorirhme d'Uzawa s'écrit alors, pour le $k$-ème itéré $({\bf u}^k,p^k))$ donné, ${\bf u}^{k+1}$ réalise le minimum de ${\bf u}\to L_r( {\bf u},p^k)$, soit : $$ \begin{cases} & -\mu \Delta {\bf u}^{k+1}-r\nabla \nabla\cdot {\bf u}^{k+1} + \nabla {p}^k =f \text{ dans } \Omega,\\ & {\bf u}^{k+1}= 0 \text{ sur } \partial \Omega \end {cases} $$ et $$ p^{k+1}=p^k-r\nabla \cdot {\bf u}^{k+1}. $$
</font>
On se donne une formulation faible vérifiant les hypothèses du théorème de Lax-Milgram:
Pour tout $v\in V$ (espace de Hilbert), on cherche $u\in V$, $$ a(u,v)=l(v). $$ Soit $(v_n)_{n\in \mathbb N}$ une base de $V$. On note $V_m=Vect(v_0,v_1,\cdots, v_{m-1})$.
On peut définir la formulation faible approchée: Pour tout $v\in V_m$,on cherche $u_m\in V_m$, $$ a(u_m,v)=l(v). $$ L'unique solution $u_m\in V_m$ est une approximation de $u$ grâce au lemme de Céa : il existe $C>0$ tel que : $$ \|u-u_m\|_V\le C\| u - P_{V_m}u\|_V, $$ où $P_{V_m}$ est la projection orthogonale sur $V_m$.
La formulation faible approchée se résume à chercher $u_m=\sum_{0\le j<N} \alpha_j v_j$ solution de $$ a\left(\sum_{0\le j<N} \alpha_j v_j,v_i\right )=l(v_i),\quad 0\le i<N, $$ qui n'est autre, par linéarité de $a(.,.)$, que le système linéaire $Ax=b$ avec $$ \begin {cases} & A_{i,j}=a(v_j,v_i)\\ & x_j=\alpha_j\\ &b_i=l(v_i). \end {cases} $$ Les méthodes éléments finis conformes rentrent dans ce cadre. Le choix d'une base de fonctions trigonométriques s'appuyant sur la transformée de Fourier est un autre choix de méthode de Galerkin. Une base construite sur les vecteurs propres d'un opérateur (laplacien par exemple) intervenant dans une EDP est encore un choix pertinent. </font>
Si on choisit comme base, pour la méthode de Galerkin, les fonctions trigonométriques introduites dans la transformée de Fourier discrète, on dispose d'une méthode efficace pour une classe de problème bien précise. La base est tronquée pour les hautes fréquences. Cette base est constituée de fonctions périodiques, c'est une première limitation.
L'avantage de cette base est qu'il est peu coûteux de projeter sur cette base grâce à la transformée de Fourier rapide (FFT) dont la complexité est en $N\log(N)$ où $N$ est la dimension de la base mais aussi la taille du vecteur de discrétisation de la fonction considérée sur une grille uniforme. La FFT s'applique sur ce vecteur.
L'énorme avantage de ce choix de base est que la transformée de Fourier transforme les dérivations en un multiplateur par une fonction polynôme en Fourier. C'est très pratique pour obtenir des opérateurs de dérivation devenus diagonaux en Fourier et ainsi résoudre trivialement des EDP linéaires à coefficients constants pour des conditions limites périodiques (sur un domaine cartésien).
Si $f$ est $L$-périodique et $x_j=\frac j N L$ pour $0\le j<N$, alors $f_N$ le projeté de $f$ sur cette espace tronqué des hautes fréquences se décompose comme suit : $$ f_N(x)=\sum_{k=0}^{N-1} S_k \exp(2i\pi k\frac x L) $$ et le coefficient de Fourier $S_k$ est donné par $$ S_k=\sum_{j=0}^{N-1} f(x_j) \exp(-2i\pi k \frac {x_j} L) $$ La transformée de fourier rapide (FFT) qui calcule $S_k$ s'appuie sur la formule de Moivre pour faire en sorte qu'un terme de la somme serve en cascade à plusieurs coefficients $S_k$. Ainsi, la complexité en $O(N^2)$ pour le calcul des $N$ coefficients est réduite à $O(N\log(N))$.
Si on veut résoudre en périodique $-u''(x)=f(x)$ ($f$ à moyenne nulle), on résoudra $u_N''(x)=f_N(x)$, $u_N(x)=\sum_{k=0}^{N-1} u_k \exp(2i\pi k\frac x L)$ et donc $u''_N(x)=\sum_{k=0}^{N-1} (2i\pi \frac k L)^2 u_k\exp(2i\pi k\frac {x} L) =-f_N(x)$. Il ne reste plus qu'à identifier les coefficients de fourier pour en déduire $u_k$, $0\le k <N$. En Fourier, le système est ainsi diagonal, chaque coefficient se déduit indépendament des autres (coût $O(N)$).
La reconstruction de $u_N$ ne se fait pas en calculant $u_N(x_j)$ par la formule précédente (complexité en $O(N^2)$) mais en applicant la FFT inverse. </font>
from scipy.fftpack import fft, ifft
import time
L = 8
N = 2**15 # nombre de points d'échantillonnage
print('N=',N)
dx = L/(N)
# construction du signal
x = linspace(0.0,L-dx,N)
s=10*x**2*(x-L/2)**3
s-=s@(0*s+1)/N
#puis le calcul de la TFD par l'algorithme de FFT :
start = time.time()
tfd = fft(s)
#plot(x,tfd,'r');
tfd[1:]=tfd[1:]/(2*pi/L*arange(1,N))**2
z=2*ifft(tfd).real#fft(tfd)/N#z=z[-1:-N-1:-1]
print("elapsed_time =", time.time() - start)
figure(figsize=(8,8))
grid(True)
plot(x,s,'r')
plot(x,z,'b')
plot(x[1:N-1],(2*z[1:N-1]-z[0:-2]-z[2:N])/dx**2,'g-.')
title('Solution EDP et terme source')
xlabel('x')
ylabel("z(x) s(x) -z''(x)");
N= 32768 elapsed_time = 0.0035390853881835938
En prenant pour base des fonctions trigonométriques périodiques réelles, on obtient le même résultat.
from scipy.fftpack import rfft, irfft
import time
# paramètres
L = 8
N = 2**14
print('N=',N)
dx = L/(N)
x = linspace(0.0,L-dx,N)
kk=zeros(N)
kk[1::2]=range(1,int(N/2)+1)
kk[0::2]=range(0,int(N/2))
s=10*x**2*(x-L/2)**3
s-=s@(0*s+1)/N
Uf=rfft(s)
Uf[1:]=Uf[1:]/(2*pi*kk[1:]/L)**2
U=irfft(Uf)
figure(figsize=(8,8))
plot(x,U);
figure(figsize=(8,8))
plot(x[1:N-1],(2*U[1:N-1]-U[0:N-2]-U[2:N])/dx**2);
plot(x[1:-1],s[1:-1],'-.');
N= 16384
En travaillant avec la base de fonction trigonométrique respectant des conditions limites de type Neumann homogène, on peut résoudre le problème de Neumann suivant.
from scipy.fftpack import dct, idct
import time
N=2**12
L=2
dx=L/N
x = linspace(dx/2,L-dx/2,N)
f=x**2#-4/3.#ones(N)#sin(pi*1.5*x/L)
#f[0:int(5*N/6)]=0.
#f=2*ones(N)
#f[0:int(N/2)]=0.
f=f-f@ones(N)/N
#f[0:int(N/10)]=0
#f[-1-int(N/10):-1]=0
start = time.time()
U=dct(f)
U[1:]=U[1:]/(pi*arange(1,N)/L)**2
U=0.5*idct(U)/N
print('simu time',time.time()-start)
figure(figsize=(8,8))
plot(x,U);
uex=-x**4/12+4/3.*x**2/2
uex=uex-uex@ones(N)/N
plot(x,uex,'-.');
print(N,max(abs(uex-U)))
simu time 0.0017199516296386719 4096 2.648359076484752e-08
On peut s'assurer que cette approximation est bien d'ordre 2. Voici ci-dessous un script python pour vérifier l'ordre de la méthode pour un problème de Neuman 1D.
L=2.
def pb_neuman(N):
dx=L/N
x = linspace(dx/2,L-dx/2,N)
f=x**2
f=f-f@ones(N)/N
start = time.time()
U=dct(f)
U[1:]=U[1:]/(pi*arange(1,N)/L)**2
U=0.5*idct(U)/N
simu_time=time.time()-start
uex=-x**4/12+4/3.*x**2/2
uex=uex-uex@ones(N)/N
return [max(abs(uex-U)),simu_time]
errtab=[]
Ntab=[]
comp=[]
for l in range(8,21):
N=2**l+10
Ntab.append(N)
errtab.append(pb_neuman(N)[0])
comp.append(pb_neuman(N)[1])
Ntab=array(Ntab)
errtab=array(errtab)
comp=array(comp)
plot(log(Ntab),log(errtab),'o')
plot(log(Ntab),polyfit(log(Ntab),log(errtab), 1)[0]*log(Ntab)+polyfit(log(Ntab),\
log(errtab), 1)[1])
print("l'ordre de la méthode est :",-polyfit(log(Ntab),log(errtab), 1)[0])
l'ordre de la méthode est : 1.9995112156059645
On peut étendre la résolution des EDPs à coefficients constants sur domaine cartésien en dimension supérieure. Pour le problème de Neumann en dimension 2, traité ci-dessous, une comparaison en temps calcul est donnée avec une approche par diférences finies ou volume fini sur grille cartésienne.
def Laplace(N,M,l,h):
D0=(2/(l**2)+2/(h**2))*ones(N*M)# diagonale principale
D1=-1/l**2*ones(N*M)# surdiagonale
D1[N::N]=0.#correction de la surdiagonale (voisin de droite n existe pas au bord droit)
D0[0::N]-=1/l**2
D0[N-1::N]-=1/l**2
DM1=-1/l**2*ones(N*M)# sousdiagonale
DM1[N-1::N]=0.#correction de la sousdiagonale (voisin de gauche n existe pas au bord gauche)
DN=-1/h**2*ones(N*M)
D0[0:N]-=1/h**2
D0[N*(M-1):]-=1/h**2
A=spdiags(D0,[0],N*M,N*M)+spdiags(D1,[1],N*M,N*M)+spdiags(DM1,[-1],N*M,N*M)
A=A+spdiags(DN,[N],N*M,N*M)+spdiags(DN,[-N],N*M,N*M)
return A.tocsr()
from scipy.fftpack import dctn, idctn
N=1400
M=1600
L=2
H=2
dx=L/N
dy=H/M
LAP=Laplace(N,M,dx,dy)
x = linspace(dx/2,L-dx/2,N)
y = linspace(dy/2,H-dy/2,M)
X, Y = meshgrid(x,y)
f=(X-0.5*L)**2+(Y-0.5*H)**2
#f=(4*pi**2/L**2+9*pi**2/H**2)*cos(2*pi*X/L)*cos(3*pi*Y/H)
f=f-reshape(f,N*M)@ones(N*M)/(N*M)
start=time.time()
U=dctn(f)
Ti,Tj=meshgrid(range(N),range(M))
Ti[0,0]=1.
Tj[0,0]=1.
U=U/(pi**2*Ti**2/L**2+pi**2*Tj**2/H**2)
U=0.25*idctn(U)/(N*M)
#U=U/sqrt(2)
print('grille',N,M,'temps simu fft',time.time()-start)
start=time.time()
U_DF=spsolve(LAP,reshape(f,N*M))
U_DF-=U_DF@ones(N*M)/(N*M)
print('grille',N,M,'temps simu difference finie',time.time()-start)
err=reshape(U,N*M)-U_DF
print(max(abs(err))/max(abs(U_DF)))
grille 1400 1600 temps simu fft 0.11378121376037598 grille 1400 1600 temps simu difference finie 68.54114484786987 2.2480726889205454e-06
from mpl_toolkits import mplot3d
fig = figure(figsize=(8,8))
ax = fig.gca(projection='3d')
ax.plot_surface(X[::10,::10], Y[::10,::10], U[::10,::10], rstride=1, cstride=1,
cmap='viridis', edgecolor='none');
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u');
ax.view_init(30, 30)