Cours Calcul Scientifique 2013

Calcul scientifique Licence de M´ Mecanique-3` e´ canique-3eme e` me ann´ annee e´ e Universit´ Universite´ d’Aix-Marseille, 2013-2014 Uwe Ehrenstein 12 septembre 2013

Table des mati` matieres e` res 1

Inter Interpol polati ation on et int´ integration e´ gration num´ nume´ rique 3 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Inte Interrpolatio tion polyn lynomiale . . . . . . . . . . . . . . . . . . . . . . 5 1.2 1.2.1 Polynˆ olynômes de Lagrange . . . . . . . . . . . . . . . . . . 8 1.3 Intégration egration numériqu eriquee : les formul formules es de Newto Newton n et Coates Coates . . . . . 9 1.4 L’erreur ’erreur dans dans les formule formuless de Newton Newton et et Coates Coates : la formule formule de Peano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.4.1 Erreurs Erreurs dans les formules formules des trapèzes ezes et de Simpson pour l’intervalle [ a, b] . . . . . . . . . . . . . . . . . . . . . . 19

2

Resolution e´ solution num´ numerique e´ rique des équations equations diff erentie e´ rentielles lles ordin ordinair aires es (EDO) (EDO) 2.1 Resultats ´ gén´ enéraux sur les EDO . . . . . . . . . . . . . . . . . . 2.1 2.1.1 Systèmes emes d’équations equations différentielle erentielless linéaires eaires a` coeffic coefficients ients constants . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. 2.1.2 2 Calc Calcul ul de l’exp l’expon onen entie tielle lle de la matr matric icee . . . . . . . . . . . 2.2 Schémas emas a` un pas pour la résolut olutio ion n d’un ’une EDO . . . . . . . . . . 2.2.1 2.2.1 Ordre Ordre d’un d’un sch´ schéma, ema, consistance, consistance, stabilit´ stabilite´ et con converge vergence nce . 2.2. 2.2.2 2 Les Les sch´ schémas de Runge-Kutta . . . . . . . . . . . . . . . .

23 27 34 36 42

Resolution e´ solution num´ numerique e´ rique directe de syst emes e` mes lin´ line´ aires 3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . LU d’un 3.2 Decomposition LU ecomposition ´ d’unee matr matric icee trid tridia iago gona nale le . . . . . . LU de matrices . . . . . . . . . . . . . . 3.3 Decomposition LU ecomposition ´ 3.3.1 Algorithme de Gauss . . . . . . . . . . . . . . . 3.3.2 Décomposition LU ecomposition LU ave avecc perm permut utat atio ions ns des des ligne ligness

49 49 52 55 55 60

3

4

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

Norm Normes es de de mat matri rice ces, s, m´ methodes e´ thodes it´ iteratives e´ ratives de r esolution e´ solution de syst emes e` mes lin´ line´ aires 4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Normes de matrices . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Applicatio Application n : conditionn conditionnemen ementt d’un système eme linéaire . . . 1

21 21

67 67 68 73

` TABLE DES MATI ERES

4.3 4.4

Conditio ition ns de convergence . . . . . . . . . . . . . . . . . . . . . 74 Methode ´ de Jacobi et de Gauss-Seidel, méthode ethode de relax relaxati ation on . . . 76 4.4. 4.4.1 1 Quel Quelqu ques es résultats esultats de convergence de méthodes ethodes itérat e rativ ives es . 78

2

Chapitre 1 Interpolation et int´ integration e´ gration num´ numerique e´ rique 1.1 1.1

Moti otivati vation on

En gén´ enéral eral il n’est guère ere possible de trouver la valeur valeur exacte d’une int egrale e´ grale

Z a

b

f ( x) dx

f , ou sauf dans le cas où on connaˆıt ıt explicitement une primitive de la l a fonction f , f permet par exemple une intégration alors lorsque l’expression de la fonction f permet egration par parties ou un changement de variable approprié. e. L’idée ee est alors d’approcher l’intégrale egrale par une somme. La méthode ethode la plus rudimentaire rudimentaire est d’utiliser d’utiliser les sommes de Riemann. Si on introduit les n + 1 points j x j = a + (b n

− a),

j = 0, 1,

· · ·, n

alors une somme de Riemann, appelée ee encore formule des rectangles à gauche, est donnée ee par 1 n−1 S g = ∑ f ( x j ). n j=0 L’expression S g correspond a` la somme des aires de tous les rectangles de base x j , x j+1 ] (de longueur h longueur h = 1/n) et de hauteur f ( x j ), qui est la valeur de f ( x j ) “à [ x gauche” du petit intervalle [ x j , x j+1 ]. Dans la formule des rectangles “à droite” on prend la valeur f ( x j+1 ) pour obtenir 1 n−1 S d d = ∑ f ( x j+1). n j=0 3

Interpolation et int´ egration num´ erique

Ces deux approximations sont illustrées ees sur la figure 1.1 par les aires hachurées. ees. On sait que, pour des fonctions f continues, ces ces sommes tendent précis´ ecisément ement vers l’intégrale, egrale, lorsque n lorsque n tend tend vers l’infini, c.-à-d. a-d. lorsqu’on opère ere une sous-division de plus en plus fine de l’intervalle d’intégration. egration. f (x)

S g

a

b

x

b

x

f (x)

S d

a

F IG . 1.1 – Formule des rectangles rectangles à gauche (haut) et des rectangles à droite (bas).

Pour construire ces sommes, on peut dire que l’on approche la fonction dans l’intervalle [ x j , x j+1 ] par une constante, égale egale a` f ( x j ) pour la somme S somme S g et egale e´ gale x j , x j+1 ] la a` f ( x j+1 ) pour la somme S d d . On peut donc dire que sur l’intervalle [ x fonction est approchée ee par un polynôme ome de degré zéro, ero, c.c.-à-d. a` -d. une constante. 4

Interpolation polynomiale

L’idée ee de l’interpolation polynomiale est précis´ ecisément ement d’approcher une fonction tion sur sur des des inte interv rval alle less donn´ données ees par par des des polynˆ polynômes o mes de degr´ degrés es plus plus ou moins moins elev´ e´ levés, es, selon le nombre de points que l’on consid considère e` re dans l’intervalle.

1.2

Interp Interpolat olation ion polynomi polynomiale ale

Tout d’abord il faut préciser eciser ce que l’on entend par interpolation. interpolation. Il s’agit de construire une fonction ayant des valeurs données ees en des points donnés. es. Plus précis´ ecisément, ement, soient n soient n + 1 données ees

( x0 , y0 ), ( x1 , y1 ),

·· · · · , ( x , y ), n

n

(1.1)

ou` x j , j = 0, , n, d esignent e´ signent par exemple les abscisses et y j , j = 0, , n les ordonnées ees pouvant etre eˆ tre les valeurs d’une fonction en x j , c.-à-d. a-d. y j = f ( x j ). Lors x j on ne conna qu’on qu’on rel` relève eve par exemp exemple le des “mesur “mesures” es” aux points points x connaˆˆıt ıt par précis´ ecisément ement la fonction f associ´ associée ee a` cette mesure, mais seulement la valeur de cette fonction (la mesure) aux points x points x j . On suppose que les abscisses sont distinctes, c.-à-d. a-d. xi = x j si i = j et on cherche donc une fonction g telle g telle que

···

···



g( xi ) = yi .



(1.2)

Evidem Evidemmen ment, t, on cherc cherche he la foncti fonction on g sous sous une une cert certai aine ne form formee et dans dans une une proc´ procédure edure d’interpolation on ecrit g e´ crit g comme comme une combinaison lin´ lin eaire e´ aire de n de n + 1 fonctions h0 ( x), h1 ( x), données. ees. On ecrit e´ crit ainsi

· · · , h ( x) n

n

g( x) =

∑ c j h j ( x)

(1.3)

j=0

et les conditions d’interpolation deviennent deviennent h0 ( xi )c0 + h1 ( xi )c1 +

· · · + h − ( x )c − + h ( x )c = y , n 1

i

n 1

n

i

n

i

i = 0, 2,

· · · , n. (1.4)

Introduisant le vecteur des inconnues de n de n + 1 composantes T

d = = ( c0 , c1 ,

···,c )

z = ( y0 , y1 ,

· · · , y )

n

et le vecteur solution n

T

on peut encore ecrire e´ crire le système eme sous forme matricielle pour une matrice A de coefficient a coefficient a i j Ad = = z, avec

ai j = h j−1 ( xi−1 ), i = 1, 5

· · · , n + 1, j = 1, · · · , n + 1.

(1.5)


d et Si la matrice A matrice A est est inversible, alors le système eme ci-dessus admet une solution d et une seule en fonction du vecteur z vecteur z.. L’interpolation polynomiale consiste a` prendre comme fonctions h j ( x) des monômes omes h j ( x) = x j , j = 0, , n, (1.6)

···

n que l’on note pn et c.-à-d. a-d. on cherche g cherche g sous sous la forme d’un polynôme ome de degré n que n

∑ c j x j .

pn ( x) =

(1.7)

j=0

A de (1.5) sont Dans ce cas les coefficients de la matrice A de

−

j 1

ai j = xi−1 , et la matrice A matrice A qui qui en résulte esulte est appelée ee matrice de Vandermonde de Vandermonde

A =

  

1 1 .. .

x20 x21 .. .

x0 x1 .. .

1 xn−1 x2n−1 x2n 1 xn

   − 

··· ···

xn0 xn1 .. .

··· ···

xnn 1 xnn

.

(1.8)

On peut montrer montrer (voir des ouvrages ouvrages d’algèbre ebre linéaire) eaire) que le d eterminant e´ terminant de cette matrice est det( A) = ∏ ( x j xi ) 0 i< j n

≤ ≤

−

et et donc ce déterminant eterminant est non nul si x i = x j quand i quand i = j. Le système eme (1.5) est = (c0 , , cn )T quel donc inver inversible sible,, c’est-` c’est-à-dire a-dire qu’il existe existe une solution solution unique d = que soit z soit z = ( y0 , , yn )T . On peut donc enoncer e´ noncer le théor` eorème eme suivant.





···

···

Th´ Theor` e´ oreme e` me 1 Soient ( (n + 1) points xi , i = 0, , n, distincts deux a` deux et ( (n + 1) valeurs y0 , y1 , , yn . Alors il existe un unique polyn ome degr e´ n, not e´ ˆ de degr ´ ´ pn , tel que pn ( xi ) = yi , i = 0, , n.

···

···

≤

···

Supposons maintenant que les valeurs valeurs y j correspondent a` la valeur valeur d’une d’une foncfonction f aux aux points x points x j et soit l’unique polynôme ome de degré n tel que

≤

f ( x j ) = p n ( x j ), j = 0,

· · · , n. (1.9) On dira que p interpole f aux aux points x points x , j = 0, · · · , n. Par conséquent, equent, f et p co¨ıncident ınci dent aux a ux point po intss x , appelés es points d’interpolation ; mais qu’en est-il est-il en un point x point x = e´ noncer le théor` eorème eme suivant.  x ? On peut enoncer n

j

n

j

j

6

Interpolation polynomiale

Th´ Theor` e´ oreme e` me 2 Supposons que les points x j , j = 0, , n sont dans un intervalle [a, b] et que f est ( ( n + 1) fois continument d erivable dans l’intervalle [ a, b]. Soit ˆ ´ ´ x [a, b] et on introduit la fonction

···

∈

φ( x) = ( x x0 )( x x1 )

− · · · ( x − x ).

−

n

(1.10)

On note Rn ( x) = f ( x)

− p ( x)

n

(1.11)

ξ x dans le plus petit intervalle qui l’erreur d’interpolation. Alors il existe un point ξ contient x, x0 , xn (donc ξ x [min( x, x0 , , xn ), max( x, x0 , , xn )]) tel que

···

∈

···

Rn ( x) =

···

φ( x) (n+1) f (ξ x) (n + 1)!

(1.12)

ξ x . avec f (n+1) (ξ x ) la d eriv´ eriv eme de f au point ξ ´ ´ ee ´ n + 1 ème La formule (1.12) est evidente e´ vidente si x si x = x j , car dans ce cas R n ( x j ) = φ ( x j ) = 0. Pour démontrer emontrer le résultat esultat lorsque x lorsque x = x j , on introduit une fonction



F (t ) = Rn (t )φ( x) Rn( x)φ(t ).

− La fonction s’annule aux points x points x , j = 0, · · · n, car R car R ( x ) = 0 et φ ( x ) = 0, mais j

n

j

j

eegalement ´ galement par construction au point x. Donc, F (t ) possède ede (n + 2) zéros. eros. Or, d’après es le théor` eorème eme de Rolle, si une fonction dérivable erivable s’annule en deux points, il y au moins un point entre ces deux zéros eros o` ou` la dériv´ erivée ee de la fonction s’annule. Ici F (t ) a n a n + 2 zéros, eros, donc il y a au moins n moins n + 1 points ou` la dériv´ erivée F ee F ′ s’annule. On peut ensuite appliquer le théor` eorème eme de Rolle a` F ′ , ensuite a` F ′′ etc. On en déduit eduit qu’il existe au moins un point ξ x tel que la dériv´ erivée n ee n + 1 ème eme de F de F (t ) (n+1)

F (n+1) (t ) = Rn

(t )φ( x) Rn ( x)φ(n+1) (t )

−

(n+1)

s’annule. Mais R Mais R n (t ) = f (t ) pn(t ) et R et R n erivée n ee n + 1 (t ) = f (n+1) (t ), car la dériv´ eme e` me de p de p n (t ) est identiquement egale e´ gale a` zéro). ero). Par ailleurs, il est facile de constater n+1) ( que φ (t ) = (n + 1)! et l’expression de l’erreur d’interpolation (1.12) s’ensuit.

−

Evidemment, sauf dans des cas particuliers le point ξ x , qui dépend epend de x de x pour pour des points d’interpolation x d’interpolation x j , j = 0, , n, donnés, es, n’est pas connu explicitement. De (1.12) on peut par exemple déduire eduire la majoration

···

....( x − x )| | R R ( x)| ≤ C + |( x −(n x+)1....( )! n

n 1

0

n

avec 7

C n+1 = max f (n+1) ( x) (1.13) x [a,b]

∈

|

|


....( x xn ) = φ( x) . La majoration de l’erreur est donc fonction de ( x x0 )....( On peut essayer essayer de trouver trouver une majoration majoration de cette quantit´ quantité. e. On suppose avoir ordonné les points dans l’ordre croissant

| −

x0 <

· · · < x < x + < · · · < x i 1

i

− | |

|

n

et que h que h soit soit la distance maximale entre deux points successifs. On suppose que x est x est tel que x h2 /4. que x i < x < x i+1 alors on peut affirmer que ( x xi )( x xi+1 ) x xi−k (k + 1)h, k = 1, , i et x x xi+k kh, k = 2, n i. On peut Ensuite, x en déduire eduire la majoration 2 n−1 h φ( x) n!h 4 et substituant cette expression dans (1.13) on trouve la majoration

| − |≤

| − | − |≤

···

|

| R R ( x)| ≤ n

− |≤ ··· −

|≤

C n+1 hn+1 avec 4(n + 1)

C n+1 = max f (n+1) ( x) x [a,b]

∈

|

|

(1.14)

avec h avec h la la plus grande distance entre deux points d’interpolation voisins.

1.2. 1.2.1 1

ômes Polyn olyn ˆ omes de Lagrange

Une fac f aç on commode com mode de d e d´ determiner e´ terminer le polynôme ome d’interpolation pn ( x) qui interpole une fonction f ( x) aux points distincts x0 , x1 , x2 , , xn est d’utiliser les polynômes omes de Lagrange.

···

Definition e´ finition 1 Soient donn´ 1 points distincts x 0 , x1 , , xn ; les polyn omes donnes ´ n + 1 ˆ de Lagrange L 0 , L1 , , Ln associ´ associes ´ a` ces points sont des polyn omes ˆ de degr ´ degr e´ n d efinis efin d e fac f aç on a` ce que ´ ´ is de

···

···

L j ( xk ) = pour j = 0, 1,



1 si 0 si

j = k j = k

k = = 0, 1,



· · · , n,

(1.15)

· · · , n.

Soit donc L j ( x) qui par définition efinition s’annule en n points x points x k , k = j et il s’écrit ecrit par conséquent equent



n

L j ( x) = a ∏ ( x xk ) k =0

−

k = j



a . On en déduit et la condition L condition L j ( x j ) = 1 fournit la constante a. eduit que le j eme e` me polynôme ome de Lagrange Lagrange s’´ s’ ecrit e´ crit n

L j ( x) =

−

x xk

∏ x j − xk , k =0

k = j



8

j = 0,

· · · n.

(1.16)

Int´ egration num´ erique : les formules de Newton et Coates

On peut alors aisément ement construire l’unique polynôme ome d’interpolation pn ( x) de degré n tel que pn ( xk ) = f ( xk ), k = = 0, , n.

≤

···

En effet, il peut s’écrire ecrire sous la forme n

pn ( x) =

∑ f ( x j ) L j ( x).

(1.17)

j=0

En effet, l’expression l’expression ci-dessus est bien un polyn ome o ˆ me de degré

≤ n et

n

pn ( xk ) =

∑ f ( x j ) L j ( xk ) = f ( xk ),

j=0

d’après es la définition efinition (1.15) des polynômes omes de Lagrange. C’est précis´ ecisément ement l’interpolation polynomiale qui permet de construire des formules d’intégration. egration.

1 .3

Integration e´ gration num´ numerique e´ rique : les formules formules de Newton et Coates

On suppose donné un intervalle [ c, d ] et on cherche a` evaluer e´ valuer I = =

Z

d

f ( x)dx

c

f . L’idée f par un popour une fonction (continue) f . ee est d’approcher la fonction f par lynôme ome de degré l qui interpole f en des points discrets dans l’intervalle l’intervalle [ c, d ]. Soient donc une sous-division de l de l + 1 points de l’intervalle, c.-à-d. a-d.

≤

x j = c + jh j h, j = 0,

· · ·, l

h =

et

d

− c. l

D’après es l’expression (1.17) le polynôme ome pl qui interpole f en f en ces points peut s’écrire ecrire a` l’aide des polynômes omes de Lagrange et l

pl ( x) =

∑ f ( x j ) L j ( x).

j=0

Une formule d’intégration egration numérique erique est obtenue par la somme

Z c

d

pl ( x) =

l

∑ f ( x j )

j=0

9

Z c

d

L j ( x)dx

(1.18)


et il faut alors evaluer e´ valuer les int´ integrales e´ grales des polyn polynômes oˆ mes de Lagrange. Faisons le changement de variable l

L j (c + ht ) = φ j,l (t ) = ∏

x = c + ht , donc

k =0

t j

k = j



− k − k

(1.19)

(on ecrit e´ crit φ j,l (t ) car ces fonctions dépendent ependent bien sûr ur de l ). On peut alors ecrire e´ crire hdt ) (etant e´ tant donné que d que dxx = hdt )

Z c

d

L j ( x)dx = hα j,l avec

α j,l =

Z

l

0

φ j,l (t )dt .

(1.20)

Donc, la formule d’intégration egration s’écrit ecrit

Z

l

d

pl ( x) = h ∑ α j,l f ( x j ).

c

(1.21)

j=0

Exemples : 1. Pour l Pour l = 1, il y a dans ce cas 2 points dans l’intervalle,

Z

α0,1 =

α1,1 =

Z

Z t − 1 1 (t )dt = = dt = = 1

φ0,1

0

et

Par conséquent equent

1

Z

0

1

0

φ1,1 (t )dt = =

d

p1 ( x)dx =

( 1)

Z 0

−

1

t dt = =

2

1 . 2

h

( f ( x0 ) + f ( x1 )) 2 2. Pour l Pour l = 2, donc avec 3 points dans l’intervalle, c

α0,2 =

Z

2

0,2

0

2

α1,2 =

et

(1.22)

Z (t − 1)(t − 2) 1 φ (t )dt = = dt = = , 2 3 Z Z t (t − 2) 4

2

0

2

φ1,2 (t )dt = =

= − dt = 3 Z Z t (t − 1) 1 α , = φ , (t )dt = = dt = = . 0

0

2

22

0

( 1)

2

22

0

2

3

Dans ce cas, on obtient la formule

Z c

d

h p2 ( x)dx = ( f ( x0 ) + 4 f ( x1 ) + f ( x2 )) . 3 10

(1.23)

Int´ egration num´ erique : les formules de Newton et Coates

Considérons erons maintenant un intervalle [ a, b] et une fonction f ( x) continue sur cet intervalle. On sous-divise l’intervalle en N en N + 1 points xi = a + ih, i = 0, 1,

· · · , N , avec

h =

−

b a . N

L’idée ee est de considérer erer des sous-intervalles a` l’intérieur erieur de [ a, b] avec l avec l + 1 points > et d’interpoler sur ces sous-intervalles f par par des polynômes omes de degré l , pour N pour N > l (et en gén´ N grand devant l ). Plus précis´ N est un enéral eral N grand ecisément, ement, supposons que N est = lM l M et multiple de l de l,, c’est-à-dire N a-dire N = et on définit efinit les M les M sous-intervalles sous-intervalles

[ x xil , x(i+1)l ], i = 0,

· · · , M − 1

(1.24)

dont chacun contient l contient l + 1 points. Ces intervalles intervalles jouent le rôle ole de l’interv l’intervalle alle [c, d ] ci-dessus : interpolant f sur sur cet intervalle par le polynôme ome d’interpolation pi,l ( x) de degré l, e´ crire par (1.21) pour les points x points x il , xil +1 , , x(i+1)l l , on peut ecrire

Z

···

x(i+1)l

xil

l

pi,l ( x)dx = h ∑ α j,l f ( xil + j ).

(1.25)

j=0

Evidemment, ces sommes sont des approximations approximations de la vraie int egrale e´ grale

Z

x(i+1)l

xil

f ( x)dx

et l’analyse de l’erreur fait l’objet du paragraphe suivant. Raccordant toutes ces formules on obtient une formule d’int d’intégration e´ gration qui est une approximation de

Z

b

a

f ( x)dx .

On note cette formule d’intégration I egration I N ,l ( f ) et elle fait intervenir N intervenir N + 1 points, avec N = = lM l M . Cette formule est appelée ee de Newton Newton et Coates ; on l’obtient en sommant les expressions (1.25) et donc

−

M 1

I N ,l ( f ) = h



l

∑ ∑ α j,l f ( xil + j )

i=0

j=0



,

N = = lM l M .

(1.26)

Exemples : 1. Prenons l = 1 dans la formule ci-dessus : alors il a et´ e´ té montré plus haut que α0,1 = α1,1 = 1/2 et on obtient la formule bien connue des trapèzes ezes

−

N 1

I N ,1 ( f ) = h

1 ( f ( xi + f ( xi+1 )) 2 i=0

∑

11


que l’on peut encore écrire ecrire N −1 h I N ,1 ( f ) = ( f ( x0 ) + f ( x N )) + h ∑ f ( xi ). 2 i=1

(1.27)

Bien Bien sur, uˆ r, le nom vient du fait que sur chaque intervalle [ x xi, xi+1 ] on approche la fonction par un polynôme ome de degré 1, donc une droite, et l’aire obtenue est celle du trapèze eze qui en résulte esulte (cf. figure 1.2).

f (x)

a

b

x

FIG . 1.2 – Schéma ema illustrant la formule des trapèzes. ezes.

2. Pour l Pour l = 2 nous avons avons montre montre plus haut que α0,2 = 1/3, α1,2 = 4/3 et α2,2 = 1/3 et la formule correspondante correspondante s’écrit ecrit

−

M 1

I N ,2 ( f ) = h

1 ( f ( x2i ) + 4 f ( x2i+1 ) + f ( x2i+2 )) , 3 i=0

∑

N = = 2 M ,

ou encore h 4h M −1 2h M −1 I N ,2 ( f ) = ( f ( x0 ) + f ( x N ))+ )) + = 2 M . (1.28) ∑ f ( x2i+1) + 3 ∑ f ( x2i ), N = 3 3 i=0 i=1 Cette formule est appelée ee la formule la formule de Simpson. Simpson . 12

L’erreur ’erreur dans les formules formules de Newton Newton et Coates : la formule de Peano

1.4

L’err erreur dans les formu rmules les de Newt Newton on et Coate oatess : la formule de Peano

Pour evaluer e´ valuer l’erreur que l’on fait en approchant l’int egrale e´ grale par une somme du type (1.26), il convient de considérer erer d’abord l’intervalle [ xil , x(i+1)l ] avec l avec l + 1 points que l’on note pour pour simplifier a` nouveau [c, d ]. Sur de tels intervalles on note egrale et la formule d’intégration egration et R( f ) l’erreur entre l’intégrale R( f ) =

Z

l

d

c

f ( x)dx

− h ∑ α f ( x ) j

j

(1.29)

j=0

ou` pour simplifier on omet l’indice l l’indice l pour pour les coefficient α j de la formule. Tout d’abord il faut remarquer que si la fonction f elle-mˆ elle-même eme est un polynˆ polynome oˆ me de degré inférieur erieur ou egal e´ gal a` l , alors R( f ) = 0. En effet, dans ce cas f est est identique a` son polyn polynome oˆ me d’interpolation pl . En effet, l , la fonction r si f est est un polynôme ome de degré inférieur erieur ou egal e´ gal a` l, fonction r ( x) = f ( x) pl ( x) (qui est un polynôme ome de degré inférieur erieur ou egal e´ gal a` l ) s’annule en les l les l + 1 points x j , j = 0, , l ; or, un polynˆ l a au plus l polynome oˆ me non nul de degré l a plus l z zéros e´ ros réels eels et on en déduit eduit que r que r ( x) = 0. Soit maintenant f quelconque et et afin de trouver une formule gén´ enérale erale de l’erreur d’intégration, egration, on suppose que f est l + 1 fois continûment ument dérivable erivable dans l +1 [c, d ]. On ecrira e´ crira f C [c, d ]. On rappelle la formule de Taylor avec reste sous forme intégrale egrale

−

···

∈ ∈

f ( x) = (on ecrit e´ crit

( x − c)l 1 ′ (l ) f (c) + f (c)( x − c) + · · · + f (c) +

f ( j)

l!

l!

Z

x

c

f (l +1) (t )( )( x

l

− t ) dt

(1.30)

f ). On introduit la fonction pour la dériv´ erivée ee j ème eme de f ). ql ( x, t ) =



