4.2 Équation de Poisson

4.2.1 Problème physique

On considère une membrane carrée de coté $a$ qui se déforme sous l'effet d'une charge surfacique $f(x,y)$ . La membrane est sous tension et fixée sur les bords. On suppose qu'en chacun des points la tension $T$ est constante et tangente à la membrane (on néglige les forces élastiques dues à la déformation de la membrane). On note $u(x,y)$ la déformée. Les forces exercées sur un élément de membrane $dx$ $dy$ sont:

Figure 4.1: Équilibre statique d'une membrane dans le plan (x,z)

1. des forces de tension exercées sur les cotés de l'élément. Pour le coté de longueur $dy$ et d'abcisse $x$ , cette force est perpendiculaire au coté $dy$ et tangente à la surface $u(x,y)$ . Elle est située dans le plan $(x,z)$ et a pour composantes:
   
   $$Tdy\left[\begin{array}{c} -\cos\theta_{x}\\ 0\\ \sin\theta_{x}\end{array}\right]$$
   en notant $\theta_{x}$ l'angle de la surface $u(x,y)$ avec l'horizontal. De même pour le coté d'abscisse $y$ , la force de tension s'écrit:
   
   $$Tdx\left[\begin{array}{c} 0\\ -\cos\theta_{y}\\ \sin\theta_{y}\end{array}\right]$$

2. des forces de chargement verticales:
   
   $$f(x,y)dxdy$$

L'équilibre statique conduit donc aux équations suivantes:

$$Tdy(\cos\theta_{x+dx}-\cos\theta_{x}) =$$ 0

$$Tdx(\cos\theta_{y+dy}-\cos\theta_{y}) =$$ 0

$$Tdy(\sin\theta_{x}-\sin\theta_{x+dx})+Tdx(\sin\theta_{y}-\sin\theta_{y+dy}) = f(x,y)dxdy$$

avec $\tan\theta_{x}=\frac{\partial u}{\partial x}$ et $\tan\theta_{y}=\frac{\partial u}{\partial y}$ . En supposant que les angles $\theta$ sont petits, i.e:

$$\cos\theta_{x}\approx1, \sin\theta_{x}\appro... ...s\theta_{y}\approx1, \sin\theta_{y}\approx\frac{\partial u}{\partial y}$$

On effectue des développements limités à l'ordre 1 dans les équations précédentes, qui conduisent à l'équation d'équilibre de la membrane:

$$-\frac{\partial}{\partial x}\left(T\frac{\partial u}{\partial x}\... ...)-\frac{\partial}{\partial y}\left(T\frac{\partial u}{\partial y}\right)=f(x,y)$$ (4.1)

auquel on ajoute les conditions aux limites:

$$u(0,y)=u(a,y)=u(x,0)=u(x,a)=0$$

En effectuant un changement de variables, on obtiens le problème modèle suivant, qui est une équation de Poisson:

