Analyse Numérique I - Cours (Licence Mathématiques - PDF)

Université Aix Marseille Licence de mathématiques Cours d’Analyse numérique Raphaèle Herbin 10 septembre 2015 Table des matières 1 Systèmes linéaires 5 1.1...

Université Aix Marseille Licence de mathématiques Cours d’Analyse numérique Raphaèle Herbin 10 septembre 2015 Table des matières 1 Systèmes linéaires 5 1.1 Objectifs................................................ 5 1.2 Pourquoi et comment ?........................................ 5 1.2.1 Quelques rappels d’algèbre linéaire............................. 5 1.2.2 Discrétisation de l’équation de la chaleur.......................... 11 1.2.3 Exercices........................................... 15 1.2.4 Suggestions pour les exercices................................ 21 1.2.5 Corrigés des exercices.................................... 21 1.3 Les méthodes directes......................................... 28 1.3.1 Définition........................................... 28 1.3.2 Méthode de Gauss, méthode LU............................... 28 1.3.3 Méthode de Choleski..................................... 36 1.3.4 Quelques propriétés..................................... 42 1.3.5 Exercices........................................... 44 1.3.6 Suggestions.......................................... 48 1.3.7 Corrigés............................................ 48 1.4 Normes et conditionnement d’une matrice.............................. 60 1.4.1 Normes, rayon spectral.................................... 61 1.4.2 Le problème des erreurs d’arrondis............................. 66 1.4.3 Conditionnement et majoration de l’erreur d’arrondi.................... 66 1.4.4 Discrétisation d’équations différentielles, conditionnement “efficace"........... 70 1.4.5 Exercices........................................... 71 1.4.6 Suggestions pour les exercices................................ 78 1.4.7 Corrigés............................................ 79 1.5 Méthodes itératives.......................................... 95 1.5.1 Définition et propriétés.................................... 95 1.5.2 Quelques exemples de méthodes itératives.......................... 97 1.5.3 Les méthodes par blocs.................................... 103 1.5.4 Exercices, énoncés...................................... 106 1.5.5 Exercices, suggestions.................................... 113 1.5.6 Exercices, corrigés...................................... 114 1.6 Valeurs propres et vecteurs propres.................................. 133 1.6.1 Méthode de la puissance et de la puissance inverse..................... 133 1.6.2 Méthode QR......................................... 134 1.6.3 Exercices........................................... 135 1.6.4 Suggestions.......................................... 139 1.6.5 Corrigés............................................ 139 1 TABLE DES MATIÈRES TABLE DES MATIÈRES 2 Systèmes non linéaires 144 2.1 Rappels et notations de calcul différentiel.............................. 144 2.2 Les méthodes de point fixe...................................... 147 2.2.1 Point fixe de contraction................................... 147 2.2.2 Point fixe de monotonie................................... 151 2.2.3 Vitesse de convergence.................................... 153 2.2.4 Méthode de Newton dans IR................................. 155 2.2.5 Exercices........................................... 156 2.3 Méthode de Newton dans IR n.................................... 164 2.3.1 Construction et convergence de la méthode......................... 164 2.3.2 Variantes de la méthode de Newton............................. 166 2.3.3 Exercices........................................... 169 3 Optimisation 204 3.1 Définitions et rappels......................................... 204 3.1.1 Extrema, points critiques et points selle............................ 204 3.1.2 Convexité........................................... 206 3.1.3 Exercices........................................... 208 3.2 Optimisation sans contrainte..................................... 210 3.2.1 Définition et condition d’optimalité............................. 210 3.2.2 Résultats d’existence et d’unicité............................... 211 3.2.3 Exercices........................................... 214 3.3 Algorithmes d’optimisation sans contrainte............................. 220 3.3.1 Méthodes de descente.................................... 220 3.3.2 Algorithme du gradient conjugué.............................. 224 3.3.3 Méthodes de Newton et Quasi–Newton........................... 227 3.3.4 Résumé sur les méthodes d’optimisation........................... 230 3.3.5 Exercices........................................... 231 3.4 Optimisation sous contraintes..................................... 258 3.4.1 Définitions.......................................... 258 3.4.2 Existence – Unicité – Conditions d’optimalité simple.................... 258 3.4.3 Conditions d’optimalité dans le cas de contraintes égalité.................. 259 3.4.4 Contraintes inégalités..................................... 262 3.4.5 Exercices........................................... 263 3.5 Algorithmes d’optimisation sous contraintes............................. 267 3.5.1 Méthodes de gradient avec projection............................ 267 3.5.2 Méthodes de dualité..................................... 269 3.5.3 Exercices........................................... 272 4 Equations différentielles 275 4.1 Introduction.............................................. 275 4.2 Consistance, stabilité et convergence................................. 278 4.3 Théorème général de convergence.................................. 280 4.4 Exemples............................................... 282 4.5 Explicite ou implicite ?........................................ 283 4.5.1 L’implicite gagne......................................... 283 4.5.2 L’implicite perd.......................................... 284 4.5.3 Match nul........................................... 284 4.6 Etude du schéma d’Euler implicite.................................. 284 4.7 Exercices............................................... 286 4.8 Corrigés................................................ 294 Analyse numérique I, télé-enseignement, L3 2 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 Introduction L’objet de l’analyse numérique est de concevoir et d’étudier des méthodes de résolution de certains problèmes mathématiques, en général issus de la modélisation de problèmes “réels", et dont on cherche à calculer la solution à l’aide d’un ordinateur. Le cours est structuré en quatre grands chapitres : – Systèmes linéaires – Systèmes non linéaires – Optimisation – Equations différentielles. On pourra consulter les ouvrages suivants pour ces différentes parties (ceci est une liste non exhaustive !) : – A. Quarteroni, R. Sacco et F. Saleri, Méthodes Numériques : Algorithmes, Analyse et Applications, Springer 2006. – P.G. Ciarlet, Introduction à l’analyse numérique et à l’optimisation, Masson, 1982, (pour les chapitre 1 à 3 de ce polycopié). – M. Crouzeix, A.L. Mignot, Analyse numérique des équations différentielles, Collection mathématiques ap- pliquées pour la maitrise, Masson, (pour le chapitre 4 de ce polycopié). – J.P. Demailly, Analyse numérique et équations différentielles Collection Grenoble sciences Presses Univer- sitaires de Grenoble – L. Dumas, Modélisation à l’oral de l’agrégation, calcul scientifique, Collection CAPES/Agrégation, Ellipses, 1999. – E. Hairer, polycopié du cours "Analyse Numérique", http ://www.unige.ch/ hairer/polycop.html – J. Hubbard, B. West, Equations différentielles et systèmes dynamiques, Cassini. – J. Hubbard et F. Hubert, Calcul Scientifique, Vuibert. – P. Lascaux et R. Théodor, Analyse numérique matricielle appliquée à l’art de l’ingénieur, tomes 1 et 2, Masson, 1987 – L. Sainsaulieu, Calcul scientifique cours et exercices corrigés pour le 2ème cycle et les éécoles d’ingénieurs, Enseignement des mathématiques, Masson, 1996. – M. Schatzman, Analyse numérique, cours et exercices, (chapitres 1,2 et 4). – D. Serre, Les matrices, Masson, (2000). (chapitres 1,2 et 4). – P. Lascaux et R. Theodor, Analyse numérique sappliquée aux sciences de l’ingénieur, Paris, (1994) – R. Temam, Analyse numérique, Collection SUP le mathématicien, Presses Universitaires de France, 1970. Et pour les anglophiles... – M. Braun, Differential Equations and their applications, Springer, New York, 1984 (chapitre 4). – G. Dahlquist and A. Björck, Numerical Methods, Prentice Hall, Series in Automatic Computation, 1974, Englewood Cliffs, NJ. 3 TABLE DES MATIÈRES TABLE DES MATIÈRES – R. Fletcher, Practical methods of optimization, J. Wiley, New York, 1980 (chapitre 3). – G. Golub and C. Van Loan, Matrix computations, The John Hopkins University Press, Baltimore (chapitre 1). – R.S. Varga, Matrix iterative analysis, Prentice Hall, Englewood Cliffs, NJ 1962. Pour des rappels d’algègre linéaire : – Poly d’algèbre linéaire de première année, P. Bousquet, R. Herbin et F. Hubert, http ://www.cmi.univ- mrs.fr/ herbin/PUBLI/L1alg.pdf – Introduction to linear algebra, Gilbert Strang, Wellesley Cambridge Press, 2008 Analyse numérique I, télé-enseignement, L3 4 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 Chapitre 1 Systèmes linéaires 1.1 Objectifs On note Mn (IR) l’ensemble des matrices carrées d’ordre n. Soit A ∈ Mn (IR) une matrice inversible, et b ∈ IR n , on a comme objectif de résoudre le système linéaire Ax = b, c’est–à– dire de trouver x solution de :  x ∈ IR n (1.1) Ax = b Comme A est inversible, il existe un unique vecteur x ∈ IR n solution de (1.1). Nous allons étudier dans les deux chapitres suivants des méthodes de calcul de ce vecteur x : la première partie de ce chapitre sera consacrée aux méthodes “directes" et la deuxième aux méthodes “itératives". Nous aborderons ensuite en troisième partie les méthodes de résolution de problèmes aux valeurs propres. Un des points essentiels dans l’efficacité des méthodes envisagées concerne la taille des systèmes à résoudre. La taille de la mémoire des ordinateurs a augmenté de façon drastique de 1980 à nos jours. Le développement des méthodes de résolution de systèmes linéaires est liée à l’évolution des machines infor- matiques. C’est un domaine de recherche très actif que de concevoir des méthodes qui permettent de profiter au mieux de l’architecture des machines (méthodes de décomposition en sous domaines pour profiter des architectures parallèles, par exemple). Dans la suite de ce chapitre, nous verrons deux types de méthodes pour résoudre les systèmes linéaires : les méthodes directes et les méthodes itératives. Pour faciliter la compréhension de leur étude, nous commençons par quelques rappels d’algèbre linéaire. 1.2 Pourquoi et comment ? Nous donnons dans ce paragraphe un exemple de problème dont la résolution numérique recquiert la résolution d’un système linéaire, et qui nous permet d’introduire des matrices que nous allons beaucoup étudier par la suite. Nous commençons par donner ci-après après quelques rappels succincts d’algèbre linéaire, outil fondamental pour la résolution de ces systèmes linéaires. 1.2.1 Quelques rappels d’algèbre linéaire Quelques notions de base Ce paragraphe rappelle des notions fondamentales que vous devriez connaître à l’issue du cours d’algèbre linéaire de première année. On va commencer par revisiter le produit matriciel, dont la vision combinaison linéaire de lignes est fondamentale pour bien comprendre la forme matricielle de la procédure d’élimination de Gauss. 5 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES Soient A et B deux matrices carrées d’ordre n, et M = AB. Prenons comme exemple d’illustration       1 2 −1 0 5 4 A= ,B = et M = 0 1 3 2 3 2 On note ai,j bi,j et mi,j , i, j = 1,... n les coefficients respectifs de A, B et M. Vous savez bien sûr que n X mi,j = ai,k bk,j. (1.2) k=1 Si on écrit les matrices A et B sous forme de lignes (notées `i ) et colonnes (notées cj ) :   `1 (A)   A = ...  et B = c1 (B)... `n (B) `n (A) Dans nos exemples, on a donc         −1 0 `1 (A) = 1 2 , `2 (A) = 0 1 , c1 (B) = c2 (B) =. 3 2 L’expression (1.2) s’écrit encore mi,j = `i (A)cj (B), qui est le produit d’une matrice 1 × n par une matrice n × 1, qu’on peut aussi écrire sous forme d’un produit scalaire : mi,j = (`i (A))t · cj (B) où (`i (A))t désigne la matrice transposée, qui est donc maintenant une matrice n × 1 qu’on peut identifier à un vecteur de IR n. C’est la technique “habituelle” de calcul du produit de deux matrices. On a dans notre exemple :     0 m1,2 = `1 (A) c2 (B) = 1 2. 2     1 0 = (`i (A))t · cj (B) = · 2 2 = 4. Mais de l’expression (1.2), on peut aussi avoir l’expression des lignes et des colonnes de M = AB en fonction des lignes de B ou des colonnes de A : n X `i (AB) = ai,k `k (B) (1.3) k=1 Xn cj (AB) = bk,j ck (A) (1.4) k=1 Dans notre exemple, on a donc :       `1 (AB) = −1 0 +2 3 2 = 5 4 ce qui montre que la ligne 1 de AB est combinaison linéaire des lignes de B. Le colonnes de AB, par contre, sont des combinaisons linéaires de colonnes de A. Par exemple :       1 2 4 c2 (AB) = 0 +2 = 0 1 2 Il faut donc retenir que dans un produit matriciel AB, Analyse numérique I, télé-enseignement, L3 6 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES les colonnes de AB sont des combinaisons linéaires des colonnes de A les lignes de AB sont des combinaisons linéaires des lignes de B. Cette remarque est très importante pour la représentation matricielle de l’élimination de Gauss : lorqu’on calcule des systèmes équivalents, on effectue des combinaisons linéaires de lignes, et donc on multiplie à gauche par une matrice d’élimination. Le tableau ci-dessous est la traduction littérale de “Linear algebra in a nutshell”, par Gilbert Strang 1 Pour une matrice carrée A, on donne les caractérisations du fait qu’elle est inversible ou non. A inversible A non inversible Les vecteurs colonne sont indépendants Les vecteurs colonne sont liés Les vecteurs ligne sont indépendants Les vecteurs ligne sont liés Le déterminant est non nul Le déterminant est nul Ax = 0 a une unique solution x = 0 Ax = 0 a une infinité de solutions. Le noyau de A est réduit à {0} Le noyau de A contient au moins un vecteur non nul. Ax = b a une solution unique x = A−1 b Ax = b a soit aucune solution, soit une infinité. A a n (nonzero) pivots A a r < n pivots A est de rang maximal : rg(A) = n. rg(A) = r < n La forme totatement échelonnée R de A est la matrice identité R a au moins une ligne de zéros. L’image de A est tout IR n. L’image de A est strictement incluse dans IR n. L’espace L(A) engendré par les lignes de A est tout IR n. L(A) est de dimension r < n Toutes les valeurs propres de A sont non nulles Zero est valeur propre de A. At A is symétrique définie positive 2 At A n’est que semi- définie. TABLE 1.1: Extrait de “Linear algebra in a nutshell”, G. Strang On rappelle pour une bonne lecture de ce tableau les quelques définitions suivantes : Définition 1.1 (Pivot). Soit A ∈ Mn (IR) une matrice carrée d’ordre n. On appelle pivot de A le premier élément non nul de chaque ligne dans la forme échelonnée de A obtenue par élimination de Gauss. Si la matrice est inversible, elle a donc n pivots (non nuls). Définition 1.2 (Valeurs propres). Soit A ∈ Mn (IR) une matrice carrée d’ordre n. On appelle valeur propre de A tout λ ∈ Cl tel qu’il existe x ∈ Cln , x 6= 0 tel que Ax = λx. L’élément x est appelé vecteur propre de A associé à λ. Définition 1.3 (Déterminant). Il existe une unique application, notée det de Mn (IR) dans IR qui vérifie les pro- priétés suivantes (D1) Le déterminant de la matrice identité est égal à 1. (D2) Si la matrice à est obtenue à partir de A par échange de deux lignes, alors detà = −detA. 1. Voir la page web de Strang www.mit.edu/~gs pour une foule d’informations et de cours sur l’algèbre linéaire. Analyse numérique I, télé-enseignement, L3 7 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES (D3) Le déterminant est une fonction linéaire de chacune des lignes de la matrice A. (D3a) (multiplication par un scalaire) si à est obtenue à partir de A en multipliant tous les coefficients d’une ligne par λ ∈ IR, alors det(Ã) = λdet(A).       `1 (A) `1 (A) `1 (A) ..  ..  ..  .  .  .        (D3b) (addition) si A = `k (A), à = `k (A) et B = `k (A) + `k (A)    ˜   ˜ , alors .  .  .  ..  ..  ..  `n (A) `n (A) `n (A) det(B) = det(A) + det(Ã). On peut déduire de ces trois propriétés fondamentales un grand nombre de propriétés importantes, en particulier le fait que det(AB) = detA detB et que le déterminant d’une matrice inversible est le produit des pivots : c’est de cette manière qu’on le calcule sur les ordinateurs. En particulier on n’utilise jamais la formule de Cramer, beaucoup trop coûteuse en termes de nombre d’opérations. On rappelle que si A ∈ Mn (IR) une matrice carrée d’ordre n, les valeurs propres sont les racines du polynôme caractéristique PA de degré n, qui s’écrit : PA (λ)) = det(A − λI). Matrices diagonalisables Un point important de l’algèbre linéaire, appelé “réduction des endomorphismes” dans les programmes français, consiste à se demander s’il existe une base de l’espace dans laquelle la matrice de l’application linéaire est diago- nale ou tout au moins triangulaire (on dit aussi trigonale). Définition 1.4 (Matrice diagonalisable dans IR). Soit A une matrice réelle carrée d’ordre n. On dit que A est diagonalisable dans IR s’il existe une base (u1 ,... , un ) de IR n et des réels λ1 ,... , λn (pas forcément distincts) tels que Aui = λi ui pour i = 1,... , n. Les réels λ1 ,... , λn sont les valeurs propres de A, et les vecteurs u1 ,... , un sont les vecteurs propres associés. Vous connaissez sûrement aussi la diagonalisation dans Cl : une matrice réelle carrée d’ordre n admet toujours n valeurs propres dans Cl, qui ne sont pas forcément identiques. Une matrice est diagonalisable dans Cl s’il existe une base (u1 ,... , un ) de Cln et des nombres complexes λ1 ,... , λn (pas forcément distincts) tels que Aui = λi ui pour i = 1,... , n. Ceci est vérifié si la dimension de chaque sous espace propre Ei = Ker(A − λId) (appelée multiplicité géométrique) est égale a la multiplicité algébrique de λi , c.à.d. son ordre de multiplicité en tant que racine du polynôme caractéristique.  0 0 Par exemple la matrice A = n’est pas diagonalisable dans Cl (ni évidemment, dans IR). Le polynôme 1 0 caractéristique de A est PA (λ) = λ2 , l’unique valeur propre est donc 0, qui est de multiplicité algébrique 2, et de multiplicité géométrique 1, car le sous espace propre associé à la valeur propre nulle est F = {x ∈ IR 2 ; Ax = 0} = {x = (0, t), t ∈ IR}, qui est de dimension 1. Ici et dans toute la suite, comme on résout des systèmes linéaires réels, on préfère travailler avec la diagonalisation dans IR ; cependant il y a des cas où la diagonalisation dans Cl est utile et même nécessaire (étude de stabilité des Analyse numérique I, télé-enseignement, L3 8 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES systèmes diférentiels, par exemple). Par souci de clarté, nous préciserons toujours si la diagonalisation considérée est dans IR ou dans Cl. Lemme 1.5. Soit A une matrice réelle carrée d’ordre n, diagonalisable dans IR. Alors A = P diag(λ1 ,... , λn )P −1 , ‘ où P est la matrice dont les vecteurs colonnes sont égaux aux vecteurs propres u1 ,... , un associées au valeurs propres λ1 ,... , λn D ÉMONSTRATION – Par définition d’un vecteur propre, on a Aui = λi ui pour i = 1,... N , et donc, en notant P la matrice dont les colonnes sont les vecteurs propres ui ,     Au1... Aun = A u1... un = AP et donc   λ1 0... 0 ....      0 λ2..  AP = λ1 u1... λn un = u1... un  .  = P diag(λ1 ,... , λn ). ........ ...  0... 0 λn Notons que dans ce calcul, on a fortement utilisé la multiplication des matrices par colonnes, c.à.d. X ci (AB) = ai,j cj (B). j=1,n Remarquons que P lest aussi la matrice définie (de manière unique) par P ei = ui , où (ei )i=1,...,n est la base canonique de IR n , c’est-à-dire que (ei )j = δi,j. La matrice P est appelée matrice de passage de la base (ei )i=1,...,n à la base (ui )i=1,...,n ; (il est bien clair que la i-ème colonne de P est constituée des composantes de ui dans la base canonique (e1 ,... , en ). La matrice P est inversible car les vecteurs propres forment une base, et on peut donc aussi écrire : P −1 AP = diag(λ1 ,... , λn ) ou A = P diag(λ1 ,... , λn )P −1. La diagonalisation des matrices réelles symétriques est un outil qu’on utilisera souvent dans la suite, en particulier dans les exercices. Il s’agit d’un résultat extrêmement important. Lemme 1.6 (Une matrice symétrique est diagonalisable dans IR). Soit E un espace vectoriel sur IR de dimension finie : dimE = n, n ∈ IN∗ , muni d’un produit scalaire i.e. d’une application E × E → IR, (x, y) → (x | y)E , qui vérifie : ∀x ∈ E, (x | x)E ≥ 0 et (x | x)E = 0 ⇔ x = 0, ∀(x, y) ∈ E 2 , (x | y)E = (y | x)E , ∀y ∈ E, l’application de E dans IR, définie par x → (x | y)E est linéaire. p Ce produit scalaire induit une norme sur E, kxk = (x | x)E. Soit T une application linéaire de E dans E. On suppose que T est symétrique, c.à.d. que (T (x) | y)E = (x | T (y))E , ∀(x, y) ∈ E 2. Alors il existe une base orthonormée (f 1 ,... , f n ) de E (c.à.d. telle que (f i , f j )E = δi,j ) et (λ1... λn ) ∈ IR n tels que T (fi ) = λi f i pour tout i ∈ {1... n}. Analyse numérique I, télé-enseignement, L3 9 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES n Conséquence immédiate : Dans le cas où E = IR Pn, le produit scalaire canonique de x = (x1 ,... , xn ) et t y = (y1 ,... , yn ) est défini par (x | y)E = x · y = i=1 xi yi. Si A ∈ Mn (IR) est une matrice symétrique, alors t l’application T définie de E dans E par : T (x) = Ax est linéaire, et : (T x|y) = Ax · y = x · At y = x · Ay = (x | T y). Donc T est linéaire symétrique. Par le lemme précédent, il existe (f 1 ,... , f n ) et (λ1... λn ) ∈ IR tels que T f i = Af i = λi f i ∀ i ∈ {1,... , n} et fi · f j = δi,j , ∀ (i, j) ∈ {1,... , n}2. Interprétation algébrique : Il existe une matrice de passage P de (e1 ,... , en ) base canonique dans (f 1 ,... , f n ) dont la première colonne de P est constituée des coordonnées de fi dans (e1... en ). On a : P ei = fi. On a alors P −1 AP ei = P −1 Af i = P −1 (λi fi ) = λi ei = diag(λ1 ,... , λn )ei , où diag(λ1 ,... , λn ) désigne la matrice diagonale de coefficients diagonaux λ1 ,... , λn. On a donc :   λi 0 ..  P −1 AP = .  = D. 0 λn De plus P est orthogonale, i.e. P −1 = P t. En effet, P t P ei · ej = P ei · P ej = (fi |fj ) = δi,j ∀i, j ∈ {1... n}, et donc (P t P ei − ei ) · ej = 0 ∀j ∈ {1... n} ∀i ∈ {1,... n}. On en déduit P t P ei = ei pour tout i = 1,... n, i.e. P t P = P P t = Id. D ÉMONSTRATION du lemme 1.6 Cette démonstration se fait par récurrence sur la dimension de E. 1ère étape. e On suppose dimE = 1. Soit e ∈ E, e 6= 0, alors E = IRe = f1 avec f 1 =. Soit T : E → E linéaire symétrique, kek on a : T f 1 ∈ IRf 1 donc il existe λ1 ∈ IR tel que T f 1 = λ1 f 1. 2ème étape. On suppose le lemme vrai si dimE < n. On montre alors le lemme si dimE = n. Soit E un espace vectoriel normé sur IR tel que dimE = n et T : E → E linéaire symétrique. Soit ϕ l’application définie par : ϕ: E → IR x → (T x|x). L’application ϕ est continue sur la sphère unité S1 = {x ∈ E| kxk = 1} qui est compacte car dimE < +∞ ; il existe donc e ∈ S1 tel que ϕ(x) ≤ ϕ(e) = (T e | e) = λ pour tout x ∈ E. Soit y ∈ E \ {0}, et soit t ∈]0, kyk 1 [ alors e + ty 6= 0. On en déduit que :   e + ty (e + ty) e + ty ∈ S1 et donc ϕ(e) = λ ≥ T | ke + tyk ke + tyk ke + tyk E donc λ(e + ty | e + ty)E ≥ (T (e + ty) | e + ty). En développant on obtient : λ[2t(e | y) + t2 (y | y)E ] ≥ 2t(T (e) | y) + t2 (T (y) | y)E. Comme t > 0, ceci donne : λ[2(e | y) + t(y | y)E ] ≥ 2(T (e) | y) + t(T (y) | y)E. En faisant tendre t vers 0+ , on obtient 2λ(e | y)E ≥ 2(T (e) | y), Soit 0 ≥ (T (e) − λe | y) pour tout y ∈ E \ {0}. De même pour z = −y on a 0 ≥ (T (e) − λe|z) donc (T (e) − λe | y) ≥ 0. D’où (T (e) − λe | y) = 0 pour tout y ∈ E. On en déduit que T (e) = λe. On pose f n = e et λn = λ. L Soit F = {x ∈ E; (x | e) = 0}, on a donc F 6= E, et E = F IRe : on peut décomposer x ∈ E comme (x = x − (x | e)e + (x | e)e). L’application S = T |F est linéaire symétrique et on a dimF = n − 1. et S(F ) ⊂ F. On peut donc utiliser l’hypothèse de récurrence : ∃(λ1... λn−1 ) ∈ IR n et ∃(f 1... f n−1 ) ∈ E n tels que ∀ i ∈ {1... n − 1}, Sf i = T f i = λi f i , et ∀i, j ∈ {1... n − 1}, f i · f j = δi,j. Et donc (λ1... λn ) et (f 1 ,... , f n ) conviennent. Analyse numérique I, télé-enseignement, L3 10 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES x x x ui x x u(x) x x x x x x0 = 0 x1 ··· xi = ih ··· xN +1 = 1 F IGURE 1.1: Solution exacte et approchée de −u00 = f 1.2.2 Discrétisation de l’équation de la chaleur Dans ce paragraphe, nous prenons un exemple très simple pour obtenir un système linéaire à partir de la discréti- sation d’un problème continu. L’équation de la chaleur unidimensionnelle Discrétisation par différences finies de −u00 = f Soit f ∈ C([0, 1], IR). On cherche u tel que − u00 (x) = f (x) (1.5a) u(0) = u(1) = 0. (1.5b) Remarque 1.7 (Problèmes aux limites, problèmes à conditions initiales). L’équation différentielle −u00 = f admet une infinité de solutions. Pour avoir existence et unicité, il est nécessaire d’avoir des conditions supplémentaires. Si l’on considère deux conditions en 0 (ou en 1, l’origine importe peu) on a ce qu’on appelle un problème de Cauchy, ou problème à conditions initiales. Le problème (1.5) est lui un problème aux limites : il y a une condition pour chaque bord du domaine. En dimension supérieure, le problème −∆u = f nécessite une condition sur au moins “un bout” de frontière pour être bien posé : voir le cours d’équations aux dérivées partielles de master pour plus de détails à ce propos. On peut montrer (on l’admettra ici) qu’il existe une unique solution u ∈ C 2 ([0, 1], IR). On cherche à calculer u de manière approchée. On va pour cela introduire la méthode de discrétisation dite par différences finies. Soit n ∈ IN∗ , on définit h = 1/(n + 1) le pas de discrétisation, c.à.d. la distance entre deux points de discrétisation, et pour i = 0,... , n + 1 on définit les points de discrétisation xi = ih (voir Figure 1.1), qui sont les points où l’on va écrire l’équation −u00 = f en vue de se ramener à un système discret, c.à.d. à un système avec un nombre fini d’inconnues u1 ,... , un. Remarquons que x0 = 0 et xn+1 = 1, et qu’en ces points, u est spécifiée par les conditions limites (1.5b). Soit u(xi ) la valeur exacte de u en xi. On écrit la première équation de (1.5a) en chaque point xi , pour i = 1... n. −u00 (xi ) = f (xi ) = bi ∀i ∈ {1... n}. (1.6) Supposons que u ∈ C 4 ([0, 1], IR) (ce qui est vrai si f ∈ C 2 ). Par développement de Taylor, on a : h2 00 h3 000 h4 (4) u(xi+1 ) = u(xi ) + hu0 (xi ) + u (xi ) + u (xi ) + u (ξi ), 2 6 24 h2 h3 000 h4 (4) u(xi−1 ) = u(xi ) − hu0 (xi ) + u00 (xi ) − u (xi ) + u (ηi ), 2 6 24 Analyse numérique I, télé-enseignement, L3 11 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES avec ξi ∈ (xi , xi+1 ) et ηi ∈ (xi , xi+1 ). En sommant ces deux égalités, on en déduit que : h4 (4) h4 u(xi+1 ) + u(xi−1 ) = 2u(xi ) + h2 u00 (xi ) + u (ξi ) + u(4) (ηi ). 24 24 On définit l’erreur de consistance, qui mesure la manière dont on a approché −u00 (xi ) ; l’erreur de consistance Ri au point xi est définie par u(xi+1 ) + u(xi−1 ) − 2u(xi ) Ri = −u00 (xi ) −. (1.7) h2 On a donc : u(xi+1 ) + u(xi−1 ) − 2u(xi ) |Ri | = − + u00 (xi ) h2 h4 (4) h4 ≤ u (ξi ) + u(4) (ηi ) 24 24 2 h ≤ ku(4) k∞. (1.8) 12 où ku(4) k∞ = supx∈]0,1[ |u(4) (x)|. Cette majoration nous montre que l’erreur de consistance tend vers 0 comme h2 : on dit que le schéma est consistant d’ordre 2. On introduit alors les inconnues (ui )i=1,...,n qu’on espère être des valeurs approchées de u aux points xi et qui sont les composantes de la solution (si elle existe) du système suivant ( ui+1 + ui−1 − 2ui − = bi , ∀i ∈ J; 1, nK, h2 (1.9) u0 = un+1 = 0.   u1 ..  On cherche donc u = .  ∈ IR n solution de (1.9). Ce système peut s’écrire sous forme matricielle :Kn u = b un   b1 où b = ...  et Kn est la matrice carrée d’ordre n de coefficients (ki,j )i,j=1,n définis par :   bn  2  , ∀ i = 1,... , n,  ki,i  = h2 1 (1.10)   ki,j = − 2 , ∀ i = 1,... , n, j = i ± 1,  h ki,j = 0, ∀ i = 1,... , n, |i − j| > 1. On remarque immédiatement que Kn est tridiagonale. On peut montrer que Kn est symétrique définie positive (voir exercice 10 page 19), et elle est donc inversible Le système Kn u = b admet donc une unique solution. C’est bien, mais encore faut il que cette solution soit ce qu’on espérait, c.à.d. que chaque valeur ui soit une approximation pas trop mauvaise de u(xi ). On appelle erreur de discrétisation en xi la différence de ces deux valeurs : ei = u(xi ) − ui , i = 1,... , n. (1.11) mSi on appelle e le vecteur de composantes ei , on déduit de la définition 1.11 de l’erreur de consistance et des équations (exactes) 1.6 que Kn e = R et donc e = Kn−1 R. (1.12) Le fait que le schéma soit consistant est une bonne chose, mais cela ne suffit pas à montrer que le schéma est convergent, c.à.d. que l’erreur entre maxi=1,...,n ei tend vers 0 lorsque h tend vers 0, parce que A dépend de Analyse numérique I, télé-enseignement, L3 12 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES h ! Pour cela, il faut de plus que le schéma soit stable, au sens où l’on puisse montrer que kKn−1 k est borné indépendamment de h, ce qui revient à trouver une estimation sur les valeurs approchées ui indépendante de h. La stabilité et la convergence font l’objet de l’exercice 47, où l’on montre que le schéma est convergent, et qu’on a l’estimation d’erreur suivante : h2 (4) max {|ui − u(xi )|} ≤ ku k∞. i=1...n 96 Cette inégalité donne la précision de la méthode (c’est une méthode dite d’ordre 2). On remarque en particulier que si on raffine la discrétisation, c’est–à–dire si on augmente le nombre de points n ou, ce qui revient au même, si on diminue le pas de discrétisation h, on augmente la précision avec laquelle on calcule la solution approchée. L’équation de la chaleur bidimensionnelle Prenons maintenant le cas d’une discrétisation du Laplacien sur un carré par différences finies. Si u est une fonction de deux variables x et y à valeurs dans IR, et si u admet des dérivées partielles d’ordre 2 en x et y, l’opérateur laplacien est défini par ∆u = ∂xx u + ∂yy u. L’équation de la chaleur bidimensionnelle s’écrit avec cet opérateur. On cherche à résoudre le problème : −∆u = f sur Ω =]0, 1[×]0, 1[, (1.13) u = 0 sur ∂Ω, On rappelle que l’opérateur Laplacien est défini pour u ∈ C 2 (Ω), où Ω est un ouvert de IR 2 , par ∂2u ∂2u ∆u = + 2. ∂x2 ∂y Définissons une discrétisation uniforme du carré par les points (xi , yj ), pour i = 1,... , M et j = 1,... , M avec xi = ih, yj = jh et h = 1/(M + 1), representée en figure 1.2 pour M = 6. On peut alors approcher les dérivées secondes par des quotients différentiels comme dans le cas unidimensionnel (voir page 11), pour obtenir un système linéaire : Au = b où A ∈ Mn (IR) et b ∈ IR n avec n = M 2. Utilisons l’ordre“lexicographique" pour numéroter les inconnues, c.à.d. de bas en haut et de gauche à droite : les inconnues sont alors numérotées de 1 à n = M 2 et le second membre s’écrit b = (b1 ,... , bn )t. Les composantes b1 ,... , bn sont définies par :pour i, j = 1,... , M , on pose k = j + (i − 1)M et bk = f (xi , yj ). Les coefficients de A = (ak,` )k,l=1,n peuvent être calculés de la manière suivante :   Pour i, j = 1,... , M, on pose k = j + (i − 1)M,   4     ak,k = 2,   (h   1    − 2 si j 6= M,  ak,k+1 = h   0 sinon,   (   1   − 2 si j 6= 1,   ak,k−1 = h   0 sinon, (  1   − 2 si i < M,   ak,k+M = h   0 sinon,   (   1    − 2 si i > 1,  ak,k−M = h   0 sinon,          Pour  k = 1,... , n, et ` = 1,... , n; ak,` = 0, ∀ k = 1,... , n, 1 < |k − `| < n ou |k − `| > n. Analyse numérique I, télé-enseignement, L3 13 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES y 31 32 33 34 35 36 25 26 27 28 29 30 19 20 21 22 23 24 13 14 15 16 17 18 7 8 9 10 11 12 i=1 1 2 3 4 5 6 x j=1 F IGURE 1.2: Ordre lexicographique des inconnues, exemple dans le cas M = 6 La matrice est donc tridiagonale par blocs, plus précisément si on note   4 −1 0...... 0 −1 4 −1 0... 0    ..........  0 .....   D=. , ..    ......  0... −1 0... 0 −1 4 les blocs diagonaux (qui sont des matrices de dimension M × M ), on a :   D −Id 0...... 0 −Id D −Id 0... 0    0 −Id D −Id ··· 0   A=.......... ,.. (1.14) . ....  . ..   0. −Id D −Id 0... 0 −Id D où Id désigne la matrice identité d’ordre M , et 0 la matrice nulle d’ordre M. Matrices monotones, ou à inverse positive Une propriété qui revient souvent dans l’étude des matrices issues de la discrétisation d’équations différentielles est le fait que si leur action sur un vecteur u donne un vecteur positif v (composante par composante) alors le vecteur u de départ doit être positif (composante par composante) ; on dit souvent que la matrice est “monotone”, ce qui n’est pas un terme très évocateur... Dans ce cours, on lui préfèrera le terme “à inverse positive” ; en effet, on montre à la proposition qu’une matrice A est monotone si et seulement si elle inversible et à inverse positive. Analyse numérique I, télé-enseignement, L3 14 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES Définition 1.8 (IP-matrice ou matrice monotone). Si x ∈ IR n , on dit que x ≥ 0 [resp. x > 0] si toutes les composantes de x sont positives [resp. strictement positives]. Soit A ∈ Mn (IR), on dit que A est une matrice monotone si elle vérifie la propriété suivante : Si x ∈ IR n est tel que Ax ≥ 0, alors x ≥ 0, ce qui peut encore s’écrire : {x ∈ IR n t.q. Ax ≥ 0} ⊂ {x ∈ IR n t.q. x ≥ 0}. Proposition 1.9 (Caractérisation des matrices monotones). Une matrice A est monotone si et seulement si elle inversible et à inverse positive (c.à.d. dont tous les coefficients sont positifs). La démonstration de ce résultat est l’objet de l’exercice 8. Retenez que toute matrice monotone est inversible et d’inverse positive que cette propriété de monotonie est utilisée pour établir une borne de kA−1 k pour la matrice de discrétisation du Laplacien, dont on a besoin pour montrer la convergence du schéma. C’est donc une propriété qui est importante au niveau de l’analyse numérique. 1.2.3 Exercices Exercice 1 (Théorème du rang). Soit A ∈ Mn,p (IR) (n, p ≥ 1). On rappelle que Ker(A) = {x ∈ IR p ; Ax = 0}, Im(A) = {Ax, x ∈ IR p } et rang(A) = dim(Im(A)). Noter que Ker(A) ⊂ IR p et Im(A) ⊂ IR n. Soit f1 ,... , fr une base de Im(A) (donc r ≤ n) et, pour i ∈ {1,... , r}, ai tel que Aai = fi. 1. Montrer que la famille a1 ,... , ar est une famille libre de IR p (et donc r ≤ p). 2. On note G le sous espace vectoriel de IR p engendré par a1 ,... , ar. Montrer que IR p = G ⊕ Ker(A). En déduire que (théorème du rang) p = dim(Ker(A)) + dim(Im(A)). 3. On suppose ici que n = p. Montrer que l’application x 7→ Ax (de IR n dans IR n ) est injective si et seulement si elle est surjective. Exercice 2 (rang(A)=rang(At )). Soit A ∈ Mn,p (IR) (n, p ≥ 1). 1. Soient P une matrice inversible de Mn (IR) et Q une matrice inversible de Mp (IR). Montrer que dim(Im(P A)) = dim(Im(AQ)) = dim(Im(A)). Montrer aussi que les matrices P t et Qt sont inversibles. Soit f1 ,... , fr une base de Im(A) (donc r ≤ p) et, pour i ∈ {1,... , r}, ai tel que Aai = fi. Soit ar+1 ,... , ap une base de Ker(A) (si Ker(A) 6= {0}). La famille a1 ,... , an est une base de IR p (voir l’exercice 1). De même, on complète (si r < n) f1 ,... , fr par fr+1 ,... , fn de manière à avoir une base f1 ,... , fn de IR n. Enfin, on note P et Q les matrices telles que P ei = ai (pour tout i = 1,... , p) et Qfj = ēj (pour tout j = 1,... , n) ou e1 ,... , ep est la base canonique de IR p et ē1 ,... , ēn est la base canonique de IR n. On pose J = QAP. 2. Montrer que P et Q sont inversibles. 3. calculer les colonnes de J et de J t et en déduire que les matrices J et J t sont de même rang. 4. Montrer que A et At sont de même rang. Exercice 3 (Vrai ou faux ? Motiver les réponses... ). Corrigé en page 21 On suppose dans toutes les questions suivantes que n ≥ 2. 1. Soit Z ∈ IR n un vecteur non nul. La matrice ZZ t est inversible. 2. La matrice inverse d’une matrice triangulaire inférieure est triangulaire supérieure. Analyse numérique I, télé-enseignement, L3 15 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES 3. Les valeurs propres sont les racines du polynôme caractéristique. 4. Toute matrice inversible est diagonalisable dans IR. 5. Toute matrice inversible est diagonalisable dans Cl. 6. Le déterminant d’une matrice A est égal au produit de ses valeurs propres. 7. Soit A une matrice carrée telle que Ax = 0 =⇒ x = 0, alors A est inversible. 8. Soit A une matrice carrée telle que Ax ≥ 0 =⇒ x ≥ 0, alors A est inversible. 9. Une matrice symétrique est inversible. 10. Une matrice symétrique définie positive est inversible. Le système linéaire n+1 X ai,j xj = 0 pour tout i = 1,... , n j=1 admet toujours une solution non nulle. Exercice 4 (Sur quelques notions connues). Corrigé en page 21 1. Soit A une matrice carrée d’ordre n et b ∈ IR 3. Peut il exister exactement deux solutions distinctes au système Ax = b ? 2. Soient A, B et C de dimensions telles que AB et BC existent. Montrer que si AB = Id et BC = Id, alors A = C. 3. Combien y a -t-il de matrices carrées d’ordre 2 ne comportant que des 1 ou des 0 comme coefficients ? Combien d’entre elles sont inversibles ?   3 2 4. Soit B =. Montrer que B 1024 = Id. −5 −3 Exercice 5 (La matrice K3 ). Suggestions en page 21. Corrigé en page 22 Soit f ∈ C([0, 1], IR). On cherche u tel que − u00 (x) = f (x), ∀x ∈ (0, 1), (1.15a) u(0) = u(1) = 0. (1.15b) 1. Calculer la solution exacte u(x) du problèmes lorsque f est la fonction identiquement égale à 1 (on admettra que cette solution est unique), et vérifier que u(x) ≥ 0 pour tout x ∈ [0, 1]. On discrétise le problème suivant par différences finies, avec un pas h = 1 4 avec la technique vue en cours. 2. A l’aide de dévloppements de Taylor, écrire l’approximation de u (xi ) au deuxième ordre en fonction de 00 u(xi ), u(xi−1 ) et u(xi+1 ). En déduire le schéma aux différences finies pour l’approximation de (1.15), qu’on écrira sous la forme : K3 u = b, (1.16)       u1 b1 f (x1 ) où K3 est la matrice de discrétisation qu’on explicitera, u = u2  et b = b2  = f (x2 ). u3 b3 f (x3 ) 3. Résoudre le système linéaire (1.16) par la méthode de Gauss. Comparer ui et u(xi ) pour i = 1, 2, 3, et expliquer pourquoi l’erreur de discrétisation u(xi ) − ui est nulle. 4. Reprendre les questions précédentes en remplaçant les conditions limites (1.15b) par : u(0) = 0, u0 (1) = 0. (1.17) Analyse numérique I, télé-enseignement, L3 16 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES 5. Soit c ∈ IR. On considère maintenant le problème suivant : − u00 (x) = c, ∀x ∈ (0, 1), (1.18a) 0 0 u (0) = u (1) = 0, (1.18b) (a) Montrer que le problème (1.18) admet soit une infinité de solutions, soit pas de solution. e3u = e (b) Ecrire la discrétisation du problème (1.18), toujours avec h = 41 , sous la forme K b en explicitant e e K3 et b. (c) Montrer que la matrice K e 3 n’est pas inversible : on part d’un problème continu mal posé, et on obtient par discrétisation un problème discret mal posé... Exercice 6 (Matrices symétriques définies positives). Suggestions en page 21, corrigé en page 23. On rappelle que toute matrice A ∈ Mn (IR) symétrique est diagonalisable dans IR (cf. lemme 1.6 page 9). Plus précisément, on a montré en cours que, si A ∈ Mn (IR) est une matrice symétrique, il existe une base de IR n , notée {f 1 ,... , f n }, et il existe λ1 ,... , λn ∈ IR t.q. Af i = λi f i , pour tout i ∈ {1,... , n}, et f i · f j = δi,j pour tout i, j ∈ {1,... , n} (x · y désigne le produit scalaire de x avec y dans IR n ). 1. Soit A ∈ Mn (IR). On suppose que A est symétrique définie positive, montrer que les éléments diagonaux de A sont strictements positifs. 2. Soit A ∈ Mn (IR) une matrice symétrique. Montrer que A est symétrique définie positive si et seulement si toutes les valeurs propres de A sont strictement positives. 3. Soit A ∈ Mn (IR). On suppose que A est symétrique définie positive. Montrer qu’on peut définir une unique 1 matrice B ∈ Mn (IR), symétrique définie positive t.q. B 2 = A (on note B = A 2 ). Exercice 7 (Diagonalisation dans IR ). Corrigé en page 24. Soit E un espace vectoriel réel de dimension n ∈ IN muni d’un produit scalaire, noté (·, ·). Soient T et S deux applications linéaires symétriques de E dans E (T symétrique signifie (T x, y) = (x, T y) pour tous x, y ∈ E). On suppose que T est définie positive (c’est-à-dire (T x, x) > 0 pour tout x ∈ E \ {0}). 1. Montrer que T est inversible. Pour x, y ∈ E, on pose (x, y)T = (T x, y). Montrer que l’application (x, y) 7→ (x, y)T définit un nouveau produit scalaire sur E. 2. Montrer que T −1 S est symétrique pour le produit scalaire défini à la question précédente. En déduire, avec le lemme 1.6 page 9, qu’il existe une base de E, notée {f 1 ,... , f n } et une famille {λ1 ,... , λn } ⊂ IR telles que T −1 Sf i = λi f i pour tout i ∈ {1,... , n} et t.q. (T f i , f j ) = δi,j pour tout i, j ∈ {1,... , n}. Exercice 8 (IP-matrice). Corrigé en page 25 Soit n ∈ IN? , on note Mn (IR) l’ensemble des matrices de n lignes et n colonnes et à coefficients réels. Si x ∈ IR n , on dit que x ≥ 0 [resp. x > 0] si toutes les composantes de x sont positives [resp. strictement positives]. Soit A ∈ Mn (IR), on dit que A est une IP-matrice si elle vérifie la propriété suivante : Si x ∈ IR n est tel que Ax ≥ 0, alors x ≥ 0, ce qui peut encore s’écrire : {x ∈ IR n t.q. Ax ≥ 0} ⊂ {x ∈ IR n t.q. x ≥ 0}. 1. Soit A = (ai,j )i,j=1,...,n ∈ Mn (IR). Montrer que A est une IP-matrice si et seulement si A est inversible et A−1 ≥ 0 (c’est-à-dire que tous les coefficients de A−1 sont positifs).   a b 2. Soit A = une matrice réelle d’ordre 2. Montrer que A est une IP-matrice si et seulement si : c d    ad < bc,  ad > bc, a < 0, d < 0 ou a > 0, d > 0, (1.19)   b ≥ 0, c ≥ 0 b ≤ 0, c ≤ 0. Analyse numérique I, télé-enseignement, L3 17 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES 3. Montrer que si A ∈ Mn (IR) est une IP-matrice alors At (la transposée de A) est une IP-matrice. 4. Montrer que si A est telle que n X ai,j ≤ 0, pour tout i, j = 1,... , n, i 6= j, et ai,i > |ai,j |, pour tout i = 1,... , n, (1.20) j=1 j6=i alors A est une IP-matrice. 5. En déduire que si A est telle que n X ai,j ≤ 0, pour tout i, j = 1,... , n, i 6= j, et ai,i > |aj,i |, pour tout i = 1,... , n, (1.21) j=1 j6=i alors A est une IP-matrice. 6. Montrer que si A ∈ Mn (IR) est une IP-matrice et si x ∈ IR n alors : Ax > 0 ⇒ x > 0. c’est-à-dire que {x ∈ IR n t.q. Ax > 0} ⊂ {x ∈ IR n t.q. x > 0}. 7. Montrer, en donnant un exemple, qu’une matrice A de Mn (IR) peut vérifier {x ∈ IR n t.q. Ax > 0} ⊂ {x ∈ IR n t.q. x > 0} et ne pas être une IP-matrice. 8. On suppose dans cette question que A ∈ Mn (IR) est inversible et que {x ∈ IR n t.q. Ax > 0} ⊂ {x ∈ IR n t.q. x > 0}. Montrer que A est une IP-matrice. 9. (Question plus difficile) Soit E l’espace des fonctions continues sur IR et admettant la même limite finie en +∞ et −∞. Soit L(E) l’ensemble des applications linéaires continues de E dans E. Pour f ∈ E, on dit que f > 0 (resp. f ≥ 0) si f (x) > 0 (resp. f (x) ≥ 0) pour tout x ∈ IR. Montrer qu’il existe T ∈ L(E) tel que T f ≥ 0 =⇒ f ≥ 0, et g ∈ E tel que T g > 0 et g 6> 0 (ceci démontre que le raisonnement utilisé en 2 (b) ne marche pas en dimension infinie). Exercice 9 (M-matrice). Corrigé en page 26 Dans ce qui suit, toutes les inégalités écrites sur des vecteurs ou des matrices sont à entendre au sens composante par composante. Soit A = (ai,j )i,j=1,...,n une matrice carrée d’ordre n. On dit que A est une M -matrice si A est une IP-matrice (A est inversible et A−1 ≥ 0, voir exercice 8) qui vérifie de plus (a) ai,i > 0 pour i = 1,... , n ; (b) ai,j ≤ 0 pour i, j = 1,... , n, i 6= j ; 1. Soit A uneIP-matrice  ; montrer que A est une M -matrice si et seulement si la propriété (b) est vérifiée. a b 2. Soit A = une matrice réelle d’ordre 2. Montrer que A est une M-matrice si et seulement si : c d   ad > bc, a > 0, d > 0, (1.22)  b ≤ 0, c ≤ 0.     −1 2 2 −1 3. Les matrices A = et B = sont-elles des IP-matrices ? des M-matrices ? 2 −1 −1 2 4. Soit A la matrice carrée d’ordre 3 définie par :  1  2 −1 2 A= 0 1 −1  −1 0 1 Analyse numérique I, télé-enseignement, L3 18 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES Montrer que A−1 ≥ 0 mais que A n’est pas une M –matrice. 5. Soit A une matrice carrée d’ordre n = m + p, avec m, p ∈ IN tels que m ≥ 1 et p ≥ 1, vérifiant :  ai,i ≥ 0,   ai,j ≤ 0, pour j = 1,... , n, j 6= i,    Xn pour i = 1,... , m, (1.23) ai,i + ai,j = 0     j=1  j6=i  ai,i = 1, pour i = m + 1,... , n. (1.24) ai,j = 0, pour j = 1,... , n, j 6= i, ∀i ≤ m, ∃(k` )`=1,...,Li ; k1 = i, kLi > m, et ak` ,k`+1 < 0, ∀ ` = 1,... , Li. (1.25) Soit b ∈ IR n tel que bi = 0 pour i = 1,... , m. On considère le système linéaire Au = b (1.26) 5.1 Montrer que le système (1.26) admet une et une seule solution. 5.2 Montrer que u est tel que mink=m+1,n bk ≤ ui ≤ maxk=m+1,n bk. (On pourra pour simplifier supposer que les équations sont numérotées de telle sorte que mink=m+1,n bk = bm+2 ≤ b2 ≤... ≤ bn = maxk=m+1,n bk.) 6. On considère le problème de Dirichlet suivant : −u00 = 0 sur [0, 1] (1.27a) u(0) = −1 (1.27b) u(1) = 1. (1.27c) 6.1 Calculer la solution exacte u de ce problème et vérifier qu’elle reste comprise entre -1 et 1. 6.2 Soit m > 1 et soient A et b et la matrice et le second membre du système linéaire d’ordre n = m+2 obtenu par discrétisation par différences finies avec un pas uniforme h = m 1 du problème (1.27) (en écrivant les conditions aux limites dans le système). Montrer que la solution u = (u1 ,... , un )t ∈ IR n du système Au = b vérifie −1 ≤ ui ≤ 1. Exercice 10 (Matrice du Laplacien discret 1D). Corrigé détaillé en page 23. Soit f ∈ C([0, 1]). Soit n ∈ IN? , n impair. On pose h = 1/(n + 1). Soit Kn la matrice définie par (1.10) page 12, issue d’une discrétisation par différences finies avec pas constant du problème (1.5a) page 11. Montrer que Kn est symétrique définie positive. Exercice 11 (Pas non constant). Reprendre la discrétisation vue en cours avec un pas hi = xi+1 −xi non constant, et montrer que dans ce cas,le schéma est consistant d’ordre 1 seulement. Exercice 12 (Réaction diffusion 1d.). Corrigé détaillé en page 24. On s’intéresse à la discrétisation par Différen- ces Finies du problème aux limites suivant : −u00 (x) + u(x) = f (x), x ∈]0, 1[, (1.28) u(0) = u(1) = 0. Soit n ∈ IN?. On note U = (uj )j=1,...,n une “valeur approchée” de la solution u du problème (1.28) aux points   j n+1. Donner la discrétisation par différences finies de ce problème sous la forme AU = b. j=1,...,n Exercice 13 (Discrétisation). Corrigé en page 27 Analyse numérique I, télé-enseignement, L3 19 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.2. POURQUOI ET COMMENT ? CHAPITRE 1. SYSTÈMES LINÉAIRES On considère la discrétisation à pas constant par le schéma aux différences finies symétrique à trois points du problème (1.5a) page 11, avec f ∈ C([0, 1]). Soit n ∈ IN? , n impair. On pose h = 1/(n + 1). On note u est la solution exacte, xi = ih, pour i = 1,... , n les points de discrétisation, et (ui )i=1,...,n la solution du système discrétisé (1.9). 1. Montrer que si u ∈ C 4 ([0, 1], alors la propriété (1.7) est vérifiée, c.à.d. : u(xi+1 ) + u(xi−1 ) − 2u(xi ) 00 h2 (4) − = −u (x i ) + R i avec |R i | ≤ ku k∞. h2 12 2. Montrer que si f est constante, alors max |ui − u(xi )| = 0. 1≤i≤n 3. Soit n fixé, et max |ui − u(xi )| = 0. A-t-on forcément que f est constante sur [0, 1] ? 1≤i≤n Analyse numérique I, télé-enseignement, L3 20 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.3. LES MÉTHODES DIRECTES CHAPITRE 1. SYSTÈMES LINÉAIRES 1.3 Les méthodes directes 1.3.1 Définition Définition 1.10 (Méthode directe). On appelle méthode directe de résolution de (1.1) une méthode qui donne exactement x (A et b étant connus) solution de (1.1) après un nombre fini d’opérations élémentaires (+, −, ×, /). Parmi les méthodes de résolution du système (1.1), la plus connue est la méthode de Gauss (avec pivot), encore appelée méthode d’échelonnement ou méthode LU dans sa forme matricielle. Nous rappelons la méthode de Gauss et sa réécriture matricielle qui donne la méthode LU et nous étudierons plus en détails la méthode de Choleski, qui est adaptée aux matrices symétriques. 1.3.2 Méthode de Gauss, méthode LU Soit A ∈ Mn (IR) une matrice inversible, et b ∈ IR n. On cherche à calculer x ∈ IR n tel que Ax = b. Le principe de la méthode de Gauss est de se ramener, par des opérations simples (combinaisons linéaires), à un système triangulaire équivalent, qui sera donc facile à inverser. Commençons par un exemple pour une matrice 3 × 3. Nous donnerons ensuite la méthode pour une matrice n × n. Un exemple 3 × 3 On considère le système Ax = b, avec     1 0 1 2 A= 0 2 −1  b =  1 . −1 1 −2 −2 On écrit la matrice augmentée, constituée de la matrice A et du second membre b.     1 0 1 2 à = A b =  0 2 −1 1 . −1 1 −2 −2 Gauss et opérations matricielles Allons y pour Gauss : La première ligne a un 1 en première position (en gras dans la matrice), ce coefficient est non nul, et c’est un pivot. On va pouvoir diviser toute la première ligne par ce nombre pour en soustraire un multiple à toutes les lignes d’après, dans le but de faire apparaître des 0 dans tout le bas de la colonne. La deuxième équation a déjà un 0 dessous, donc on n’a rien besoin de faire. On veut ensuite annuler le premier coefficient de la troisième ligne. On retranche donc (-1) fois la première ligne à la troisième 3 :     1 0 1 2 1 0 1 2  0 2 −1 1  `3 ←to` +` −→3 1  0 2 −1 1 −1 1 −2 −2 0 1 −1 0   1 0 0 Ceci revient à multiplier à à gauche par la matrice E1 = 0 1 0. 1 0 1 La deuxième ligne a un terme non nul en deuxième position (2) : c’est un pivot. On va maintenant annuler le deuxième terme de la troisième ligne ; pour cela, on retranche 1/2 fois la ligne 2 à la ligne 3 : 3. Bien sûr, ceci revient à ajouter la première ligne ! Il est cependant préférable de parler systématiquement de “retrancher” quitte à utiliser un coefficient négatif, car c’est ce qu’on fait conceptuellement : pour l’élimination on enlève un multiple de la ligne du pivot à la ligne courante. Analyse numérique I, télé-enseignement, L3 28 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.3. LES MÉTHODES DIRECTES CHAPITRE 1. SYSTÈMES LINÉAIRES     1 1 1 2 1 0 1 2  0 2 −1 1 `3 ←`−→ 3 −1/2`2  0 2 −1 1 . 0 1 −1 0 0 0 − 21 − 21   1 0 0 Ceci revient à multiplier la matrice précédente à gauche par la matrice E2 = 0 1 0. On a ici obtenu une 0 − 12 1 matrice sous forme triangulaire supérieure à trois pivots : on peut donc faire la remontée pour obtenir la solution du système, et on obtient (en notant xi les composantes de x) : x3 = 1 puis x2 = 1 et enfin x1 = 1. On a ainsi résolu le système linéaire. On rappelle que le fait de travailler sur la matrice augmentée est extrêmement pratique car il permet de travailler simultanément sur les coefficients du système linéaire et sur le second membre. Finalement, au moyen des opérations décrites ci-dessus, on a transformé le système linéaire Ax = b en U x = E2 E1 b, où U = E2 E1 A est une matrice triangulaire supérieure. Factorisation LU Tout va donc très bien pour ce système, mais supposons maintenant qu’on ait à résoudre 3089 systèmes, avec la même matrice A mais 3089 seconds membres b différents 4. Il serait un peu dommage de recommencer les opérations ci-dessus 3089 fois, alors qu’on peut en éviter une bonne partie. Comment faire ? L’idée est de “factoriser” la matrice A, c.à.d de l’écrire comme un produit A = LU , où L est triangulaire inférieure (lower triangular) et U triangulaire supérieure (upper triangular). On reformule alors le système Ax = b sous la forme LU x = b et on résout maintenant deux systèmes faciles à résoudre car triangulaires : Ly = b et U x = y. La factorisation LU de la matrice découle immédiatement de l’algorithme de Gauss. Voyons comment sur l’exemple précédent. 1/ On remarque que U = E2 E1 A peut aussi s’écrire A = LU , avec L = (E2 E1 )−1. 2/ On sait que (E2 E1 )−1 = (E1 )−1 (E2 )−1. 3/ Les matrices inverses E1−1 et E2−1 sont faciles à déterminer : comme E2 consiste à retrancher 1/2 fois la ligne 2 à la ligne 3, l’opération inverse consiste à ajouter 1/2 fois la ligne 2 à la ligne 3, et donc   1 0 0 E2−1 = 0 1 0. 0 21 1     1 0 0 1 0 0 Il est facile de voir que E1−1 =  0 1 0 et donc L = E1−1 E2−1 =  0 1 0. −1 0 1 −1 12 1 La matrice L est une matrice triangulaire inférieure (et c’est d’ailleurs pour cela qu’on l’appelle L, pour “lower” in English...) dont les coefficients sont particulièrement simples à trouver : les termes diagonaux sont tous égaux à un, et chaque terme non nul sous-diagonal `i,j est égal au coefficient par lequel on a multiplié la ligne pivot i avant de la retrancher à la ligne j. 4/ On a bien donc A = LU avec L triangulaire inférieure (lower triangular) et U triangulaire supérieure (upper triangular). La procédure qu’on vient d’expliquer s’appelle méthode LU pour la résolution des systèmes linéaires, et elle est d’une importance considérable dans les sciences de l’ingénieur, puisqu’elle est utilisée dans les programmes informatiques pour la résolution des systèmes linéaires. Dans l’exemple que nous avons étudié, tout se passait très bien car nous n’avons pas eu de zéro en position pivotale. Si on a un zéro en position pivotale, la factorisation peut quand même se faire, mais au prix d’une permutation. 4. Ceci est courant dans les applications. Par exemple on peut vouloir calculer la réponse d’une structure de génie civil à 3089 chargements différents. Analyse numérique I, télé-enseignement, L3 29 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.3. LES MÉTHODES DIRECTES CHAPITRE 1. SYSTÈMES LINÉAIRES (1) (1) a1,1 a1,N 0 (i+1) A(i+1) = 0 ai+1,i+1 (i+1) ai+2,i+1 (i+1) (i+1) 0 0 aN,i+1 aN,N F IGURE 1.3: Allure de la matrice de Gauss à l’étape i + 1 Le résultat général que l’on peut démontrer est que si la matrice A est inversible, alors il existe une matrice de permutation P , une matrice triangulaire inférieure L et une matrice triangulaire supérieure U telles que P A = LU : voir le théorème 1.19. Le cas général d’une matrice n × n De manière plus générale, pour une matrice A carrée d’ordre n, la méthode de Gauss s’écrit : On pose A(1) = A et b(1) = b. Pour i = 1,... , n − 1, on cherche à calculer A(i+1) et b(i+1) tels que les systèmes A(i) x = b(i) et A(i+1) x = b(i+1) soient équivalents, où A(i+1) est une matrice dont les coefficients sous-diagonaux des colonnes 1 à i sont tous nuls, voir figure 1.3 : Une fois la matrice A(n) (triangulaire supérieure) et le vecteur b(n) calculés, il sera facile de résoudre le système A(n) x = b(n). Le calcul de A(n) est l’étape de “factorisation", le calcul de b(n) l’étape de “descente", et le calcul de x l’étape de “remontée". Donnons les détails de ces trois étapes. Etape de factorisation et descente Pour passer de la matrice A(i) à la matrice A(i+1) , on va effectuer des combinaisons linéaires entre lignes qui permettront d’annuler les coefficients de la i-ème colonne situés en dessous de la ligne i (dans le but de se rapprocher d’une matrice triangulaire supérieure). Evidemment, lorsqu’on fait ceci, il faut également modifier le second membre b en conséquence. L’étape de factorisation et descente s’écrit donc : (i+1) (i) (i+1) (i) 1. Pour k ≤ i et pour j = 1,... , n, on pose ak,j = ak,j et bk = bk. (i) 2. Pour k > i, si ai,i 6= 0, on pose : (i) (i+1) (i) ak,i (i) ak,j = ak,j − (i) ai,j , pour k = j,... , n, (1.33) ai,i (i) (i+1) (i) ak,i (i) bk = bk − b. (i) i (1.34) ai,i La matrice A(i+1) est de la forme donnée sur la figure 1.3. Remarquons que le système A(i+1) x = b(i+1) est bien équivalent au système A(i) x = b(i). (i) Si la condition ai,i 6= 0 est vérifiée pour i = 1 à n, on obtient par le procédé de calcul ci-dessus un système linéaire A(n) x = b(n) équivalent au système Ax = b, avec une matrice A(n) triangulaire supérieure facile à inverser. On (i) verra un peu plus loin les techniques de pivot qui permettent de régler le cas où la condition ai,i 6= 0 n’est pas vérifiée. Analyse numérique I, télé-enseignement, L3 30 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.3. LES MÉTHODES DIRECTES CHAPITRE 1. SYSTÈMES LINÉAIRES Etape de remontée Il reste à résoudre le système A(n) x = b(n). Ceci est une étape facile. Comme A(n) est une (i) matrice inversible, on a ai,i 6= 0 pour tout i = 1,... , n, et comme A(n) est une matrice triangulaire supérieure, on peut donc calculer les composantes de x en “remontant", c’est–à–dire de la composante xn à la composante x1 : (n) bn xn = (n) , an,n   1  (i) X (n) xi = (n) b − ai,j xj  , i = n − 1,... , 1. ai,i j=i+1,n Il est important de savoir mettre sous forme algorithmique les opérations que nous venons de décrire : c’est l’étape clef avant l’écriture d’un programme informatique qui nous permettra de faire faire le boulot par l’ordinateur ! Algorithme 1.11 (Gauss sans permutation). 1. (Factorisation et descente) Pour i allant de 1 à n, on effectue les calculs suivants : (a) On ne change pas la i-ème ligne (qui est la ligne du pivot) ui,j = ai,j pour j = i,... , n, y i = bi (b) On calcule les lignes i + 1 à n de U et le second membre y en utilisant la ligne i. Pour k allant de i + 1 à n : a `k,i = ak,i i,i (si ai,i = 0, prendre la méthode avec pivot partiel) pour j allant de i + 1 à n, uk,j = ak,j − `k,i ui,j (noter que ak,i = 0) Fin pour yk = bk − `k,i yi Fin pour 2. (Remontée) On calcule x : yn xn = un,n Pour i allant de n − 1 à 1, xi = yi Pour j allant de i + 1 à n, xi = xi − ui,j xj Fin pour 1 xi = xi ui,i Fin pour Coût de la méthode de Gauss (nombre d’opérations) On peut montrer (on fera le calcul de manière détaillée pour la méthode de Choleski dans la section suivante, le calcul pour Gauss est similaire) que le nombre d’opérations nécessaires nG pour effectuer les étapes de factorisation, descente et remontée est 32 n3 + O(n2 ) ; on rappelle qu’une fonction f de IN dans IN est O(n2 ) veut dire qu’i existe un réel constant C tel que f (n) ≤ Cn2. On a donc limn→+∞ nnG3 = 23 : lorsque n est grand, le nombre d’opérations se comporte comme n3. En ce qui concerne la place mémoire, on peut très bien stocker les itérés A(i) dans la matrice A de départ, ce qu’on n’a pas voulu faire dans le calcul précédent, par souci de clarté. Analyse numérique I, télé-enseignement, L3 31 Université d’Aix-Marseille, R. Herbin, 10 septembre 2015 1.3. LES MÉTHODES DIRECTES CHAPITRE 1. SYSTÈMES LINÉAIRES Décomposition LU Si le système Ax = b doit être résolu pour plusieurs second membres b, on a déjà dit qu’on a intérêt à ne faire l’étape de factorisation (i.e. le calcul de A(n) ), qu’une seule fois, alors que les étapes de descente et remontée (i.e. le calcul de b(n) et x) seront faits pour chaque vecteur b. L’étape de factorisation peut se faire en décomposant la matrice A sous la forme LU. Supposons toujours pour l’instant que lors de l’algorithme de Gauss,