( x

l

− t ) , 0,

ce qui permet d’écrire ecrire

Z c

x

(l +1)

f

(t )( )( x

l

− t ) dt = =

Z c

d

≥

si x si x t si x si x < t

(1.31)

f (l +1) (t )ql ( x, t )dt .

(1.32)

Ecrivant le reste sous forme intégrale egrale avec la fonction ql ( x, t ) fait que la borne supérieure erieure de l’intégrale egrale est d est d et et non pas x pas x.. La formule de Taylor devient donc 1 f ( x) = p ( x) + l!

Z

d

c

13

f (l +1) (t )ql ( x, t )dt

(1.33)


avec

( x − c)l ′ (l ) f (c) + f (c)( x − c) + · · · + f (c)

p( x) =

l! qui est un polynôme ome de degré inférieur erieur egal e´ gal a` l . Donc, d’après es ce qui préc` ecède, ede, R( p) =

Z

l

d

p( x)dx

c

− h ∑ α p( x ) = 0. j

j

j=0

Il s’ensuit que 1 R( f ) = l!

Z Z d

c

d

(l +1)

f

c



(t )ql ( x, t )dt dx

−

1 l h ∑ α j l ! j=0

Z

d

c

f (l +1) (t )ql ( x j , t )dt .

Or, un peut intervertir l’ordre d’intégration egration dans l’int´ l’i ntégrale egrale double et l’int´ l’integrale e´grale d’une somme etant e´ tant la somme des intégrales, egrales, on obtient le théor` eorème eme

Th´ Theor` e´ oreme e` me 3 Soit Soit f C l +1 [c, d ] ; alors l’erreur (1.29) commise en approchant l’int egrale par la formule d’int ´ d’int egration de Newton et Coates est ´ ´ ´

∈ ∈

R( f ) =

Z

d

c

f (l +1) (t )K l (t )dt ,

(1.34)

K l (t ) ´ etant la fonction dite de Peano dont l’expression est K l (t ) =

1 l!

Z

l

d

c

ql ( x, t )dx

− h ∑ α q ( x , t ) j l

j

j=0



(1.35)

avec ql ( x, t ) la fonction donn ee ´ par (1.31).

Exemples de fonctions de Peano : 1. Considérons erons d’abord le cas l cas l = 1 et la formule des trapèzes ezes (1.22). Dans ce cas [c, d ] contient contient deux points points et d et d c = h. Il convient alors de considérer erer un intervalle intervalle type de longueur h longueur h,, par exemple [0, h]. Il suffit de déterminer eterminer la fonction de Peano car cette cette foncti fonction on dans dans tout tout autre autre interv intervall allee de longue longueur ur h peut etre eˆ tre K 1 (t ) pour [0, h] car h peut obtenue par translation de la variable. Pour l = 1 la fonction q 1 ( x, t ) est d’après es (1.31) si x t ( x t )1, si x q1 ( x, t ) = 0, si x si x < t

−



Alors pour t pour t

∈ [0, h] on aura

Z 0

h

q1 ( x, t )dx =

−

Z t

≥

h

( x 14

− t )dx =

(h

− t ) 2

2


t de la deuxième (la borne inférieure erieure t de eme intégrale egrale ci-dessus etant e´ tant due au fait que q1 ( x, t ) = 0 si x si x < t ). ). Si t Si t [0, h], on déduit eduit de l’expres l ’expression sion de q de q 1 que q que q1 (0, t ) = 0 et q et q 1 (h, t ) = h t . Les coefficients de la formule sont α0 = α1 = 1/2 et on obtient finalement par (1.35)

∈

−

K 1 (t ) =

(h

2

− t ) − h (h − t ) = (h − t )()(−t ) , t ∈ [0, h]. 2 2 2

(1.36)

On observe que cette fonction, représent´ esentée ee sur la figure 1.3, est de signe constant négatif egatif sur l’intervalle. l ’intervalle. 0

−0.02

−0.04

−0.06

−0.08

−0.1

−0.12

−0.14 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

FIG . 1.3 – Fonction de Peano K 1 (t ) associ associee e´ e a` la formule formule des trap` trapèzes ezes (tracée ee pour h = 1).

2. Le calcul calcul pour pour l = 2 et pour pour la form formul ulee de Simp Simpso son, n, est est plus plus comp compliq liqu´ ué ; pren prenon onss dans l’expression (1.23) l’intervalle de longueur 2 h centré en 0, a` savoir [ h, h]. La formule formu le est construite construit e de faç on a` ce que si p( x) est un polynôme ome de degré au plus 2, alors h h p( x)dx ( p( h) + 4 p(0) + p(h)) = 0. 3 −h

−

Z

−

−

p avec p de degré 2 au plus. Soit s( x) un polynôme ome de degré 3, alors s = ax 3 + p avec Bien sûr, ur, 3h (s( h) + 4s(0) + s(h)) = 3h ( p( h) + 4 p(0) + p(h)) et on aura

−

Or,

Z

Z

−

h

−h

p( x)dx

h

−h

s( x)dx =

− h3 (s(−h) + 4s(0) + s(h)) = 0.

Z

h

−h

p( x)dx 15

car

Z

h

−h

ax3 dx = 0


Z

et finalement

h

−h

s( x)dx

− h3 (s(−h) + 4s(0) + s(h)) = 0

pour tout polyn polynome s oˆ me s de de degré inférieur erieur ou egal e´ gal a` 3. D’une manière ere gén´ enérale, erale, si une formule d’intégration egration de Newton et Coates est exacte pour des polynômes omes de degré l , avec l entier pair, alors elle est exacte pour des polynômes omes de degré l + 1. La démonstration emonstration dans le cas gén´ enéral eral se fait aisément ement en s’inspirant de la démonstra emonstration tion pour le cas l cas l = 2 ci-dessus. ci-dessus. Le fait d’avoir d’avoir démontr´ emontré le r´ esultat pour l’intervalle [ h, h] n’enlève eve rien a` la gén´ enéralit´ eralité. e. En effet, soit x soit x [c, d ] avec les 3 points c points c,, c + h et d et d = variable y = c + 2h, alors par la translation y = x c h la variable y e´ grale car d car dxx = dy [ h, h]. Or, une telle translation ne change ni la nature de l’int egrale d y, ni les degrés es des polynômes. omes. Donc, pour la formule de Simpson avec l = 2 la fonction de Peano peut être etre prise avec l avec l + 1 = 3 et d’après es (1.35)

−

∈ − −

−

1 K 3 (t ) = 3! avec

Soit t Soit t

Z

h

∈ [−h, h], alors

Z

−

q3 ( x, t ) =



h

−h

tandis que

q3 ( x, t )dx

−h

q3 ( x, t )dx =

q3 (0, t ) = q3 (h, t ) =



h (q3 ( h, t ) + 4q3(0, t ) + q3(h, t )) )) 3

−

( x

3

− t ) ,

≥

si x si x t si x si x < t

0,

Z

h

( x

t

3

− t ) dx =

q3 ( h, t ) = 0, si

−

∈

 − − ≤ −≤≤ h

t

t 3 si h 0 si 0

(h

− t )

4

4

h

< 0 t < t h

≤ ≤ (h − t ) , si − h ≤ t ≤ h. 3

ecède ede ≤ t ≤ h, d’aprèses ce qui préc` 1 (h − t ) h(h − t ) − K (t ) = 6 4 3

(1.37)

Pour 0



4

3



− 172 (h − t ) (h + 3t ), 0 ≤ t ≤ h. Le calcul, qu’on ne détaille etaille pas ici, pour −h ≤ t ≤ 0 permet de se convaincre que K (−t ) = K (t ), c’est-à-dire a-dire la fonction K fonction K (t ), représent´ esentée ee sur la figure 1.4 pour h = 1, est egalement e´ galement de signe constant constant et n´ n egative e´ gative pour −h ≤ t ≤ h. 3

3

=

3

3

3

Lorsque la fonction de Peano est de signe constant on peut démontrer démontrer le résultat esultat suivant, suivant, a` partir de la formule (1.34). 16


0

−0.002

−0.004

−0.006

−0.008

−0.01

−0.012

−0.014 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

K 3 (t ) associée FIG . 1.4 – Fonction de Peano Peano K ee a` la formule de Simpson (tracée ee pour h = 1).

Proposition Proposition 1 Si K l (t ) est de signe constant dans [c, d ] , , et si f C l +1 [c, d ] , , alors il existe un point ξ ]c, d [ tel que l’erreur d’int ´ d’int egration donn´ donnee ´ ´ par (1.34) peut s’´ s’ecrire ´ 1 R( f ) = f (l +1) (ξ) R(gl ) (1.38) (l + 1)!

∈ ∈

∈

avec R(gl ) l’erreur d’int egration pour la fonction g l ( x) = xl +1 . ´ ´ Pour la preuve, on fait appel a` une variante du théor` eorème eme de la moyenne. Soit donc K donc K l (t ) de signe constant, par exemple K exemple K l (t ) 0, c 0, c t d (le (le cas K cas K l (t ) 0 se traite de manière ere analogue). On peut alors écrire ecrire a` partir de (1.34) que (l +1)

min f

≤≤

c t d

(t )

Z c

≥

d

K l (t )dt

≤ ≤ ≤

(l +1)

≤ R( f ) ≤ max f ≤≤ c t d

Z

≤

d

(t )( )(t )

c

K l (t )dt .

Or, d’après es hypothèse ese f (l +1) (t ) est continue car f C l +1 [c, d ] et par conséquent, equent, par le théor` eorème eme de la valeur intermédiaire, ediaire, il existe ξ ]c, d [ tel que (l +1)

R( f ) = f

(ξ)

Z

∈ ∈

∈

d

c

K l (t )dt .

Pour prouver la formule (1.38), il reste l’intégrale egrale de la fonction de Peano à calculer. culer. Pour ce faire, il s’avère ere commode de calculer l’erreur pour la fonction parl +1 ticulière ere f ( x) = g l ( x) = x : en effet, la dériv´ erivée ee (l + 1) eme e` me de cette fonction etant e´ tant ( l + 1)!, on aura R(gl ) = ( l + 1)! 17

Z c

d

K l (t )dt


et le résultat esultat (1.38) s’ensuit. A nouveau nous traitons le cas particulier pour lesquels les fonctions de Peano ont et´ e´ té calculées ees plus haut. 1. Nous avons vu que pour l = 1 et la formule des trapèzes ezes la fonction K fonction K 1 (t ) est de signe constant. On prend à nouveau l’intervalle [ 0, h] et d’après es (1.38), si egration f C 2 [0, h], l’erreur d’intégration

∈ ∈

R( f ) =

Z

h

f ( x)dx

0

peut s’écrire ecrire

f ′′ (ξ)

R( f ) =

2

avec g avec g 1 ( x) = x2 . Or,

Z

R(g1 ) =

h

− h2 ( f (0) + f (h))

R(g1 ),

2

x dx

0

h

ξ ]0, h[

∈

h3

2

− 2h = − 6

et par conséquent equent h3

− 12 f ′′(ξ),

R( f ) =

ξ ]0, h[.

∈

(1.39)

2. On suppose maintenant l = 2 et on considère ere la formule de Simpson (1.23) : nous avons montré plus haut qu’alors la formule est aussi exacte pour des polynômes omes de degré 3. Soit donc f C 4 [ h, h]. La fonction fonction K 3 (t ) etant e´ tant de signe constant pour la formule de Simpson,

∈ ∈

R( f ) =

−

f (4) (ξ) 24

R(g3 ),

ξ ]

∈ − h, h[,

etant l’erreur d’intégration egration R( f ) étant R( f ) =

Z

h

−h

f ( x)dx

Pour f = g3 = x4 on obtient R(g3 ) =

Z

h

− h3 ( f (−h) + 4 f (0) + f (h)). )). 4

x dx

−h

h

4

− 3 (2h ) = −

4h5 15

et par conséquent equent R( f ) =

−

h5 (4) f (ξ), 90 18

ξ ]

∈ − h, h[.

(1.40)

Erreurs dans les formules des trap` ezes et de Simpson pour l’intervalle [ [a, b]

1.4. 1.4.1 1

Erre Erreur urss dans dans les form formul ules es des trap trapezes e` zes et de Simpson pour l’intervalle [ a, b]

Comme il a et´ e´ té dit plus haut, l’intérˆ erêt et d’une formule d’intégration egration est de pouvoir approcher l’intégrale egrale sur un intervalle [a, b] avec un grand nombre de points ezes (1.27), alors pour chaque sous N = = lM l M . Prenons d’abord la formule des trapèzes intervalle [ xi , xi+1 ] on aura bien entendu une erreur de la forme (1.39) et Ri ( f ) =

Z

xi+1

xi

f ( x)dx

−

h ( f ( xi ) + f ( xi+1 )) = 2

−

h3 ′′ f (ξi ), 12

ξi ] x xi , xi+1 [.

∈

(1.41) Dans la formule des trapèzes ezes (1.27), les intéegrales grales sur les sous-intervalles sont sommées ees et par conséquent equent l’erreur pour [ a, b] est R N ,1 ( f ) =

Z

f ( x)dx I N ,1 =

−

a

−

N 1

b

∑ Ri ( f ).

i=0

Or, par la formule de la moyenne on pourra écrire ecrire N 1

−

h3 N −1

i=0

i=0

∑ Ri( f ) = − 12

pour une valeur η

h3 ′′ ∑ f (ξi) = − N f ′′(η) 12

∈]a, b[. Or, ici h =

−

b a N

et on obtient le résultat esultat pour la formule des trapèzes ezes R N ,1 ( f ) =

−

h2 (b 12

− a) f ′′(η), avec

η ]a, b[.

∈

(1.42)

Pour aboutir aboutir a` une une expr expres essi sion on de l’err l’erreu eurr d’in d’int´ tégra e gratio tion n pour la form formul ulee de Simp Simpso son n = 2 M et (1.28) on procède ede de manière ere analogue. Dans ce cas N = et R N ,2 ( f ) = avec Ri ( f ) =

=

Z

x2i+2

x2i

−

Z a

−

M 1

b

f ( x)dx I N ,2 =

f ( x)dx

−

− h3 ( f ( x

h5 (4) f (ξi ), ξ i 90

∑ Ri ( f )

i=0

2i ) + 4 f ( x2i+1 ) + f ( x2i+2 ))

∈] x x

2i , x2i+2 [.

19


Or, a` nouveau par par le th eor` e´ orème eme de la valeur intermédiaire ediaire il existe η ]a, b[ tel que

∈

−

M 1

∑ Ri ( f ) = −

i=0

h5 M f (4) (η). 90

= (b Etant donné que h = (b a)/ N = (b a)/(2 M ), on aura M = finalement on obtient la formule d’erreur

−

R N ,2 ( f ) =

−

h4

− 180 (b − a) f ( )(η), avec 4

η ]a, b[.

∈

− a)/(2h) et

(1.43)

Comparant R N ,1 avec R avec R N ,2 , on observe que la formule de Simpson est en O(h4 ) tandis que la formule des trapèzes ezes n’est qu’en O (h2 ). De mettre en œuvre la formule de Simpson est a peine plus complexe que d’utiliser la formule des trapèzes, ezes, ce qui fait que la formule de Simpson est largement utilis ee. e´ e. Il faut cependant être etre conscient, que ce résultat esultat d’erreur pour la méthode ethode de Simpson est obtenu pour des fonctions 4 fois continûment ument d´ derivables e´ rivables dans [ a, b].

20

Chapitre 2 Resolution e´ solution num´ numerique e´ rique des equations e´ quations diff erentielles e´ rentielles ordinaires (EDO) Avant vant d’abor d’aborder der quelqu quelques es méthodes ethodes d’approx d’approximatio imation n d’équations equations différentielles erentielles ordin ordinair aires es (on utilise utilisera ra l’abr´ l’abréviat eviation ion EDO), EDO), il con convient vient de passe passerr en revue revue quelqu quelques es résultats esultats gén´ enéraux. eraux.

2 .1

Resultats e´ sultats g´ gen´ e´ neraux e´ raux sur les EDO

Une equation e´ quation différentielle erentielle est une relation de la forme du (t ) = f (t , u(t )) )) dt

(2.1)

ou` u(t ) est l’inconnue, la solution de l’équation equation différentielle, erentielle, qu’il convient de déterminer eterminer tandis que f (t , u) est une fonction donnée. ee. La solution u(t ) d epend e´ pend de t de t R (t (t est est par exemple exemple le temps) temps) ; u(t ) peut etre eˆ tre une fonction scalaire mais u 1 (t ), u2(t ), , un (t ) et dans ce cas la foncaussi vectorielle avec n avec n composantes composantes u tion f a egalement n e´ galement n composantes composantes f 1 (t , u), f 2(t , u), f n (t , u) qui peuvent etre eˆ tre des u . Donc, f (t , u) est une apfonctions non linéaires eaires de toutes les composantes de u. plication de R Rn a` valeurs dans Rn. Le problème eme a` valeur initiale consiste à adjoindre a` l’équation equation (2.1) une condition dite initiale en t en t 0 avec

∈ ∈

··· ···

×

u(t 0 ) = u0

(2.2)

(on prendra souvent t 0 = 0). Ici nous considérons erons des equations e´ quations différentielles erentielles d’ordre 1, c.-à-d. a-d. seule la dériv´ erivée ee première ere de la fonction u fonction u (t ) intervient. En fait, une equation e´ quation d’ordre supérieur erieur avec avec des conditions initiales peut toujours etre eˆ tre reformul formuler er comme comme un syst` système eme d’ODE. d’ODE. Preno Prenons ns par exemp exemple le le probl` problème e me du pend pendul ulee

θ′′ (t ) =

− sin(θ(t )), )),

θ(0) = θ0 , 21

θ′(0) = θ1 .

R´ esolution num´ erique des equations ´ diff´ erentielles ordinaires (EDO)

Alors, ecrivant u e´ crivant u 1 = θ et u et u2 = θ′ on aura u aura u 1′ = u2 et le système eme s’écrit ecrit d dt

   u1 u2

=

−

 

u2 sin(u1 )

,

u1 (0) u2 (0)

   =

θ0 θ1

.

Pour qu’une EDO ait une solution, il faut que la fonction f (t , u) ait quelques propriétés es de régularit´ egularité. e. Le théor` eorème eme fondamental quant a` l’existence et l’unicité de la solution d’une EDO avec condition initiale peut etre eˆ tre résum´ esumé comme suit.

Th´ Theor` e´ oreme e` me 4 (Th´ (Theor e` me de Cauchy Lipschitz) ´ eme Soit f (t , u) est une application de R Rn a` valeurs dans Rn . Soient u0 Rn et t 0 R donn´ donnes ´ : on suppose qu’il existe un domaine D Rn contenant u0 et un intervalle [t 0 , t 1] ainsi qu’une constante L > 0 , tels que quels que soient v, w D et quel que soit t [t 0 , t 1] , ,

×

∈

∈

∈ ∈

∈

∈

|| f (t , v) − f (t , w)|| ≤ L||v − w||

(2.3)

(o` (ou` ˙ d esigne une norme de Rn , par exemple exemple la norme euclidienne). On dira ´ ´ que f (t , u) est lipschitzienne par rapport a` u de constante de Lipschitz L. Alors il existe un plus grand temps T [t 0 , t 1] , , tel que l’´ l’equation ´ diff erentielle ´ ´

||||

∈

du (t ) = f (t , u(t )), )), dt poss` possede e` de une solution et une seule pour t

u(t 0 ) = u0 ,

∈ [t , T ]. 0

Afin d’interpréter eter ce théor` eorème, eme, prenons la cas d’une EDO scalaire : il est alors u , alors f (t , u) est lipschitfacile de voir que si f (t , u) est dérivable erivable par rapport a` u, zienne. En effet, effet, par le th eor` e´ orème eme des accroissements finis

 | ≤ ∈    ∈

| f (t , v) − f (t , w)

| − | 

∂ f max (t , u) v u D ∂u

w

∈ ∈ [t , t ]

et on pourra prendre comme constante de Lipschitz le maximum pour t de la quantité ∂ f (t , u) . max u D ∂u

0

1

Pourquoi est-il est-il n´ necessaire e´ cessaire de supposer que f (t , u) est lipschitzienne lipschitzienne et que signifie signifie l’existence d’un plus grand T grand T tel que la solution existe pour t pour t [t 0, T ] ? Prenons par exemple exemple l’´ l’equation e´ quation différentielle erentielle du (t ) = (u(t )) ))2 dt 22

∈

Systèmes emes d’equations ´ diff´ erentielles lin´ eaires a` coefficients coefficients constants

a . Ici la fonction f (u) = u 2 ne dépend avec u(0) = a > 0, donc t 0 = 0 et u0 = a. epend pas de t de t et et la fonction est lipschitzienne pour tout domaine contenant a et pour tout intervalle de temps [0, t 1]. Il est facile de voir que la solution de l’équation equation différentielle erentielle est 1 u(t ) = −1 a t

−

et par conséquent equent la solution tend vers l’infini quand t 1/a. On peut donc déduire eduire que le plus grand temps pour l’existence de la solution est T = 1/a. Prenons maintenant comme exemple du (t ) = dt

√ ≥



u(t ),

→ →

u(0) = 0.

√

Ici, f (u) = u, u 0 et f (u) / u = 1/ u tend vers l’infini quand u 0. Par f n’est pas lipschitzienne au voisinage de u = 0. L’ODE ne vérifie conséquent, equent, f n’est erifie pas les conditions du théor` eorème. eme. On observe que cette équation equation différentielle erentielle n’a pas de solution unique : en effet

|

|||

u(t )

→

≡0

est solution, mais aussi

