Cours et exercices corrigés de L3 : Analyse numérique matricielle : Conditionnement, factorisation LU/Cholesky et gradient

1. Systèmes linéaires $Ax = b$ : conditionnement et méthodes directes

• Conditionnement : pour $A \in \text{GL}_n(\mathbb{R})$ et une norme subordonnée, $\kappa(A) = \|A\| \, \|A^{-1}\| \ge 1$. Si $Ax = b$ et $A(x + \delta x) = b + \delta b$ :

$$\frac{\|\delta x\|}{\|x\|} \le \kappa(A) \frac{\|\delta b\|}{\|b\|}$$
Pour la norme 2 et $A$ symétrique définie positive, $\kappa_2(A) = \frac{\lambda_{\max}}{\lambda_{\min}}$.

• Décomposition LU : pour $A$ inversible, il existe une unique écriture $A = LU$ ($L$ triangulaire inférieure à diagonale unité, $U$ triangulaire supérieure) si et seulement si tous les mineurs principaux dominants de $A$ sont non nuls. Sinon, on pivote : $PA = LU$ avec $P$ matrice de permutation (existe toujours). Coût : environ $\frac{2}{3}n^3$ opérations.
• Cholesky : si $A$ est symétrique définie positive, il existe une unique $L$ triangulaire inférieure à diagonale strictement positive telle que $A = L \, {}^tL$ (coût environ $\frac{1}{3}n^3$).
• QR : toute $A$ inversible s'écrit $A = QR$ avec $Q$ orthogonale et $R$ triangulaire supérieure (Gram-Schmidt, ou réflexions de Householder, numériquement plus stables).

2. Méthodes itératives

• Principe : on écrit $A = M - N$ avec $M$ facile à inverser et l'on itère $x_{k+1} = M^{-1}N x_k + M^{-1}b$. L'erreur $e_k = x_k - x$ vérifie $e_{k+1} = B e_k$ avec $B = M^{-1}N$.
• Critère : la méthode converge pour tout $x_0$ si et seulement si le rayon spectral $\rho(B) < 1$. (Si $\rho(B) < 1$, il existe une norme subordonnée telle que $\|B\| < 1$ ; si $\rho(B) \ge 1$, un vecteur propre associé fournit une erreur qui ne tend pas vers 0.)
• Jacobi : $M = D$ (diagonale de $A$). Gauss-Seidel : $M = D - E$ (partie triangulaire inférieure de $A$, avec $A = D - E - F$).
• Conditions suffisantes : si $A$ est à diagonale strictement dominante ($|a_{ii}| > \sum_{j \neq i} |a_{ij}|$), Jacobi et Gauss-Seidel convergent ; si $A$ est symétrique définie positive, Gauss-Seidel converge.
• Vitesse : asymptotiquement, l'erreur est multipliée par environ $\rho(B)$ à chaque itération.

3. Interpolation de Lagrange

• Existence et unicité : pour $n+1$ abscisses distinctes $x_0, \dots, x_n$ et des valeurs $y_0, \dots, y_n$, il existe un unique $P \in \mathbb{R}_n[X]$ tel que $P(x_i) = y_i$ pour tout $i$ :

$$P = \sum_{i=0}^n y_i L_i, \qquad L_i(x) = \prod_{j \neq i}\frac{x - x_j}{x_i - x_j}$$
(Unicité : la différence de deux solutions est de degré $\le n$ et a $n+1$ racines.)

• Forme de Newton : $P(x) = f[x_0] + f[x_0, x_1](x - x_0) + \dots + f[x_0, \dots, x_n]\prod_{i<n}(x - x_i)$ avec les différences divisées $f[x_i, \dots, x_{i+k}] = \frac{f[x_{i+1}, \dots, x_{i+k}] - f[x_i, \dots, x_{i+k-1}]}{x_{i+k} - x_i}$ ; ajouter un point ne demande qu'un terme de plus.
• Erreur : si $f \in \mathcal{C}^{n+1}([a ; b])$ et les $x_i \in [a ; b]$, pour tout $x \in [a ; b]$, il existe $\xi_x$ dans le plus petit intervalle contenant $x, x_0, \dots, x_n$ tel que

$$f(x) - P(x) = \frac{f^{(n+1)}(\xi_x)}{(n+1)!}\prod_{i=0}^n (x - x_i)$$

• Runge et Tchebychev : avec des nœuds équirépartis, l'erreur peut tendre vers l'infini quand $n$ croît (exemple $f(x) = \frac{1}{1 + 25x^2}$ sur $[-1 ; 1]$). Les nœuds de Tchebychev $x_k = \cos\left(\frac{(2k+1)\pi}{2n+2}\right)$ minimisent $\max_{[-1;1]}\left|\prod (x - x_i)\right|$, qui vaut alors $2^{-n}$.

