suivant: 4.3 Équation des ondes monter: 4.

Schémas différences
finies précédent: 4.1 Introduction Table des matières

Sous-sections
 4.2.1 Problème physique
 4.2.2 Étude de la solution analytique

o 4.2.2.1 cas d'un chargement constant

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

o 4.2.3.1 Précision et erreur de troncature

 4.2.4 Expérimentation numérique avec Matlab

4.2 Équation de Poisson

4.2.1 Problème physique

On considère une membrane carrée de coté qui se déforme sous l'effet d'une charge

surfacique . La membrane est sous tension et fixée sur les bords. On suppose qu'en
chacun des points la tension est constante et tangente à la membrane (on néglige les forces

élastiques dues à la déformation de la membrane). On note la déformée. Les forces

exercées sur un élément de membrane 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 et d'abcisse , cette force est

 perpendiculaire au coté et tangente à la surface .

Elle est située dans le plan et a pour composantes:

en notant l'angle de la surface avec l'horizontal. De

même pour le coté d'abscisse , la force de tension s'écrit:

2. des forces de chargement verticales:

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

avec et . En supposant que les angles sont petits, i.e:

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:

(4.1
)

auquel on ajoute les conditions aux limites:

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

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

en série de Fourier vérifiant les conditions aux limites:

(4.3
)

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

d'où les valeurs de , en multipliant cette relation par . En

intégrant sur le domaine et en utilisant l'orthogonalité des fonctions , il vient:

 (4.4
)

4.2.2.1 cas d'un chargement constant

Dans le cas d'un chargement constant , la valeur du coefficient de Fourier se

calcule simplement avec Maple et on trouve:

Ce coefficient est non nul si et seulement si et sont tous les

deux impaires. La solution exacte s'écrit donc:

avec

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

 Figure 4.2: solution exacte de (4.2)

pour

La valeur maximale de la déformation se trouve au centre et a pour expression:

On peut calculer une valeur approchée très précise de cette série avec Maple, et on trouve

(pour ):
(4.5
)

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

Sur un maillage de points suivant et points suivant (figure 4.3), la discrétisation

par différences finies centrées de l'équation (4.2) s'écrit:

(4.6
)
avec les conditions aux limites:

et

Les pas de discrétisation en espace sont équidistants et vérifient et

.

Figure 4.3: discrétisation différences finies du

laplacien

Ce schéma conduit à un système matriciel de inconnues . Pour écrire ce

système sous la forme matricielle , on doit transformer la matrice des

inconnues en un vecteur inconnu . Pour cela on numérote les inconnues par lignes,

i.e. on effectue la transformation d'indice vers le mono-indice .

Avec ce changement d'indice, l'équation aux différences (4.6) s'écrit:

pour tous les noeuds internes avec

et en notant , et .
Les conditions aux limites s'écrivent pour les noeuds

frontières , avec et

, avec .
La matrice 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 , les 2 co-diagonales adjacentes de coefficients et les 2 co-

diagonales distantes de de la diagonale de coefficients . On constate aussi que la

matrice est non symétrique, à cause de la façon d'appliquer les conditions aux limites. En

effet pour un noeud sur la frontière, on applique la condition aux limites dans la

ligne de la matrice en annulant la ligne et en mettant sur la diagonale. On ne tiens pas

compte de cette condition aux limites dans les équations où intervient la valeur de , i.e.

dans les lignes de ayant un coefficient non nul dans la colonne . 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 (figure 4.4b). Dans le cas , il faut en outre retrancher

la colonne du second membre .

[matrice initiale] [matrice symétrisée]

 Figure 4.4: matrice du laplacien pour et

On note aussi que le nombre de coefficients non nuls de la matrice est de l'ordre

de , ce qui est beaucoup plus petit que le nombre de

coefficients de la matrice . 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 et en , donc

l'erreur de troncature du schéma est d'ordre . Elle s'écrit:

La précision du schéma (4.6) est donc d'ordre 2 en espace, i.e. en .

4.2.4 Expérimentation numérique avec Matlab

La fonction Matlab laplace2d (4.2.4) calcule la matrice et le second membre sur un

maillage différences finis de points pour une fonction définie aux noeuds du
maillage et des conditions aux limites homogènes. Compte tenu des remarques sur la structure
de la matrice , 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 dans un vecteur Ac, ainsi que leurs indices (i,j) dans deux autres
vecteurs I et J. Pour la matrice ci dessous:
on utilise les 3 tableaux suivants:

Pour utiliser cette structure de données avec Matlab, on

initialise avec la fonction Matlab spalloc (ligne 18), qui permet

de créer une matrice creuse de éléments, au lieu d'utiliser la

fonction zeros(N,N) qui crée une matrice carrée de éléments.

Pour , on a et le stockage de sous forme

de matrice creuse nécessite alors kilo-octets de mémoire au

lieu des méga-octets nécessaire au stockage de tous les

coefficients de , 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 est supérieur à réels (soit

kilo-octets) puisque pour chaque valeur non nulle on stocke

aussi les indices et 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 on transforme la matrice des valeurs aux noeuds du maillage

en un vecteur de colonne de dimension avec la fonction Matlab reshape. La boucle

d'assemblage de 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 et (lignes 29 à 32) et et (lignes 34 à 37), en annulant la

ligne et la colonne de ainsi que le second membre puis en imposant .

Le script Matlab (4.2.4) résout numériquement le problème (4.2) dans le cas d'une

fonction (ligne 8). On utilise la fonction Laplace2d précédente pour calculer la

matrice et le second membre du problème. On résout le système (ligne 11) avec

l'opérateur standard 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 puisque la largeur de bande vaut .

En utilisant la structure particulière de la matrice (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 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%:

où on a noté la valeur exacte (4.5) et est la valeur approchée.

 Figure 4.5: solution numérique ( ) 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 et ).

Chaque maillage est caractérisée par un pas de discrétisation . Sur la figure

(4.6), on a tracé l'évolution de l'erreur relative en fonction de , et on trouve sur une

échelle logarithmique une droite de pente 2. Cela montre que l'erreur relative est en ,

ce qui confirme le calcul de l'erreur de troncature en .

Figure 4.6: Erreur relative en fonction du pas

du maillage

suivant: 4.3 Équation des ondes monter: 4. Schémas différences

finies précédent: 4.1 Introduction Table des matières

Pr. Marc BUFFAT

[email protected]
2008-04-07