MACS 2 EDP

Sup'Galilée Année 2018-2019

1
Travaux pratiques - Méthodes des différences finies 2D

TP 6 : Matrices creuses, Laplacien 2D et plus si affinité

Table des matières

1 Les matrices creuses 2
1.1 Stockage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.1 Occupation mémoire . . . . . . . . . . . . . . . . . . . . 3
1.1.2 Insertion d'éléments . . . . . . . . . . . . . . . . . . . . 3
1.1.3 Construction ecace de matrices creuses . . . . . . . . . 4
1.2 Matrices diagonales et tridiagonales . . . . . . . . . . . . . . . 4
1.3 Produit de Kronecker de deux matrices . . . . . . . . . . . . . 5

2 Notations et approximation des dérivées secondes 5

3 Equation de Poisson avec conditions de Dirichlet 6
3.1 Assemblage de la matrice sans conditions aux limites . . . . . . 7
3.2 Quelques indices . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.2.1 Indices des bords . . . . . . . . . . . . . . . . . . . . . . 8
3.2.2 Evaluation de fonctions sur la grille . . . . . . . . . . . 9
3.2.3 Evaluation de fonctions sur une partie de la grille . . . . 11
3.3 Assemblage du second membre sans conditions aux limites . . 11
3.4 Prise en compte des conditions aux limites . . . . . . . . . . . . 12
3.5 Résolution du système . . . . . . . . . . . . . . . . . . . . . . . 13
3.5.1 Méthode 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.5.2 Méthode 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 13
En résolvant numériquement des EDPs par des méthodes de type diérences nies, éléments nis et/ou
volumes nis, on abouti souvent à la résolution de grands systèmes linéaires creux (i.e. les matrices des systèmes
linéaires contiennent un très grand nombre de zéros). Au niveau mémoire (de l'ordinateur), il est alors recom-
mandé (voir necessaire) d'utiliser un stockage particulier pour ces matrices : CSC, CSR, COO, LIL, DOK, ...
Ces diérentes formats permettent de ne stocker en mémoire que les éléments non nuls des matrices.
Sous Matlab/Octave, le format utilisé est le format CSC (Compressed Sparse Column). Il est très ecace
pour les operations arithmétiques et les produits matrices vecteurs. Nous donnons deux utilisations (parmis
d'autres) de la fonction sparse de Matlab/Octave :
‚ M= sparse(m,n);
Cette commande crée la matrice creuse nulle M de dimension mˆn
‚ M = sparse(I,J,K,m,n);
Cette commande crée la matrice creuse M de dimension mˆn telle que
M(Ig(k),Jg(k)) Ð M(Ig(k),Jg(k)) + Kg(k).
Les vecteurs Ig, Jg et Kg ont la même longueur. Les éléments nuls de Kg ne sont pas pris en compte et
les éléments de Kg ayant les mêmes indices dans Ig et Jg sont sommés.
On donne maintenant dans plusieurs langages interprétés l'analogue de ces deux commandes :
‚ Python ([Link] module) :
 M=sparse.<format>_matrix((m,n))
 M=sparse.<format>_matrix((Kg,(Ig,Jg)),shape=(m,n))
où <format> est le format de stockage de la matrice creuse (i.e. csc, csr, lil, ...),
‚ R : format CSC (seulement ?),
 M=sparse(i=,j=,dims=c(m,n))
1. En cours de rédaction, version du 19 décembre 2018 à 13:41

1
  M=sparse(i=Ig,j=Jg,x=Kg,dims=c(m,n))
‚ Scilab : format row-by-row( ?)
 M=sparse([],[],[m,n])
 M=sparse([Ig,Jg],Kg,[m,n])
Dans tous ces langages, de très nombreuses fonctions existent permettant d'utiliser les matrices creuses
comme des matrices classiques. Il est donc possible d'eectuer des opérations aussi élémentaires que des sommes,
produits de matrices et vecteurs mais aussi des opérations plus complexes comme des factorisations (LU, Cho-
lesky,QR, ...), résolution de systèmes linéaires (méthodes directes ou itératives), calcul de valeurs propres et
vecteurs propres, ...
Pour les langages de programmation usuels (Fortan, C, C++, ...) des libraires existent permettant de réaliser
ces mêmes opérations sur les matrices creuses :
‚ C/C++ avec les librairies SuiteSparse de T. Davis [?] ainsi que BLAS (Basic Linear Algebra Subroutines),
LAPACK (Linear Algebra PACKage) et leurs dérivées.
‚ CUDA avec la librairie cuSparse [?] de Nvidia.

1 Les matrices creuses

Nous avons commencé à utiliser les matrices creuses (sparse matrix ) dans les TPs précédents. Par exemple,
pour générer la matrice K P Md pRq dénie par
¨ ˛
´2 1 0 ... 0
.. .. ‹
. . ‹
˚
˚ 1 ´2 1
K “ ˚ 0 ... .. ..
˚ ‹
. . (1.1)
˚ ‹
0‹
˚ . .. .. ..
˚ ‹
˝ .. . . .
‹
1‚
0 ... 0 1 ´2

nous avons écrit une fonction Lap1D pouvant ressembler à ceci :

Listing 1  Fonction Lap1D

function K=Lap1D ( d )
K=sparse ( d , d ) ;
K ( 1 , 1 ) =´2; K ( 1 , 2 ) =1;
for i =2: d´1
K ( i , i ) =´2;
K ( i , i ´ 1)=1;
K ( i , i +1)=1;
end
K ( d , d ) =´2; K ( d , d ´ 1)=1;
end

Toutefois cette manière de construire une matrice creuse est très loin d'être la plus ecace.
Dans un premier temps, une rapide explication du pourquoi vous est proposée et ensuite nous verrons
comment améliorer celà en utilisant judicieusement la fonction sparse . Il est aussi possible d'utiliser conjoin-
tement les fonctions speye et spdiags pour créer la matrice K mais ces commandes ne seront plus utilisables
pour assembler les matrices provenant des méthodes d'éléments nis et volumes nis, c'est pourquoi nous nous
focaliserons sur la fonction sparse .

1.1 Stockage

Sous Matlab et Octave, une matrice creuse (ou sparse matrix ), A P MM,N pRq, est stockée sous le format
CSC (Compressed Sparse Column), utilisant trois tableaux :

iap1 : nnzq, jap1 : N ` 1q et aap1 : nnzq,

2
où nnz est le nombre d'éléments non nuls de la matrice A. Ces tableaux (cachés mais contenus dans l'object
sparse Matlab/Octave) sont dénis par
‚ aa : contient l'ensemble des nnz éléments non nuls de A stockés par colonnes.
‚ ia : contient les numéros de ligne des éléments stockés dans le tableau aa.
‚ ja : permet de retrouver les éléments d'une colonne de A, sachant que le premier élément non nul de la
colonne k de A est en position japkq dans le tableau aa. On a jap1q “ 1 et japN ` 1q “ nnz ` 1.
Pour illustrer, voici un exemple simple avec la matrice
¨ ˛
1. 0. 0. 6.
A “ ˝0. 5. 0. 4.‚
0. 1. 2. 0.
on a M “ 3, N “ 4, nnz “ 6 et

aa 1. 5. 1. 2. 6. 4.

ia 1 2 3 3 1 2

ja 1 2 4 5 7

Le premier élément non nul de la colonne k “ 3 de A est 2, ce nombre est en position 4 dans aa, donc jap3q “ 4.

1.1.1 Occupation mémoire

Pour illustrer le gain mémoire, nous allons calculer la taille mémoire occupée par la matrice tridiagonale K
dénie en (1.1) dans le cas d'un stockage plein (full matrix ) ou creux (sparse matrix ). Cette matrice contient
d ` 2pd ´ 1q “ 3d ´ 2 éléments non nuls.
Sous Matlab/Octave le stockage creux de cette matrice va utiliser p3d ´ 2q double (tableau aa), p3d ´ 2q
int (tableau ia) et pd ` 1q int (tableau ja) pour un total de
p3d ´ 2q double et p4d ´ 1q int.

Or un double occupe 8 Octets (en anglais bytes ) et un int 4 Octets. La place mémoire occupée par la matrice
stockée sous forme sparse est alors de

p3d ´ 2q ˆ 8 ` p4d ´ 1q ˆ 4 “ p40 d ´ 20q Octets.

Par contre, si la matrice était stockée classiquement (full), alors la place mémoire serait de

d2 Octets.

On rappelle que 1 Mo (Méga-octet) correspond à 106 Octets, 1 Go (Giga-octet) correspond à 109 Octets et 1 To
(Tera-octet) correspond à 1012 Octets. On donne maintenant la place mémoire occupée par la matrice suivant
sa dimension et son stockage (full ou sparse) :
d full matrix sparse matrix
103 1 Mo 40 Ko
104 100 Mo 400 Ko
105 10 Go 4 Mo
106 1000 Go 40 Mo
107 100 To 400 Mo

1.1.2 Insertion d'éléments

Regardons les opérations à eectuer sur les tableaux aa, ia et ja si l'on modie la matrice A par la commande

A (1,2)=8;
La matrice devient alors ¨
8.
˛
1. 0. 6.
A “ ˝0. 5. 0. 4.‚.
0. 1. 2. 0.

3
Dans cette opération un élément nul de A a été remplacé par une valeur non nulle 8 : il va donc falloir la stocker
dans les tableaux sachant qu'aucune place n'est prévue. On suppose que les tableaux sont susamment grands
(pas de souci mémoire), il faut alors décaler d'une case l'ensemble des valeurs des tableaux aa et ia à partir de
la 3ème position puis copier la valeur 8 en aap2q et le numéro de ligne 1 en iap2q :
aa 1. 8. 5. 1. 2. 6. 4.

ia 1 1 2 3 3 1 2
Pour le tableau ja, il faut à partir du numéro de colonne 2 plus un, incrémenter de `1 :
ja 1 2 5 6 8
La répétition de ces opérations peut devenir très coûteuse en temps CPU lors de l'assemblage de matrices
utilisant des méthodes similaires à celle utilisée en Listing 1 (et on ne parle pas ici des problèmes de réallocation
dynamique qui peuvent arriver !).

1.1.3 Construction ecace de matrices creuses

Pour générer ecacement une matrice creuse, il faut donc éviter d'insérer de manière répétitive des éléments
dans celle-ci. Pour celà, on va utiliser l'appel suivant de la fonction sparse :
M = sparse(I,J,K,m,n);
Cette commande génère une matrice creuse m par n telle que M(I(k), J(k)) = K(k) . Les vecteurs I , J et
K ont la même longueur. Il faut noter que tous les éléments nuls de K sont ignorés et que tous les éléments
de K ayant les mêmes indices dans I et J sont sommés. Cette sommation est très utile lors de l'assemblage
de matrices obtenues par une méthode de type éléments nis.
Par exemple, voici deux commandes permettant de générer la matrice A précédente :
A = sparse([1 1 2 3 3 1 2],[1 2 2 2 3 4 4],[1 8 5 1 2 6 4],3,4) ;
A = sparse([1 2 1 1 2 3 3],[4 4 1 2 2 2 3],[6 4 1 8 5 1 2],3,4) ;
Nous allons maintenant appliquer ceci dans un cadre un peu général : les matrices tridiagonales.

1.2 Matrices diagonales et tridiagonales

Q. 1 (Matlab) Ecrire la fonction spMatDiag permettant à partir d'un vecteur v P Rd de retourner la matrice
diagonale (creuse) D P Md pRq de diagonale v (i.e. Dpi, iq “ v piq, @i P v1, dw) en créant (sans boucle) les tableaux
I, J et K puis en générant la matrice D à l'aide de la commande D = sparse(I,J,K,d,d); ‚
Soit A P Md pRq la matrice tridiagonale dénit par
¨ ˛
v 1 w1 0 ¨¨¨ ¨¨¨ 0
. ..
w2 . . .
˚ ‹
˚u1 v2 ‹
..
˚ 0 ... .. .. ..
˚ ‹
. . . .
˚ ‹
(1.2)
‹
A“˚ ˚. . .. .. ..
‹
˚ .. .. . . .
‹
0 ‹
˚.
˚ ‹
˝ .. .. .. ‹
. . vd´1 wd´1 ‚
0 ¨¨¨ ¨¨¨ 0 ud´1 vd
avec u P Rd´1 , v P Rd et w P Rd´1 .
Q. 2 (Matlab) 1. Ecrire la fonction spMatTriDiagSca permettant à partir des vecteurs u , v et w de
retourner la matrice creuse A dénie en (1.2). Pour celà on va tout d'abord créer une matrice creuse de
Md pRq puis on va la remplir à l'aide d'une boucle for sur le principe du Listing 1.
2. Ecrire la fonction spMatTriDiagVec permettant à partir des vecteurs u , v et w de retourner la matrice
creuse A en créant (sans boucle) les tableaux I, J et K puis en générant la matrice A à l'aide de la
commande A = sparse(I,J,K,d,d);
3. Ecrire un programme permettant de mesurer (tic-toc) puis de représenter graphiquement les performances
de spMatTriDiagSca et spMatTriDiagVec en fonction de d.
(Attention à la capacité mémoire de la machine !) ‚

4
1.3 Produit de Kronecker de deux matrices

Pour générer ecacement et simplement les matrices creuses obtenus par des méthodes de diérences nies
nous utiliserons le produit de Kronecker de deux matrices.
Soient A P Mm,n pRq et B P Mp,q pRq. Le produit tensoriel de Kronecker de A par B, noté A b B, est la
matrice de Mmp,nq pRq dénie avec des blocs de dimension p ˆ q par
¨
A
1,1 B ¨¨¨ 1,n B
˛
A
A b B “ ˝ ... ..
.
.. ‹ (1.3)
˚
. ‚
Am,1 B ¨¨¨ Am,n B
Le produit de Kronecker est bilinéaire et associatif : Si les dimensions des matrices A, B et C sont compatibles
on a @λ P K
A b pB ` λ ¨ Cq “ pA b Bq ` λpA b Cq
pA ` λ ¨ Bq b C “ pA b Cq ` λpB b Cq
A b pB b Cq “ pA b Bq b C
Par contre, il n'est pas commutatif.
Nous allons écrire une fonction permettant d'eectuer ce produit. Mais auparavant, quelques petits rappels
(ou non) des possibilités oertes par le langage Matlab.
Il est facile sous Malab/Octave de modier une matrice (creuse ou non). Par exemple,
‚ A(I , J)=M
Les tableaux d'indices I et J sont de dimension respective n et m. La matrice M est de dimension
n-par-m. La matrice A est modiée de telle sorte que

A(I(i), J(j))=M(i,j), @i P v1, nw, @j P v1, mw

‚ A(I , J)=A(I,J)+M
Même chose en sommant.
‚ A(i ,:) =a
Ici a est un scalaire. Tous les éléments de la ligne i de la matrice A valent a .
‚ A (:, j)=a
Ici a est un scalaire. Tous les éléments de la colonne j de la matrice A valent a .
‚ ...
Sous Matlab/Octave, il est possible de récupérer tous les éléments non nuls d'une matrice creuse A ainsi
que leurs indices de ligne et de colonne à l'aide de la commande

[ I , J , K]=nd(A);

Ici les trois tableaux I, J et K ont même longueur et on a A(I(k), J(k)) == K(k) pour tout k inférieur ou égal à
length(K).
Q. 3 (Matlab) Ecrire la fonction spMatKron permettant à partir deux matrices creuses carrées A P Mn pRq
et B P Mm pRq de retourner la matrice creuse C “ A b B. Pour celà, on initialisera la matrice retournée C à
l'aide de la commande C=sparse(n*m,n*m) puis on utilisera uniquement une boucle sur les éléments non-nuls
de la matrice A pour remplir la matrice C . ‚

On pourra comparer cette fonction à la fonction kron de Matlab/Octave.

2 Notations et approximation des dérivées secondes

Soient Ω “sa, brˆsc, drĂ R2 et Γ “ BΩ la frontière du domaine Ω. On note ΓN , ΓS , ΓO et ΓE respectivement
les frontières nord, sud, ouest et est. on a

Γ “ ΓN Y ΓS Y ΓO Y ΓE .

Soit v : Ω ÝÑ R susament régulière.

5
Q. 4 (sur feuille) Montrer que l'on a
B2 v vpx ` h, yq ´ 2vpx, yq ` vpx ´ h, yq
px, yq “ ` Oph2 q (2.1)
Bx2 h2
B2 v vpx, y ` hq ´ 2vpx, yq ` vpx, y ´ hq
px, yq “ ` Oph2 q. (2.2)
By 2 h2
‚
N
i“0 et pyj qj“0 les discrétisation régulières, respectivement, des intervalles ra, bs et rc, ds dénes par
On note pxi qN x y

xi “ a ` ihx , @i P v0, Nx w et yj “ c ` jhy , @j P v0, Ny w (2.3)

avec hx “ pb ´ aq{Nx et hy “ pd ´ cq{Ny . On note aussi

nx “ Nx ` 1, ny “ Ny ` 1 et N “ nx ˆ ny (2.4)

Q. 5 (sur feuille) Montrer que l'on a @i Pw0, Nx v, @j Pw0, Ny v

def B2 v B2 v
∆vpxi , yj q “p ` qpxi , yj q
Bx2 By 2
vpxi`1 , yj q ´ 2vpxi , yj q ` vpxi´1 , yj q
“
h2x
vpxi , yj`1 q ´ 2vpxi , yj q ` vpxi , yj´1 q
` ` Oph2x q ` Oph2y q (2.5)
h2y
‚

3 Equation de Poisson avec conditions de Dirichlet

Soient f : Ω ÝÑ R et g : Γ ÝÑ R deux fonctions données. On veut résoudre le problème suivant

´∆u “ f, dans Ω (3.1)

u “ g, sur Γ (3.2)
N
On utilise par la suite les discrétisations pxi qN
i“0 et pyj qj“0 dénies en section 2.
x y

Q. 6 (sur feuille) 1. Montrer que ce problème peut s'écrire après discrétisation sous la forme
Ui`1,j ´ 2Ui,j ` Ui´1,j Ui,j`1 ´ 2Ui,j ` Ui,j´1
´ 2
´ “ f pxi , yj q, @pi, jq Pw0, Nx vˆw0, Ny v, (3.3)
hx h2y
U0,j “ gpa, yj q, @j P v0, Ny w, (3.4)
UNx ,j “ gpb, yj q, @j P v0, Ny w, (3.5)
Ui,0 “ gpxi , cq, @i P v0, Nx w, (3.6)
Ui,Ny “ gpxi , dq, @i P v0, Nx w. (3.7)

avec Ui,j « upxi , yj q.

2. Les inconnues du problème discrétisé sont les Ui,j , pour i P v0, Nx w et j P v0, Ny w. Compter le nombre
d'équations distinctes du problème discrétisé (bord compris). Que peut-on en conclure ? ‚

En notant βx “ ´ h12 , βy “ ´ h12 et µ “ 2p h12 ` 1

h2y q, les équations (3.3) peuvent aussi s'écrire sous la forme
x y x

βx pUi`1,j ` Ui´1,j q ` µUi,j ` βy pUi,j`1 ` Ui,j´1 q “ f pxi , yj q, @pi, jq Pw0, Nx vˆw0, Ny v (3.8)

Pour tout j P v0, Ny w, on note U:,j le vecteur de Rnx dénit par

¨ ˛
U0,j
U:,j “ ˝ ... ‚.
˚ ‹

UNx ,j

6
On note V P RN le vecteur bloc ¨ ˛
U:,0
˚
˚ U:,1 ‹
‹
V “˚ .. ‹
.
˚ ‹
˝ ‚
U:,Ny

Q. 7 (sur feuille) Explicitez la bijection F : v0, Nx w ˆ v0, Ny w ÝÑ v1, N w telle que

@pi, jq P v0, Nx w ˆ v0, Ny w, Vk “ Ui,j , avec k “ Fpi, jq.

Dans le cas de la numérotation en pi, jq P v0, Nx w ˆ v0, Ny w on parlera de numérotation 2D et pour la

numérotation en k P v1, N w on parlera de numérotation globale.

Q. 8 (Matlab) 1. Ecrire la fonction k=bijF(i,j,nx) correspondant à la bijection F ( numerotation

2D vers numerotation globale).
2. Ecrire la fonction réciproque [ i , j]=bijRecF(k,nx) correspondant à F -1 ( numerotation globale vers
numerotation 2D). ‚

Chacune des équations du problème discret (3.3)-(3.7) correspond à une discrétisation en un point pxi , yj q.
Nous choisissons d'écrire ces équations en utilisant la même numérotation que lors de la construction du vecteur
V : l'équation écrite au point pxi , yj q sera écrite en ligne k “ Fpi, jq du système.
Q. 9 (sur feuille) Etablir que le problème discret (3.3)-(3.7) peut s'écrire sous la forme du système linéaire
bloc ¨ ˛
E O ¨¨¨ ¨¨¨ ¨¨¨ O O
˚ M D M O ¨¨¨ O O ‹ ¨
B:,0
˛
.. .. . . .. .. ‹
˚ ‹
˚ B:,1 ‹
. . . . . ‹
˚ ‹
˚ O M
.
˚ ˚ ‹
.. .. .. .. .. ‹V “ ˚ .. ‹
˚ ‹
(3.9)
˚ ‹
. . . . O . ‹
˚ ˚
O ‹
.
˚
.. .. .. ‚
˚ ‹
.. .. ..
˚ ‹ ˚ ‹
. . . . . M O ‹
˚ ‹ ˝
˚
B:,Ny
˚ ‹
˝ O O ¨¨¨ O M D M ‚
O O ¨¨¨ ¨¨¨ ¨¨¨ O E
où chaque bloc de la matrice est une matrice pNx ` 1q par pN x ` 1q. La matrice O P MNx `1 pRq est la matrice
nulle. Les matrices creuses D, M et E ainsi que les vecteurs B:,j P RN x`1 , pour tout j P v0, Ny w, devront être
donnés explicitement. On notera B le vecteur second membre de RN et A P MN pRq la matrice du système. ‚

3.1 Assemblage de la matrice sans conditions aux limites

On note In la matrice identité de Mn pRq et Jn la matrice de Mn pRq dénie par

¨ ˛
0 0 ¨¨¨ ¨¨¨ 0
.. .. ‹
. .‹
˚
˚0 1
Jn “ ˚ ... . . . .. .. . ‹ P M pRq.
˚ ‹
. . .. ‹
˚
n
˚. ..
˚ ‹
˝ .. . 1 0‚
‹

0 ... ... 0 0

Dans cette partie nous allons générer/assembler la matrice du système (3.9) sans tenir compte des
conditions aux limites (i.e. les lignes du système correspondant à des points du bord sont nulles ... pour
l'instant).

7
 Dans ce cas, on note ´Axy la matrice ainsi obtenue. C'est une matrice bloc de MN pRq avec ny lignes bloc
composées de blocs carrés de dimension nx qui s'écrit sous la forme :
¨ ˛ ¨ ˛
O O ¨¨¨ ¨¨¨ ¨¨¨ O O O O ¨¨¨ ¨¨¨ ¨¨¨ O O
˚ O Tx O O ¨¨¨ O O ‹ ˚ S Ty Sy O ¨ ¨ ¨ O O ‹
‹ ˚ y
.. .. ‹ ˚ .. .. ‹
˚ ‹
˚ .. .. ‹ ˚ .. ..
˚ O O T
˚ x . . . . ‹ ˚ O Sy Ty . . . . ‹‹
˚ . . . .
Axy “ ˚ .. O . .. . .. . . . O .. ‹ ` ˚ .. . .. . .. . .. O .. ‹ (3.10)
‹ ˚ ‹
˚ ‹ ˚
O ‹
˚ .. .. . . .. ‹ ˚ .. .. . . ..
˚ ‹ ˚ ‹
˚ . . . . Tx O O ‹ ‹ ˚ . . . . Ty Sy O ‹
‹
˚
˚ ‹
˝ O O ¨¨¨ O O Tx O ‚ ˝ O O ¨ ¨ ¨ O Sy Ty Sy ‚
O O ¨¨¨ ¨¨¨ ¨¨¨ O O O O ¨¨¨ ¨¨¨ ¨¨¨ O O

où ¨ ˛
0 0 0 ¨¨¨ ¨¨¨ 0
.. .. ‹
. .‹
˚
˚1 ´2 1
.. .. .. .. .. ‹
˚ ‹
. . . . .‹
˚
1 2 1 ˚0
Sy “ 2 Jnx , Ty “ ´ 2 Jnx et Tx “ 2 ˚
˚ ... .. .. .. ..
‹
hy hy hx ˚ . . . .
‹
0‹
˚.
˚ ‹
˝ .. .. ‹
. 1 ´2 1‚
0 ... ... 0 0 0
On peut noter que les matrices Tx , Ty et Sy sont des matrices de Mnx pRq.
En notant Ax P Mnx pRq et Ay P Mny pRq les matrices dénies par
¨ ˛ ¨ ˛
0 0 0 ¨¨¨ ¨¨¨ 0 0 0 0 ¨¨¨ ¨¨¨ 0
.. .. ‹ .. .. ‹
. .‹ . .‹
˚ ˚
˚1 ´2 1 ˚1 ´2 1
.. .. .. .. .. .. .. .. .. .. ‹
˚ ‹ ˚ ‹
. . . . .‹ . . . . .‹
˚ ‹ ˚
1 ˚0 1 ˚0
‹ et Ay “ h2 ˚ .
‹
Ax “ 2 ˚
˚ ... .. .. .. .. .. .. .. ..
˚ ‹
hx ˚ . . . . ˚ .. . . . .
‹
0‹ ‹ y 0‹
˚. ˚.
˚ ˚ ‹
˝ .. .. ‹
˝ .. .. ‹
. 1 ´2 1 ‚ . 1 ´2 1‚
0 ... ... 0 0 0 0 ... ... 0 0 0

on déduit de (3.10)
Axy “ Iny b Ax ` Ay b Inx . (3.11)

Q. 10 (Matlab) 1. Ecrire la fonction Lap2DAssembling retournant la matrice bloc creuse A en utilisant

(3.11) et les fonctions déjà écrites.
2. Proposer au moins un programme permettant de tester/valider la matrice ainsi obtenue. ‚

Sur le même principe, il est possible de généraliser en dimension quelconque la construction de la matrice
du Laplacien obtenue par diérences nies sans prise en compte des conditions aux limites. Par exemple en
dimension 3, cette matrice s'écrit
` ˘
Axyz “ Jnz b Axy ` Az b pJny b Jnx q
“ Jnz b Jny b Ax ` Jnz b Ay b Jnx ` Az b Jny b Jnx

3.2 Quelques indices

3.2.1 Indices des bords

La frontière Γ du domaine peut se décomposer en 4 sous-ensembles :

Γ “ ΓN Y ΓS Y ΓE Y ΓO .

Q. 11 (sur feuille) 1. Pour ΓS , donner le sous-ensemble DS d'indices pi, jq dans v0, Nx w ˆ v0, Ny w tel que
pxpiq, ypjqq P ΓS .

8
 2. En déduire le sous-ensemble IS (ordonné suivant DS ) correspondant dans la numérotation globale.
3. Même chose pour ΓN , ΓE et ΓO . ‚

Q. 12 (matlab) Ecrire la fonction BordDomaine retournant la structure de données (help struct) contenant
les champs suivants :
‚ Sud ayant pour valeur IS ,
‚ Nord ayant pour valeur IN ,
‚ Est ayant pour valeur IE ,
‚ Ouest ayant pour valeur IO .
Chacun des tableaux Ix sera stocké sous forme ligne. ‚

Une fois cette fonction écrite, il est aisé en numérotation globale de récupérer l'ensemble idxBD des indices des
points du bords ainsi que l'ensemble idxBDc des indices des points strictement intérieurs (complémentaire du
précédant) :

BD=BordDomaine ( . . . ) ; % . . . ´> a remplacer

idxBD=unique ( [ BD . Sud , BD . Nord , BD . Est , BD . Ouest ] ) ;
idxBDc=setdiff ( 1 : N , idxBD ) ; % N=nx * ny

3.2.2 Evaluation de fonctions sur la grille

Soit F P RN le vecteur (bloc) dénit par
¨ ˛
F:,0 ¨ ˛
˚ F ‹ f px0 , yj q
˚ :,1 ‹ ..
F “˚ ˚ .. ‹ avec F:,j “ ˝ ‚ P RN x`1 , @j P v0, Ny w.
˚ ‹
.
‹
˝ . ‚
f pxNx , yj q
F:,Ny

En Algorithme 1, est présenté une version non vectorisée de la construction du vecteur F . Dans cet algorithme
l'ordre des boucles est primordial pour respecter le choix de la numérotation globale.

Algorithme 1 Algorithme non vectorisé de calcul du vecteur F

kÐ1
for j Ð 1 to ny do
for i Ð 1 to nx do
xpiq, y pjqqq
F pkq Ð f px
k Ðk`1
end for
end for
Q. 13 (Matlab) Ecrire la fonction EvalFun2DSca permettant à partir d'une fonction f donnée et des dis-
crétisations x et y de retourner le vecteur F associé. On pourra utiliser des boucles for . ‚

Pour le cas où la fonction f est vectorisée sous Matlab/Octave, il est alors possible de calculer le vecteur F
sans boucle. Pour celà, nous allons construire les deux vecteurs blocs X et Y de Rn dénit par
¨ ˛
X:,0 ¨ ˛
˚ X ‹ x0
:,1 ‹
. ‹
‹ avec X:,j “ ˚˝ .. ‚ P RN x`1 , @j P v0, Ny w
˚
X “˚ ..
.
˚ ‹
xNx
˝ ‚
X:,Ny

9
et ¨ ˛
Y:,0 ¨ ˛
˚ Y:,1 ‹ yj
.‹
˝ .. ‚ P RN x`1 , @j P v0, Ny w.
‹ avec Y:,j “ ˚
˚ ‹
Y “˚ ..
.
˚ ‹
yj
˝ ‚
Y:,Ny
Avec ces deux vecteurs, le calcul du vecteur F pourra être vectorisé (sous la condition que f le soit) :

F=f(X,Y);

Il nous faut donc écrire une version vectorisée du calcul des vecteurs X et Y .

Algorithme 2 Calcul de F avec construction de Algorithme 3 Calcul de F avec construction de

X et Y (version 1) X et Y (version 2)
kÐ1 1: k Ð 1

for j Ð 1 to ny do 2: for j Ð 1 to ny do

for i Ð 1 to nx do 3: for i Ð 1 to nx do
X pkq Ð x piq 4: X pkq Ð x piq
Y pkq Ð y pjq 5: Y pkq Ð y pjq
X pkq, Y pkqqq
F pkq Ð f pX 6: k Ðk`1
k Ðk`1 7: end for
end for 8: end for
end for 9: X,Y q
F Ð f pX

On peut noter que l'ordre des boucles dans l'Algorithme 3 permet d'armer que la valeur de k en ligne 4 est
donnée par Fpi ´ 1, j ´ 1q.

Algorithme 4 Calcul de F avec construction de

X et Y (version 3)
for j Ð 1 to ny do
for i Ð 1 to nx do
k Ð bijFpi ´ 1, j ´ 1, nx q
X pkq Ð x piq
Y pkq Ð y pjq
end for
end for
X,Y q
F Ð f pX

Grâce à l'application réciproque F -1 il est alors possible de transformer la double boucle en j et i en une seule
boucle sur k : c'est l'objet de l'Algorithme 5. On en déduit alors l'Algorithme 6 totalement vectorisé.

Algorithme 5 Calcul de F avec construction de Algorithme 6 Calcul de F avec construction de

X et Y (version 4) X et Y (version 5 : vectorisée)
for k Ð 1 to N do rI, Js Ð bijRecFp1 : N, nx q
ri, js Ð bijRecFpk, nx q X Ð x pI ` 1q
X pkq Ð x pi ` 1q Y Ð y pJ ` 1q
Y pkq Ð y pj ` 1q F Ð f pXX,Y q
end for
X,Y q
F Ð f pX

Un programme Matlab/Octave complet basé sur ce dernier algorithme est proposé en Listing 3.

10
 Listing 2  Calcul et représentation d'une fonction sur la grille 2D
clear a l l
close a l l
a=´pi ; b=pi ; c =0; d=pi ;
f=@ ( x , y ) 4 * cos ( x + y ) . * sin ( x ´ y) ;
Nx =133; Ny =222;
x=linspace ( a , b , Nx +1) ;
y=linspace ( c , d , Ny +1) ;
nx=Nx +1; ny=Ny +1; N=nx * ny ;

[ I , J ]= bijRecF ( 1 : N , nx ) ;
X=x ( I +1) ; Y=y ( J +1) ; % A t t e n t i o n X e t Y v e c t e u r s l i g n e
F=f ( X , Y ) ; % e t donc F a u s s i !

surf ( x , y , reshape ( F , nx , ny ) ' )

shading interp

Il faut noter qu'il existe sous Matlab/Octave deux fonctions meshgrid et ndgrid qui peuvent être utilées
pour le calcul des vecteurs X et Y . Ces deux fonctions n'utilisent pas la fonction bijRecF .

Q. 14 (Matlab) Ecrire la fonction EvalFun2DVec permettant à partir d'une fonction f et des discrétisations
x et y de retourner le vecteur F associé sans utiliser de boucle for . ‚

3.2.3 Evaluation de fonctions sur une partie de la grille

Soit g : ΓN ÝÑ R. On suppose cette fonction vectorisée sous Matlab/Octave. Pour évaluer cette fonction
aux points de la grille 2D appartenant à ΓN , on utilise la structure de données retournée par la fonction
BordDomaine et plus particulièrement le champ Nord de cette structure. Voici une manière de créer un
vecteur G valant G(k)=g(x(i),y(j)) sur un point du bord Nord avec k “ Fpi, jq et 0 sinon.

Listing 3  Calcul d'une fonction sur le bord Nord de la grille 2D

BD=BordDomaine ( . . . ) ; % . . . ´> a remplacer
[ I , J ]= bijRecF ( 1 : N , nx ) ;
X=x ( I +1) ; Y=y ( J +1) ; % coordonnees des p o i n t s de l a g r i l l e
% en numerotation g l o b a l e
g=@ ( x , y ) sin ( x+y ) . * x+y . ^ 2 ;
Gloc=g ( X ( BD . Nord ) , Y ( BD . Nord ) ) ;
G=zeros ( N , 1 ) ;
G ( BD . Nord )=Gloc ;

Pour évaluer une fonction sur une autre partie de la grille, il sut de récupérer l'ensemble des indices dans la
numérotation globale des points de appartenant à cette partie et de remplacer [Link] par cet ensemble.

3.3 Assemblage du second membre sans conditions aux limites

Nous allons générer/assembler le vecteur second membre du système (3.9) sans tenir compte des condi-
tions aux limites. Avec le choix de la numérotation globale, ce vecteur B de RN est déni par
f pxpiq, ypjqq si pxpiq, ypjqq R Γ et k “ Fpi, jq
"
@k P v1, N w, Bk “
0 sinon
En utilisant les codes donnés (et expliqués) en section 3.2, on a immédiatement un code vectorisé permettant
d'initialiser ce vecteur : il est donné en Listing 4.

11
 Listing 4  Calcul du vecteur second membre sans conditions aux limites
BD=BordDomaine ( . . . ) ; % . . . ´> a remplacer
[ I , J ]= bijRecF ( 1 : N , nx ) ;
X=x ( I +1) ; Y=y ( J +1) ; % coordonnees des p o i n t s de l a g r i l l e
% en numerotation g l o b a l e
f=@ ( x , y ) sin ( x+y ) . * cos ( x´y ) ;
idxBD=unique ( [ BD . Sud , BD . Nord , BD . Est , BD . Ouest ] ) ;
idxBDc=setdiff ( 1 : N , idxBD ) ; % N=nx * ny

B=zeros ( N , 1 ) ;
B ( idxBDc )=f ( X ( idxBDc ) , Y ( idxBDc ) ) ;

3.4 Prise en compte des conditions aux limites

Sans tenir compte des conditions aux limites, nous avons décrit le calcul de la matrice Axy P MN pRq
(section 3.1) et du vecteur B P RN (section 3.3).
En reprenant les résultats et notations de Q.9 (page 7), on en déduit que
A “ ´Axy ` AΓ et B “ B ` B Γ
où AΓ et B Γ sont les contributions des conditions aux limites, c'est à dire que le système
AΓU “ B Γ (3.12)
contient uniquement les équations (3.4) à (3.7). Pour k “ Fpi, jq P v1, N w, si pxpiq, ypjqq P Γ alors la ligne k du
système (3.12) correspond à
U pkq “ gpxpiq, ypjqq.
Toutes les autres lignes du système (3.12) sont nulles. On a donc pour tout k P v1, N w
" t
e k si pxpiq, ypjqq P Γ avec pi, jq “ F -1 pkq
AΓ pk, :q “
0 sinon
où e tk est le k ième vecteur de la base canonique de RN , et
gpxpiq, ypjqq si pxpiq, ypjqq P Γ avec pi, jq “ F -1 pkq
"
B Γ pkq “
0 sinon

En utilisant les codes donnés (et expliqués) en section 3.2, on calcule dans le Listing 5, le vecteur Bcl corres-
pondant à B Γ et la matrice Acl correspondant à AΓ .

Listing 5  Calcul des contributions du bord

BD=BordDomaine ( . . . ) ; % . . . ´> a remplacer
[ I , J ]= bijRecF ( 1 : N , nx ) ;
X=x ( I +1) ; Y=y ( J +1) ; % coordonnees des p o i n t s de l a g r i l l e
% en numerotation g l o b a l e
g=@ ( x , y ) sin ( x+y ) . * cos ( x´y ) ;
idxBD=unique ( [ BD . Sud , BD . Nord , BD . Est , BD . Ouest ] ) ;
idxBDc=setdiff ( 1 : N , idxBD ) ; % N=nx * ny

Bcl=zeros ( N , 1 ) ;
Bcl ( idxBD )=g ( X ( idxBD ) , Y ( idxBD ) ) ;
Acl=sparse ( idxBD , idxBD , ones ( 1 , length ( idxBD ) ) , N , N ) ;

12
3.5 Résolution du système

On doit résoudre le système AU U “ B dénit en Q.9 (page 7). Or on a construit les matrices Axy , AΓ et les
vecteurs B , B Γ de telle sorte que
A “ ´Axy ` AΓ et B “ B ` B Γ

3.5.1 Méthode 1
Dans cette section nous allons directement résoudre le système AU
U “ B.
Q. 15 1. Ecrire le programme EDP2D01 permettant de résoudre numériquement l'EDP (3.1)-(3.2) à l'aide du
schéma discrétisé (3.3) à (3.7). On choisira judicieusement les données permettant sur une même gure
de représenter de la solution numérique et de la solution exacte, et sur une autre gure de représenter
l'erreur entre les deux solutions.
2. Ecrire le programme ordreEDP2D01 permettant de retouver numériquement l'ordre du schéma utilisé. ‚

3.5.2 Méthode 2
Dans cette section nous proposons une autre méthode permettant de résoudre le système AU
U “ B en
éléminant les équations portants sur les points Dirichlet.
Pour celà, on note
ID “ k P v1, N w; pxpiq, ypjqq P Γ avec pi, jq “ F -1 pkq
(

et son complémentaire
c
ID “ v1, N wzID .
Le système linéaire à résoudre AU
U “ B peut se décomposer de manière équivalente en
(3.13a)
"
U “ B pID q
ApID , :qU
U “B ô
AU c
ApID c
U “ B pID
, :qU q (3.13b)
Nous allons réécrire les deux systèmes (3.13a) et (3.13b).
‚ Pour (3.13a), nous avons par construction on a ApID , :q “ AΓ pID , :q et donc

U “ AΓ pID , :qU
ApID , :qU U
c c
“ AΓ pID , ID qU
U pID q ` AΓ pID , ID U pID
qU q

Comme AΓ pID , ID
c
q “ O et AΓ pID , ID q “ I on obtient

U “ U pID q.
ApID , :qU

Donc (3.13a) est équivalent à

U pID q “ B pID q. (3.14)
‚ Pour (3.13b), nous avons par construction ApID
c c
, :q “ ´Ax,y pID , :q et donc
c c
ApID U “ ´Ax,y pID
, :qU U
, :qU
c c c c
“ ´Ax,y pID , ID qU
U pID q ´ Ax,y pID , ID U pID
qU q

En utilisant (3.14), le système (3.13b) s'écrit

c
´ Ax,y pID c
, ID c
U pID
qU c
q “ B pID c
q ` Ax,y pID , ID qB
B pID q. (3.15)

Résoudre le système linéaire AU U “ B est alors équivalent à calculer U pID q par (3.14) et déterminer U pID
c
q
en résolvant le système linéaire (3.15).
Q. 16 1. Ecrire le programme EDP2D02 permettant de résoudre numériquement l'EDP (3.1)-(3.2) à l'aide
du schéma discrétisé (3.3) à (3.7) que l'on résoudra par (3.14) et (3.15). On choisira judicieusement les
données permettant sur une même gure de représenter de la solution numérique et de la solution exacte,
et sur une autre gure de représenter l'erreur entre les deux solutions.
2. Ecrire le programme ordreEDP2D02 permettant de retouver numériquement l'ordre du schéma utilisé. ‚

13