M1 Mathématiques Fondamentales et Appliquées (MFA) Université Paris–Sud Grands Systèmes Linéaires Centre d’Orsay Année 2009 − 2010

Examen - Lundi 18 Janvier 2010 - Durée : 3 heures

Les notes de cours ainsi que les scripts de TP sont autorisés.

Les résultats des questions théoriques devront être fournis sur papier et pour les résultats des questions numériques, on rendra un répertoire contenant tous les listings des programmes. Les figures pourront également être imprimées sur papier.

Exercice - 1 Méthode du gradient conjugué projeté

On s’intéresse dans cet exercice à la recherche du minimum dans un espace affine E de la fonction

J (x) = \frac{1}{2}(Ax, x) − (b, x)

où A est une matrice n × n symétrique positive (i.e (A, x, x) ≥ 0 ∀x ∈ ℝ^{n}). Et b est un vecteur de ℝ^{n}. A est supposée définie positive sur l’espace vectorielle associé à E.

On supposera que E est défini par E = {x ∈ ℝ^{n}, Bx + c = 0} où c ∈ R^{m} et B est une matrice m × n telle que B^{T} est de rang m < n.

On désigne par P l’opérateur de projection orthogonal sur l’espace vectoriel associé à E.

Q-1 : Opérateur de projection sur l’espace vectoriel associé à E

Q-1-1 : On suppose que les vecteurs lignes de B sont linéairement indépendants. Montrer que la matrice BB^{T} est alors régulière.

Q-1-2 : On pose Q = B^{T} (BB^{T})^{−1}B.

Montrer que QQ = Q, Q^{T} = Q, Q(I − Q) = (I − Q)Q = 0.

Montrer que Q est une opérateur de projection orthogonal sur Im(B^{T}).

En déduire que P = I − Q.

Q-2 : Construction de l’algorithme

Soit x_{0} un point qui vérifie Bx_{0} + c = 0. On introduit la nouvelle variable y par la formule

x = x_{0} + P y. On pose G(y) = J (x_{0} + P y)

Q-2-1 : Soit ȳ un point de minimum de la fonction G. Montrer que le point correspondant

x̄ = x_{0} + P ȳ réalise le minimum de J sous la contrainte Bx + c = 0.

Q-2-2 : Montrer que l’algorithme du gradient conjugué appliqué à la recherche du minimum G est

Initialisation : y_{0} = 0, r_{0} = P (b − Ax_{0}), d_{0} = r_{0} Itérations : k = 0, 1, …

\begin{aligned}α_{k} & = & ∥r_{k}∥^{2}/(d_{k}, (P AP)d_{k}) \\ y_{k+1} & = & y_{k} + α_{k}d_{k} \\ r_{k+1} & = & r_{k} − α_{k}(P AP)d_{k} \\ β_{k+1} & = & ∥r_{k+1}∥^{2}/∥r_{k}∥^{2} \\ d_{k+1} & = & r_{k+1} + β_{k+1}d_{k} \end{aligned}

Q-3 : Simplification de l’algorithme

On souhaite réécrire l’algorithme précédent de sorte à faire intervenir directement x Q-3-1 : Montrer que P d_{k} = d_{k} (On pourra le faire par récurrence.) Q-3-2 : En déduire que l’algorithme peut s’écrire

Initialisation : x_{0} tel que Bx_{0} + c = 0, r_{0} = b − Ax_{0}, z_{0} = P r_{0}, d_{0} = z_{0} Itérations : k = 0, 1, …

\begin{aligned}α_{k} & = & (z_{k}, z_{k})/(d_{k}, Ad_{k}) \\ x_{k+1} & = & x_{k} + α_{k}d_{k} \\ r_{k+1} & = & r_{k} − α_{k}Ad_{k} \\ z_{k+1} & = & P r_{k+1} \\ β_{k+1} & = & ∥z_{k+1}∥^{2}/∥z_{k}∥^{2} \\ d_{k+1} & = & z_{k+1} + β_{k+1}d_{k} \end{aligned}

Q-4 : Soit m_{0} l’ordre de multiplicité de la valeur propre nulle de P AP. Montrer que l’algorithme précédent converge en au plus n − m_{0} itérations. Q-5 : Implémentation

Programmer l’algorithme à l’aide d’une fonction Matlab dont on précisera les arguments, prévoir de passer la matrice de projection P en paramètre.

Exercice - 2 Applications

On considère un problème du type −εΔu = f, dans Ω, \frac{∂u}{∂n} = 0 sur ∂Ω. On s’intéresse pour simplifier à la version mono-dimensionnelle donnée par

