- la première partie s’intéresse à une méthode de factorisation de matrice symétrique appelée décomposition de Cholesky;
- la seconde partie construit un algorithme appelé « méthode du gradient conjugué » pour résoudre un système linéaire.
- enfin, la dernière partie utilise les algorithmes définis dans les deux premières parties afin de résoudre une équation aux dérivées partielles, l’équation de la chaleur.
Lors des tests, il est possible de comparer les algorithmes développés avec la méthode fournie par Scipy (linalg.solve
), qui permet de comparer les résultats obtenus.
Décomposition de Cholesky



![\left\{\begin{array}{rl} t_{i,i}^2 &= \displaystyle a_{i,i} - \sum_{k=1}^{i-1} t_{i,k}^2 \notag \\[8mm] t_{j,i} &= \displaystyle \frac{a_{i,j}- \sum_{k=1}^{i-1} t_{i,k}t_{j,k}}{t_{i,i}} \qquad j\geq i\notag \end{array}\right.](wp-content/cours/IS104/images/left{begin{array}{rl}t_{ii}^2&=displaystylea_{ii}-sum_{k=1}^{i-1}t_{ik}^2notag[8mm]t_{ji}&=displaystylefrac{a_{ij}-sum_{k=1}^{i-1}t_{ik}t_{jk}}{t_{ii}}qquadjgeqinotagend{array}right..png)
Le résultat de la factorisation est la matrice


- Écrire l’algorithme de factorisation dense de Cholesky et le tester. Quelle est sa complexité ?
- En utilisant cet algorithme, combien coûte une résolution de système linéaire dense
?
La factorisation dense de Cholesky, bien qu’un algorithme privilégié pour la résolution de systèmes linéaires symétriques (du fait de sa stabilité numérique) est encore trop coûteuse lorsque la matrice est creuse (i.e ne contient que peu de valeurs non nulles). Il existe un algorithme de factorisation de Cholesky dans le cas creux, mais il est complexe à implémenter. Alternativement, la factorisation de Cholesky incomplète consiste à ne calculer qu’un sous-ensemble des éléments de la matrice .
Dans ce projet, la factorisation incomplète correspond à la factorisation classique, dans laquelle on ne calcule pas les éléments lorsque
est nul. Ces éléments restent alors à zéro et n’influent pas sur la suite du calcul de la factorisation.
- Écrire un algorithme permettant de générer des matrices symétriques définies positives creuses avec un nombre de termes extra-diagonaux non nuls réglable.
- Écrire l’algorithme de factorisation de Cholesky incomplète et le tester (quels tests réaliser ?).
Quelle est sa complexité ? Penser à choisir de manière appropriée en fonction de quels paramètres établir cette complexité et expliquer en quoi elle est meilleure que celle de l’algorithme précédent.
Un préconditionneur est une matrice facile à inverser telle que
soit inférieur à
, où
est le conditionnement de la matrice
. En général, on choisit
comme un inverse approché de
(
) relativement facile à calculer.
- Évaluer si la matrice
obtenue dans les deux cas précédents est un préconditionneur de bonne ou de mauvaise qualité (penser à utiliser
np.linalg.cond
).
Méthode du gradient conjugué




function [x] = conjgrad(A,b,x) r=b-A*x; p=r; rsold=r'*r; for i=1:10^(6) Ap=A*p; alpha=rsold/(p'*Ap); x=x+alpha*p; r=r-alpha*Ap; rsnew=r'*r; if sqrt(rsnew) < 1e-10 break; end p=r+rsnew/rsold*p; rsold=rsnew; end
Le paramètre sert d’initialisation à l’algorithme, mais peut être choisi comme égal au vecteur nul.
- Par quels aspects cette implémentation ne respecte t’elle pas des standards de codage très sains ?
- Par quel aspect cette méthode est-elle différente de la méthode décrite mathématiquement juste avant ? Quel est le gain de complexité occasionné ? Est-il systématique ?
- Implémenter la méthode du gradient conjugué et la tester.
Il est possible d’utiliser cette méthode avec un préconditionneur, comme évoqué dans le paragraphe précédent. La même page propose un algorithme avec préconditionneur (mais sans pseudo-code).
- Implémenter la méthode du gradient conjugué avec préconditionneur et la tester.
Application à l’équation de la chaleur
![[0;1]\times[0;1]](wp-content/cours/IS104/images/[01]times[01].png)
.png)
.png)
![\left\{\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} = f(x,y)\right. \qquad \forall (x,y) \in [0;1]\times[0;1]](wp-content/cours/IS104/images/left{frac{partial^2T}{partialx^2}+frac{partial^2T}{partialy^2}=f(xy)right.qquadforall(xy)in[01]times[01].png)
Notre but ici consiste à déterminer, étant donnée une fonction


![[0;1]\times[0;1]](wp-content/cours/IS104/images/[01]times[01].png)

(1,1) \end{pspicture}} \raise 9mm \hbox{$\qquad\longrightarrow\qquad$} {\begin{pspicture}(-0.1,-0.6)(2,2.1) \psset{unit=2} \psframe[linecolor=notwhite](-0.01,-0.01)(1,1) \psaxes(0,0)(0,0)(1,1) \multido{\nn=0.05+0.10}{10}{ \multido{\nnnn=0.05+0.10}{10}{\pscircle(\nn,\nnnn){1pt}}} \end{pspicture}} \raise 9mm \hbox{$\qquad [x_i,y_j] = [\frac{i}{N+1};\frac{j}{N+1}]$}}](wp-content/cours/IS104/images/p2_heat_square.png)
Le domaine est alors vu comme un ensemble fini de points organisés sur un carré, et les fonctions et
sont des matrices carrées associant à un point la valeur correspondant à ce point. De plus, les dérivées partielles sont approximées par les équations linéaires suivantes :
_{ij}approxfrac{T(x_i+hy_j)+T(x_i-hy_j)-2T(x_iy_j)}{h^2}=frac{t_{i+1j}+t_{i-1j}-2t_{ij}}{h^2}.png)
_{ij}approxfrac{T(x_iy_j+h)+T(x_iy_j-h)-2T(x_iy_j)}{h^2}=frac{t_{ij+1}+t_{ij-1}-2t_{ij}}{h^2}.png)
où est la distance entre deux points consécutifs sur la même ligne ou la même colonne, c’est-à-dire ici
.
Conditions aux bords de Dirichlet : Pour bien définir le problème précédent, il est nécessaire de se donner des conditions aux limites sur les bords du domaine. Les conditions les plus simples sont les suivantes :
![\left\{\begin{array}{c} T(0,y) = T(1,y) = 0 \\ T(x,0) = T(x,1) = 0\end{array}\right. \qquad \forall (x,y) \in [0;1] \times [0;1]](wp-content/cours/IS104/images/left{begin{array}{c}T(0y)=T(1y)=0T(x0)=T(x1)=0end{array}right.qquadforall(xy)in[01]times[01].png)
- Vérifier que le problème peut se mettre alors sous la forme d’un système linéaire
où la matrice
est une matrice tridiagonale par blocs de taille
et de la forme suivante :
Noter que la solution
est un vecteur de
lignes, qu’il s’agira de retransformer en une matrice
telle que
.
L’entier
compte le nombre de points par lignes et par colonnes dans le domaine. Écrire un algorithme permettant de générer la matrice
et le vecteur
en fonction du nombre
.
- Résoudre le système dans le cas d’un radiateur placé au centre du carré (quel est le vecteur image ?).
Tracer l’image du résultat. - Résoudre le système dans le cas d’un mur chaud placé au nord. du carré (quel est le vecteur image ?).
Tracer l’image du résultat.
Pour exemple, voilà le type d’image que l’on peut obtenir en résolvant certains de ces systèmes:
