Universit´e de Toulon SeaTech - 1`ere ann´ee Analyse Num´erique Travaux Pratiques V R´esolution de syst`emes tri-diagonaux et applications Le but de ce TP est de mettre en oeuvre quelques m´ ethodes pour la r´ esolution de syst` emes lin´ eaires tridiagonales issus de probl` emes physiques. Exercice 1 (M´ethode directe : algorithme de Thomas). Soient A ∈ Mn (R) une matrice tridiagonale telle que  i = j + 1 et 1 6 j 6 n − 1 ,   αj si  βj si i = j et 1 6 j 6 n , A = (ai,j )16i,j6n = γ si i = j − 1 et 2 6 j 6 n ,  j   0 sinon et b = (bi )16i6n ∈ Rn . Dans la suite on note α = (αj )16j6n−1 , β = (βj )1 6 j 6 n et γ = (γj )26j6n . ´ 1. Ecrire une routine thomas(n,α,β,γ,b,x) permettant de calculer la solution x du syst`eme Ax = b en proc´edant `a (a) la factorisation de Gauss sous forme LU o` u L est une matrice triangulaire (inf´erieure avec Li,i = 1, 1 6 i 6 n) et U une matrice triangulaire sup´erieure. On pose ici l = (Lj+1,j )16j6n−1 , d = (Uj,j )16j6n−1 et u = (Uj−1,j )26j6n , (b) la r´esolution du syst`eme triangulaire inf´erieure puis sup´erieure. 2. Avec une matrice 3 × 3 tridiagonale de votre choix, valider votre programme. Exercice 2 (M´ethode it´erative de relaxation par points). Avec les notations de l’exercice pr´ec´edent, on d´ecompose la matrice A sous la forme A = D − E − F o` u βj si i = j et 1 6 j 6 n , D= , 0 sinon −αj si i = j + 1 et 1 6 j 6 n − 1 , E= 0 sinon et F = −γj 0 si i = j − 1 et 2 6 j 6 n , sinon Nous supposerons de plus que la matrice A est diagonale strictement dominante, i.e. |β1 | > |γ2 |, et |βj − 1| > |αj | + |γj + 1|, 2 6 j 6 n − 1 et |βn−1 | > |γn |. Pour r´esoudre le syst`eme Ax = b, on propose la m´ethode de relaxation par points D 1−ω (k+1) −E x = D + F x(k) + b, ω ∈ R∗ . ω ω ´ 1. Ecrire une routine relaxation(n,ω,x(0) ,α,β,γ,b,x) permettant de calculer la solution x du syst`eme Ax = b `a ε > 0 pr`es. 2. Soient a ∈ R et la matrice  1 A= a 0 a 1 a    0 1 a  et b =  2  . 1 3 (a) Pour quelles valeurs de a la matrice est-elle `a diagonale strictement dominante ? On note I cette intervalle. (b) Soit a ∈ I. Pour quelles valeurs de ω, la m´ethode converge t-elle ? Lorsqu’elle converge, identifier sa limite ? ´ (c) V´erifier num´eriquement vos r´esultats. Etudier l’influence du param`etre ω en fonction du nombre d’it´erations k n´ecessaires pour la convergence `a ε = 10−12 fix´e. Tracer cette fonction. Conclusions ? Indications : Consid´erer le nombre d’it´erations kmax = 106 . Exercice 3 (Applications : fl´echissement d’une poutre). On consid`ere une poutre de longueur L, ´etir´ee selon son axe, soumise `a une force transversale f (x) dx par unit´e de longueur dx et fix´ee `a ses extr´emit´es. Le d´eplacement u(x) du point d’abscisse x est solution du probl`eme aux limites :  d2 u(x)   − + c(x)u(x) = f (x), x ∈ (0, L), dx2 u(0) = 0,   u(L) = 0 o` u c(x) est une fonction positive d´ecrivant les caract´eristiques de la poutre. L . n+1 (a) On pose u0 = u(0) et un+1 = u(xn+1 ). En notant ui une approximation u(xi ) au sens des diff´erences finies, montrer que le probl`eme aux limites s’´ecrit sous la forme An Un = Bn o` u A ∈ Mn (R) est une matrice tridiagonale Un = (ui )16i6n , Bn = (f (xi ))16i6n . 1. On consid`ere la subdivision r´eguli`ere xi = iδx, 0 6 i 6 n + 1 de pas δx = (b) La matrice A est-elle inversible ? (c) On fixe L = 1 dans la suite. Calculer la solution analytique lorsque c(x) ≡ 0 et f ≡ 1. (d) R´esoudre le syst`eme lin´eaire avec la routine thomas pour diff´erentes de valeurs de n. Conclusions ? (e) R´esoudre le syst`eme lin´eaire avec la routine relaxation pour diff´erentes de valeurs ´ de n. Etudier l’influence des param`etres ε et ω sur la solution. 2 (f) Comparer les deux m´ethodes. 2. Applications (normalement, `a ce stade, vous avez choisi votre m´ethode “favorite” . . . ). On suppose que la poutre ` a un comportement non lin´eaire de la forme c(x) = x(L − x) et est soumise ` a une force f (x) = x(x − 1). Calculer la solution num´erique. Exercice 4 (Applications : diffusion de la chaleur dans une tige). On consid`ere une tige m´etallique de longueur L avec un coefficient de dissipation de chaleur ν (non n´ecessairement homog`ene en espace x et en temps t). On soumet cette tige `a une source de chaleur f (t, x), x ∈ (0, L) et on maintient la temp´erature `a ses extr´emit´es `a 0 degr´e. L’´evolution de la chaleur dans cette tige est solution du IBVP (Initial Boundary Value Problems) :  ∂   u(t, x) = F (t, x, u(t, x)), t ∈ (0, T ], x ∈ (0, L),   ∂t u(t, 0) = 0,   u(t, L) = 0,   u(0, x) = 0. o` u F (t, x, u(t, x)) = ν(t, x) ∂2 u(t, x) + f (t, x) ∂x2 1 si t < T /2 , 0 sinon Pour simplifier le probl`eme nous supposerons ν constant. On consid`ere la subdivision espace-temps r´eguli`ere suivante xi = iδx, 0 6 i 6 N + 1 de L δx2 1 pas δx = et tn+1 = tn + δt avec δt = α , α ∈ 0, . N +1 ν 2 avec f (t, x) = ´ 1. Dans un premier temps, on fixe la variable x = xi . Ecrire le sch´ema num´erique d’Euler explicite. 2. Dans un second temps, on fixe la variable t = tn au niveau du second membre F . ´ Ecrire un sch´ema au sens des diff´erences finies. 3. Montrer que les deux op´erations pr´ec´edentes conduisent au sch´ema num´erique suivant : n ui−1 − 2uni + uni+1 n un+1 = u + νδt + δtf (tn , xi ) . i i δx2 ´ 4. On fixe T = 10, L = 1. Etudier l’influence des param`etres α, n, ν . Que se passe t-il 1 si α > ? Faire le lien avec l’exercice 1 du TP pr´ec´edent. 2 5. Reprendre le probl`eme pr´ec´edent avec la m´ethode d’Euler implicite. Notamment montrer que la version implicite du sch´ema pr´ec´edent (3.) n´ecessite la r´esolution d’un n syst`eme lin´eaire ` a chaque temps tn de la forme (I + δtAN )U n+1 = U n + BN o` u AN n et BN sont similaires ` a la matrice An et le vecteur Bn de l’exercice pr´ec´edent. En utilisant votre m´ethode favorite de r´esolution, r´epondre `a la question 4. 6. Comparer les deux approches graphiquement ? Quelles sont les avantages et inconv´enients ? 3