1 u(t ) = t 2. 4 Avant d’aborder des méthodes ethodes num numeriques e´ riques de résolution, esolution, nous allons allons passer passer en revue vue quel quelqu ques es résulta esultats ts gén´ enérau eraux x conce concern rnant ant les syst` systèmes e mes d’équations equations différentielles erentielles linéaires. eaires.

2.1. .1.1

Syst ystemes e` mes d’´ d’equations e´ quations diff erentielles e´ rentielles lin´ lineaires ` e´ aires à coefficie coefficients nts constants

×

Soit A Soit A une une matrice n matrice n n à coefficients coefficients r´ reels e´ els constant constants. s. On consid` considère ere le syst` système eme d’équations equations différentielles erentielles u ′ (t ) = Au(t ),

u(t 0) = u0 ,

(2.4)

avec u0 Rn donné. e . Tout out d’ab d’abor ord d on cons consta tate te que que ce syst` système eme poss` possède ede une solutio solution n n Au et pour toute norme de R unique. En effet, ici f (u) = Au et

∈

Av − Aw|| = || A A(v − w)||. || f (v) − f (w)|| = || Av A|| et on justifiera des On verra au chapitre 4.2 la notion de norme de matrice || A majorations de la forme

|| A A(v − w)||≤|| A A||||v − w||. 23


Au est donc L = A A la norm La foncti fonction on f (u) = Au est donc lipschitzie lipschitzienne nne avec avec L normee de A de A comme comme consta constante nte de Lipsc Lipschitz hitz.. Avant vant de pours poursui uivre vre,, il convie convient nt de pr´ préciser eciser de quelle quelle façon ¸on une equation e´ quation différentielle erentielle d’ordre n peut etre eˆ tre ecrite e´ crite sous forme d’un système eme j) ( d’ordre 1. On note v note v (t ) la dériv´ erivée ee j eme e` me d’une fonction scalaire v(t ) R et on considère ere l’équation equation différentielle erentielle d’ordre n d’ordre n

|| ||

∈

v(n) (t ) + an−1 v(n−1) (t ) +

· · · + a v ′(t ) + a v(t ) = 0 1

(2.5)

0

avec les conditions intiales v(t 0 ) = v0 , v ′ (t 0 ) = v1 , (v0 , v1 ,

· · · , v( − )(t ) = v − n 1

0

n 1

(2.6)

es). Donc, si l’on note · · · , v − donnés). · · · , u = v( − ) u = v, u = v ′ , ·· (2.7) on obtie obtient nt pour pour le vect vecteu eurr u = (u , u , · · · , u ) un syst` système e me d’équations equations différentielles erentielles n 1

1

2

1

avec A = u ′ (t ) = Au(t ), avec A

n 1

n

2

  

n

0

T

1 0

1 .. .

..

. 0

1 an−1

−a −a ··· ··· − 0

1

  

. (2.8)

En effet effet,, ui′ = ui+1 , i = 1, n 1, les les ui étan e tantt défini e finiss par par (2.7 (2.7), ), et un′ = ∑ jn=1 a j−1 u j par (2.5). Au système eme (2.8) est bien entendu associée ee la condition initiale

··· −

u(t 0) =



−

T



v(t 0 ), v ′ (t 0 ), · · · , v(n−1) (t 0 )

.

Revenons a` (2.4) et considérons erons d’abord le cas scalaire n = 1 et la solution de l’´ l’equation e´ quation (2.4) avec A avec A = a R s’écrit ecrit bien entendu

∈

u(t ) = u0 e a(t −t 0 ) .

×

n n, c’es On cherche cherche a` gén´ enéral e ralis iser er ce résulta esultatt pour pour des des matric matrices es n c’est-` t-à-dire a-dire on cher cherche che a` définir efinir l’exponentielle d’une matrice. On se rappelle que l’exponentielle d’un nombre nombre r´ reel e´ el est donnée ee par une série erie ce qui conduit conduit par analogie a` la définition efinition suivante

Definition e´ finition 2 Soit A une matrice n trice A, not ´ not ee ´ e A , par la s erie ´ A

e =

× n, alors on d efinit l’exponentielle de la ma´ ´

∞

1 j 1 2 A I A A + = + + ∑ j! j ! 2 j=0

· · · + j! j1! A + · · · j

avec la convention que A 0 = I, I , I etant ´ la matrice identit e´ ´ n 24

× n.

(2.9)

Systèmes emes d’equations ´ diff´ erentielles lin´ eaires a` coefficients coefficients constants

Donc, e Donc, e A est une série erie de de matrices matrices et et on constate constate que cette cette série erie converge. En effet, soit une somme partielle de cette série erie N

1

A j , ∑ j! j !

S N =

j=0

alors

N

||S || ≤ N

N

1 A A j ∑ j! j j=0 !

|| || ≤ ∑ j! j1! || A A|| =

j

j 0

d’après es la définition efinition et les propriét´ etés es des normes de matrices qui seront abordées ees ultérieurement erieurement au chapitre 4.2. On note a = A A et donc

|| || ∞

1 || || || || ≤ a ∑ →∞ ! j! j = A

lim S N = e

n

j

= ea

j 0

et on déduit eduit de l’inégalit´ egalité ci-dessus la convergence convergence de la série. erie. Avec cette définition, efinition, on observe que si l’on note 0 la matrice n les coefficients sont nuls, alors e0 = I .

× n dont tous

Pour des nombres réels, eels, on a la relation bien connue ea+b = e a eb : cette relation B , sauf dans le cas particulier n’est en gén´ enéral eral pas vérifi´ erifiée ee pour deux matrices A, B, (on admettra ce résultat) esultat) : si les matrices A matrices A, B, commutent ( AB = BA), alors e alors e A+ B = e A e B .

(2.10)

On déduit eduit par exemple de cette propriét´ eté que e A− A = e0 = I = = e A e− A et donc l’inverse de la matrice e matrice e A est simplement

 − e A

1

= e− A

(2.11)

La définition efinition de la norme d’une matrice nous permet permet d’´ d’ enoncer e´ noncer le résultat esultat fondamental suivant.

Th´ Theor` e´ oreme e` me 5 La solution du syst ` syst eme e` me d’´ d’equations diff erentielles ´ ´ ´ u ′ (t ) = Au(t ), s’´ s’ecrit ´

u(t 0) = u0

u(t ) = e(t −t 0 ) A u0 . 25

(2.12)


Pour démontrer emontrer ce résultat, esultat, observons d’abord que u(t 0) = e0 u0 = Iu I u0 = u0 donc le vecteur u(t ) v erifie e´ rifie la condition initiale. Ecrivons maintenant le terme exponentiel sous sous forme de s erie e´ rie (t t 0 ) A

e

−

∞

=

∑

j=0

(t t 0 ) j

−

j! j!

A j .

On peut facilement justifier pour cette série erie que la dériv´ erivée ee par rapport au temps de la somme de la série erie est egale e´ gale a` la somme obtenue en dérivant erivant terme par terme de la série, erie, a` savoir d (t −t 0 ) A e dt

 

∞

∞ j (t t 0 ) j−1 j (t t 0 ) j−1 j = ∑ A =∑ A j! j ( j ) ! 1 ! j=1 j=1

−

∞

− −

∞ (t t 0 ) j−1 j−1 (t t 0 )k k = A∑ A = A ∑ A = Ae(t −t 0 ) A ( j 1)! k ! j=1 k =0

− −

Par conséquent, equent, si

−

u(t ) = e(t −t 0 ) A u0 ,

alors

d (t −t 0 ) A e u0 = Ae(t −t 0 ) A u0 = Au(t ), dt ce qui achève eve la démonstration emonstration du théor` eorème. eme. u ′(t ) =

 

Traitons brièvement evement un système eme d’équations equations différentielles erentielles non homogène, ene, c’esta-dire a` -dire de la forme (2.13) u ′ (t ) = Au(t ) + g(t ), u(t 0 ) = u0 avec une fonction g(t ) donnée. ee. Faisons d’abord abstraction de la condition au temps initiale : alors on cherche u(t ) comme une somme d’une solution du probl` bleme e` me homogène ene v ′ (t ) = Av(t ) et d’une solution particulière ere de w ′ (t ) = Aw(t ) + g(t ).

(2.14)

Il est clair que v(t ) = etA v0 , quel que soit le vecteur v0 . Pour la solution particonstante : on ecrit culière, ere, on applique la méthode ethode dite de la de la variation de la constante : e´ crit w(t ) = etA c(t ) 26

Calcul de l’exponentielle de la matrice

pour c pour c (t ) une fonction vectorielle dépendant ependant du temps. Dérivant w erivant w (t ) on obtient w ′ (t ) =

d tA e c(t ) + etA c ′ (t ) = AetA c(t ) + etA c ′ (t ) = Aw(t ) + etA c ′ (t ) dt



et injectant cette expression dans (2.14), on obtient etA c ′ (t ) = g(t ). De (2.11) on déduit eduit que c ′ (t ) = e−tA g(t ) dont la primitive s’annulant en t en t 0 s’´ s’ecrit e´ crit

Z

t

c(t ) =

t 0

e−sA g(s)ds ,

l’intégration egration etant e´ tant a` opérer erer composante par composante du vecteur e vecteur e −sA g(s). On en déduit eduit que la solution particulière ere qui s’annule en t en t 0 s’´ s’ecrit e´ crit w(t ) = etA

Z

t

t 0

e−sA g(s)ds =

Z

t

t 0

e(t −s) A g(s)ds .

(2.15)

Reste a` déterminer eterminer pour la solution gén´ enérale erale v(t ) = e tA v0 le vecteur v0 tel que v(t 0 ) = u0 : or v(t 0 ) = et 0 A v0 = u0 implique

v0 = e−t 0 A u0 .

La solution de (2.13) s’écrit ecrit en additionnant v additionnant v(t ) et w et w (t ) et donc u(t ) = e

(t t 0 ) A

−

Z

t

u0 +

t 0

e(t −s) A g(s)ds .

(2.16)

On voit que l’exponentielle de la matrice A est la quantité-clef e-clef qui permet de résoud esoudre re les syst` systèmes e mes d’équations equations différenti erentiell elles es linéaires, eaires, a` coeffic coefficients ients constants constants.. Au paragraphe suivant nous allons donner un mode d’emploi pour le calcul de l’exponentielle d’une matrice, avant de l’illustrer par l’exemple d’un système eme simple 2 2.

×

2.1.2 2.1.2

Calcul Calcul de l’expone l’exponentie ntielle lle de la la matri matrice ce

On supposera pour simplifier que t 0 = 0 et donc la solution du système eme (2.4) s’écrit ecrit u(t ) = etA u0 . 27


On peut remarquer que de choisir t choisir t 0 = 0 n’enlève eve rien a` la gén´ enéralit´ eralité. e. En effet, la (t −t 0 ) A u0 qui vérifie u solution u solution u (t ) = e erifie u (t 0 ) = u0 est identique a` la solution u solution u (t ) = tA t 0 A t −t 0 ) A − ( e u˜0 qui vérifie u u0 car nous avons vu que e = etA e−t 0 A . erifie u (0) = u˜0 = e

Digression : valeurs et vecteurs propres d’une matrice. Ici il convient convient de rappele rappelerr la notion notion de valeur valeur propre propre d’une matrice A matrice A.. Il s’agit d’un nombre λ , complexe en gén´ enéral, eral, tel qu’il existe un vecteur x = 0, appelé vecteur propre associ associe´ a` la valeur propre λ , tel que



A x = λ x. Cette relation peut encore s’écrire ecrire ( A λ I ) x = 0 pour un vecteur x = 0 avec I la matrice identité n n. Par conséquent, equent, pour cette valeur λ la matrice (carrée) ee) A λ I n’est n’est pas inversible (on dit que cette matrice a un noyau non nul) et donc det( A λ I ) = 0. A partir partir de la définition efinition du déterminant eterminant on peut se convaincre que (pour tout nombre λ) p(λ) = det( A λ I ) est en fait un polynôme ome en λ de degré n. Un polynôme n a précis´ n zéros ome de degré n a ecisément ement n z e´ ros (complexes ou réels, eels, compt´ comptés es avec leurs multiplicités, es, c’est-à-dire a-dire un zéro ero double compte deux fois etc.). Donc toute matrice n n possède ede n valeurs propres, complexes ou réelles, eelles, compt´ comptées ees avec leurs multiplicités es eventuelles, e´ ventuelles, qui qui sont les z´ z eros e´ ros de p(λk ) = det( A λk I ) = = 1, 2, , n. A chaque valeur propre est associé un vecteur propre xk tel que 0, k = Axk = λk xk . Supposons que les n valeurs n valeurs propres λ1 , λ2 , , λn de A de A sont sont distinctes et formons la matrice P dont les colonnes sont précis´ ecisément ement les vecteurs propres x1 , x1 , , xn associés. es. On forme la matrice diagonale D diagonale D

−

−

×

−



−

×

−

···

···

···

D =

  

λ1

(0) λ2 ..

.

λn

(0)

  

avec sur la diagonale les valeurs propres de A. Les règles egles de multplication de = 1, , n sous forme matrices matrices permettent permettent d’écrire ecrire les n egalit´ e´ galités es Axk = λ k xk , k = matricielle AP = PD.

···

En effet, la kème eme colonne de AP de AP est est précis´ ecisément Ax ement Ax k , si la kème eme colonne de P de P est est le vecteur x vecteur xk , et la kème eme colonne de PD de PD est est λk xk . On peut montrer que les vecteurs propres asociés es a` des valeurs propres distinctes sont linéairement eairement indépendants ependants et il s’ensuit que la matrice P est inversible. Multipliant les deux membres de l’´ l’egalite e´ galite ci-dessus par l’inverse P l’inverse P −1 de P de P,, on obtient l’égalit´ egalité A = PDP−1 . 28


A sont distinctes, on peut utiliser cette relaDans ce cas où les valeurs propres de A sont tion pour trouver une expression relativement relativement simple des puissances pu issances successives successives de A. En effet, A2 = PDP −1 PDP−1 = PD 2 P−1 et en itérant erant on trouve bien sûr ur j j −1 A = PD P de fac faç on a` ce que etA =

∞



j

t

∞

j

t

PD j P−1 = P ∑ D j ∑ j! j! j! j!

j=0

j=0

− P

1

= Pe PetD P−1 .

On peut se convaincre aisément ement que

t j j D = j! j !

  

t j j λ j! j! 1

(0) j

t j λ j! j! 2

..

. t j j λ j! j! n

(0)

  

et par conséquent equent l’él´ elément ement de la kème eme position sur la diagonale de la matrice (diagonale) e (diagonale) etD est l’exponentielle e l’exponentielle et λk et

etD =

  

t λ1

(0)

e

t λ2

e

..

. et λn

(0)

  

.

Donc, par etA = Pe tD P−1 , on peut affirmer que chaque él´ elément ement de etA est une = 1, 2, , n. Bien sûr, combinais combinaison on linéaire eaire des et λk , k = ur, si on est en mesure de − 1 déterminer eterminer explicitement P ainsi que P , ce qui souvent pour n 3 n’est pas aisé, e, on trouve une expression explicite de etA et donc aussi de u (t ) = etA u0 . On peut donc enoncer e´ noncer le résultat esultat suivant.

···

≥

Proposition Proposition 2 Chaque composante u j (t ) du vecteur solution u(t ) = ( u1 (t ), u2(t ),

· · · , u (t )))) n

de (2.4) est de la forme n

u j (t ) =

∑ b jk eλ t k

k =1

pour des coefficients coefficients b jk . 29

T


Il convient de remarquer ici que les valeurs valeurs propres d’une ma trice A trice A ` à coefficients réels eels ne sont pas forcément ement réelles. eelles. Cependant, s’il existe une valeurs propre λ k = αk + iβk avec β k = 0, la valeur conjuguée ee complexe ¯λk = αk iβk est egalement e´ galement valeur propre, la matrice A étant etant r´ reelle. e´ elle. Donc, une combinaison linéaire eaire a` valeur ¯ réelle eelle faisant intervenir eλk t = e αk t eiβk t et eλk t = e αk t e−iβk t peut toujours toujours s’ ecrire e´ crire sous la forme eαk t (c cos (βk t ) + d sin sin (βk t )



−

d . avec des coefficients réels c eels c et et d Pour le système eme (2.8) equivalent e´ quivalent a` l’équation equation scalaire (2.5) d’ordre n d’ordre n,, u 1 (t ) = v(t ) et la solution de (2.5) s’écrit ecrit n

v(t ) =

∑ bk eλ t

(2.17)

k

k =1

A de (2.8). Cette fonction est solution de avec λk les valeurs propres de la matrice A de l’´ l’equation e´ quation différentielle erentielle (2.5) quels que soient les coefficients b coefficients b k , si on ne précise ecise t λ pas les conditions initiales. En particulier les fonctions e fonctions e k sont solutions et injectant cette expression dans (2.5), on trouve

λnk + an−1 λnk −1 +

· · · + a λ + a = 0. 1 k

0

On en déduit eduit le resultat e´ sultat que les valeurs propres propres de la matri ce A ce A dans dans (2.8) associ associee e´ e a` l’équation equation (2.5) d’ordre n d’ordre n sont sont les zéros eros du polynôme ome caract caracteristique e´ ristique p(λ) = λn + an−1 λn−1 +

···+a λ +a . 1

0

(2.18)

Exemple : On considère ere l’équation equation du pendule avec un coefficient de frottement α v ′′ =

−v + αv ′

et d’après es ce qui préc` ecède, ede, en notant u1 = v et u2 = v ′ , ce système eme peut encore s’écrire ecrire d u1 u1 0 1 = . (2.19) u2 1 α dt u2

  

 

−

A du système Les valeurs propres de la matrice A du eme sont

λ1/2 =

α

±

√ α 2

2

−4.

On suppose que α = complexes avec

(2.20)

 ±2. Si |α| < 2, alors les valeurs propres sont conjuguées ees λ1/2 = α/2

 ± − i

30

1

(α/2)2


et la solution u solution u1 (t ) = v(t ) s’écrit ecrit t (α/2)

v(t ) = e

(a cos(ωt ) + b sin(ωt )) )) avec ω =

 − 1

(α/2)2 .

Si α < 0, alors la solution est amortie car alors et (α/2) décroˆ ecr oˆıt ıt pour pou r t > 0. Si ce > 0. Si α > 2, alors les deux pendant α > 0, alors la solution est amplifiée ee pour t pour t > solution sont réelles eelles et la solution s’écrit ecrit

||

v(t ) = aet λ1 + bet λ1 avec λ1/2 donnés es par (2.20). Considérons erons maintenant maintenant le cas o u` la matrice A matrice A dans dans (2.4) est de la forme

A =

On peut ecrire A e´ crire A sous sous la forme

  

λ

1 .. .

et N = =

On peut se convaincre que

2

N =

  

0

0 .. .

1 .. . 0

(0)

  

. 0 0

1 0 0

.

(2.21)

I matrice identité n

0

1 .. .

(0) ..

(0)

  

1

λ

(0)

(0) ..

..

λ

A = λ I + N avec

  

(0)

,

···,

. 0

1 0

  

N n−1 =

×n

.

  

0

(0)

0 0

(2.22)

··· ···

0 1 0 0 . . .. . .. .. 0 0 0

  

et qu’ainsi N n = 0 est la matrice identiquement égale egale a` zéro. ero. Donc, lorsqu’on forme les puissances successives de N les les el´ e´ léments ements non nuls “remontent”. Les matrices λ I et N et N commutent commutent et d’après es ce qui préc` ecède ede etA = et (λ I + N ) = et λ I etN . 31


Or, etant e´ tant donné que N que N n = 0, e 0, etN est en fait une somme finie et tN

e

t 2 2 = I + t N + N + 2

···

t n−1 + N n−1 . (n 1)!

−

D’après es les produits successifs de N ci-dessus ci-dessus

etN =

   

1

t 2 2

t

··· t 2 2

t

1

..

..

.

t n−2 (n 2)!

.

1

(0)

−

··· ..

t n−1 (n 1)! t n−2 (n 2)!

   

− − .. .

.

t 2 2

t 1

t 1

.

La matrice e matrice et λ I est diagonale avec eλt sur la diagonale et donc

etA =

   

eλt

t eλt eλt

t 2 λt 2e λt

te ..

.

···

2

t λt e 2

..

.

eλt

(0)

t n−2 λt (n 2)! e

−

··· ..

t n−1 λt (n 1)! e t n−2 λt (n 2)! e

− −

.

.. .

t 2 λt 2e λt

t eλt eλt

te eλt

   

Chaque el´ e´ lément ement de etA est ainsi de la forme p(t )eλt avec p(t ) un polynôme ome de degré n 1. On peut donc conclure, que si dans le système eme (2.4) la matrice A matrice A est de la forme (2.21), alors les composantes ui (t ) de la solution u (t ) sont de la forme ui (t ) = p i (t )eλt

≤ −

≤ − 1. On observe que λ est valeur propre

pour des polynômes omes pi (t ) de degrés es n de A de A donn´ donnée ee par (2.21) de multiplicité n. n .

×

A matrice n n n quelconque : si ses n valeurs n valeurs propres sont disRevenons alors a` A matrice tinctes, alors A est diagonalisable. Supposons cependant que certaines des valeurs propres sont sont multiples, plus pr´ precis´ e´ cisément ement supposons que A que A poss` possède m ede m valeurs valeurs propres distinctes de multiplicit´ multiplicites e´ s respectivement l respectivement l 1 , l2 , , lm (avec l (avec l1 + l2 + + lm = n). n). Alors A A peut-être Alors A n’est n’est plus forcément ement diagonalisable : cependant cependant A peut-ˆ etre mise 1 − Jordan , a` savoir il existe des matrices P P telles que sous forme dite de dite de Jordan, matrices P et et P

···

A = PJP PJ P−1 32

···

(2.23)


avec J avec J une une matrice par blocs

             ×  ∗    ∗  ∗  J 1

0

J 2

J = =

..

.

0

(2.24)

J m

l k chaque bloc J bloc J k etant etant une sous-matrice sous-matrice l k ´

λk

J k k =

,

lk de la forme

(0)

..

.

..

.

,

λk

= 0 ou 1 .

(2.25)

λk

(0)

Donc, chaque sous-matrice J sous-matrice J k e de la valeur propre k est de la taille de la multiplicit´ λk , qui se trouve sur la diagonale de J k k , les coefficients sur la quasi-diagonale immédiatement ediatement au-dessus de la diagonale sont égaux egaux a` 1 ou egaux égaux a` 0, selon des cas particuliers, les autres coefficients de J k e´ tant nuls. Bien-sûr, ur, si toutes k etant les valeurs propres sont simples, les J k 1 donc des k sont en fait des matrices 1 scalaires et on retrouve une matrice J diagonale. diagonale. La mise sous forme de Jordan est assez complexe complexe en g´ g en´ e´ néral eral et la m´ methode e´ thode est donnée ee dans des ouvrages d’algèbre ebre linéaire. eaire. Ce qu’il convient de retenir est que chaque bloc est de la forme

×

J k k = λk I k k + N k k

(2.26)

×

avec I k e lk lk et N k ` (2.22), mais ici de k matrice identit´ k une matrice similaire a taille l taille l k lk (et tous les el´ e´ léments ements juste au-dessus ne sont pas forcément ement egaux e´ gaux a` n ème 1). Nous avons vu que la n ` eme puissance de la matrice (2.22) est identiquement i dentiquement egale e´ gale a` 0 et de manière ere analogue on peut affirmer que

×

l

N k k = 0 Par (2.23), on aura à nouveau A j = PJ j P−1 et donc (analogue au cas diagonalisable) etA = Pe PetJ P−1 . Les multiplications de J de J avec avec elle-même eme se faisant par blocs, il est facile de voir que 0 etJ 1 tJ 2 e etJ = (2.27) .. .

        

etJ m

0

33

  


a egalement e´ galement une structure par blocs. La décomposition ecomposition (2.26) (2.26) est similaire a` celle tJ k pour la matrice (2.21) : donc, donc, chaque chaque él´ elément ement de e de e est de la forme p k (t )eλk t avec pk (t ) polynôme ome (en fait un monôme) ome) de degré lk 1. On peut ainsi enoncer e´ noncer le résultat esultat gén´ enéral. eral.

≤ −

Proposition Proposition 3 Soit le syst eme e` me (2.4) tel que A poss ede e` de m valeurs propres λ k , k = = 1, , m distinctes de multiplicit es ´ ´ respectivement l 1 , lm (avec l 1 + lm = n). Alors les composantes u j (t ) de la solution u (t ) de (2.4) sont de la forme

···

···

···

m

u j (t ) =

∑ p jk (t )eλ t

(2.28)

k

k =1

avec p jk (t ) des polynomes de degr ´ degr es ˆ ´

≤ l − 1. k

On peut peut remar remarque querr que si l si lk = 1, les p les p jk (t ) sont sont des des polynˆ polynômes omes de degr´ degré zéro, ero, donc des constantes. Par Par cons´ cons equent, e´ quent, on retrouve évidemment evidemment le résultat esultat de la proposi A sont simples. tion 5 lorsque les valeurs propres de A sont

×

Exemple : On reprend le système 2 2 (2.19), les valeurs propres de la matrice en question etant e´ tant (2.20). Pour α = 2, λ1 = 1 est double et d’après es ce qui préc` ecède ede la solution v solution v(t ) = u1 (t ) s’écrit ecrit

±

±

v(t ) = ( a + bt )e±t , a b étant pour le polynôme ome p1 (t ) = a + bt de de degré 1, les l es coefficients coefficients a etant d´ determin´ e´ terminés es par la condition initiale.

2.2

Schemas e´ mas a` un pas pour la r´ resolution e´ solution d’une EDO

Sauf dans des cas particuliers, par exemple pour des systèmes emes d’équations equations différentielles erentielles linéaires eaires a` coefficients coefficients constants, il n’est en g´ g en´ e´ néral eral pas possible de trouver des solutions explicites et analytiques d’une équation equation différentielle erentielle u ′ (t ) = f (t , u(t )), )),

u(t 0) = u0 ,

(2.29)

a` partir du moment moment o` ou` la fonction f (t , u) avec u avec u Rn est quelque peu complexe et que n 2. La variable t sera sera appelée e e désormais esormais le temps et on cherche à trouver une solution approchée ee de la solution pour t > t 0 . On supposera toujours que f (t , u) satisfait a` la conditions de Cauchy-Lipschitz du théor` eorème eme 4 pour tout a-dire l’´ l’equation e´ quation différentielle erentielle possède ede une solution unique dans t [t 0, T ], c’est-à-dire [t 0 , T ] avec u avec u (t 0 ) = u0 . Une solution numérique erique sera sera d´ determin´ e´ terminée ee en des temps discrets

∈

≥

∈

t k k +1 = t k k + h, 34

k = = 0, 1,

···

Sch´ emas a` un pas pour la r´ esolution d’une EDO

et on appellera h le pas de temps. On notera U k temps t k k l’approximation au temps t k de la vraie solution u solution u (t k ema ema a` un k ). Lorsqu’on met en œuvre ce qui est appelé un sch´ pas, U pas, U k ecurrence k est solution de la récurrence U 0 = u0 U k k +1 = U k k + hF (t k k , U k k , h),