(1)  \{ \begin{matrix}−εu′′(x) & = & f(x) & dans ]0, 1[, \\ u′(1) & = & 0 & , \\ u′(0) & = & 0 \end{matrix}

où f est une fonction continue sur l’intervalle [0, 1] et ε > 0 est un réel. Q-1 : Montrer que si ∫_{0}^{1} f(x) dx ≠ 0 alors le problème (1) n’admettra pas de solution.

Dans la suite on suppose que ∫_{0}^{1} f(x) dx = 0. On pourra par exemple prendre f(x) = \frac{1}{2} − x. Q-2 : Montrer que si une fonction u est solution de (1), il en sera de même de u + c, ∀c ∈ ℝ Q-3 : Montrer que (1) a une unique solution si l’on impose à cette solution d’être à moyenne nulle, c’est-à-dire

(2)  ∫_{0}^{1} u(x) dx = 0.

(On pourra déterminer analytiquement cette solution en résolvant le problème.)

Ainsi, le problème que l’on considère dans la suite est le problème (1)-(2) qu’on choisit de discrétiser par une méthode de différences finies.

On commence par recouvrir [0, 1] d’une grille uniforme de pas h = \frac{1}{n−1}, pour un n > 1 ∈ ℕ. Ceci génère les noeuds

X_{n} = {x_{i}, x_{i} = (i − 1) × h, i ∈ {1, …, n}} .

On écrit ensuite que l’équation est vérifiée en chaque noeud interne, soit

(3)  \{ \begin{matrix}−εu′′(x_{i}) & = & f(x_{i}) & & ∀i ∈ \{2, …, n − 1\}, \\ u′(x_{1}) & = & 0, \\ u′(x_{n}) & = & 0, \\ ∫_{0}^{1} u(x) dx & = & 0 & . \end{matrix}

On remplace ensuite les opérateurs différentiels par des différences divisées et on approche l’intégrale par une formula de quadrature. On choisit la formule des trapèzes pour approcher l’intégrale sur chaque intervalle (i.e ∫_{x_{i}}^{x_{i}+h} u(x) dx = \frac{h}{2}(u(x_{i}) + u(x_{i} + h)) + O(h^{3}) )

Q-3-1 : Montrer que u solution de (3) vérifie :

(4)  \{ \begin{matrix}ε\frac{−u(x_{i−1}) + 2u(x_{i}) − u(x_{i+1})}{h^{2}} + O(h^{2}) & = & f(x_{i}), ∀i = 2, …, n − 1 \\ \frac{u(x_{2}) − u(x_{1})}{h} + O(h) & = & 0, \\ \frac{u(x_{n}) − u(x_{n−1})}{h} + O(h) & = & 0, \\ h (\frac{1}{2}u(x_{1}) + ∑_{i=2}^{n−1} u(x_{i}) + \frac{1}{2}u(x_{n})) + O(h^{2}) & = & 0. \end{matrix}

Pour tout i = 1, …, n, on désigne par u_{i} la valeur approchée de u(x_{i}) lorsqu’on se débarrasse des restes dans les formules (4). De même on définit f_{i} = f(x_{i}), ∀i = 2, …, n − 1 et f_{1} = 0, f_{n} = 0.

Q-3-2 : Montrer qu’on obtient le problème

(5)  \{ \begin{matrix}u_{1} − u_{2} & = & 0, \\ −u_{i−1} + 2u_{i} − u_{i+1} & = & \frac{h^{2}}{ε} f_{i}, ∀i = 2, …, n − 1 \\ −u_{n−1} + u_{n} & = & 0, \\ u_{1} + (∑_{i=2}^{n−1} 2u_{i}) + u_{n} & = & 0. \end{matrix}

Q-3-3 : Montrer que l’équation (5) se met sous la forme \{ \begin{matrix}Ax & = & b, \\ Bx & = & 0, \end{matrix}

où A est une matrice carrée de taille n et B est une matrice 1 × n, x et b sont des vecteurs de taille n, tous donnés par :

A = ( \begin{matrix}1 & −1 & 0 & ⋯ & 0 & 0 \\ −1 & 2 & −1 & 0 & & 0 \\ 0 & −1 & 2 & −1 & ... & ⋮ \\ ⋮ & ... & ... & ... & ... & 0 \\ 0 & & ... & −1 & 2 & −1 \\ 0 & 0 & ⋯ & 0 & −1 & 1 \end{matrix} ) , x = ( u_{1} ; u_{2} ; ⋮ ; u_{n−1} ; u_{n} ) , b = \frac{h^{2}}{ε} ( f_{1} ; f_{2} ; ⋮ ; f_{n−1} ; f_{n} ) , B^{T} = ( 1 2 ⋮ 2 1 )

On rappelle qu’on a posé f_{1} = f_{n} = 0. Vérifier qu’on a effectivement Bb = 0.

Q-4 : Montrer que ce problème rentre dans le formalisme de l’exercice précédent.

Q-5 : Implémentation : On prendra dans les tests f(x) = \frac{1}{2} − x, ε = 1

Q-5-1 : Écrire un script Matlab qui construit les matrices A, B et le vecteur b et qui teste l’algorithme du gradient conjugué projeté de l’exercice précédent.

Q-5-2 : Comparer la solution obtenue avec celle obtenue en résolvant le système linéaire Ax = b par la fonction minres de Matlab .

Exercice - 3 Accélération d’une méthode itérative

Soit à résoudre le système Ax = b. On utilise une méthode itérative qui produit une suite (x_{k}) d’itérés que l’on supposera converger vers x. Pour accélérer la convergence de (x_{k}), on la transforme en une nouvelle suite (y_{k}) donnée par

y_{k} = x_{k} − Z_{k}Λ_{k}, k = 0, 1 …

où Z_{k} est une matrice n × n_{k} de rang n_{k} (n est la dimension du système) et Λ_{k} ∈ ℝ^{n_{k}}. On pose

e_{k} = x − x_{k}
ε_{k} = x − y_{k}
r_{k} = b − Ax_{k} = Ae_{k}
ρ_{k} = b − Ay_{k} = Aε_{k}

Soit G une matrice symétrique définie positive. Q-1 : Quel est le vecteur Λ_{k} qui minimise ∥ε_{k}∥_{G}^{2} = (ε_{k}, Gε_{k}). Q-2 : Montrer que, pour ce vecteur Λ_{k}, on a ∥ε_{k}∥_{G}^{2} = ∥e_{k}∥_{G}^{2} − e_{k}^{T} GZ_{k}(Z_{k}^{T} GZ_{k})^{−1}(GZ_{k})^{T} e_{k}.

Q-3 : Montrer que ∥ε_{k}∥_{G} ≤ ∥e_{k}∥_{G} et que \{ \begin{matrix}ε_{k} & ⊥_{G} & Z_{k} \\ ε_{k} & ⊥_{G} & y_{k} − x_{k} \end{matrix} où u ⊥_{G} v signifie que u^{T} Gv = 0.

Q-4 : Montrer que

∥y_{k} − x_{k}∥_{G} ≤ ∥e_{k}∥_{G}.

Q-5 : Pour le choix G = A^{T} A, montrer que ∥e_{k}∥_{G} = ∥r_{k}∥, ∥ε_{k}∥_{G} = ∥ρ_{k}∥ et que

∥r_{k}∥^{2} − ∥ρ_{k}∥^{2} = r_{k}^{T} AZ_{k} ((AZ_{k})^{T} AZ_{k})^{−1} (AZ_{k})^{T} r_{k}.

Montrer que ρ_{k} = r_{k} − P_{k}r_{k} avec P_{k} = U_{k}(U_{k}^{T} U_{k})^{−1}U_{k}^{T} et U_{k} = AZ_{k}.

Que représente la matrice P_{k}? Q-6 : Pour le choix G = I, on pose Z_{k}^{T} = Y_{k}^{T} A. Montrer que

∥e_{k}∥^{2} − ∥ε_{k}∥^{2} = r_{k}^{T} Y_{k} ((A^{T} Y_{k})^{T} A^{T} Y_{k})^{−1} Y_{k}^{T} r_{k}.

Montrer que ε_{k} = e_{k} − P_{k}e_{k} avec P_{k} = Z_{k}(Z_{k}^{T} Z_{k})^{−1}Z_{k}^{T} .

Que représente la matrice P_{k}? Q-7 : Lorsque la matrice A est symétrique définie positive, on choisit G = A. Montrer que

∥e_{k}∥_{A}^{2} − ∥ε_{k}∥_{A}^{2} = r_{k}^{T} Z_{k} (Z_{k}^{T} AZ_{k})^{−1} Z_{k}^{T} r_{k}.

Montrer que ρ_{k} = r_{k} − P_{k}r_{k} avec P_{k} = AZ_{k}(Z_{k}^{T} AZ_{k})^{−1}Z_{k}^{T} .

Que représente la matrice P_{k}?

On pose A = B^{T} B. Que peut-on dire de la matrice Q_{k} = V_{k}(V_{k}^{T} V_{k})^{−1}V_{k}^{T} avec V_{k} = BZ_{k}? Q-8 : Évaluez ces accélérateurs sous Matlab pour k = 2, 3, 4, sur une méthode de votre choix. La matrice A pourra être une matrice quelconque. Commentez les résultats.