4. Démonstration : formule de l'erreur d'interpolation

Théorème. Soit $f \in \mathcal{C}^{n+1}([a ; b])$, $x_0, \dots, x_n$ distincts dans $[a ; b]$ et $P$ le polynôme d'interpolation de $f$ en ces points. Pour tout $x \in [a ; b]$, il existe $\xi_x$ dans le plus petit intervalle $I_x$ contenant $x, x_0, \dots, x_n$ tel que $f(x) - P(x) = \frac{f^{(n+1)}(\xi_x)}{(n+1)!}\pi(x)$, avec $\pi(t) = \prod_{i=0}^n (t - x_i)$.

Démonstration.
1. Si $x$ est l'un des $x_i$, les deux membres sont nuls ; tout $\xi_x$ convient. Supposons désormais $x \notin \{x_0, \dots, x_n\}$, de sorte que $\pi(x) \neq 0$.
2. Fonction auxiliaire. Posons $K = \frac{f(x) - P(x)}{\pi(x)}$ et $\phi(t) = f(t) - P(t) - K\pi(t)$ pour $t \in [a ; b]$. Alors $\phi \in \mathcal{C}^{n+1}$.
3. Zéros. $\phi(x_i) = f(x_i) - P(x_i) - 0 = 0$ pour $0 \le i \le n$, et $\phi(x) = 0$ par choix de $K$ : $\phi$ s'annule en $n+2$ points distincts de $I_x$.
4. Rolle itéré. Entre deux zéros consécutifs de $\phi$, le théorème de Rolle fournit un zéro de $\phi'$ : $\phi'$ a au moins $n+1$ zéros distincts dans $I_x$. En itérant, $\phi^{(j)}$ a au moins $n+2-j$ zéros dans $I_x$ ; donc $\phi^{(n+1)}$ s'annule en un point $\xi_x \in I_x$.
5. Calcul. $P$ est de degré $\le n$, donc $P^{(n+1)} = 0$ ; $\pi$ est unitaire de degré $n+1$, donc $\pi^{(n+1)} = (n+1)!$. Ainsi $0 = \phi^{(n+1)}(\xi_x) = f^{(n+1)}(\xi_x) - K(n+1)!$.
6. D'où $K = \frac{f^{(n+1)}(\xi_x)}{(n+1)!}$, c'est-à-dire la formule annoncée. $\square$

Corollaire. $\|f - P\|_{\infty, [a;b]} \le \frac{\|f^{(n+1)}\|_\infty}{(n+1)!}\|\pi\|_{\infty, [a;b]}$.

5. Intégration numérique

Une formule de quadrature $\int_a^b f \approx \sum w_i f(x_i)$ est de degré d'exactitude $m$ si elle est exacte pour tous les polynômes de degré $\le m$ mais pas pour $X^{m+1}$.

• Trapèzes : $\int_a^b f \approx \frac{b-a}{2}(f(a) + f(b))$ ; degré 1 ; erreur $-\frac{(b-a)^3}{12}f''(\xi)$ pour $f \in \mathcal{C}^2$.
• Simpson : $\int_a^b f \approx \frac{b-a}{6}\left(f(a) + 4f\left(\frac{a+b}{2}\right) + f(b)\right)$ ; degré 3 ; erreur $-\frac{(b-a)^5}{2880}f^{(4)}(\xi)$ pour $f \in \mathcal{C}^4$.
• Méthodes composites : en découpant $[a ; b]$ en $N$ intervalles de longueur $H = \frac{b-a}{N}$, l'erreur devient $-\frac{(b-a)H^2}{12}f''(\xi)$ (trapèzes, ordre 2) et $-\frac{(b-a)H^4}{2880}f^{(4)}(\xi)$ (Simpson, ordre 4).
• Gauss-Legendre : avec $n+1$ nœuds choisis comme racines du polynôme de Legendre de degré $n+1$ (et des poids adaptés), la formule a le degré d'exactitude maximal $2n+1$.
• Newton pour $f(x) = 0$ : $x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$. Si $f \in \mathcal{C}^2$ et $r$ est une racine simple ($f'(r) \neq 0$), la convergence est locale et quadratique : $|x_{k+1} - r| \le C|x_k - r|^2$ pour $x_0$ assez proche de $r$.
• EDO $y' = f(t, y)$, $f$ lipschitzienne en $y$ : Euler explicite $y_{n+1} = y_n + hf(t_n, y_n)$ est d'ordre 1 (erreur globale $O(h)$) ; Runge-Kutta classique (RK4) est d'ordre 4 (erreur globale $O(h^4)$). Les schémas explicites demandent un pas petit sur les problèmes raides.