t k k +1 = t k k + h,

k = = 0, 1, 2,

···

(2.30)

F ` pour pour une foncti fonction on F apr´ a` précise eciserr. Donc, Donc, au temps temps initial initial t 0 , U 0 est egale e´ gale a` la cond condit itio ion n initiale u initiale u 0 de l’EDO. Connaissant Connaissant U 0 on peut déterminer U eterminer U 1 par (2.30) et ainsi de suite pour obtenir l’approximation U l’approximation U k k de la solution pour des temps discrets successifs. Pour simplifier l’exposé de la méthode, ethode, on supposera que u (t ) R. Dans le cas vectoriel il conviendra simplement d’appliquer l’approche donnée ee pour le cas scalaire composante par composante au cas vectoriel.

∈

Si l’on prend par exemple dans (2.30) F (t k k , U k k , h) = f (t k k , U k k ) schema avec f la fonction fonction de l’EDO (2.29), on obtient obtient la récurrence ecurrence appelée le ee le sch´ ´ d’Euler explicite U k k = = 0, 1, (2.31) k +1 = U k k + h f (t k k , U k k ),

···

Si la soluti solution on exac exacte te u(t k etait etait aussi aussi solutio solution n de (2.31) (2.31),, alors alors la soluti solution on numérique erique k ) ´ et la solution sol ution exacte co¨ıncideraient, ınciderai ent, ce qui bien bi en sˆ s ur u ˆ r est faux en gén´ enéral. eral. En fait, la solution exacte ne vérifie erifie (2.31) qu’à une erreur près. es. En effet, par la formule de Taylor on voit aisément ement que 2 u(t k k +1 ) = u(t k k ) + h f (t k k , u(t k k )) + O(h )

′ k ) = f (t k k , u(t k k )). Il convient car si u si u (t k k ) est solution exacte, par (2.29) on aura u (t k de remarquer ici que par la suite il sera toujours supposé que f (t , u), et ainsi que ordres de dérivabilit´ erivabilité suffisants pour opérer erer des développements eveloppements de u(t ) ont des ordres Taylor nécessaires. ecessaires. D’autr D’autres es méthode ethodess d’appr d’approxi oximat mation ion sont sont définies efinies a` part partir ir des des form formul ules es d’in d’int´ tégration. egration. A partir de (2.29) on peut écrire ecrire u(t k k +1 ) = u(t k k ) +

Z

t k +1

t k

f (t , u(t )) ))dt

(2.32)

et prenons par exemple exemple la formule des trap` trap ezes e` zes pour approcher approcher l’int´ l’int egrale. e´ grale. Il a et´ e´ té montré au chapitre 3 que

Z

t k +1

t k

f (t , u(t )) ))dt = =

h 2

[ f (t k k , u(t k k )) + f (t k k +1 , u(t k k +1))]+ O(h3 ) 35


et une solution approchée ee de l’EDO est donnée ee par le schéma ema U k k +1 = U k k + hF (t k k , U k k , h),

F (t k k , U k k , h) =

1 [ f (t k k , U k k ) + f (t k k +1 , U k k+ )]. (2.33) 1 )]. 2

La foncti fonction on F (t k efinit efinit le le schéma ema (2.33 (2.33)) fait fait interv interveni enirr U k qui dépend epend k , U k k , h) qui d´ k+ 1 qui implicitement de U k eterminer U eterminer U k esoudre esoudre une k : pour d´ k +1 en fonction de U k k il faut r´ implicite, implicite, tandis e´ quation qui est non linéaire equation eaire en gén´ enéral. eral. On parle d’une m d’une m ethode ´ explicite, dans la mesure où U k que le schéma ema Euler (2.31) est explicite, k +1 est obtenue par U k une simple relation algébrique ebrique en fonction de de U ema ema (2.33) est k . Cependant le sch´ plus précis ecis que le schéma ema d’Euler explicite dans un sens qui sera précis´ ecisé ci-après. es.

2.2.1 2.2.1

Ordr Ordree d’un d’un sch´ schema, e´ ma, consistance, stabilit´ stabilite´ et convergence

Supp Suppos oson onss que que pour pour le temp tempss t k solutio tion n num´ numérique erique co¨ co¨ıncide ıncide avec avec la solution solution k la solu exacte, c’est-à-dire U a-dire U k e´ valuer l’erreur entre la solution k = u (t k k ) et on cherche à evaluer exac exacte te et la solu solutio tion n num´ numériqu eriquee au temps temps t k aleurr U k etant etant donn´ donnée ee k +1 = t k k + h. La valeu k +1 ´ par (2.30), on aura u(t k k +1 )

− U + = u(t + ) − u(t ) − hF (t , u(t ), h) k k 1

k k 1

k k

k k

k k

etant e´ tant donné que U k ese. ese. On note k = u(t k k ) par hypoth` e(t k k , h) = u(t k k +1 )

− U +

k k 1

l’erreur appel appelee e´ e erreur locale et e(t k k , h) = u(t k k +1)

− u(t ) − hF (t , u(t ), h). k k

k k

k k

(2.34)

Il s’agit s’agi t donc de l’erreur l ’erreur que l’on l’ on fait en avanç ant d’un pas avec le sch´ sch ema, éma, partant de U de U k ema d’Euler (2.31) par exemple, avec F (t k ema k = u (t k k ). Pour le sch´ k , u(t k k ), h) = ′ etant la solution exacte, u(t k f (t k k , u(t k k )) = u (t k k ), u(t k k ) étant k ) + hF (t k k , u(t k k ), h) est sur la droite passant par u par u (t k ` la courbe u courbe u (t ) comme illustré sur la figure k ) et tangente a 2.1. On note t note t k et on introduit la quantité τ(t , h) = e(t , h)/h, c’est-à-dire a-dire k = t et

τ(t , h) =

u (t + + h) h

− u(t ) − F (t , u(t ), h)

(2.35)

erreur de discr etisation qui qui est est l’err l’erreu eurr loca locale le divi divis´ sée e e par par h : cette quantit´ quantité, e , appe appel´ lée ee erreur ´ ´ locale, locale, est obtenue en injectant la solution exacte u(t k ema ema (2.30). k ) dans le sch´ La solution approchée ee est solution du schéma ema tandis que la solution exacte est solution a` l’erreur τ (t , h) près. es. 36

Ordre d’un sch´ ema, consistance, stabilit´ e et convergence

e

u(t)

2

e

e1

3

e0

t0

t1

t2

t3

t4

t

FIG . 2.1 – Illustration de l’erreur de discrétisation etisation locale ek = e (t k k , h) pour le schéma ema d’Euler.

Definition e´ finition 3 Le sch´ sch ema 0 , s’il existe une ´ (2.30) est d’ordre p pour un entier p ∗ constante constante C > , 0 < h < > 0 et h > 0 , tels que pour tout t ]t 0 , T [ , , et pour tout h > 0 , 0 ∗ h , l’erreur de discr etisation locale τ(t , h) d efinie par (2.35) v erifie ´ ´ ´ ´ ´

≥

∈ p

|τ(t , h)| ≤ Ch .

On ´ On ´ ecrira que τ(t , h) = O(h p). Le sch´ schema ´ est dit consistant consistant si p

(2.36)

≥ 1.

Calcul pratique de l’ordre d’un sch ema e´ ma Pour déterminer eterminer l’ordre d’un schéma, ema, il convient d’opérer erer un développement eveloppement h de la quantité τ(t , h). Considérons de Taylor par rapport a` h de erons d’abord u(t + h) h

− u(t )

et par un dévelo eveloppe ppemen mentt de Taylor aylor (en suppo supposa sant nt que u que u(t ) est p est p + 1 fois fois continˆ continûment ument dérivable) erivable) u(t + + h) h

− u(t ) = u ′(t ) + h u ′′(t ) + · · · + h − u( )(t ) + O(h ). p 1

p! p!

2

Or, u Or, u (t ) est solution de l’EDO et donc u ′ (t ) = f (t , u(t )) )) 37

p

p

(2.37)


et ensuite u ′′ (t ) =

∂ f ∂ f ∂ f ∂ f (t , u(t )) )) + (t , u(t )) )) u ′ (t ) = (t , u(t )) )) + (t , u(t )) )) f (t , u(t )). )). ∂t ∂u ∂t ∂u

Evidemment, les dériv´ erivées ees suivantes de u s’expriment egalement e´ galement en fonction des dériv´ erivées ees parti partiell elles es de f (t , u(t ), la comp comple lexi xit´ té de l’expr l’expres essio sion n augme augmenta ntant nt consi consid´ dérablement erablement avec l’ordre l’ordre de la d eriv´ e´ rivée. ee. S’arrêtant etant a` l’ordre 2, on peut ecrire e´ crire u(t + h) h

− u(t ) = f (t , u(t )))) + h

2





∂ f ∂ f (t , u(t )) )) + (t , u(t )) )) f (t , u(t )) )) + O(h2 ). ∂t ∂u

(2.38) h et bien sûr La fonction F fonction F (t , u(t ), h) peut egalement e´ galement etre eˆ tre développ´ eveloppée ee par rapport rapport a` h et ur

∂F h2 ∂2 F F (t , u(t ), h) = F (t , u(t ), 0) + h (t , u(t ), 0) + (t , u(t ), 0) + ∂h 2 ∂h2 h p−1 ∂ p−1 F + (t , u(t ), 0) + O(h p) (2.39) p 1 − ( p 1)! ∂h

···

···

−

L’erreur de discrétisation etisation locale τ(t , h) définie efinie par (2.35) etant e´ tant la différence erence entre (2.38) et (2.39), on obtient jusqu’à l’ordre 2 en h en h

τ(t , h) = ( f (t , u(t )) )) F (t , u(t ), 0)) ∂ f 1 ∂ f + h (t , u(t )) )) + (t , u(t )) )) f (t , u(t )) )) ∂u 2 ∂t



−

+ O(h2 ).

−

∂F (t , u(t ), 0) ∂h



(2.40)

On peut bien sûr ur expliciter les termes d’ordres supérieurs erieurs de τ(t , h), au prix de calcul fastidieux. De l’expression ci-dessus ci-dessus on peut conclure que le sch´ sch éma ema est au moins d’ordre 1, c’est-à-dire a-dire τ (t , h) = O(h), si et seulement si F (t , u(t ), 0) = f (t , u(t )) ))

(2.41)

Ensuite il est au moins d’ordre 2, si en plus

∂F 1 (t , u(t ), 0) = ∂h 2





∂ f ∂ f (t , u(t )) )) + (t , u(t )) )) f (t , u(t )) )) . ∂t ∂u

(2.42)

Evaluant l’ordre des deux schémas emas (2.31) et (2.33) introduits ci-dessus. Pour le schéma ema d’Euler (notant t (notant t k ), k = t ), F (t , u(t ), h) = f (t , u(t )) )) 38


)) et le sch´ et donc donc F (t , u(t ), 0) = f (t , u(t )) schéma e ma est est d’or d’ordr dree 1, mais mais il n’es n’estt pas pas d’or d’ordr dree h . Prenons le schéma 2, car ∂ F (t , u(t ), h)/∂h = 0, F 0, F ne ne dépendant ependant pas de h. ema (2.33) associé a` la formule des trapèzes, ezes, alors (t ( t k k = t , t k k +1 = t + h) 1 F (t , u(t ), h) = ( f (t , u(t )) )) + f (t + h, u(t + h)) ) . 2 Or, par Taylor





∂ f ∂ f f (t + h, u(t + h)) = f (t , u(t )) )) + h (t , u(t )) )) + (t , u(t )) ))u ′ (t ) + 0(h2 ) ∂t ∂u

)) car pour evaluer Or, a` nouveau u ′ (t ) = f (t , u(t )) e´ valuer l’ordre on injecte la solution exacte de l’EDO dans le schéma ema et donc F (t , u(t ), h) = f (t , u(t )) )) +

h 2





∂ f ∂ f (t , u(t )) )) + (t , u(t )) )) f (t , u(t )) )) + O(h2 ) ∂t ∂u

Par conséquent, equent, dans ce cas on obtient

τ(t , h) = O(h2 ), en comparant l’expression ci-dessus avec (2.38). Nous allons préciser eciser maintenant de quelle manière ere l’ordre d’un schéma ema est lié a` l’erreur entre la solution exacte et la solution approchée. ee.

Stabilit´ Stabilite´ et convergence On reprend le schéma ema a` un pas U k k +1 = U k k + hF (t k k , U k k , h),

U 0 = u0 ,

k = = 0, 1, 2,...,

(2.43)

et on suppose qu’à chaque pas en temps une erreur εk s’introduit. Il peut s’agir par exemple d’erreurs d’erreurs d’arrondis des nombres r´ r eels e´ els lors d’une mise en œuvre sur ordinateur. ordinateur. Donc, au lieu l ieu du sch ema e´ ma exacte, on aura aura plut ot oˆ t V k k +1 = V k k + hF (t k k , V k k , h) + εk ,

V 0 = u0 ,

k = = 0, 1, 2,...,

(2.44)

Definition e´ finition 4 Le sch´ sch ema (2.30) est dit stable, s’il existe une constante M > 0 et ´ ∗ h > 0 tels h, 0 < h < h∗ et pour toute suite εk on ait 0 tels que pour tout h, 0

−

N 1

|U − V | ≤ M ∑ |ε | N

N

k =0

pour tout N tel que t N N = t 0 + Nh

∈]t , T [. 0

39

k


Un schéma ema est donc stable si les erreurs introduites à chaque itération eration du schéma ema ne font au pire que s’additionner. Dans quelles conditions un schéma ema (2.30) (2.30) est-il stable ?

Th´ Theor` e´ oreme e` me 6 Pou Pourr que que le sch schema em (2.30) 0) soit soit stab stable le,, il suffi suffitt qu’il qu’il exis existe te une une cons consta tant ntee ´ a (2.3 Λ telle que pour tout t [t 0 , T ] , , pour tout v, w R et pour tout h [0, h∗ ] on ait

∈ ∈ |F (t , v, h) − F (t , w, h)| ≤ Λ|v − w|

∈

(2.45)

Pour démontrer emontrer ce théor` eorème, eme, formant a` partir de (2.43) et (2.44))

|U − V + h(F (t , U , h) − F (t , V , h)) − ε | ≤ |U − V | + h|F (t , U , h) − F (t , V , h)| + |ε | ≤ (1 + hΛ) |U − V | + |ε | =

|U + − V + | k k 1

k k 1

k k

k k

k k

k k

k k

k k

k k

k k

k k

k k

k k

k k

k k

k k

k

k

k

si F (t , u, h) v´ verifie e´ rifie l’inegalit´ e´ galité (2.45). Mais on peut ecrire e´ crire l’in´ l’inegalit´ e´ galité c-dessus en rempl rem plac aç ant k ant k par k par k 1 etc. et c. de d e faç on a` obtenir

−

k

k +1

j

|U + − V + | ≤ (1 + hΛ) |U − V | + ∑ (1 + hΛ) |ε − |. k k 1

k k 1

0

0

k j

j=0

Or U Or U 0 = V 0 d’après es hypothèse ese et tenant compte que (1 + hΛ) j = k + 1) 0 j k , on peut ecrire écrire (posant N (posant N =

≤ ≤

N

− 

N

k 1

si

N 1

|U − V | ≤ (1 + hΛ) ∑ |ε | = N

≤ (1 + hΛ) +

j

.

j 0

Or, de l’inégalit´ egalité bien connue

(1 + hΛ) on déduit eduit

(1 + hΛ) N

≤e Λ h

≤ e Λ ≤ e( − hN

T t 0 )Λ

pour pour tout tout N tel tel que que t N = t 0 + Nh ]t 0 , T [. Par cons´ conséquen equent, t, la consta constante nte M M = = e(T −t 0 )Λ dans la définition efinition de stabilité. e. On remarque ici que la condition (2.45) est une condition de type Lipschitz pour F qui définit la fonction fonction F qui efinit le schéma. ema. Nous avons supposé que pour la fonction f qui f qui définit efinit l’EDO, la condition de Cauchy-Lipschitz (2.3) est vérifi´ erifiée. ee. Donc, le schéma ema d’Euler explicite (2.31) et le schéma ema des trapèzes ezes (2.33) sont stables, etant e´ tant donn´ données ees les expres expressions sions de F de F (t , u(t ), h) dans ces cas en fonctio fonction n de f (t , u(t )) )).

∈

40


Par la suite nous allons aborder la question de convergence, c’est-à-dire nous allons evaluer e´ valuer l’erreur entre la solution exacte et la solution numérique erique pour un temps t temps t ]t 0 , T [. Il convient alors d’introduire un pas de temps h N de faç on a` ce que t t 0 . t N = t 0 + Nh N h N = t , c.-à-d. a-d. h N = (2.46) N N pas avec Dans ce cas U cas U N , la solution numérique erique obtenue en mettant en œuvre N pas le schéma ema numérique, erique, est une approximation de la vraie solution sol ution au temps t temps t et et on note E (t , h N ) = U N u(t ) (2.47)

∈ ∈

−

−

l’erreur. l’erreur. On suppose que le sch ema e´ ma (2.43) est d’ordre p et d’après es la définition efinition de l’ordre du schéma, ema, on peut ecrire e´ crire pour la solution exacte u(t k k +1) = u(t k k ) + h N F (t k k , u(t k k ), h N ) + e(t k k , h N )

(2.48)

avec e(t k k , h N ) = h N τ(t k k , h N )

(2.49)

p

l’erreur locale. On rappelle que τ (t k ema ema est d’ordre p. k , h N ) = O(h N ) si le sch´

Th´ Theor` e´ oreme e` me 7 Si le sch ema (2.43) est stable stable et d’ordre d’ordre p avec avec p entier positif, alors ´ (2.43) ¯ > 0 et il existe une constante K > , 0 < h N < ¯h 0 et h > 0 tels 0 tels que pour tout h N , 0 p N

| E E (t , h ) ≤ K h − u(t ). En particulier, si p ≥ 1 , donc si le sch´ schema ´ est consis N

avec E (t , h N ) = U N tant, alors

lim E (t , h N ) = 0.

→0

h N

On note que h N

a` N → ∞ par (2.46). → 0 est ´ equivalent ´

Avant de démontrer emontrer ce r esultat e´ sultat fondamental, il convient d’observer que l’erreur sera sera d’autan d’autantt plus plus petite que l’ordre l’ordre du sch´ schéma ema sera elev´ e´ levé, e, ci qui démontr emontree l’int´ l’intérˆ erêt et d’employer précis´ ecisément ement des schémas emas d’un ordre p 2. L’ordre L’ordre 1, c’est-` c’est-a-dire a` -dire la consistance, consistance, est n´ n ecessaire e´ cessaire pour la convergence. Pour la démonstration, emonstration, on utilise utilise la definition e´ finition 4 de stabilité, e, en prenant

≥

V k k = u(t k k )

et

εk = e(t k k , h N )

car alors (2.48) et (2.44) sont identiques. Si le schéma ema est stable, on aura d’après es la définition efinition 4

−

N 1

−

N 1

|U − u(t )| ≤ M ∑ |e(t , h )| = Mh ∑ |τ(t , h )|, = = N

N

k k N

k 1

N

k 0

41

k k N


> O, O, et ˜h > 0 tel que si 0 < h N < h∗ . Si le schéma ema est d’ordre p, alors il existe C existe C > si 0 < h N < ˜h p τ(t k k , h N ) Ch N

|

|≤

p ∗˜ ¯ car τ(t k k , h N ) = O(h N ). Donc, prenons h = min(h , h), alors si 0 p N

p N

≤ h ≤ h¯ N

p N

K h |U − u(t )| ≤ MN h Ch ≤ M (T − t )Ch = Kh = t + Nh ∈]t , T [, ce qui car N car N h ≤ (t − t ) ≤ (T − t ) étant etant donné que t que t = t = N

0

N

N

N

0

0

0

N N

0

N

achève eve la démonstration emonstration du du th´ theor` e´ orème. eme.

2.2. 2.2.2 2

Les Les sch schemas e´ mas de Runge-Kutta Runge-Kutta

On peut construire des schémas emas a` un pas d’ordres d’ordres de plus en plus elev´ e´ levés, es, en emboˆıtant ıtant en quelque q uelque sorte des formules d’int d’ intégration e´ gration numérique. erique. Ces schémas, emas, schemas de Runge-K Runge-Kutta utta,, se construisent de façon connus connus sous sous le nom de de sch´ ¸on suivante. suivante. ´ Un schéma ema a` un pas est une règle egle qui permet permet de d´ d eterminer e´ terminer la solution numérique erique au temps t temps t k erique au temps t temps t k k +1 = t k k + h en fonction de la solution numérique k . On introduit q duit q temps temps intermédiaires ediaires t k k , j = t k k + c j h,

0

≤ c ≤ 1, j

j = 1,

· · · q.

(2.50)

La solution exacte de l’EDO est telle que u(t k k +1 ) = u(t k k ) +

Z

t k +1

t k

f (t , u(t )) ))dt

et approchant approchant l’int´ l’integrale e´ grale par une somme (une formule formule d’int egration) e´ gration) on ecrit e´ crit pour la solution approchée ee q