$$-\Delta U=F \Omega=]0,1[*]0,1[$$ (4.2)

$$U(0,y)=U(1,y)=U(x,0)=U(x,1)=0$$

4.2.2 Étude de la solution analytique

Pour déterminer la solution générale de l'équation de Poisson (4.2), on décompose $U(x,y)$ en série de Fourier vérifiant les conditions aux limites:

$$U(x,y)=\sum_{k=1}^{\infty}\sum_{p=1}^{\infty}U_{kp}\sin\left(k\pi x\right)\sin\left(p\pi y\right)$$ (4.3)

En remplaçant dans (4.2), on obtiens:

$$\sum_{k=1}^{\infty}\sum_{p=1}^{\infty}(k^{2}+p^{2})\pi^{2}U_{kp}\sin\left(k\pi x\right)\sin\left(p\pi y\right)=F(x,y)$$

d'où les valeurs de $U_{kp}$ , en multipliant cette relation par $\sin\left(k\pi x\right)\sin\left(p\pi x\right)$ . En intégrant sur le domaine $\Omega$ et en utilisant l'orthogonalité des fonctions $\sin(k\pi x)$ , il vient:

$$U_{kp}=\frac{4}{(k^{2}+p^{2})\pi^{2}}\int_{0}^{1}\int_{0}^{1}F(x,y)\sin\left(k\pi x\right)\sin\left(p\pi y\right)dxdy$$ (4.4)

4.2.2.1 cas d'un chargement constant $F=-1$

Dans le cas d'un chargement constant $F=-1$ , la valeur du coefficient de Fourier $U_{kp}$ se calcule simplement avec Maple et on trouve:

$$U_{kp}=\frac{-4}{(k^{2}+p^{2})\pi^{4}kp}\left(1-(-1)^{k}\right)\left(1-(-1)^{p}\right)$$

Ce coefficient est non nul si et seulement si $k$ et $p$ sont tous les deux impaires. La solution exacte s'écrit donc:

$$U(x,y)=\sum_{l=1}^{\infty}\sum_{m=1}^{\infty}U_{lm}\sin\left((2l-1)\pi x\right)\sin\left((2m-1)\pi y\right)$$

$$U_{lm}=\frac{-16}{((2l-1)^{2}+(2m-1)^{2})\pi^{4}(2l-1)(2m-1)}$$

L'allure de la solution est donnée sur la figure (4.2).

Figure 4.2: solution exacte de (4.2) pour $F=-1$

La valeur maximale $U_{max}$ de la déformation se trouve au centre et a pour expression:

$$U_{max}=\sum_{l=1}^{\infty}\sum_{m=1}^{\infty}\frac{16(-1)^{l+1}(-1)^{m+1}}{((2l-1)^{2}+(2m-1)^{2})\pi^{4}(2l-1)(2m-1)}$$

On peut calculer une valeur approchée très précise de cette série avec Maple, et on trouve (pour $m=l=200$ ):

$$U_{max}=-0.07367135123$$ (4.5)

4.2.3 Schéma aux différences finies pour le laplacien

Sur un maillage de $N_{x}$ points suivant $x$ et $N_{y}$ points suivant $y$ (figure 4.3), la discrétisation par différences finies centrées de l'équation (4.2) s'écrit:

$$\frac{-U_{i+1,j}+2U_{i,j}-U_{i-1,j}}{dx^{2}}+\frac{-U_{i,j+1}+2U_{i,j}-U_{i,j-1}}{dy^{2}} = F_{i,j}$$ (4.6)

$$\forall i=2,N_{x}-1 \forall j=2,N_{y}-1$$

avec les conditions aux limites:

$$U_{1,j}=U_{N_{x},j}=0 \forall j=1,N_{y}$$   et  $$U_{i,1}=U_{i,N_{y}}=0 \forall i=1,N_{x}$$

Les pas de discrétisation en espace sont équidistants et vérifient $dx=\frac{1}{N_{x}-1}$ et $dy=\frac{1}{N_{y}-1}$ .

Figure 4.3: discrétisation différences finies du laplacien

Ce schéma conduit à un système matriciel de $N=N_{x}N_{y}$ inconnues $U_{i,j}$ . Pour écrire ce système sous la forme matricielle $\mathcal{A}X=\mathcal{B}$ , on doit transformer la matrice des inconnues $U_{i,j}$ en un vecteur inconnu $X_{k}$ . Pour cela on numérote les inconnues par lignes, i.e. on effectue la transformation d'indice $(i,j)$ vers le mono-indice $k=i+(j-1)N_{x}$ . Avec ce changement d'indice, l'équation aux différences (4.6) s'écrit:

$$cU_{k-N_{x}-1}+bU_{k-1}+aU_{k}+bU_{k+1}+cU_{k+N_{x}+1}=F_{k}$$

pour tous les noeuds internes $k=i+(j-1)N_{x}$ avec $1<i<N_{x}$ et $1<j<N_{y}$ en notant $a=\frac{2}{dx^{2}}+\frac{2}{dy^{2}}$ , $b=\frac{-1}{dx^{2}}$ et $c=\frac{-1}{dy^{2}}$ .

Les conditions aux limites s'écrivent $U_{k}=0$ pour les noeuds frontières $k=1+(j-1)N_{x}$ , $k=N_{x}+(j-1)N_{x}$ avec $1<j<N_{y}$ et $k=i$ , $k=i+(N_{y}-1)N_{x}$ avec $1<i<N_{x}$ .

La matrice $\mathcal{A}$ est une matrice penta-diagonale dont la forme est donnée sur la figure (4.4a). On vérifie que la matrice possède bien au maximum 5 coefficients non nuls répartis sur la diagonale de coefficients $a$ , les 2 co-diagonales adjacentes de coefficients $b$ et les 2 co-diagonales distantes de $N_{x}-1$ de la diagonale de coefficients $c$ . On constate aussi que la matrice $\mathcal{A}$ est non symétrique, à cause de la façon d'appliquer les conditions aux limites. En effet pour un noeud $k$ sur la frontière, on applique la condition aux limites $X_{k}=0$ dans la ligne $k$ de la matrice en annulant la ligne et en mettant $1$ sur la diagonale. On ne tiens pas compte de cette condition aux limites dans les équations où intervient la valeur de $X_{k}$ , i.e. dans les lignes de $\mathcal{A}$ ayant un coefficient non nul dans la colonne $k$ . Pour conserver la symétrie de la matrice, qui traduit la symétrie du problème physique, il faut aussi annuler les coefficients de la colonne $k$ (figure 4.4b). Dans le cas $X_{k}=X_{0}$ , il faut en outre retrancher la colonne $X_{0}\mathcal{A}_{k,j}$ du second membre $\mathcal{B}$ .

Figure 4.4: matrice du laplacien pour $N_{x}=7$ et $N_{y}=7$

[matrice initiale]

[matrice symétrisée]

On note aussi que le nombre de coefficients non nuls de la matrice $\mathcal{A}$ est de l'ordre de $5N_{x}N_{y}\approx5N^{2}$ , ce qui est beaucoup plus petit que le nombre de coefficients $N_{x}^{2}N_{y}^{2}\approx N^{4}$ de la matrice $\mathcal{A}$ . On tiendra compte de ces propriétés lors de la résolution du système matriciel.

4.2.3.1 Précision et erreur de troncature

On utilise une discrétisation centrée et d'ordre 2 de la dérivée seconde en $x$ et en $y$ , donc l'erreur de troncature du schéma est d'ordre $O(dx^{2},dy^{2})$ . Elle s'écrit:

$$ErrT=\frac{1}{12}\frac{\partial^{4}U}{\partia... ...rac{1}{12}\frac{\partial^{4}U}{\partial y^{4}} dy^{2}+O(dx^{4},dy^{4})$$

La précision du schéma (4.6) est donc d'ordre 2 en espace, i.e. en $O(dx^{2},dy^{2})$ .

4.2.4 Expérimentation numérique avec Matlab

La fonction Matlab laplace2d (4.2.4) calcule la matrice $\mathcal{A}$ et le second membre $\mathcal{B}$ sur un maillage différences finis de $N_{x}N_{y}$ points pour une fonction $F$ définie aux noeuds $(i,j)$ du maillage et des conditions aux limites homogènes. Compte tenu des remarques sur la structure de la matrice $\mathcal{A}$ , on utilise une structure de donnée de matrice creuse (sparse matrix en anglais) qui permet de ne stocker que les éléments non nuls. Pour cela on stocke les coefficients non nuls de $\mathcal{A}$ dans un vecteur Ac, ainsi que leurs indices (i,j) dans deux autres vecteurs I et J. Pour la matrice $\mathcal{A}$ ci dessous:

$$\mathcal{A}=\left[\begin{array}{ccccc} 4 & -1... ...4 & -1 & 0\\ 0 & 0 & -1 & 4 & 0\\ -1 & 0 & 0 & 0 & 4\end{array}\right]$$

on utilise les 3 tableaux suivants:

$$Ac = \left[\begin{array}{ccccccccccc} 4 & -1 & -1 & -1 & 4 & 4 & -1 & -1 & 4 & -1 & 4\end{array}\right]$$

$$I = \left[\begin{array}{ccccccccccc} 1 & 1 & 1 & 2 & 2 & 3 & 3 & 4 & 4 & 5 & 5\end{array}\right]$$

$$J = \left[\begin{array}{ccccccccccc} 1 & 2 & 5 & 1 & 2 & 3 & 4 & 3 & 4 & 1 & 5\end{array}\right]$$

Pour utiliser cette structure de données avec Matlab, on initialise $\mathcal{A}$ avec la fonction Matlab spalloc (ligne 18), qui permet de créer une matrice creuse de $5N$ éléments, au lieu d'utiliser la fonction zeros(N,N) qui crée une matrice carrée de $N^{2}$ éléments. Pour $N_{x}=N_{y}=100$ , on a $N=10^{4}$ et le stockage de $A$ sous forme de matrice creuse nécessite alors $640$ kilo-octets de mémoire au lieu des $100$ méga-octets nécessaire au stockage de tous les coefficients de $\mathcal{A}$ , ce qui rend la résolution du problème possible sur un ordinateur de bureau. On remarque aussi que la taille nécessaire au stockage de $\mathcal{A}$ est supérieur à $5N$ réels (soit $400$ kilo-octets) puisque pour chaque valeur non nulle $\mathcal{A}_{pq}$ on stocke aussi les indices $p$ et $q$ correspondants.

programme Matlab 4.2.4: Fonction laplace2d: calcul de la matrice du laplacien

function [A,B]=laplace2d(F,nx,ny)
% entree:
% matrice F du second membre Fij valeur de F au noeud (i,j)
% nx,ny  nombre de points en x et en y
% sortie:
% matrice A et second membre B
% en utilisant un stocage creux 
% pble  -lap(U)=f avec des C.L. homogene
dx=1/(nx-1); dy=1/(ny-1);
% coefficiant du shema pour les nds
% (i,j-1) (i-1,j) (i,j) (i+1,j) (i,j+1)
coeff=[-1/dy^2,-1/dx^2,2/dx^2+2/dy^2,-1/dx^2,-1/dy^2];
% decalage / au neud (i,j) dans la numérotation
num = [  -nx, -1, 0, 1, nx];
% assemblage de la matrice pour un stocage par ligne
% i.e le noeud (i,j) a pour adresse k=(j-1)*nx+i
N=nx*ny; % dimension
A=spalloc(N,N,5*N); % matrice creuse de 5 elts maxi / ligne
B=reshape(F,N,1);
for i=2:nx-1
    for j=2:ny-1
        k=(j-1)*nx+i;
        A(k,k+num)=coeff;
    end
end
% conditions aux limites
% ======================
% C.L. sur les frontieres i=1,i=nx (Dirichlet)
for j=1:ny
    k=(j-1)*nx+1; A(k,:)=0; A(:,k);A(k,k)=1.0; B(k)=0;
    k=(j-1)*nx+nx; A(k,:)=0; A(:,k)=0; A(k,k)=1.0; B(k)=0;
end;
% C.L. sur les frontieres j=1,j=ny (Dirichlet)
for i=1:nx
    k=i; A(k,:)=0; A(:,k)=0; A(k,k)=1.0; B(k)=0;
    k=(ny-1)*nx+i; A(k,:)=0; A(:,k)=0; A(k,k)=1.0; B(k)=0;
end;
% fin

Pour le second membre $\mathcal{B}$ on transforme la matrice $F$ des valeurs aux noeuds du maillage en un vecteur de colonne de dimension $N$ avec la fonction Matlab reshape. La boucle d'assemblage de $\mathcal{A}$ correspond aux lignes 20 à 25, et se fait ligne par ligne en utilisant la structure matricielle de Matlab (ligne 23). Les conditions aux limites sont imposées sur les cotés $x=0$ et $x=1$ (lignes 29 à 32) et $y=0$ et $y=1$ (lignes 34 à 37), en annulant la ligne et la colonne $k$ de $\mathcal{A}$ ainsi que le second membre $\mathcal{B}_{k}$ puis en imposant $\mathcal{A}_{k,k}=1$ .

Le script Matlab (4.2.4) résout numériquement le problème (4.2) dans le cas d'une fonction $F=-1$ (ligne 8). On utilise la fonction Laplace2d précédente pour calculer la matrice $\mathcal{A}$ et le second membre $\mathcal{B}$ du problème. On résout le système (ligne 11) avec l'opérateur standard $\setminus$ de Matlab, qui pour des matrices creuses utilisent un algorithme de Gauss par bande. C'est la méthode la plus efficace sous Matlab, même si elle nécessite un stockage temporaire important, de l'ordre de $2N_{x}^{2}N_{y}$ puisque la largeur de bande vaut $N_{x}$ . En utilisant la structure particulière de la matrice $\mathcal{A}$ (tri-diagonale par blocs), on pourrait utiliser un algorithme très efficace, qui est l'extension de l'algorithme de Thomas. Son principe est décrit dans l'annexe 3DBloc, mais il n'est pas implémenté sous Matlab.

La fin du script permet la visualisation en 3D de la solution calculée.

programme Matlab 4.2.4: Résolution du problème (4.2)

% resolution du laplacien
clear
% dimension
nx=21; ny=21;
dx=1/(nx-1); dy=1/(ny-1);
X=[0:dx:1]; Y=[0:dy:1];
% assemblage matrice et 2nd membre
F=-ones(nx,ny);
[A,B]=laplace2d(F,nx,ny);
% resolution
U=A\B;
% transformation de U en matrice pour visualisation
U1=reshape(U,nx,ny);
% visualisation
surfc(X,Y,U1); title('deformee'); shading interp;

On a tracé le résultat obtenu pour $N_{x}=N_{y}=21$ sur la figure (4.5). En comparant avec la solution exacte (figure 4.2), on constate une bonne concordance. L'erreur relative sur la valeur maximale de la déformée est inférieure à 3%:

$$Err=\left\vert\frac{U_{max}-U_{max}^{h}}{U_{max}}\right\vert\approx0.027$$

où on a noté $U_{max}$ la valeur exacte (4.5) et $U_{max}^{h}$ est la valeur approchée.

Figure 4.5: solution numérique ( $N_{x}=N_{y}=11$ ) du problème (4.2)

Pour terminer cette étude, nous avons effectué une étude de précision en calculant cette erreur relative pour différents maillages (avec un nombre de points identique suivant $x$ et $y$ ). Chaque maillage est caractérisée par un pas de discrétisation $h=dx=dy$ . Sur la figure (4.6), on a tracé l'évolution de l'erreur relative $Err$ en fonction de $h$ , et on trouve sur une échelle logarithmique une droite de pente 2. Cela montre que l'erreur relative est en $O(h^{2})$ , ce qui confirme le calcul de l'erreur de troncature en $O(dx^{2},dy^{2})$ .

Figure 4.6: Erreur relative en fonction du pas $h$ du maillage