U k k +1 = U k k + h ∑ bi f (t k k + ci h, U k k ,i ) i=1

pour des coefficients bi a` préciser eciser et U et U k e´ signe la la solution approchée ee de la k, i d esigne solution exacte u exacte u (t k k + ci h). Or, u(t k k + ci h) = u(t k k ) +

Z

t k +ci h

t k

f (t , u(t )) ))dt

et en employant a` nouveau des formules d’intégration egration numérique, erique, on d etermine e´ termine , q par des expressions U k k ,i , i = 1,

···

q

U k k ,i = U k k + h ∑ ai j f (t k k + c j h, U k k , j ) j=1

42

Les sch´ emas de Runge-Kutta Runge-Kutta

pour des coefficients ai j a` préciser. eciser. En résum´ esumé, e, un schéma ema de Runge-Kutta se présente esente de la l a faç on suivante. On pose po se U 0 = u 0 la condition initiale de l’EDO et = 0, 1, 2, ensuite pour k pour k =

···

q

U k k ,1 = U k k + h ∑ a1 j f (t k k + c j h, U k k , j ) j=1 q

U k k ,2 = U k k + h ∑ a2 j f (t k k + c j h, U k k , j ) j=1

. ..

(2.51) q

U k k ,q = U k k + h ∑ aq j f (t k k + c j h, U k k , j ) j=1

et finalement

q

U k k +1 = U k k + h ∑ bi f (t k k + ci h, U k k ,i )

(2.52)

i=1

, q, bi , i = est la solution approchée ee au temps t k k +1 . Les coefficients c j , j = 1, 1, , q ainsi que ai j , i = 1, , q, j = 1, , q sont a` déterminer eterminer en fonction de l’ordre du schéma ema souhaité. e. Observant d’abord que

···

···

···

···

q

U k k, i = U k k + h ∑ ai j f (t k k + c j h, U k k , j ) j=1

est censée ee etre eˆ tre l’approximation de la solution au temps t k notant t k k + ci h ; notant t k = t et injectant la solution exacte on a u(t + + ci h) h

− u(t ) =

q c j h, u(t + c j h)) + τi (t , h). ∑ ai j f (t k k +

j=1

Pour Pour que cette cette expr express ession ion soit soit consis consistan tante, te, c’est-` c’est-à-dire a-dire que que l’err l’erreu eurr de discr´ discrétisation etisation locale partielle τi (t , h) 0 quand h quand h 0, il suffit que

→

→

q

)) = ci f (t , u(t )). )). ∑ ai j f (t , u(t ))

j=1

En effet,

u(t + + ci h) h h→0 lim

et

q

− u(t ) = c u ′(t ) = c f (t , u(t )))) i

i

q

)) quand ∑ ai j f (t + c j h, u(t + c j h)) → ∑ ai j f (t , u(t ))

j=1

j=1

43

h

→ 0.


U k De même, eme, injectant la solution exacte dans (2.52) ((U e´ tant la solution apk +1 etant prochée ee au temps t temps t k notant t k ), k + h et notant t k = t ), u(t + + h) h

− u(t ) =

q

+ ci h, u(t + ci h)) + τ(t , h) ∑ bi f (t +

i=1

)) quand h et le memb membre re a` gauch gauchee tenda tendant nt vers vers f (t , u(t )) locale τ(t , h) 0 quand h quand h 0 si

→

→

→ 0, l’err l’erreu eurr de disc discr´ rétisation etisation

q

f (t , u(t )) )) =

)). ∑ bi f (t , u(t )).

i=1

Le fait que les erreurs de discrétisation etisation locales tendent tendent vers z ero e´ ro quand h quand h tend tend vers zéro ero exprime la consistance du schéma ema : par conséquent, equent, d’après es ce qui préc` ecède, ede, le schéma ema est au moins d’ordre 1, si q

q

∑ bi = 1

et

ci =

i=1

∑ ai j , i = 1, · · · , q.

(2.53)

j=1

Essayons d’ d’ecrire e´ crire les les conditions conditions pour que le sch´ schéma ema soit au moins mo ins d’ordre 2. L’expression (2.52) peut peut s’´ s’ ecrire e´ crire U k k +1 = U k k + hF (t k k , U k k , h) avec

q

F (t k k , U k k , h) =

ci h, U k k ,i ). ∑ bi f (t k k +

i=1

avec U k es es par (2.51). On pose t pose t k et il faut pouvoir vérifier erifier les conditions k ,i donn´ k = t et (2.42) pour que le schéma ema soit d’ordre 2. Soit donc (en injectant injec tant la solution exacte dans l’expression ci-dessus) q

F (t , u(t ), h) = ∑ bi f (t + + ci h, u(t + ci h)), )), i=1

avec d’après es (2.51) q

u(t + ci h) = u(t ) + h ∑ ai j f (t + + c j h, u(t + c j h)), )), j=1

ce qui donne lieu a` l’expression q



q



F (t , u(t ), h) = ∑ bi f t + + ci h , u(t ) + h ∑ ai j f (t + + c j h, u(t + c j h)) . i=1

j=1

44


Par conséquent, equent, q

∂F (t , u(t ), 0) = ∂h

∑ bici

i=1 q

∑ bi

+

i=1

∂ f (t , u(t )) )) ∂t





q

∂ f (t , u(t )) )) ∑ ai j f (t , u(t )) )) . ∂u j=1

Les relations (2.53) doivent etre eˆ tre vérifi´ erifiées, ees, car pour etre eˆ tre d’ordre 2 il faut etre eˆ tre q d’ordre 1, donc ∑ j=1 ai j = c i et par conséquent equent le schéma ema est d’ordre au moins 2, c’est-à-dire a-dire la condition (2.42) est vérifi´ erifiée, ee, si en plus de (2.53) on ait q

1

∑ bici = 2 .

(2.54)

i=1

Evidemment, d’écrire ecrire les relations pour des ordres supérieurs erieurs est de plus en plus fastidieux et on se contente d’avoir démontr´ emontré les conditions d’ordre 2. On peut montrer par exemple que le schéma ema est au moins d’ordre 3 si en plus des relations (2.53) et (2.54) les conditions suivantes entre entre les coefficients coefficients q

∑

i=1

bi c2i

q

1 = 3

∑

et

q

∑ bi ai j c j =

i=1 j=1

1 6

(2.55)

sont vérifi´ erifiées. ees. On représente esente en gén´ enéral eral les coefficients coefficients d’un sch´ sch ema e´ ma de Runge-Kutta sous forme d’un tableau c1 c2 .. .

a11 a21 .. .

a12 a22 .. .

··· ···

a1q a2q .. .

cq

aq1 b1

aq2 b2

··· ···

aqq bq

Si dans les formules (2.51) les coefficient a coefficient a i j = 0, i

≤ j ≤ q, alors

U k k ,1 = U k k et pour i pour i

≥2

−

i 1

U k k, i = U k k + h ∑ ai j f (t k k + c j h, U k k , j ) j=1

, i 1 et et par conséquent equent les U k ependent ependent que des U k k ,i ne d´ k , j , pour j = 1 , peuvent peuvent par conséquent equent etre eˆ tre détermin´ eterminés es aisément ement au fur et a` mesure. On parle

··· −

45


ema explicite. alors d’un sch d’un sch´ ´ explicite.

Exemples : On cherche a` déterminer eterminer tous les schémas emas de Runge-Kutta explicites d’ordre au moins 2 avec q = 2. Le tableau ci-dessus devient alors c1 c2

0 a21 b1

0 0 b2

D’après es les conditions (2.53) et (2.54), le schéma ema est au moins d’ordre 2, si c1 = 0, a21 = c2 , b1 + b2 = 1, b2 c2 =

1 . 2

Prenons par exemple b exemple b 1 = 0, c 0, c2 = 21 et b et b 2 = 1 (c (c1 = 0 et a et a 21 = 21 ). Alors U k k ,1 = U k k ,

h U k f (t k k ,2 = U k k + k , U k k ,1) 2

et on trouve le schéma ema



1 U k h, U k k +1 = U k k + h f t k k + k ,2 2



et donc U k k +1 = U k k + hF (t k k , U k k , h) avec



h 1 F (t k h, U k f (t k k , U k k , h) = f t k k + k + k , U k k ) 2 2



modifie´ . Pour c qui est le schéma ema dit d’Euler dit d’Euler modifi´ Pour c 2 = 1, a 1, a 21 = 1 et b 1 = b 2 = 1/2 on Heun avec trouve le schéma ema dit de dit de Heun avec 1 1 F (t k )). k , U k k , h) = f (t k k , U k k ) + f (t k k + h, U k k + h f (t k k , U k k )). 2 2 Un schéma ema classique de Runge-Kutta explicite avec q avec q = 4, dont on peut montrer qu’il est d’ordre 4, est donné par le tableau 0 1/2 1/2 1

0 1/2 0 0 1/6

0 0 1/2 0 1/3 46

0 0 0 1 1/3

0 0 0 0 1/6


ce qui donne lieu au schéma ema U k k ,1 = U k k h

             

U k f t k k ,2 = U k k + k , U k k ,1 2 h h , U k k, 2 U k = U + f t + , k 3 k k k k 2 2 h , U k k ,3 U k k ,4 = U k k + h f t k k + 2 h h U k f t k , U k k ,2 k+ 1 = U k k + k , U k k ,1 + 2 f t k k + 6 2 h , U k k ,3 + f t k k + + 2 f t k k + h, U k k ,4 2

Pour conclure ce chapitre, il convient d’observer, que pour les schémas emas de Runge F est vérifi´ Kutta la condition (2.45) de Lipschitz pour F est erifiée, ee, dans la mesure où F est donnée, ee, par des combinaisons plus ou moins complexes, en fonction de f (à condition bien sur uˆ r que f satisfait aux conditions de Cauchy-Lipschitz Cauchy-Lipschitz garantissant l’existence et et l’unicit´ l’unicite´ de la solution de l’EDO). Ces schémas emas sont donc stables et les résultats esultats de convergence ci-dessus s’appliquent.

47


48

Chapitre 3 Resolution e´ solution num´ numerique e´ rique directe de syst` systemes e` mes lin´ lineaires e´ aires 3.1 3.1


Supposons donné le problème eme suivant : on souhaite souhai te connaˆıtre ıtre la l a temp erature e´ rature d’une barre métallique etallique chauffée ee a` ses deux extrémit´ emités es et plongée ee dans une pièce ece elle-même eme a` une une temp´ température erature donnée. e e. La barr barree est est assi assimil´ milée ee a` un segm segmen entt de droi droite te b.. On suppose dans cette simplification que la dont les extrémit´ emités sont notées a ees a et et b température erature de la barre ne dépend epend que de x et T et et on cherche T ( x), a < x < b, b , avec T avec T (a) = T a et T (b) = T b (o u` T a et T b sont les températures eratures aux extrémit´ emités). es). La température erature environnante est désign´ esignée ee par T e . Il y a perte de chaleur due a` la convection de l’air que l’on modélise elise par une fonction r ( x). On note κ le le coefficient de diffusion thermique que l’on suppose constant. Le probl eme e` me est régi egi par une equation e´ quation différentielle erentielle Te

Ta

Tb b

a

x

Te

F IG . 3.1 – Schéma ema de principe d’une barre mince chauff chauff ee e´ e aux extrémit´ emités. es.

d 2 T ( x) κ + r ( x)(T ( x) dx 2

−

− T ) = 0, a < x < b, T (a) = T , T (b) = T . e

a

b

(3.1)

Sauf dans des cas particuliers, il n’est en gén´ enéral eral guère ere possible de déterminer eterminer une solution analytique de cette équation. equation. Donc, même eme pour ce cas relativement 49

R´ esolution num´ erique directe de syst` emes lin´ eaires

simple une approche “numérique” erique” s’avère ere nécessaire ecessaire afin de trouver une solution approchée ee de la température erature T ( x) le long de la barre. Une méthode ethode classique, er finies , consiste a` appro conn connue ue sous sous le nom nom des des diff erence approche cherr les dériv´ erivées ees par des ´ ´ encess finies quotients aux différences. erences. La technique utilise la formule de Taylor. aylor. Prenons par exemple une fonction f ( x) “suffisamment dérivable”, erivable”, alors la formule de Taylor donne f ( x

± h) = f ( x) ±

d f h 2 d 2 f h ( x) + ( x) dx 2 dx 2

±

h3 d 3 f ( x) + O(h4 ). (3.2) 3 6 dx

f soit 4 fois continûment Pour que cette expression ait un sens, il suffit que f soit ument dérivable. erivable. Dans l’expres l ’expression sion ci-dessus la l a notion “grand O” est employée ee : d’une p manière ere gén´ enérale O erale O (h ) est une expression t. q. 0

≤  

O(h p) h p

 ≤ 

M , quand

h

→ 0,

> 0 indépendant pour M pour M > ependant de h de h et et p. En d’autres termes, O termes, O (h p) est une expression expression p de l’ordre de grandeur h grandeur h . De l’expression (3.2) en déduit eduit par exemple d f f ( x + h) ( x) = dx h

− f ( x) + O(h)

et l’erreur O(h) dans l’expression est d’autant plus faible que h sera petit. En combinant les expressions expressions f ( x + h), f ( x) et f ( x + h) on montre que d 2 f f ( x = dx 2

− h) − 2 f ( x) + f ( x + h) + O(h ). 2

h2

(3.3)

On vérifiera erifiera cette expression à titre d’exercice. La méthode ethode des différences erences finies consiste consiste a` faire abstraction de l’erreur l ’erreur dans ces expressions. expressions. On ecrira e´ crira par exemple à la place de la “vraie” dériv´ erivée ee seconde le quotient aux différences erences f ( x

− h) − 2 f ( x) + f ( x + h) . h2

Revenons Revenons a` notre problème eme de départ epart et a` la détermination etermination de T de T ( x). Dans les formules de Taylor, un pas discret h discret h intervient intervient et l’idée ee fondamental est de diviser −a l’intervalle [ a, b] en sous-intervalles précis´ ecisément ement de longueur h longueur h.. On pose h pose h = nb+ 1 et on définit efinit les n les n + 2 points x j = a + jh j h,

≤ j ≤ n + 1. a et x x + = b et b et il y a n x , 1 ≤ j ≤ n à l’intérieur Ainsi, x Ainsi, x = a et a n points points x erieur de l’intervalle 0

0

n 1

j

ethode d’approximation d’approximation de la d eriv´ e´ rivée ee seconde permet alors de trouver [a, b]. La méthode 50

Motivation

j n, de la température une approximation de T ( x j ), 1 erature au points x j (sachant que T que T ( x0 ) = T a , T ( xn+1 ) = T b et les températures eratures aux bords sont des paramètres etres du problème). eme). Utilisant la relation (3.3) on peut donc ecrire e´ crire l’approximation

≤ ≤

d 2 T ( x j ) dx 2

≈ T ( x − ) − 2T h( x ) + T ( x + ) , j 1

j 1

j

2

car x car x j±1 = x j h. On note T note T j j l’approximation de T ( x j ) et on obtient le système eme suivant, suivant, a` partir de (3.1) et en écrivant ecrivant le quotient aux diff erences e´ rences a` la place de la dériv´ erivée ee seconde pour tous les points x j à l’intérieur erieur de l’intervalle [ a, b]

±

−κ T − − 2T + T + h 2



j 1

j

j 1



+ r ( x j )T j = r ( x j )T e ,

j = 1,

· · · , n.

(3.4)

La fonction r fonction r ( x j ) est donnée ee dans l’expression ci-dessus ainsi que T 0 . Les inconnues du système eme ci-dessus sont les T j j , a` savoir les approximations de la solution exacte T ( x j ) aux points x j . Evaluant l’équation equation pour j = 1, on voit apparaˆıtre ıtre erature imposée ee en x en x = a, ere analogue T 0 = T ( x0 ) = T a, qui est la température a, et de manière b . Les inconpour j = n la température erature T n+1 = T ( xn+1 ) = T b imposée ee en x = b. nues sont ainsi les approximations T j des valeurs de la température erature aux points x j , j = 1, n, a` l’intérieur erieur du domaine [a, b]. Les T j j peuvent etre eˆ tre considér´ erées ees T . T . A partir de (3.4), et en tecomme les composantes d’un vecteur que l’on note  nant compte que pour j = 1 et j = n + 1 les données T ees T a et T et T b interviennent, il est T est facile de voir (exercice) que le vecteur  est solution du système eme linéaire eaire

···

A T =  B

(3.5)

pour la matrice

A =

  

2 + s( x1 ) 1 1 2 + s( x2 ) .. .

−

−

(0)

−1

(0)

−

−

..

.. . . 1 2 + s( xn−1 ) 1 1 2 + s( xn )

−

et le second membre qui contient les paramètres etres du problème eme

 = B

  

s( x1 )T e + T a s( x2 )T e .. . s( xn−1 )T e s( xn )T e + T b 51

  

,

  

(3.6)

(3.7)

R´ esolution num´ erique directe de syst` emes lin´ eaires 2

avec s avec s( x j ) = hκ r ( x j ), j = 1, n. La matrice A matrice A,, avec n avec n lignes lignes et n et n colonnes, colonnes, est dite tridiagonale. En effet, seuls les el´ e´ léments ements sur la diagonale ainsi que leurs voisins immédiats ediats a` droite et a` gauche sont non nuls. Plus n est grand, plus la solution approchée ee sera précise ecise et etant e´ tant donnée ee la taille du système eme linéaire eaire il faut imaginer une une m´ methode e´ thode numérique erique en vue de sa résolution. esolution. Les paramètres etres du problème eme B. B. L’idée (T a , T b , T e) n’interviennent que dans l’expression du second membre  ee est alors de ne pas résoudre esoudre directement le système eme (3.5) pour un second membre donné, e, mais plut plutot oˆ t d’opérer erer une décomposition ecomposition de A une fois pour toutes, afin de pouvoir aisément ement résoudre esoudre les systèmes emes successifs lorsqu’on fait varier les paramètres etres du problème. eme.

···

Decomposition e´ composition LU d’une matrice tridiagonale tridiagonale

3.2

On ecrit e´ crit d’une manière ere gén´ enérale erale une matrice tridiagonale n tridiagonale n réels eels ou complexes) sous la forme

A =

  

a1 b2

c1 a2 .. .

(0)

× n (à coefficients

   −

(0)

c2 .. .

..

.

(3.8)

bn−1 an−1 cn 1 bn an

L’idée ee est d’écrire ecrire A sous la forme du produit de deux matrices LU . U . La matrice L aura L aura la particularité d’être etre triangulaire inférieure, erieure, c’est-à-dire a-dire ses el´ e´ léments ements au U sera triangulaire supérieure, dessus de la diagonale sont nuls, tandis que U erieure, c.a-d. a` -d. ses el´ e´ léments ements en dessous de la diagonale sont nuls. La notation L et U est U de “upper”. On gén´ enéralement eralement utilisée, L ee, L étant etant la premi` première ere lettre de “lower” et U de suppose que ai = 0, i = 1,



· · · , n, c = 0, i = 1, · · · , n − 1, b = 0, i = 2, · · · , n. i

i

(3.9)

On cherche L cherche L et et U sous la forme

L =

  

1

β2

(0)

1 .. .

(0) ..

.

βn−1

1

βn

1

  

, U =

52

  

α1

γ 1 α2

γ 2 ..

(0)

.

   − 

(0) ..

.

.

αn−1 γ n 1 αn

(3.10)

D´ ecomposition LU d’une d’une matrice tridiagonale

On opérant erant le produit matriciel LU matriciel LU on on obtient (exercice)

LU LU =

  

α1 γ 1 α1 β2 α2 + γ 1 β2 ..

γ 2 ..

.

(0) ..

.

.

αn−2 βn−1 αn−1 + γ n−2 βn αn−1 βn

(0)

  

−1

γ n−1 αn + γ n−1 βn

.

(3.11) A on a bien En identifiant les el´ e´ léments ements de la matrice LU avec avec ceux de la matrice A on A = LU si si

βi = βn =

bi

αi−1 bn

αn−1

α1 = a1 ,

γ 1 = c1

,

αi = ai

− γ − β , i 1 i

,

αn = an

− γ − β .

γ i = ci ,

i = 2, 3,

· · ·n − 1

(3.12)

n 1 n

L et U U ,, a` condition que α i = 0. Les Ces relations déterminent eterminent les coefficients de L et relations ci-dessus s’apparentent s’apparentent a` un algorithm algorithmee dans la mesure mesure où si l’on l’o n connaˆ co nnaˆıt ıt αi−1 on peut déterminer eterminer β i dont on a besoin pour calculer α i etc. ; par ailleurs ailleurs les γ i = ci . L’égalit´ egalité α 1 = a 1 initialise le calcul et cet algorithme peut donc être etre mis en œuvre sur un ordinateur. Pour certaines classes classes de matrices il est assur´ assur e´ que lors des etapes e´ tapes successives αi = 0.





Th´ Theor` e´ oreme e` me 8 On suppose que les coefficients de la matrice A donn ee ´ par (3.8) avec (3.9) sont tels que

|a | > |b |, |a | ≥ |b + | + |c − |, i = 2, · · · , n − 1, |a | ≥ |c − |. 1

2

i

i 1

i 1

n 1

n

On parle dans ce cas d’une matrice à diagonale dominante. Alors

|β | < 1, i = 2, · · · , n,

et αi = 0, i = 1,

· · · , n. Ce résultat esultat se démontre emontre par recurren e´ currence ce ; en effet effet α = a = ese  0 par hypothèse sur les coefficients de la matrice A et par la dominance diagonale on trouve que |β | = |b |/|α | < 1. On fait alors l’hypothèse  0, i = 1, · · · , m ese de récurrence ecurrence : α = et |β | < 1, i = 1, · · · , m + 1 (avec m ≤ n − 2). De l’expression de α + donnée ee par (3.12) (3.12) on d eduit e´ duit la minoration |α + | ≥ |a + | − |γ ||β + | > |a + | − |γ | car |β + | < 1 par hypothèse ese de récurrence. ecurrence. Or, par la dominance diagonale on aura |a + | − |γ | ≥ |b + | et finalement on en déduit eduit que |α + | > |b + |. Par hypothèse b ese b + =  0 et donc α + = 0 et aussi |β + | = |b + |/|α + | < 1. Le i



1

2

2

1

1

i

m 1

i

m 1

m 1

m

m 1

m 1

m 1

m 1

m

m 2

m 2

m 1

m 1

m 2

théor` eorème eme est ainsi démontr´ emontré.

53

m 2

m 2

m 1

m


Qu’a-t-on gagné en ecrivant A e´ crivant A = LU ? ? Précis´ ecisément, ement, admettons qu’on cherche a` résoudre esoudre Ax = d avec d avec d un un vecteur second membre donné et A et A une une matrice de la forme (3.8) avec (3.9). On peut donc ecrire e´ crire LUx = d , ou encore de manière ere equivalente e´ quivalente sous forme de deux syst` syst emes e` mes linéaires eaires avec des matrices triangulaires L triangulaires L et et U U Ly = d ,

U x = y.

On note d i les coefficients du vecteur d et et les coefficients yi de y s’obtiennent aisément, L ement, L ´ étant etant sous la forme (3.10). En effet y1 = d 1 ,

yi+1 = d i+1

− β + y , i 1 i

i = 1,

· · · , n − 1.

Connaissant y, le vecteur x s’obtient facilement, a` nouveau grâce ace a` la structure particulière ere de U de U .. On obtient xn = yn /αn ,

xn−i = ( yn−i

− γ − x − + )/α − , n i n i 1

n i

i = 1,

· · · , n − 1.

Revenon Revenonss un instan instantt au syst` système e me (3.5 (3.5)) avec vec (3.6 (3.6)) du prob probl` lème e me de la barr barree métaletallique. On voit bien que cette matrice satisfait au théor` eorème, eme, donc est a` diago2 ))/κ avec nale dominante, si s( x) 0. Or, s( x) = (h r ( x))/ avec κ le le coefficient de diffusion thermique (positif par convention) et r ( x) la fonction de transfert de chaleur : on montre que l’équation equation différentielle erentielle (3.1) est est bien pos´ pos ee, e´ e, c’est-à-dire a-dire que l’´ l’equation e´ quation admet une solution unique, précis´ ecisément ement si r si r ( x) 0. Enfin, la décomposition A ecomposition A = LU permet permet aisément ement de calculer le déterminant eterminant det( A) de A de A.. En effet, le déterminant eterminant d’un produit de matrices est egal égal au produit des déterminants, eterminants, c.-à-d. a-d. det ( A) = det( LU eterminant LU ) = det( L) det(U ). Or, le déterminant d’une matrice triangulaire est le produit des él´ eléments ements sur la diagonale de la matrice ; par conséquent equent det ( L) = 1 et donc detA = det d et (U ) = α1 α2 αn .

≥

≥

· ···

U est non nul pour les matrices a` Ce produit des el´ e´ léments ements sur la diagonale de U est diagonale dominante d’après es le théor` eorème. eme. Nous avons donc montré au passage qu’une matrice n matrice n n à diagonale dominante est inversible. La décomposition LU ecomposition LU ayant ayant et´ e´ té illustrée ee pour une matrice tridiagonale, nous allons maintenant exposer la gén´ enéralisation eralisation a` des matrices n n sans structure particulière. ere.

×

×

54

D´ ecomposition LU de de matrices

3 .3 3.3.1 3.3.1

Decomposition e´ composition LU de de matrices Algor Algorit ithm hmee de Gauss Gauss

Soit A Soit A une une matrice n matrice n cherche a` résoudre esoudre

note a × n ; On note a

i j ses

el´ e´ léments ements (réels eels ou complexes) et on

Ax = b avec b avec b

n

n

eme s’écrivant ecrivant ∈ R (ou C ), ce système a x + a x + ··· ··· a x + a x + ··· ··· 11 1

12 2

21 1

22 2

.. . an1 x1 + an2 x2 +

+ a1n xn = b1 + a2n xn = b2

(3.13)

.. .. . . + ann xn = bn .

··· ···

Le but de l’algorithme dit LU dit LU est est décrire ecrire la matrice sous la forme A = LU avec U avec U matrice matrice triangulaire supérieure, erieure, obtenue par par un algorithme algorith me de Gauss, et L et L matrice triangulaire inférieure erieure telle que l que l ii = 1, i = 1, 2, , n (l (lii étant etant les el´ e´ léments ements L ). On suppose que a 11 = 0 et on retranche de la i ème sur la diagonale de L). eme ligne de la matrice A matrice A sa sa première ere ligne multipliée ee par

···



ai1 , pour a11

i = 2, 3,

· · · , n,

ce qui donne a11 x1 +

+

a12 x2

 

a22

an2

−

a21 a a11 12

−

an1 a11 a12

.. .

 

x2

··· ··· + ··· ···

x2 +

··· ···

+ + +

=

a1n xn

 

a2n

ann

−

a21 a a11 1n

−

an1 a11 a1n

.. .

 

xn = b2

xn = b2

b est devenu Suite a` ces opérations erations le système Ax eme Ax = b est A(1) x = b(1) avec les coefficients de la matrice A (1) tels que (1)

= a1 j ,

(1)

= ai j

a1 j ai j

−

j = 1, 2, , n ai1 a1 j , i = 2, 3 a11

···

55

b1

· · · , n, j = 1, 2, · · · , n.

− −

a21 b a11 1

.. .

an1 a11 b1

(3.14)


Le but de cette opération eration est bien entendu de faire apparaˆ apparaˆıtre ıtre des z eros e´ ros dans la première ere colonne, a` partir de la deuxième eme ligne, c.-à-d. a-d. (1)

ai1 = 0,

i = 2, 3,

· · · n.

De même eme le vecteur b vecteur b (1) est tel que (1)

(1)

b1 = b1 ,

− aa

i1

· · · , n. On introduit la notation suivante : soient deux vecteurs x vecteurs x , y ∈ C alors x alors x peut peut etre eˆ tre n lignes), tandis que y le considér´ eré comme une matrice n × 1 ( a` une colonne et n lignes), vecteur transposé s’apparente a` une matrice 1 × n. Le produit xy est ainsi bien défini efini et egal e´ gal a` une matrice n matrice n × n d’él´ eléments x ements x y (avec x (avec x et y et y les coefficients de x de x et y et y respectivement). respectivement). On note I note I la la matrice identit identité´ n × n et e et e le premier vecteur vecteur de bi

= bi

11

b1 , i = 2, 3,

n

T

T

i j

i

j

1

la base canonique, dont le premier coefficient est égal egal a` 1 et les autres sont nuls. Soit le vecteur (l’exposant T (l’exposant T désigne esigne le transposé, e, car par commodité d’écriture ecriture on ecrira e´ crira toujours les vecteurs vecteurs colonnes comme le transpos´ transpos e´ d’un vecteur ligne) v1 =



a21 a31 , , 0, a11 a11

···

an1 , a11



T

alors si l’on note L1 = I

T 1 1

−v e

alors (il convient de s’en convaincre en explicitant le calcul) A(1) = L1 A,

b(1) = L1 b.

Le nouveau syst syst eme A e` me A (1) x = b(1) s’écrit ecrit (1)

(1)

a11 x1 + a12 x2 +

··· ··· + ··· ···

(1)

a22 x2 .. . (1)

an2 x2 +

··· ···

(1)

(1)

+ a1n x n = b1 (1) (1) + a2n x n = b2 .. .

(1)

.. .

(3.15)

(1)

+ ann x n = bn

Mainte Maintenan nantt on consid` considère ere le sous-s sous-syst` ystème eme a` part partir ir de la deux deuxi` ième e me lign lignee et la deux deuxi` ième eme (1)

colonne et on suppose que a22 = 0. On retranche des lignes i, i = 3, 4, deuxième eme ligne ligne multip multipli´ liée e e par par

(1)

· · · , n la

(1)

ce qui fait fait appar apparaˆ aˆıtre ıtre des des zéros eros dans dans la deuxi` deuxième eme

ai2



a22

colonne a` partir de la 3ème eme ligne. Cette opération eration peut s’écrire ecrire A(2) = L2 A(1) , 56

b(2) = L2 b(1)

Algorithme de Gauss

avec L2 = I

T 2 2,

v2 =

−v e

(1)



(1)

(1)

a32 a42

0, 0, (1) , (1) , a22 a22

···,

an2

(1)

a22



T

e2 étan e tantt le deux deuxi` ième e me vect vecteu eurr de la base base cano canoni niqu que. e. La matr matric icee A(2) est est par par cons´ conséquent equent

A(2) =

   

(2)

a11

a12 (2)

0 . .. .. .

a22

0

0

··· ··· ( ) ··· a ( ) ··· a

(2)

2 23 2 33

0 .. .

(2)

a1n

(2)

a2n

(2)

a3n .. .. . .

.. .

(2)

···

an3

(2)

ann

   

.

(3.16)

Ensuite on applique l’algorithme a` A(2) en faisant fai sant apparaˆıtre ıtre des z eros e´ ros sur la 3 (2)

eme e` me colonne a` partir de la 4 eme e` me ligne en supposant que a que a 33 = 0. Pour Pour gén´ enérali eraliser ser ce proc´ procéd´ edé, e , on intr introd odui uitt le vect vecteu eurr vi , en supp suppos osan antt que que l’él´ elément ement (i−1) aii de la matrice A(i−1) (i défini efini comme

vi =



≥ 1, avec la convention que A( ) = A) A) est non nul,

  ··· 0,

0

(i 1) (i 1) ai+1,i ai+2,i , 0, (i 1) , (i 1) ,

aii

− −

aii

− −

(i 1)

− ani · · · , (i−1) a ii

T

 

(les i (les i premi` premières eres composantes du vecteur v i sont nulles). Soit alors Li = I

T i i ,

−v e

avec ei le ième eme vecteur de la base canonique (dont seul la ieme e` me composante est i) ( non nulle et egale e´ gale a` 1), alors les coefficients de la matrice A telle que A(i) = Li A(i−1) , i = 1, 2,

· · · , n − 1,

sont tels que (i)

ak , j = 0, 1 Ainsi, la matrice

≤ j ≤ i, j + 1 ≤ k ≤ n.

A(n−1) = Ln−1 A(n−2) =

· · · = L − L − · · · L L A n 1 n 2

2 1

(3.17)

est triangulaire supérieure, erieure, on la note U note U avec avec

U = = A

(n 1)

−

=

  

(n 1)

−

a11

0 .. . 0 0

(n 1)

− (n−1) a

···· · ··· · ·

··· ···

0

a12 22

..

.

..

.

0 57

··· ··· ..

.

(n 1)

−

(n 1)

− (n−1) a a1n 2n

.. .

(n 1)

− (n−1) a

an−1,n−1 an−1,n 0

nn

  

.

(3.18)


Le système eme a` résoudre esoudre est alors A(n−1) x = b(n−1) avec b(n−1) = Ln−1 Ln−2

· · · L L b. 2 1

Par l’équation equation (3.17) on peut ecrire e´ crire A = L1−1 L2−1

· · · L−− L−− A( − ) 1 1 n 2 n 1

n 1

avec L avec L i−1 l’inverse de la matrice L matrice L i . − 1 T T T T T Or, L Or, L i = I + vi ei . En effet, ( I + vi eT i ) Li = I vi ei vi ei = I vi (ei vi )ei = I car on peut se convaincre aisément ement que le scalaire eT eres eres i vi = 0 car les i premi` composantes de v de v i sont nulles. Donc,

−

−

n−1 − − 1 −1 1 −1 L1 L2 · · · Ln−2 Ln−1 = ∏ ( I + vi eT i ). i=1

On observe que T vi eT i v j e j = 0,

j

si

≥i

car alors e alors e T equent i v j = 0 et par conséquent n−1 1 −1 1 −1 − − L = L1 L2 · · · Ln−2 Ln−1 = I + ∑ vi eT i . i=1

Explicitons la matrice L : on peut voir qu’elle est triangulaire inférieure erieure et elle s’écrit ecrit 1 0 0 0 0

L =

     

(0)

a21

(0)

a11

(0)

1 a32

a11

a22

.. .

(0)

an−1,1 (0)

a11

(0)

(1)

.. .

1

0

..

..

.

(1)

an−1,2 (1)

a22

(1)

an1

an2

a11

a22

(0)

0

(1)

a31

(0)

0

(1)

.

··· ··· ···

0 0

.

.. .

1

0

..

(n 3)

··· ···

On peut donc enoncer e´ noncer le théor` eorème eme suivant. 58

− (n−3) an−2,n−2 (n−3) an,n−2 (n−3) an−2,n−2 an−1,n−2

(n 2)

− 1 (n−2) an−1,n−1 an,n−1

     

(3.19)

Algorithme de Gauss (i 1)

−

Th´ Theor` e´ oreme e` me 9 Si dans dans l’algo l’algorith rithme me de Gauss Gauss les ´ les ´ el´ element em pivots aii = 0, i = ´ entss dits pivots sup erieure erieur e dont les ´ les ´ el´ elements 1, n 1 , alors il existe une matrice U triangulaire sup´ ´ ´ sont donn´ donnes inf erieure donn´ donnee ´ par (3.18) et une matrice triangulaire inf ´ ´ ´ par (3.19) telles que A = LU .

··· −



b. Au lieu Revenons a` la résolution esolution du système Ax eme Ax = b. l ieu d’appliquer l’algorithme de Gauss directement directement a` ce système, eme, il est en gén´ enéral eral préf´ eférable erable de l’´ l’ecrire e´ crire de mani maniere e` re e´ quivalente sous la forme equivalente LUx = b. (3.20) En effet, dans de nombreux problèmes emes de discrétisation, etisation, l’op´ l’operateur A e´ rateur A est est donné une fois pour toutes et le second membre b est variable. On résout esout le système eme b (ce qui est aisé car L est (3.20) de d e la façon ¸o n suivante : d’abord on o n r´ resout e´ sout Ly = b (ce triangulaire) et une fois y obtenu on détermine etermine x en résolvant esolvant U x = y avec U triangulaire. Evalu Evaluan antt le nombr nombree d’op´ d’opérati erations ons néces e cessa sair ires es pour pour la décomposition LU ecomposition LU . L’op´ ’opération eration (i 1)

−

(par a ii ) et ( n i)2 multiplications et au Li A(i−1) correspond a` n i divisions (par a tant d’additions, pour i pour i = 1, 2, , n 1. Par conséquent equent il faut

−

−

n 1

2 ∑ j j=1

2

··· − (n − 1)n(2n − 1) 2n ∼ 2n = 2 6

−

3

(pour n pour n grand grand)

3

multiplications et additions ainsi que

−

n 1

∑ j =

n(n

− 1) ∼ n

2

j=1

2

2

divisions. Il est indispensable de faire cette décomposition ecomposition une fois pour toutes lorsqu’on doit résoudre esoudre plusieurs fois la solution d’un système, eme, pour n grand, dont la matrice est A. De résoudre esoudre LU x = b n ecessite e´ cessite en effet beaucoup moins d’opérations. erations. Considérons Ly erons Ly = b ce eme b ce qui donne lieu au système

l21 y 1 ln1 y 1 + ln2 y 2 +

···

y1 + y2

= b1 = b2

···

+ yn = bn .

Par conséquent y equent y 1 = b1 et

−

i 1

yi = bi

−∑l

i j y j ,

j=1

59

i = 2,

· · · , n.


On voit que pour calculer y i il faut ( i d’où au total

−

n 1

2 ∑ j = 2 j

− 1) additions et autant de multiplications,

(n

− 1)n = n(n − 1) 2

y donne lieu au système opérations. erations. La résolution esolution de U de U x = y donne eme (en ( en commen co mmencç ant par p ar la composante x composante xn , U etant e´ tant triangulaire supérieure) erieure)

un−1,n−1 x n−1 u11 x 1 +

···

+

u1,n−1 x n−1

unn x n = yn + un−1,n x n = yn−1

···

+

u1,n x n

=

, i = 1, n

− 1.

y1 .

Maintenant x Maintenant xn = yn /unn et

xn−i =

 −− yn

i

∑ jn=n i+1 un i, j x j

− − un−i,n−i



Le nombre total de multiplications et additions est a` nouveau ( n 1)n et il faut ajouter n divisions, ce qui donne n2 opérations. erations. Donc, pour n grand le nombre total des opérations erations en vue de la résolution esolution des deux systèmes emes triangulaires est 2 de l’ordre de 2n 2n . En conclusion, il y a un facteur n/3, pour n grand, entre le nombre d’opérations erations pour la décomposition LU ecomposition LU (ou (ou de manière ere equivalente e´ quivalente pour laméthode ethode d’élimi e limina natio tion n de Gaus Gauss) s) et le nomb nombre re d’op´ d’opératio erations ns perme permetta ttant nt de résoudre esoudre de L et et U LUx = b disposant b disposant de L U ..

−

3.3.2

Decomposition e´ composition LU avec avec permutations des lignes

Dans l’algorithme ci-dessus il a et´ e´ té supposé que les el´ e´ léments ements appelés pivots es pivots (i 1)

−

aii = 0, i = 1, 2, , n 1. Or, même eme pour une matrice inversible ces coefficients peuvent etre eˆ tre nuls et il convient convient d’y remédier, edier, en permutant pe rmutant des lignes dans l’algorithme de Gauss. D’une manière ere gén´ enérale, erale, et pour des raisons de stabilité numérique, erique, on cherchera toujours pour chaque pas de l’algorithme de Gauss de mettre mettre l’´ l’el´ e´lément ement le plus grand en module en position diagonale.



··· −

60

D´ ecomposition LU avec permutations permutations des lignes

Definition e´ finition 5 Soit une matrice Pi j de la forme

Pi j =

      

1 ..

.

(0) 1

··· ··· ···

0 .. . . . . .. . .. . 1

1

1 .. . .. . . . .. . . 0

··· ··· ···

1 ..

(0)

. 1

   ←    ← 

ligne i (3.21) ligne j

Alors Pi j A = Aˆ avec Aˆ matrice obtenue a` partir de A en permutant les lignes i et j. De meme, ˆ on montre que APi j = A˜ avec A˜ obtenue a` partir de A en permutant les colonnes i et j. On a les relations suivantes Pi2j = I , c.-à-d. a-d.

Pi−j 1 = Pi j ,

Pi j = PiT j ,

i et j deux fois laisse bien la matrice invariante. car en effet de permuter des lignes i et Nous allons modifier l’algorithme du chapitre 1.2.1 en y ajoutant la possibilit e´ de A et on cherche l’entier k permuter des lignes a` chaque etape. étape. On note A note A (0) = A et l’entier k 1 t.q. ( ) |a( ), | = max | a |, ≤≤ 0 k 1 1

1 k n

0 k 1

c.-à-d. a-d. on cherche l’él´ elément ement le plus grand en module dans la première ere colonne. On opère ere la permutation entre ligne 1 et la ligne k ligne k 1 et on obtient Aˆ (0) = Pk 1 ,1 A(0) et ensuite on applique l’algorithme de Gauss à Aˆ (0) pour trouver A(1) = L1 Pk 1 ,1 A(0) . On considère ere la deuxième eme colonne de A(1) et on veut mettre l’él´ elément ement le plus grand en position de pivot, c.-à-d. a-d. on cherche l’entier k l’entier k 2 tel que ( ) |a( ), | = max | a |. ≤≤ 1 k 2 2

2 k n

61

1 k 2


On opère ere la permutation entre ligne 2 et k 2 pour A pour A (1) d’où Aˆ (1) = Pk 2 ,2 A(1) = Pk 2 ,2 L1 Pk 1 ,1 A(0) . A l’étape i etape i on on obtient donc une matrice A(i) = Li Pi Li−1 Pi−1

· · · P L P A( ) 2 1 1

0

ou` l’on note P j = Pk j j , j avec k avec k j j t.q. ( − ) |a( ,− )| = max | |. a ≤≤ j 1 k j j j

j k n

j 1 kj

et a` nouveau on aboutit aboutit a` une matrice triangulaire supérieure erieure A(n−1) = Ln−1 Pn−1 Ln−2 Pn−1

· · · P L P A( ). 2 1 1

0

(3.22)

Afin Afin de comp compre rend ndre re comm commen entt on peut peut inte interp rpr´ réter ter le résultat esultat des produits produits successif successifss des matrices intervenant dans cette expression, on prend le cas particulier n = 4 et donc A(3) = L3 P3 L2 P2 L1 P1 A0) . Par P Par P j2 = I on on peut ecrire e´ crire A(3) = L3 P3 L2 P3 P3 P2 L1 P2 P3 P3 P2 P1 A(0)

(3.23)

I . Chaque matrice L car on n’intercale que des matrices identité I . matrice L j est de la forme T L j = I v j e j . Donc

−

Pk L j Pk = I

T k j j k

− P v e P = I − (P v )(P e ) k j

k j

T

= I (Pk v j )e jT

−

> j, car alors Pk e j = e j etant si k si k > e´ tant donné que seule la jème eme composante de e j est > j. Donc par exemple non nulle et P et Pk permute des el´ e´ léments ements d’indices > j si k si k > P3 P2 L1 P2 P3 = I Dans le cas gén´ enéral L eral L k = I

T 1

− (P P v )e 3 2 1

T k k

− v e et on note L˜ = I − (P − P − ....P + v )e k

L1 . = ˜

n 1 n 2

k 1 k

T k

= I v˜k eT k .

−

(3.24)

Lk est alors obtenue en opérant La matrice ˜ erant les permutations successives à partir de l’´ l’etape k e´ tape k + 1 au vecteur v vecteur vk qui définit L efinit Lk (voir chapitre 1.2.1). Avec Avec cette d´ definition e´ finition (3.23) devient (3.25) A(3) = ˜ L3 L˜ 2 L˜ 1 P3 P2 P1 A(0) 62


L3 = L3 . En extrapolant au cas gén´ avec la convention que ˜ enéral eral on trouve A(n−1) = ˜ Ln−1 L˜ n−2

· · · L˜ L˜ P − P − · · · P P A( ) 2 1 n 1 n 2

0

2 1

Les structures des matrices L˜ k sont analogues aux structures des matrices L matrices L k . On note P note P = Pn−1 Pn−2 P2 P1 et on trouve

···

˜ PA = LU

(3.26)

= A(n−1) matrice triangulaire supérieure avec U = erieure et ˜ −1 L˜ −1 L˜ == L 1 2

· · · L˜ −− L˜ −− . 1 1 n 2 n 1

L˜ est trianguComme pour les inverses des matrice L k on peut se convaincre que L est laire inférieure erieure de la forme

L˜ =

   

1 v˜12 v˜13 .. .

0 1 v˜22 .. .

0 0 1 .. .

v˜1,n−1 v˜2n−1 v˜1n v˜2n

··· ···

··· ··· ···

0 0 0 .. .

0 0 0 .. .

..

. 1

vñ−2,n−1 0 vñ−2,n vñ−1,n 1

   

(3.27)

avec v˜k ,i la i eme e` me composante du vecteur v˜ k introduit dans (3.24). Remarque : si e´ lément ement non nul dans les différentes erentes etapes e´ tapes de A est A est inversible, il y a toujours un el´ l’algorithme a` mettre en position de pivot. On peut donc enoncer e´ noncer :

Th´ Theor` e´ oreme e` me 10 Si la matrice A est inversible, alors il existe une matrice de permutation P, une matrice U triangulaire sup erieure et une matrice L˜ triangulaire ´ inf erieure avec 1 sur la diagonale, telles que ´ ´ PA = ˜ LU LU .

(3.28)

En pratique on ne garde pas les matrices P1 , P2 , , Pn−1 pour former P mais on garde garde l’effet l’effet des permutat permutations ions sur les indices. indices. Pour Pour ce faire, on introduit introduit un tableau tableau d’entiers ( p1 , p2 ,..., pn) t.q. au départ p epart p i = i. Opérant erant les permutations successives sur les el´ e´ léments ements de ce tableau, l’effet du produit des permutations P pour tout vecteur y vecteur y sera sera ( Py)i = y pi . Donc, si on veut résoudre esoudre

···

Ax = b ˜ Ly = ˜b on forme PAx forme PAx = LUx esout successivement ˜ = Pb = ˜b avec ˜bi = b pi et on résout et ensuite U x = y. y. 63


Exemple : il convient de traiter un petit exemple afin d’illustrer la mi se en œuvre de la méthode LU ethode LU . Prenons la matrice 3 3 A =

 

×

 

3 1 6 2 1 3 1 1 1

;

On cherche a` décomposer A ecomposer A en en un produit LU produit LU , en opérant erant eventuellement e´ ventuellement des permutations de lignes au cours cours des etapes e´ tapes successives de la mise sous forme triangu U obtenues laire. Les matrices L matrices L et et U obtenues lors de la décomposition ecomposition etant e´ tant respectivement triangulaire inférieure erieure et triangulaire supérieure, erieure, il s’avère ere pratique d’avoir une L et des el´ représentation esentation compacte compacte a` la fois des él´ eléments ements successifs de L et e´ léments ements de U . On rappelle rappel le que l’on l ’on connaˆıt ıt les el´ e´ léments ements sur la diagonale de L de L qui qui sont egaux e´ gaux a` 1. On note ( i, j ) la position a` l’intersection de la i la i ème eme ligne et de la j eme e` me co A : l’él´ lonne de A de A.. Prenons la première ere colonne de A : elément ement le plus grand en valeur absolue est en position (1, 1) et lors de la première ere etape e´ tape aucune permutation n’est nécessaire. ecessaire. On soustraira de la deuxi` deuxi eme e` me ligne la première ere multipliée ee par 2 /3 afin de faire apparaˆıtre ıtre un 0 en position posit ion (2, 1). Or, ce nombre 2 /3 est bien l’él´ elément ement en position (2, 1) de L de L.. Ensuite on soustraira de la troisième eme ligne la première ere multipliée ee par 1 /3, afin de faire apparaˆıtre ıtre un u n 0 en position positi on (3, 1). A nouveau le nombre 1/3 est l’él´ elément ement en position ( 3, 1) de la matrice L matrice L.. Il est alors commode de faire apparaˆ appar aˆıtre ıtr e les l es el´ e´ léments ements successifs de la matrice L et dans un même eme tableau de L et U U dans la fac fa ç on suivante sui vante (avec (ave c les el´ e´ léments ements correspondant correspondant a` L entre eses) L entre parenthèses)

 

3 1 (2/3) 1/3 (1/3) 2/3

6 1 1

 

− . − L’étape etape suivante porte sur la sous-matrice 2 × 2 a` partir de la deuxième eme ligne et

deuxième eme colonne du tableau ci-dessus : on voit qu’il convient de permuter les lignes 2 et 3 afin de mettre dans la deuxi eme e` me colonne l’él´ elément ement maximal en position (2, 2). D’où après es permutation on obtient

 

3 1 (1/3) 2/3 (2/3) 1/3

6 1 1

− −

 

.

Des el´ e´ léments ements de la troisième eme ligne, hormis la première ere colonne qui correspond aux el´ e´ léments ements de L de L,, on soustrait soustrait la deuxième eme ligne multipliée ee par 1 /2, pour faire apparaˆ appar aˆıtre ıtr e le 0 en posit p osition ion ( 3, 2). On met précis´ ecisément ement la valeur 1 /2, l’él´ elément ement de L, L, a` cette position. On récup` ecupère ere le tableau

 

3 1 (1/3) 2/3 (2/3) (1/2) 64

6 1 1/2

− −

 

,


U ,, a` savoir d’où on tire L tire L et et U L =

On obtient donc

 

 

1 0 0 1/3 1 0 2/3 1/2 1

,

U =

 

3 1 0 2/3 0 0

6 1 1/2

− −

 

.

PA = LU avec P avec P une une matrice de permutation des lignes 2 et 3

P =

 

1 0 0 0 0 1 0 1 0

 

.

En gén´ enéral eral les matrices de permutation ne sont pas explicitées, ees, il suffit de garder le résultat esultat des permutations dans un tableau d’entiers comme expliqué ci-dessus, qui devient ici ( 1, 3, 2). Admettons qu’on cherche cherche a` résoudre esoudre

A

        x1 x2 x3

2 7 4

=

.

Alors d’appliquer P d’appliquer P revient revient a` permuter les el´ e´ léments ements 2 et 3 du second membre et A = LU on par P par PA on obtient

   

LU LU

x1 x2 x3

L

y1 y2 y3

D’abord on résout esout

   

=

=

   

2 4 7 2 4 7

   

.

ce qui est aisé etant e´ tant donnée ee la structure triangulaire de L. On trouve y trouve y 1 = 2, y2 = 10/3, y3 = 4 et la solution x solution x est est obtenue en résolvant esolvant U x = y, y, ce qui est a` nouveau aisé etant e´ tant donné que U que U est est egalement e´ galement triangulaire. On trouve en définitive efinitive x1 = 19, x2 = 7, x3 = 8.

−

−

65


66

Chapitre 4 Normes de matrices, m´ methodes e´ thodes it´ iteratives e´ ratives de r´ resolution e´ solution de syst` systemes e` mes lin´ lineaires e´ aires 4.1 4.1


La résolution esolution directe d’un système eme linéaire, eaire, d ecrite e´ crite au chapitre 1, donne lieu a` un nombre fini d’opérations erations algébriques ebriques qui peuvent peuvent etre eˆ tre exécut´ ecutées ees sur un ordinateur. Ce nombre d’opérations erations d´ depend e´ pend de la taille de la matrice. Dans certains cas, notamment pour des très es grandes matrices, il peut etre eˆ tre judicieux de faire appel a` des méthodes ethodes itératives eratives : on cherche à créer eer des séquences equences de vecteurs qui convergent vers la solution exacte. Souvent des systèmes emes linéaires eaires de très es grandes dimensions proviennent proviennent des discr´ discr etisations e´ tisations d’équations equations différentielles erentielles ou d’équations equations aux dériv´ erivées ees partielles. Or, Or, une telle discr´ disc rétisation etisation introduit naturellement des erreurs : il est donc légitime egitime de chercher la solution dans un processus it´ iteratif e´ ratif que l’on arrête ete a` partir du moment que la limite, ici la solution du système, eme, est atteinte avec une erreur définie efinie a priori. Soit donc A donc A une une matrice n n et on cherche a` r esoudre e´ soudre Ax = b ; on suppose que l’on puisse ecrire e´ crire A = M N

×

−

avec M avec M une une matrice inv i nversible. ersible. On suppose que les syst` syst èmes emes avec M avec M comme comme matrice sont facilement inversibles et on définit efinit une suite de vecteurs x(k ) par les relations de récurrence ecurrence Mx (k +1) = Nx N x(k ) + b,

k = = 0, 1, 2,

···

a` partir d’un vecteur initial donné x(0) . On voudrait que la suite de vecteurs converge quand k quand k ∞, quel que soit x soit x(0) . Ainsi, si l’on note x note x la la limite, l imite, en supposant

→

67

Normes de matrices, m´ ethodes it´ eratives de r´ esolution de systèmes lin´ eaires

qu’elle existe, on aura alors à la limite pour k pour k Mx = Nx N x + b,

→∞

Ax = ( M N ) x = b.

−

ou encore

Donc, en cas de convergence la suite tend vers la solution x. x. Evidemment, le veck ) k ) ( ( x est x est solution du système teur r teur r = x eme

−

Mr (k +1) = Nr N r (k )

(4.1)

N x(k ) + b le syst` (pour (pour le const constate aterr, il suf suffit fit de soustr soustrair airee du syst` système Mx eme Mx (k +1) = Nx système eme k +1) 1 − ( matricee M ´ etant e´ tant supposée ee inve inversi rsible ble,, on aura aura r Mx = Nx N x + b). La matric = M N r (k ) 2

et egalement e´ galement r (k ) = M −1 Nr (k −1) et par conséquent equent r (k +1) = M −1 N r (k −1) . On peut donc déduire eduire en itérant erant que

 

k

r (k ) = M −1 N r (0)

 

k

avec r avec r (0) un vecteur initial a priori arbitraire et M −1 N le produit k produit k -` -eme e` me de la matrice C = ethode itérative erative converg convergee si r si r (k ) tend vers zéro, ero, ce qui est = M −1 N . La méthode 0) ( k assuré, e, quel que soit r , si C tend vers la matrice zéro ero (c’est-à-dire a-dire la matrice dont tous les coefficients sont nuls). Il faut donc disposer d’une “mesure” de l’ordre de grandeur d’une matrice, pour par exemple evaluer e´ valuer le comportement des produits successifs C k . Pour les vecteurs, vecteurs, une telle mesure est donn´ donn ee e´ e par la norme et au paragraphe suivant nous montrons comment comment on peut g´ g en´ e´ néraliser eraliser la notion de norme a` des matrices.

4.2 4.2

 

Norm Normes es de matr matric ices es

x on Soit l’espace Rn des vecteurs de n de n composantes composantes réelles eelles ; pour tout vecteur vecteur x on considère ere les 3 normes suivantes n



n

|| x x|| = ∑ | x x |, || x x|| = ∑ x , || x x||∞ = =max | x x |. ,···, = = 1

i

2

i 1

2 i

i 1

i 1

n

i

(4.2)

×

Soit A une matrice n n ; on peut définir efinir par rapport a` une des 3 normes A . On parle alors d’une norme subordonnée vectorielles ci-dessus une norme pour A. ee x une des normes ci-dessus a` la norme vectorielle. On note de manière ere gén´ enérale erale x A qui y est associée et on définit efinit la norme A ee de la manière ere suivante.

|| ||

|| ||

Definition e´ finition 6 Pour une des normes vectorielles x x de (4.2) la norme matricielle A A associ ee comme suit : ´ est d efinie ´ ´

|| ||

|| ||

|| A A|| =

max

x Rn , x=0

∈



68

|| Ax Ax|| || x x|| .

(4.3)

Normes de matrices

De mani ere ´ e` re ´ equivalente on peut ´ peut ´ ecrire

|| A A|| =

max

x Rn , x x =1

∈ || ||

|| Ax Ax||.

(4.4)

Tout d’abord il faut montrer que les deux définitions efinitions (4.3) et (4.4) sont bien equivalentes. e´ quivalentes. Partant Partant de (4.3), on peut ecrire e´ crire

 

Ax|| 1 1 || Ax = || Ax|| = || A A || x x|| || x x|| || x x|| x || Or, Or, le vecteur vecteur y est de norm normee 1, donc donc pour pour tout tout x avec y de  0, |||| |||| = || Ay y = || || x est x = Ay|| avec y x|| = 1 alors de manière Ax|| = norme 1. Réciproquement, eciproquement, soit x soit x t. t. q. || x ere evidente e´ vidente || Ax Ax Ax x x

1 x x

|| Ax Ax|| equent, les deux définitions efinitions ci-dessus sont sont bien equivalentes. e´ quivalentes. || x x|| . Par conséquent,

|| ||

A . Plus précis´ On vérifie erifie les propriét´ etés es habituelles de normes pour A ecisément, ement, A défini soit A efini comme ci-dessus alors :

|| ||

(4.5) || A A|| ≥ 0 et || A A|| = 0, si et seulement si A = 0; ||α A|| = |α||| A A||, pour tout α ∈ R; (4.6) || A A + B|| ≤ || A A|| + || B B||, quelles que soient A, B matrices n × n. (4.7) A|| = 0 implique d’après Montrons par exemple (4.5) : || A e s la définition efinition que Ax|| = 0 pour tout vecteur x et alors que Ax = 0 quel que soit x (d’après es les || Ax

propriét´ etés es de normes vectorielles). Mais on en déduit eduit que A que A est alors la matrice identiquement zéro. ero. La relation (4.6) est evidente e´ vidente d’apr d’après e` s la définition. efinition. L’inégalit´ egalité triangulaire (4.7) se démontre emontre précis´ ecisément ement a` l’aide de l’inégalit´ egalité triangulaire pour les normes vectorielles vectorielles ; en effet

|| A A + B||

=

max ( A + B) x

|| || ≤ ||max (|| Ax Ax|| + || Bx Bx||) ||= ≤ ||max || Ax Ax|| + max || Bx Bx|| = || A A|| + || B B|| ||= || ||= || x x||=1 x x

1

x x

x x

1

1

Il reste a` etablir e´ tablir deux propriét´ etés es fondamentales des normes de matrices, à savoir

|| Ax Ax|| ≤ || A A|||| x x|| || AB AB|| ≤ || A A|||| B B||

n

pour toutes matrices n matrices n

× n A, B.

quel que soit

x

∈R ,

(4.8)

(4.9)

L’inégalit´ egalité (4.8) découle ecoule directement de la définition efinition ; on utilise cette cette inégalit´ egalité pour montrer que que

||( AB) x|| = || A A( Bx)||≤|| A A|||| Bx Bx||≤|| A A|||| B B|||| x x||. 69


Par conséquent equent pour tout x tout x = 0,



||( AB) x|| ≤ || A A|||| B B|| || x x|| et on en déduit eduit (4.9), le “max” etant e´ tant le plus petit des majorants. En fait, il est possible de caractériser eriser pour les normes vectorielles usuelles A en termes des él´ (4.2) les normes matricielles de A eléments ements de A de A que l’on note ai j . On a notamment

|| ||

n

|a |, ··· ∑ i j

(4.10)

|| A A||∞ = =max ∑ |a |. ,···, =

(4.11)

|| A A|| = 1

max

j=1, ,n i 1 = n i 1

n

ij

j 1

|| ||

A 2, l’expression est Nous allons démontrer emontrer ces relations en T.D. Pour la norme A un peu plus complexe et fait intervenir la notion de rayon spectral d’une matrice. Le rayon spectral est lié a` la notion de valeurs propres de matrices (qui a ét´ eté rappelée ee au chapitre 2.1).

Definition e´ finition 7 Soit une matrice n n que l’on note B et on note λ1 , λ2 , ,n les valeurs propres de B (complexes en g en´ en Le rayon spectral de B, not ´ not e´ ρ ( B) , ´ eral). ´ est par d efinition la plus grande des valeurs propres en module, c’est- a-dire a` -dire ´ ´

×

···

ρ( B) = max λ j . j=1, ,n

···

| |

On peut montrer la propriét´ eté suivante : soient B , C deux deux matrices n matrices n

ρ( BC ) = ρ(CB).

× n ; alors (4.12)

Ce résultat esultat assez utile en pratique se démontre emontre de la façon ¸o n suivante. suivante. Soit λ la valeurs propres de BC de BC telle telle que ρ( BC ) = λ . Multiplia Multipliant nt l’égalit´ egalité BCx = λ x, x = 0, C on aura C B(Cx) = λCx. Cx. Donc, si Cx a` gauche par C on aura CB si Cx = 0, ce vecteur est vecteur B de valeur propre λ. Si C x). Dans propre de C de CB Si Cxx = 0 alors λ = 0 (car BCx = λ x). λ car ρ(CB) est la plus grande des valeurs propres en les deux cas ρ(CB) module de CB de CB.. D’où ρ (CB) ρ( BC ). Evidemment, l’inégalit´ egalité inverse se montre de la même eme faç on, commenc comm enç ant par ρ(CB). On en déduit eduit alors l’égalit´ egalité (4.12) ci-dessus. Considérons erons mainte maintena nant nt la cas cas parti particul culier ier d’une d’une matric matricee symétrique B etrique B.. On peut facilement montrer que ses valeurs propres sont réelles. eelles. En effet, soit λ valeur propre de B, a priori complexe. Alors Bx = λ x avec x avec x Cn , x = 0. On note x¯

||

≥| |





≥

∈

70



Normes de matrices

le vecteur dont les coefficients sont les conjugués es complexes de x et on forme l’expression x¯T Bx = λ x¯T x que l’on peut encore écrire ecrire sous la forme x¯T Bx = ( BT x¯)T x = ( B x¯)T x = ¯λ x¯T x. En effet, BT = B et e´ tant une matrice réelle, eelle, de Bx = λ x on eduit B x¯ = ¯λ x. B et B etant x on déduit x¯. T T T ¯ ¯ Par conséquent equent λ x¯ x = λ x¯ x : eduit que λ = λ (car x¯ x = 0) et λ est par x : on en déduit conséquent equent un nombre réel. eel. x, y d’un y d’unee matrice On montre egalement e´ galement que deux vecteurs propres x, matrice sym´ symétrique etrique B de B de valeurs propres λ et µ µ, sont orthogonaux, c’est et µ respectivement, respectivement, avec λ = µ, T T T x. Mais egalement y x. a-dire a` -dire y x = 0. En effet, y Bx = λ y x. e´ galement y T Bx = ( By)T x = µyT x. D’où (λ µ) yT x = 0 ; or, or, λ µ = 0 par hypothèse ese et donc y donc y T x = 0. Par ailleurs, on peut toujours choisir le vecteur propre x associé a` la valeur propre λ t.q. x t.q. x T x = 1. (En effet, un vecteur propre n’est défini efini qu’à une constante multiplicative près.) es.) On peut donc choisir des vecteurs propres normés es (de norme 1) et ils sont orthogonaux orthogonaux deux à deux pour des valeurs propres distinctes lorsque B est B est symétrique. etrique. On peut en déduire eduire (la démonstration emonstration complète ete ferait plutôt ot partie d’un cours d’algèbre ebre linéaire) eaire) qu’on peut construire une base une base orthonorm ee ´ B symétrique. de vecteurs propres d’une matrice B sym´ etrique. On peut enoncer e´ noncer ces résultats esultats par le théor` eorème eme suivant.





−

− 

Th´ Theor` e´ oreme e` me 11 Soit B une matrice n n et sym´ symetrique, on note λ i , i = 1, , n ses ´ valeurs propres qui sont r ´ r eelles. Alors on peut peut construire construire une base x (i) , i = 1, , n ´ orthonorm´ orthonormee formee a` -dire ´ form´ ´ de vecteurs propres de B, c’est- a-dire

×

···

x(i) T x( j) = δi j et et

···

Bx(i) = λi x(i)

δi j = 0 si avec δi j le symbole dit de Kronecker ( δi j = 1 si 1 si i = j et δ 0 si i = j). ρ( B) le rayon spectral de B. Alors pour tout vecteur x De plus, soit ρ



T

T

| x x Bx| ≤ ρ( B) x x.

(4.13)

Il reste a` d emontrer e´ montrer l’inégalit´ egalité (4.13). En effet, les vecteurs propres forment une base orthonorm orthonormee, e´ e, on peut peut ecrire e´ crire pour tout vecteur x vecteur x n

x =

∑ αi x(i)

i=1

avec α i des nombres réels. eels. Alors, on montre a` partir des relations Bx(i) = λi x(i) et x(i) T x( j) = δi j , que (exercice) T

n

x Bx =

∑ α2i λi.

i=1

71


Prenons la valeur absolue de cette expression, on aura n

T

n

α2i

2 i

= ρ( B) xT x

| x x Bx| ≤ ∑ |λ | ≤ ρ( B) ∑ α i

i=1

i=1

car ρ( B) = maxi=1,···,n λi et on montre que x que x T x = ∑ni=1 α2i .

| |

Après es ce petit détour, etour, revenons aux normes des matrices et en particulier la norme A A 2 avec A n n. On montre que avec A matrice matrice n

|| ||

× || A A|| = 2



ρ( AT A) =



ρ( AAT )

(4.14)

A. avec ρ( AT A) le rayon spectral de la matrice A T A. Nous allons démontrer emontrer ce résultat. esultat. Les propriét´ etés es des rayons spectraux rappelées ees ci-dessus permettent d’affirmer tout d’abord que ρ( AT A) = ρ ( AAT ). En A est symétrique, suite, la matrice A matrice A T A est etrique, donc, ses valeurs propres sont réelles eelles et en plus positives ou nulles. En effet, si x = 0 est vecteur propre de valeur propre λ, x implique que xT AT Ax = λ xT x : x : or, xT AT Ax = ( Ax)T Ax 0 et alors AT Ax = λ x implique bien sur x sur xT x > 0 et on en déduit eduit que λ 0. Rappelons la définition efinition de la norme

 ≥

|| A A|| = 2

≥

max

x Rn , x=0

∈



|| Ax Ax|| || x x||

2

2

Ax 22 = ( Ax)T ( Ax) = xT AT Ax ; prena B = AT A (en ometet bien ien sûr ur Ax prenant nt dans dans (4.13) (4.13) B tant la valeur absolue), absolue), on d´ deduit e´ duit la majoration

|| ||

2 2

T

|| Ax Ax|| ≤ ρ( A A) || x x|| et donc

2 2

 || || ≤ A A

ρ( AT A).

2

A telle que ρ( AT A) = λ i et x(i) = 0 Soit maintenant λi 0 valeur propre de AT A telle vecteur propre associé : alors

≥



Ax( ) || = ( Ax( ) ) Ax( ) = x( ) A Ax( ) = λ x ( ) x( ) = ρ( A A) || x x( ) || . || Ax i

2 2

i

T

i

i T T

i

i

i T

i

T

i

2 2

On en déduit, eduit, par la définition efinition de la norme, l’inégalit´ egalité inverse

 || || ≥ A A

2

ρ( AT A),

ce qui achève eve la démonstration emonstration de (4.14). Prenons le cas particulier où la matrice A elle-même eme est symétrique. etrique. Alors AT A = A2 et etant e´ tant donné que ρ ( A2 ) = ρ( A)2 on on aura le résultat esultat

|| A A|| = ρ( A), 2

pour toute matrice A matrice A sym´ symétrique etrique. 72

(4.15)

Application : conditionnement d’un syst` eme lin´ eaire

4.2.1 4.2.1

Appli Applicati cation on : conditi conditionn onnemen ementt d’un d’un syst système e` me lin´ lineaire e´ aire

La notion de norme de matrice permet notamment d’évaluer evaluer ce qu’on appelle le conditionnement d’un d’un syst` systeme e` me linéaire. eaire. Soit Soit a` résoudre esoudre Ax = b avec A matrice n n inversible et b Rn un vecteur donné. e. Lorsqu’on evalue e´ value les coefficients réels eels de b de b à l’aide d’un ordinateur, les valeurs ne sont en gén´ enéral eral pas reproduite exactement mais plutôt ot avec des erreurs d’arrondis. Donc, au lieu de résoudre esoudre le système eme ci-dessus on aura plutôt ot un second membre b + ∆b et on cherche a` evaluer e´ valuer l’effet de ∆b sur la solution exacte x exacte x,, donc on cherche à evaluer e´ valuer x telle que l’erreur ∆ x telle A( x + ∆ x) = A(b + ∆b).

×

∈

b on peut ecrire Par Ax Par Ax = b on e´ crire

∆ x = A−1 ∆b

et par la propriét´ eté (4.8) appliquée ee a` A −1 on aura

||∆ x||≤|| A A− ||||∆b||. 1

Ax on déduit De même, eme, de b de b = Ax on eduit

||b||≤|| A A|||| x x||. On obtient ainsi la relation entre l’erreur relative sur la solution en fonction de l’erreur relative du second membre

||∆ x|| ≤ || A A|||| A A− || ||∆b|| . (4.16) || x x|| ||b|| A|||| A A− ||, notée cond La quantité || A ee cond ( A), est appelée ee conditionnement de A : 1

1

plus le conditionnement cond ( A) est grand, plus l’erreur sur la solution risque d’être etre grande, par rapport rapport a` l’erreur ∆b due au second membre. Prenons la norme 2 2

T

T

λ ( A A) || A A|| = ρ( A A) = =max ,···, i 1

n

i

A. De même, avec λi ( AT A) les valeurs propres (positives) de A T A. eme,

|| A A− || = ρ 1 2 2



( A−1 )T A−1 = ρ A−1 ( A−1 )T = ρ ( AT A)−1 .

 

 



On observe que les valeurs propres de ( AT A)−1 s’écrivent ecrivent 1/λi ( AT A). On en déduit, eduit, le rayon spectral etant e´ tant la plus grande des valeurs propres de ( AT A)−1 , que 1 . A A−1 22 = mini=1,···n λi ( AT A)

|| ||

73


Par conséquent equent

|| A A|||| A A−1|| =



maxi=1,···,n λi ( AT A) mini=1,···,n λi ( AT A)



1/2

.

(4.17)

Pour une matrice symétrique etrique cette relation devient simplement = ,···, |λ ( A)| , || A A|||| A A− || = max min = ,···, |λ ( A)| 1

i 1

n

i

i 1

n

i

si A est symétrique etrique,

(4.18)

avec λi ( A) les valeurs propres (réelles) eelles) de A de A..

4.3

Conditio Conditions ns de conver convergenc gencee

Après es ce detour, indispensable, par des normes de matrices, revenons à la résolution esolution de Ax = b par une une m´ methode e´ thode itérative erative Mx (k +1) = Nx N x(k ) + b,

k = = 0, 1, 2,

···,

A = M N .

−

(4.19)

Nous avons vu au chapitre 2.1, que si la suite de vecteurs x(k ) converge, alors la b . Aussi, d’après limite est précis´ ecisément ement la solution solution x de Ax = b. es (4.1), le vecteur ( ( k ) k ) x est x est solution du système erreur r erreur r = x eme

−

r (k ) = C k r (0) ,

C = = M −1 N ,

(4.20)

et donc la méthode ethode converge quel que soit le vecteur initial x (0) , si et seulement si C k

→ 0, quand k → ∞.

(4.21)

Nous pouvons pouvons enon e´ nonce cerr le résulta esultatt suiva suivant nt conce concerna rnant nt la conve converg rgenc encee de la méthode ethode it´ iterative e´rative (4.19).

Th´ Theor` e´ oreme e` me 12 La m´ methode eth er (4.19) con conver verge quel quel que soit soit le vecte vecteur ur initial initial ´ ode it erativ ´ ´ ativee (4.19) ( 0) x , si et seulement si ρ(C ) < 1 avec 1 avec ρ(C ) le rayon spectral de la matrice C = 1 − M N. N . Pour démontrer emontrer ce ce th´ theor` e´ orème, eme, il suffit d’après es ce qui préc ec ede e` de de démontrer emontrer que C k

→ 0, k → ∞,

si et seulement si

ρ(C ) < 1.

Pour démontrer emontrer ce résultat, esultat, on montre d’abord que si ρ(C ) 1, alors C alors C k ne tend pas vers zéro ero pour k ∞. En effet, soit λ valeur propre telle que λ = ρ (C ) et

≥

→

74

| |

Conditions de convergence

x = 0 vecteur propre associé. x et en itérant C x. Or, si λ 1, e. Alors Cx Alors Cx = λ x et erant C k x = λk x. k k ∞ (et donc C ne tend pas vers alors λ , ne tend pas vers vers zéro ero pour k k zéro). ero). Il reste a` montrer que si ρ (C ) < 1, alors C alors C 0. Nous supposons que C est est diagonalisable : c’est le cas notamment quand les valeurs propres sont distinctes. Alors on peut montrer (cf. cours d’algèbre ebre linéaire) eaire) que les n vecteurs propres sont linéairement eairement indépendants. ependants. Si l’on forme la matrice P, appelée ee matrice de passage, telle que les vecteurs vecteurs colonnes de la matrice sont les vecteurs propres propres de C , alors on peut montrer que C = = PDP−1



→

| ≥

→

avec D avec D matrice matrice diagonale avec les valeurs propres de C sur sur la diagonale. Formant 2 1 1 2 −1 − − les produits C = PDP PDP = PD P etc., alors

C k = PDk P−1 ,

avec

Dk =

 

λ1k

(0) ..

.

λnk

(0)

 

avec λ i , i = 1, , n les valeurs propres de C . Or, si λi < 1 alors λ ik 0 quand k k ∞ et donc D , et ainsi C ainsi C , tendent vers zéro ero quand k 0. Bien sûr, ur, toute k matrice n’est pas diagonalisable, donc le cas gén´ enéral eral exige plus de notions sur les matrices (il (i l faut alors connaˆıtre ıtre ce qu’on appelle la réduction eduction sous forme de Jordan Jordan des des matrices). Mais bon nombre de matrices sont diagonalisables et notamment les matrices symétriques. etriques. La cond conditi ition on nécessair ecessairee et suffisa suffisante nte de conver convergenc gencee ρ(C ) < 1 nécessite ecessite de connaˆ connaˆıtre ıtre les valeurs propres de C . On peut etablir e´ tablir une condition suffisante de convergence C . a` partir de la norme de de C

···

→

| |

→

→

Proposition Proposition 4 Si pour une norme matricielle C < 1 , alors la methode it erative ´ ´ ´ ( 0) (4.19) converge quel que soit le vecteur initial x .

|| ||||

En effet, nous avons vu (cf. (4.20) que r (k ) = C k r (0) x ( x x ´ b). Par les propriét´ avec r avec r (k ) = x(k ) x ( étant etant solution de Ax de Ax = b). etés es des normes de matrices (4.8) et (4.9) on peut écrire ecrire

−

||r ( )||≤||C |||| ||r ( )||. < 1, alors ||C || Or, si ||C || 0, quand k quand k → ∞ et donc ||r ( ) || → 0 ce qui implique | | | | → que r que r ( ) = x( ) − x → 0, quand k quand k → ∞. k

k

k

k

0

k

k

75


4.4

Methode e´ thode de Jacobi et de Gauss-Seidel, m ethode e´ thode de relaxation

b, le système On reprend, pour la résolution esolution de Ax de Ax = b, eme itératif eratif Mx (k +1) = Nx N x(k ) + b,

k = = 0, 1, 2,

· · · , avec

A = M N .

−

(4.22)

Bien Bien sur, uˆ r, l’expression de M d d efinit e´ finit une méthode ethode it´ iterative e´ rative en particulier particulier ; pour la résolution esolution de la récurrence ecurrence les syst` syst emes e` mes avec M avec M doivent doivent pouvoir etre eˆ tre résolus esolus facilement et une décomposition ecomposition possible est de prendre pour M la la matrice diagonale formée ee par la diagonale de A. Plus précis´ ecisément, ement, si l’on note a i j , i = 1, , n, j = 1, , n, les coefficients de A de A,, on introduit 3 matrices D matrices D,, E , comme suit : E , F comme

···

···

D =

E = =

  − 

  

0 ..

.

(0)

.

..

a21 .. .

..

an1

···

et alors bien sûr ur

.

an,n−1 0

a11 ..

.

(0) ..

(0)

. ann

  

, F = =

  − 

  

,

0 a12 .. .

(0)

··· .. ..

a1n .. .

.

. an−1,n 0

− − F .

A = D E

(4.23)

  

. (4.24)

(4.25)

On supposera que les él´ eléments ements sur la diagonale de A sont non nuls, c’est-à-dire a-dire aii = 0, i = 1, , n.



···

Definition e´ finition 8 Dans la m´ methode de Jacobi, on choisit M = D et N = = E + + F dans ´ la m ethode it erative (4.22) avec D , E , F les matrices donn ees ´ ´ ´ ´ par (4.23) et (4.24). (k +1)

Les composantes x i de x (k +1) s’obtiennent alors en fonction des composantes de x (k ) par l’algorithme (k +1)

xi

=

1 aii



−

i 1

bi

n

− ∑ a x( ) − ∑ a x( ) = =+ j 1

k i j j

k i j j

j i 1



, i = 1,

· · · , n.

(4.26)

Dans ce cas, la m ethode converge quel que soit le vecteur initial x (0) , si et seule´ ment si C J k 0 si k ∞ avec

→

→

C J = D−1 ( E + + F ) = D−1 ( D A) = I D−1 A.

−

76

−

(4.27)

Methode ´ de Jacobi et de Gauss-Seidel, m´ ethode de relaxation relaxation

L’expression de la matrice C J résulte esulte simplement de la décomposition ecomposition A = D E F et et le résultat esultat de convergence a ét´ eté etabli e´ tabli au paragraphe paragraphe pr ec´ e´ cédent. edent. Dans l’algo l’algorit rithme hme (4.26) (4.26) la premi` première ere des sommes sommes dispara disparaˆˆıt ıt bien sûr u r si i = 1 et la la dern derni` ière ere 0 n n et on utilisera la convention si i si i = n et convention que ∑ j=1 0 et ∑ j=n+1 0. Il est facile de = D et D et N N = E + + F donnent voir que les récurrences ecurrences (4.22) avec M = donnent effectivement lieu a` (4.26) : en effet, d’inverser la matrice diagonale D revient simplement a` diviser composante par composante par a ii = 0. Dans la méthode ethode de Gauss-Seidel, au lieu de prendre pour M seulement seulement la A , on choisit M partie diagonale de A, choisit M = D E , c’est-à-dire a-dire on choisit tout le bloc triangulaire inférieur. erieur.

− −

≡

≡

 −

Definition e´ finition 9 Dans la m ethode de Gauss-Seidel, on choisit M = D E et N = = F ´ (k +1) (k +1) dans la m´ methode it erative (4.22). Les composantes x i de x s’obtiennent ´ ´ ´ alors par l’algorithme

−

(k +1)

xi

=

1 aii



−

i 1

bi

n

− ∑ a x( + ) − ∑ a x( ) = =+ k 1 i j j

j 1

k i j j

j i 1



, i = 1,

· · · , n.

(4.28)

Dans ce cas, la l a m ethode converge quel que soit le vecteur initial x (0) , si et seule´ k ment ment si si C G 0 si k ∞ avec

→

→

C G = ( D E )−1 F = = ( D E )−1 ( D E A) = I

−

−

− −

− ( D − E )− A. 1

(4.29)

L’expression de la matrice C G résulte esulte a` nouveau de la décomposition ecomposition A = D E F . Aussi, pour résoudre esoudre le système eme avec M = D E , afin de trouver k +1) k ) ( ( en fonction de x , on exploite la structure de D E qui qui est triangulaire x

− −

−

(k +1)

−

(k +1)

inférieure. erieure. Alors connaissant x connaissant x1k +1 , , xi−1 par (4.28), on peut déterminer x eterminer xi car la première ere somme du membre de droite de (4.28) s’arrˆ s’arr ete eˆ te a` i 1. Bien Bien sˆ sur, uˆ r, il 0 n faut a` nouveau utiliser la convention ∑ j=1 0 et ∑ j=n+1 0. Une méthode ethode it´ iterative e´ rative convergera d’autant plus vite que le rayon spectral de 1 − C = sera petit. Dans le but d’accél´ elérer erer dans certains cas la convergence = M N sera de la méthode ethode itérative, erative, on peut introduire un paramètre etre ω dans la méthode ethode de Gauss-Seidel Gauss-Seidel (on parle alors d’une m ethode e´ thode de relaxation).

···

≡

≡

−

Definition e´ finition 10 Dans la m´ methode de relaxation on choisit ´ M = =

1

−

D E

ω

N = =

et

1

D D + F ,

ω

−

(k +1)

pour un param` parametre e` tre r eel ´ ´ ω. Les composantes x i par l’algorithme (k +1)

xi

=

ω aii



−

i 1

bi

n

− ∑ a x( + ) − ∑ a x( ) = =+ j 1

k 1 i j j

k i j j

j i 1

77



de x(k +1) s’obtiennent alors

(k )

+ (1 ω) xi , i = 1,

−

· · · , n. (4.30)


ω = 1 on En particulier, particulier, pour ω 1 on retrouve l’algorithme de Gauss-Seidel. On note C G (ω) =



1

 −  − 1

D E

ω

1

 −

D D + F = I

−

ω

1

 − −

D E

ω

1

A

(4.31)

et la m´ methode de relaxation converge, quel que soit le vecteur initial x (0) , si et ´ seulement seulement si C G (ω)k 0 si k ∞

→

→

La dernière ere expression dans (4.31) s’obtient en écrivant ecrivant a` nouveau F nouveau F = D E A. A.

−

4.4. 4.4.1 1

−

Quel Quelqu ques es resultats e´ sultats de convergence de m´ methodes e´ thodes it´ iteratives e´ ratives

Proposition Proposition 5 Soi Soitt A une une matr matric icee n n, de coef coeffic ficie ient ntss ai j , i = 1, , n, j = 1, avec aii = 0, i = 1, , n et a` diagonale dominante stricte, c.- a-d. a` -d.



×

···

···

· · · , n,

n

∑ |ai j | < |aii|,

i = 1,

j =1

· · · , n.

j=i



Alors la methode de Jacobi converge. ´ En effet, effet, il suffit que pour pour une norme norme matricielle C J < 1 avec C J = I D−1 A d’après es la définition efinition de la méthode ethode de Jacobi (cf. Définition efinition 8). Or, il est facile de voir que les el´ e´ léments ements de la matrice C matrice C J sont

|| ||

ci j =

−

− aa , i = j, c = 0. ij

ii

ii

Or, Or, d’apr` d’après es la condi condition tion de domina dominance nce diagon diagonale ale strict strictee et en appliq appliquan uantt la formul formulee (4.11) on obtient n ai j < 1 C J ∞ = max ∑ i=1,···,n j=1 aii

 

|| ||

 

j=i

ce qui achève eve la démonstration. emonstration.

Considéron e ronss main mainte tena nant nt la méthode ethode de relax relaxati ation on (cf. (cf. Défini e finitio tion n 6) avec vec la mamatrice C G (ω) donnée ee par (4.31), pour la décomposition ecomposition (4.25) de A de A.. Le déterminant eterminant d’un produit de matrices est est egal e´ gal au produit des déterminants eterminants de chacune des matrices qui forment le produit. Par conséquent, equent, par (4.31), det (C G (ω)) = det



1

 −  − 1

D E

ω

det

78

det ω1 D D + F . D D + F = ω det ω1 D E 1

−

  −  −

Quelques r´ esultats de convergence de m´ ethodes it´ eratives

Or, les matrices dont on prend les déterminants eterminants sont des matrices triangulaires d’après es (4.25) et les déterminants eterminants sont donc egaux e´ gaux au produit des él´ eléments ements sur la diagonale. On en déduit eduit que 1

det (C G (ω)) =

n

 − 1

ω

∏ni=1 aii

∏ni=1 aii ωn

= (1

n

− ω) .

(4.32)

Cette expression expression permet de d´ d emontrer e´ montrer le résultat esultat suivant.

Th´ Theor` e´ oreme e` me 13 Pour toute matrice A, le rayon spectral de la matrice C G (ω) associ´ sociee methode de relaxation v erifie ´ a` la m´ ´ ´

ρ (C G (ω))

≥ |ω − 1|.

(4.33)

En particulier, particulier, si la m ethode de relaxation relaxation converg converge, e, alors 0 alors 0 < ω < 2. 2. ´ La démonstration emonstration utilise le fait que le déterminant eterminant d’une matrice est egal e´ gal au produit de ses valeurs propres. D’après es (4.32), et notant λ i les valeurs propres propres de C G (ω), n

∏ λi = ( 1 − ω)n i=1

et par conséquent equent

ρ (C G (ω)) = max λi

| | ≥ |1 − ω|. Il facile de se convaincre, qui si ω ≤ 0 ou ω ≥ 2, alors |1 − ω| ≥ 1. Or, la méthode ethode i=1, ,n

···

)) < 1 et le théor` converge, si et seulement si ρ (C G (ω)) < eorème eme est ainsi démontr´ emontré.

Nous allons maintenant introduire une famille de matrices particulières eres pour laquelle on peut enoncer e´ noncer un résultat esultat gén´ enéral eral quant a` la convergence des méthodes ethodes it´ iteratives e´ ratives de relaxation.

×

Definition e´ finition 11 Soit A une matrice n n ; on dira dira que la matrice matrice poss` possede e` de la pro pri´ priet ´ e´ ´ (P), s’il existe une matrice de permutation P telle que PAPT =



D1 A12 A21 D2



,

D1 , D2

matrices diagonales.

(4.34)

De former PAP former PAPT a` partir de A revient a` faire les mêmes emes permutations sur les lignes et les colonnes telles que ces permutations donnent lieu à la structure en bloc (4.34). On peut montrer le résultat esultat suivant.

Proposition Proposition 6 Toute matrice tridiagonale poss` poss ede e` de la propri et ´ e´ ´ (P). 79


Nous admettons ce résultat esultat dans le cas gén´ enéral eral mais nous en donnons une illustration pour le cas n cas n = 3. Prenons

 

A =

a b 0 c a b 0 c a

 

Si l’on opère ere une permutation des lignes 1 et 2 suivie de la permutation des colonnes 1 et 2 on obtient a c b b a 0 c 0 a

 

 

qui est bien de la forme (4.34) avec D1 =

    a

a 0 0 a

, D2 =

etc.

Pour les matrices ayan ayantt la propri´ propriet´ e´ té (P) on peut enoncer énoncer le théor` eorème eme suivant.

Th´ Theor` e´ oreme e` me 14 Soit une matric matricee r eelle n ´ ´ aii = 0 , ayant les propri propri et ´ es ´ ´ suivantes :



telle que les ´ les ´ el´ elements sur la diagonale × n telle ´

1. A pos posssede e` de la propri et ´ e´ ´ (P). 2. La matr matric icee C J = I D−1 A a les valeurs propres propres r ´ r eelles. ´

−

3. Le rayon rayon spectr spectral al ρ(C J ) < 1 (c.-` (c.-a` -d. -d. la methode eth era Jacobi conver converge) ge).. ´ ode it erative ´ ´ tive de Jacobi Alors pour tout 0 tout 0 < ω < 2 le spectrall de la matrice matrice C G (ω) s’´ s’ecrit 2 le rayon spectra ´

ρ (C G (ω)) =

   1 4

ω

ω µ +

4(1

− 1,

−

ω) + ω2 µ2

avec µ = ρ(C J ) et

ω0 =



2

,

0 < ω < ω0

ω0 < ω < 2,

(4.35)

2

 −  −  −   −  −

1+

1 µ2

est le param` parametre e` tre de relaxation optimal. Le rayon spectral associ e´ est

ρ (C G (ω0)) = ω0

− 1 =

1

1 µ2

1+

1 µ2

µ2

=

1+

1 µ2

2

On observe que ω = 1 correspond 1 correspond en fait a` Gauss-Seidel. Alors ρ(C G (1))) = 2 2 ρ(C G) = µ = (ρ(C J )) . 80


1

ρ (CJ )

0.9

ρ (CG )

0.8

0.7

0.6

0.5

ω0

−1

0.4

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

ω0

F IG . 4.1 – Rayon spectral de la matrice ρ(C G (ω)).

µ = 0.9. Le graphe graphe du rayo rayon n spectr spectral al est est donn´ donné sur sur la figur figuree 4.1, 4.1, pour pour l’exe l’exemp mple le µ On voit ici que pour le paramètre etre de relaxation optimal ω0 le rayon spectral de C G (ω) est le plus petit : c’est alors que la convergenc convergencee est optimale.

Id´ Idee e´ e de la d emonstration e´ monstration : La démonstration emonstration de ce théor` eorème eme est assez complexe et nous allons seulement indiquer dans quelle mesure la propriét´ eté ( P) intervient. Tout d’abord, on peut se convaincre que de permuter les lignes et colonnes ne fait b en un système que transformer le système eme a` résoudre esoudre Ax = b en eme equivalent e´ quivalent : donc, bien qu’on n’opère ere pas explicitement les permutations, on peut supposer que A possède ede la structure par blocs (4.34). Evidemment, la démonstration emonstration cherche a` caractériser eriser le rayon spectral de C G (ω) et donc les valeurs propres de cette matrice. Soit donc λ valeur propre de C G (ω) et x = 0 vecteur propre associé, e, c.-àa x. D’après d. C d. C G (ω) x = λ x. es l’expression (4.31), on peut ecrire e´ crire (on multipliant par 1 ω D E )



−



1

D D + F

ω

−

 − λ

1

 −

D E

ω

x = 0

On peut multiplier cette expression expression par D par D −1 pour obtenir

− D

1

F + + λ D−1 E x =



81

λ+ω ω

− 1x


On suppose que λ est une valeur non nulle et on divise l’´ l’ egalit´ e´ galité ci-dessus par

 √ 1

√ −1 1 − D F + + λ D E



λ

x =

√ λ

:

λ+ω

√ − 1 x λω

et par conséquent equent

λ+ω

√ λω− 1

est valeur propre de

√ −1 1 − √ D F + + λ D E . 1

λ

Or, si A si A est est de la forme (4.34), de par la d´ d ecomposition A e´ composition A = D E D−1 F = =

 −

0 D1−1 A12 0 0



D−1 E = =

,

 −

− − F

0 0 1 − D2 A21 0



.

On note que d’après e s la définition efinition de la méthode ethode de Jacobi et par (4.27) C J = D−1 F + + D−1 E . Il est facile de voir qu’alors

√ −1 1 − √ D F + + λ D E = = 1

λ



I 1 0

0 λ I 2

√

  C J

I 1 0

0 λ I 2

√

−

1

avec I 1 , I 2 les matrices identités es relatives a` la décomposition ecomposition en blocs. On peut 1 1 − alors conclure que les valeurs valeurs propres de √ D F + sont celles de C de C J . + λ D−1 E sont

√

λ

Donc, les valeurs propres λ de C G (ω) sont telles que

λ+ω

√ − 1 = τ λω

avec τ valeur propre de la matrice C J associée ee a` la méthode ethode de Jacobi. On en déduit eduit que λ2 + λ(2ω 2 ω2 τ2 ) + (ω 1)2 = 0.

− −

−

On voit donc de quelle faç on les valeurs propres de C de C G (ω) dépendent ependent des valeurs propres de C J . La discussion des racines de cette équation equation (que l’on ne fait pas ici) donne lieu au résultat esultat du théor` eorème. eme.

Application a` des matrices tridiagonales Soit la matrice n matrice n

× n tridiagonale A =

  

a c

b a .. .

b .. .

(0)

c

82

(0) ..

. a c

b a

  

(4.36)


avec a avec a , b, c des nombres réels eels non nuls. On suppose également egalement que c/b > 0. Les matrices tridiagonales possèdent edent la propriét´ eté (P) du théor` eorème. eme. Pour appliquer le théor` eorème, eme, il convien convientt d’abord d’abord de calculer calculer les valeurs valeurs propres propres de A de A afin afin de connaˆıtre ıtre 1 − A. celles de C J = I D A. A , alors Soit x Soit x = 0 vecteur propre de A,

−



Ax = λ x.

(4.37)

Le vecteur a comme composantes x = ( x1 , x2 , , xn )T et si l’on rajoute formellement x0 = 0, xn+1 = 0, alors en ecrivant e´ crivant (4.37) composante par composante, on obtient les equations e´ quations

···

− λ) x + bx + = 0, j = 1, 2, · · · , n,

cx j−1 + (a

j

x0 = 0, xn+1 = 0. (4.38)

j 1

Il y a une méthode ethode gén´ enérale erale de résolution esolution d’´ d’equations e´ quations aux différences erences du type j (4.38) qui consiste a` chercher la solution sous la forme x j = αr . Injectant cette expression expression dans (4.38) on voit ais´ ais ement e´ ment que r que r est est solution de l’équation equation p(r ) = br 2 + (a

− λ)r + c = 0.

(4.39)

On admet que (4.39) a deux racines distinctes r 1 , r 2 et par linéarit´ earité de (4.38) on peut ecrire e´ crire sa solution gén´ enérale erale sous la forme j

j

x j = αr 1 + βr 2 . Par x Par x 0 = 0 on a

j

x j = α(r 1 et x et xn+1 = 0 implique α(r 1n+1 non nul) r 1 r 2

(4.40)

− r + ) = 0 et donc (α = 0 car on cherche un vecteur n 1 2

n+1



j 2

− r )

= 1, d’où

 

r 1 2ik π = exp r 2 n+1

= 1, 2, , k =

···n

(4.41)

par l’expression l’expression g´ gen´ e´ nérale erale des racines n + 1 eme ème de l’unité, e, avec ici k ici k = 0 car nous avons avons suppos´ suppos e´ que les racines sont distinctes. On peut ecrire e´ crire le polynôme ome de (4.39) sous la forme r = = b(r r 1 )(r r 2 ) d’où



−

r 1 r 2 = c/b,

r 1 + r 2 = ( λ

− a)/b

et par (4.41) c exp r 12 = exp b

 

83

2ik π n+1

.

−

(4.42)


On en déduit eduit que r 1 =

   ik π n+1

c exp exp b

r 2 =

,

  −  ik π n+1

c exp exp b

(4.43)

et par la deuxième eme relation de (4.42) les l es valeurs propres propres sont so nt

     − 

xk j = α

ik π exp n+1

ik π n+1

  

c k π . cos b n+1 (4.44) Par (4.40) la j ème eme composante x composante xk j du vecteur propre associ associ e´ x k s’´ s’ecrit e´ crit

λk = a + b

c b

+ exp

= a + 2b

     −  −      c b

j/2

exp

i jk π n+1

exp

i jk π

n+1

c b

=

j/2

sin

jk π n+1

.

(4.45) En effet, le vecteur propre etant e´ tant d efini e´ fini a` une constante constante multiplicativ multiplicativee (a ( a priori complexe) près, es, on peut choisir α = 1/(2i). Ici, les valeurs propres sont réelles, eelles, sous l’hypothèse ese que c que c/b > 0.

Donc, connaissant connaissant les valeurs propres propres de la matrice (4.36), il i l est ais e´ de calculer les valeu valeurs rs propr propres es de la matric matricee associ´ associée ee a` l a m etho e´ thode de de Jaco Jacobi bi C J = I D−1 A qui sont k π 2b c µk = k = = 1, 2, , n, cos a b n+1

−

−

  

···

et le rayon spectral est

         ≤ 

µ = ρ(C J ) = 2

b a

c cos cos b

π n+1

Pour les matrices du type (4.36) la méthode ethode de Jacobi Jacobi conv converge, erge, pour pour toute dimension n sion n de de la matrice, si b c 2 1. a b Et dans ce cas on peut utiliser la formule du théor` eorème eme pour trouver le paramètre etre de relaxation optimal, ce qui fait l’objet d’un exercice de TD.

84

Cours Calcul Scientifique 2013

Recommend Documents