Interessant [PDF]

Zitiervorschau

ENIM-Monastir 2016-2017

RECUEIL DES SUJETS DE TP D’ANALYSE NUMÉRIQUE TP Analyse numérique Génie énergétique (ENR-1)

QU TI MA É TH MA

DE ES

L’E

M NI

I DH U D FO PE UI MA AH A ÉQ M SS ed DE Im B : rA ant aie h ign e u s Zo I En t: UN n — a O I n AH ouf seig a L n R U E : MO — I ant k n fi g AN u H o sei I a D n T E AR I GU t: R H n — a K n ni T ah seig Z ma Sal En A : RE : t H ) n — e a E n ir ata af M seig c s a n En I : t(V — re) nan i g a i t se ca En (Va t n — na seig n E — DE

.

1

2

TP no 1 Initiation à Matlab

1

Initiation à Matlab

Ce premier TP a pour objectif de vous rappeler les syntaxes et structures Matlab en vue d’une utilisation intensive pour les prochaines séances TP dans le cadre du module d’Analyse Numérique. Les aspects graphiques seront également abordés lors de ce TP. Remarque : Dans le présent TP, n’hésitez pas de demander de l’aide à Matlab sur les fonctions :

faire help + nom de function puis valider. Avant de commencer la programmation, utiliser les commandes clear all, close all pour initialiser l’environnement Matlab.

Exercice 1 On note

1. 2. 3. 4.



     2 2 −4 u1 =  3  , u2 =  −1  , u3 =  −3  −1 5 6

Définir ces vecteurs sous Matlab de plusieurs manières. Calculer u1 + u2 , u1 − 3u2 , u3 /6. Calculer k u1 k2 , k u2 k1 et k u3 k∞ . Calculer le cosinus de l’angle formé par les vecteurs u1 et u2 .

Exercice 2 Soient A1 , A2 et A3 les trois matrices    4 0 1 3 , A2 =  0 2 A1 = 5 2 0 0

suivantes    2 3 0 −1 4  , A3 =  2 4 2  . 6 −1 0 3

1. Calculer les valeurs propres et les déterminants de Ai pour i = 1, .., 3 2. Pour chacune de ces matrices, trouver s’il existe P et D tel que Ai = P DP −1 . 3. Déterminer le rayon spectrale, rang et la trace de ces matrices Ai Exercice 3 On note      5 1 5/8 −1/4 1/8 A =  1/4 0 1/4  , b =  −1  , u0 =  2  , −4 1 1/8 −1/4 5/8 

et on définit, pour n ≥ 0, la suite de vecteurs un+1 = Aun + b. 3

1. Ecrire une fonction "Exercice3.m" qui permet de calculer les termes de la suite un . La syntaxe de début du programme est la suivante : function [u] = Exercice3(n) A=[5/8 -1/4 1/8; 1/4 0 1/4; 1/8 -1/4 5/8]; % on définit la matrice A b=[1 ; 1; 1]; % le vecteur b u0=[5;2;-4]; % le vecteur u0 u=... for i= ? : ? u= ... end 2. En déduire la limite de la suite un

1.1

Répresentation graphique sous Matlab

Pour construire des graphiques sous matlab on peut utilise la fonction "plot" voici quelques paramètres : Les couleurs

La représentation des coordonnées

Le type du trait

b = bleu g = vert r = rouge

m = magenta k = noir w = blanc

. point o cercle x x-marque

+ plus * étoile s carré

- solide : discontinu 1

-. discontinu 2 – discontinu 3

.

d diamant v triangle < triangle (gauche)

Exemple d’application : exécuter les commandes suivantes

4

> triangle (droite) p pentagramme h hexagramme

.

Fonction plot & comet Construire des graphes

Petit animation sous matlab

>> t=[-4*pi :0.001 :4*pi] ;

>> t=[-4*pi :0.001 :4*pi] ;

>> plot(t,sin(t),’d-’)

>> comet(t,t.*sin(t))

>> plot(cos(t), sin(t),’s :’)

>> comet(t.*cos(t), t.*sin(t))

>> plot(t,t.*sin(t),’> comet(sin(t)3 , 1 + cos(t) − cos(t)4 )

1.2

fonction subplot

Pour construire plusieurs graphes dans une seule image Fonction subplot >> >> >> >> >> >>

t = 0 :pi/10 :2*pi ; [X,Y,Z] = cylinder(4*cos(t)) ; subplot(2,2,1) ; mesh(X) ; title(’X’) ; subplot(2,2,2) ; mesh(Y) ; title(’Y’) ; subplot(2,2,3) ; mesh(Z) ; title(’Z’) ; subplot(2,2,4) ; mesh(X,Y,Z) ; title(’X,Y,Z’) ;

Exercice 4 Soient f (x) = 10x − 9e−x et g(x) =

x 2

− sin(x2 )

1) Ecrire les fonctions f et g sous Matlab, ensuite calculer f (0) , f ′ (0), g(0) et g ′ (0) 2) Construire les fonctions f et g sur [0, 1]. 3) Montrer graphiquement que f et g admettent des racines dans [0, 1]

5

TP no 2 Approximation des racines d’une fonction réelle

Objectif de TP L’objectif de ce tp est d’établir les méthodes d’approximation des racines d’une fonction réelle, c’està-dire la résolution approchée du problème suivant : étant donné f :]a, b[⊂ R → R, trouver α ∈ C tel que f (α) = 0 Commençons tout d’abord par cet exemple, considérons l’équation : 10x − 9e−x = 0 Les méthodes usuelles de transformation (transposition, utilisation de la fonction logarithme,...) ne permettent pas de résoudre explicitement cette équation. Pourtant on observe graphiquement qu’elle admet une solution unique sur [0, 1]. Pour cela, il suffit de taper sur le workspace de Matlab les commandes : x=linspace(0,1,50); y=10*x-9*exp(-x); figure; plot(x,y);xlabel(’x’);ylabel(’y’); Dans ce sujet, on va exposer les principales méthodes itératives de résolution d’une équation de la forme f (x) = 0, où f est une fonction continue définie sur un intervalle [a, b]. On se placera dans le cas où, localement, il y a une unique racine, pour en donner une algorithme d’approximation.

2

Méthode de Dichotomie

La méthode de dichotomie est fondée sur la propriété suivante : Proposition (Théorème des zéros d’une fonction continue) Soit une fonction continue f : [a, b] → R, si f (a)f (b) < 0, alors ∃α ∈]a, b[ tel que f (α) = 0. Cette méthode consiste en une succession de divisions par deux de l’intervalle pour approcher de plus en plus la racine de l’équation f (x) = 0, jusqu’à ce qu’une précision ε soit atteinte. a+b i) On partage [a, b] en deux intervalles égaux [a, a+b 2 ] et [ 2 , b]. a+b ii) Si le signe de f ( a+b 2 ) est le même que celui de f (a), la racine c appartient à l’intervalle [ 2 , b]. a+b Sinon, elle appartient à l’intervalle [a, 2 ].

iii) On réitère le procédé avec l’intervalle obtenu contenant c. vi) On arrête l’itération lorsque la longueur de l’intervalle devient inférieure à un nombre ε fixé au départ.

6

Remarque : A l’étape n, c appartient à l’intervalle de travail, qui a pour longueur b−a 2n D’une itération à la suivante, l’erreur est donc multipliée par 21 . Algorithme de dichotomie

                                          

Donner f, ε > 0, a et b n = 0;                           

tant que b − a > ε faire c=

a+b 2 ;

si f (x) = 0 alors afficher la valeur exacte de la racine de f est ”c”    sortir        si f (x)f (a) < 0 alors b = c;   si non    si non a = c; 

n = n + 1;

afficher la valeur approchée du point fixe de f est ”c” afficher le nombre d’itérations pour aboutir à ce résultat est ”n”

Une valeur approchée à ε près de la racine α est c.

2.1

Application

On reprend l’équation initiale 10x − 9e−x = 0. Posons la fonction f sur l’intervalle [0, 1] par f (x) = 10x − 9e−x , ′

est continue, dérivable et sa dérivée f vérifie ′

f (x) = 10 + 9e−x > 0, ∀x ∈ [0, 1], 7

donc f est strictement croissante sur [0, 1]. 9 De plus f (0) = −9 et f (1) = 10 − ≈ 6.69 sont de signes contraires. e On peut donc utiliser la méthode de dichotomie pour calculer à 10−6 près la solution de l’équation proposée.

Remarque : La méthode de dichotomie a l’avantage d’exiger peu d’hypothèse sur la fonction. Elle sert parfois de moyen de calcul d’une initialisation pour les algorithmes des autres méthodes. L’incovénient majeur de cette méthode est la lenteur de convergence de son algorithme.

3

Méthode du point fixe

Parmi les méthodes de résolution de l’équation f (x) = 0, la méthode dite des approximations successives (ou du point fixe) est la plus importante. Son principe est basé sur la constructiion d’une suite itérative approchant de plus en plus la racine exacte, son premier élément (appelé initialisation) pouvant être n’importe quel point de l’intervalle de travail [a, b]. La méthode du point fixe s’applique à des équations de la forme ϕ(x) = x. A titre de rappel, on peut toujours écrire l’équation f (x) = 0 sous une forme équivalente de ce type. Par exmple, l’équation 10x − 9e−x = 0 est équivalente à x=

9 −x e . 10

On se place dans le cas où la fonction ϕ : [a, b] → [a, b] vérifie les hypothèses : (H1) ϕ est continue et dérivable sur [a, b]. ′

(H2) ∃M ∈]0, 1[: ∀x ∈ [a, b] |ϕ (x)| ≤ M

On dira que ϕ est contraction stricte. La méthode du point fixe est basée sur la proposition suivante. Proposition : Lorsque ϕ vérifie les deux hypthèses (H1) et (H2), alors il existe une unique racine c de l’équation ϕ(x) = x, applée point fixe de ϕ. D’où l’algorithme de la méthode du point fixe :

8

Algorithme du point fixe                           

3.1

on choisit x0 ∈ [a, b], ε > 0 et Nmax le nombre maximal d’itérations n = 0; Répéter    xn+1 = ϕ(xn );     n = n + 1;   

Jusqu’à |xn − xn−1 | ≤ ε ou n ≥ Nmax

afficher la valeur approchée du point fixe de ϕ est ”xn ” afficher le nombre d’itérations pour aboutir à ce résultat est ”n”

Application

Pour calculer à 10−6 près la solution, dans l’intervalle [0, 1], de l’équation x=

9 −x e , 10

par la méthode du point fixe, on procède comme suite 1. On définit la fonction ϕ, telle que

9 −x e , 10 cette fonction est continue et dérivable sur [0, 1]. ϕ(x) =

(a) Vérifier les hypothèses (H1) et (H2) de la fonction ϕ (b) Ecrire une fonction matlab "PointFixe.m" pour approcher la racine α de la fonction f 2. Utiliser ce programme pour déterminer la racine de la fonction f2 (x) = x2 + ln(x) sur [1/2, 1]

4

Méthode de Newton

La méthode de Newton est une méthode de point fixe pour la fonction

9

f (x) f ′ (x) si n) xn − ff′(x (xn )

ϕ(x) = x −

xn+1 =

f ′ (x) 6= 0

Algorithme de Newton

                          

4.1

donner f, f ′ , ε > 0, x0 et Nmax le nombre maximal d’itérations n = 0; 

Répéter    f (x )  xn+1 = xn − f ′ (xnn ) ;     n = n + 1;   Jusqu’à |xn − xn−1 | ≤ ε ou n ≥ Nmax afficher la valeur approchée du point fixe de ϕ est ”xn ” afficher le nombre d’itérations pour aboutir à ce résultat est ”n”

Application 1. Ecrire une fonction "Newton.m" basé sur la méthode de Newton pour approcher la racine α de la fonction f donner une estimation de cette racine à ε = 10−6 prés ainsi que le nombre d’itérations nécessaires pour l’obtenir. 2. Dans quel intervalle doit-on choisir la valeur initiale pour asurer la convergence. 3. Ecrire une programme qui permet de comparer entre les trois méthodes

5 5.1

Quelques problème physique à résoudre Problème 1 (Équation d’état d’un gaz)

Nous voulons déterminer le volume V occupé par un gaz dont la température est T et dont la pression est p. L’équation d’état est donnée par :

10



p+a

N2  (V − N b) = kN T. V2

Où a et b sont deux coefficients qui dépendent du gaz considéré, N est le nombre de molécules contenus dans le volume V et k est la constante de Boltzmann. Nous devons donc résoudre une équation non linéaire dont la racine est V0 .

Exemple du CO2 :

— Pour CO2 (dioxyde de carbone) les coefficients a et b dans ce problème prennent les valeurs suivantes : a = 0.401P a.m6 , b = 42.7 × 10−6 m3 (P a signifie Pascal). — Question : Trouver le volume occupé par 1000 molécules de CO2 à la température T = 300k et la pression p = 3.5 × 107 P a, par les deux méthodes (dichotomie et Newton), avec une tolérance de 10−12 (la constante de Boltzmann vaut k = 1.3806503 × 10−23 Joule.K −1 ). On considére l’intervalle ]0.01, 0.06[.

5.2

Problème 2 (mécanique statique)

Considérons le système mécanique représenté par les quatre barres rigides ai de la Figure 1. Pour une valeur admissible de l’angle β, déterminons la valeur de l’angle α entre les barres a1 et a2 . Partant de la relation vectorielle

a1 − a2 − a3 − a4 = 0 et remarquant que la barre a1 est toujours alignée avec l’axe des x, on peut déduire les relations suivantes entre β et α.

a1 a2 + a22 − a23 + a24 a1 cos(β) − cos(α) − cos(β − α) = − 1 a2 a4 2a2 a4 où ai est la longueur connue de la i-ème barre. Cette égalité, appelée équation de Freudenstein, peut être récrite comme suit : f (α) = 0, où

f (x) =

a1 a1 a2 + a22 − a23 + a24 cos(β) − cos(x) − cos(β − x) + 1 a2 a4 2a2 a4 11

Figure 1 – Le système de quatre barres du problème statique Une expression explicite de la solution n’existe que pour des valeurs particulières de β. Signalons également qu’il n’y a pas existence d’une solution pour toutes les valeurs de β, et qu’une solution peut ne pas être unique. Pour résoudre cette équation pour toute valeur de β entre 0 et π, nous devrons avoir recours à des méthodes numériques

Question : Utiliser la méthode de Newton pour résoudre le problème 2 pour β ∈ [0, 2π 3 ] avec une tolérance

de 10−5 . Supposer que les longueurs des barres sont a1 = 10cm, a2 = 13cm, a3 = 8cm et a4 = 10cm. Pour chaque valeur de β, considérer deux valeurs initiales, x0 = −0.1 et x0 = 2π 3

12

TP no 3 Résolution des systèmes linéairess (Méthodes Directes)

6

Objectif de TP

Dans ce TP, on va s’intéresser à la résolution du système matricielle de la forme Ax = b où A une matrice carrée connue et b un vecteur donné et x le vecteur inconnu. Dans un premier temps, nous allons utiliser deux méthodes de résolution : l’algorithme d’élimination de Gauss et la méthode de factorisation LU. ⊗ La méthode d’élimination de Gauss (avec et sans pivot) est une méthode pour transformer un système en un autre système équivalent (ayant les mêmes solutions) qui est triangulaire et donc facile à résoudre. Les opérations autorisées pour transformer ce système sont : • échange de deux lignes. • multiplication d’une ligne par un nombre non nul. • addition d’un multiple d’une ligne à une autre ligne. On est alors ramené à la résolution d’un système linéaire T x = ˜b, où la matrice T est triangulaire et inversible. Ceci se fait très facilement par substitution récursive. On appelle ce procédé remontée dans le cas d’une matrice T triangulaire supérieure et descente dans le cas d’une matrice T triangulaire inférieure. Il faudra remarquer que l’on résoud ainsi le système T x = ˜b sans inverser la matrice T. ⊗ La deuxième méthode proposée consiste à factoriser la matrice A en un produit de deux matrices triangulaires A = LU, où L est triangulaire supérieure (L pour lower en Anglais) et U est triangulaire sup ’erieure (U pour upper en Anglais). Il s’agit en fait du même algorithme que celui de l’élimination de Gauss dans le cas particulier où on ne pivote jamais. Une fois établie la fatorisation LU de A, la résolution du système linéaire Ax = b est équivalente à la simple résolution de deux systèmes triangulaires Ly = b puis U x = y, soit LU x = b.

7 7.1

Sujet de TP Cas d’une matrice triangulaire inférieure ou supérieure 1. Ouvrir un fichier qu’on nomme trianginf.m. Dans ce fichier, on va créer la fonction Matlab qui à partir d’une matrice inférieure A et d’un second membre b, résout le système AX = b. avec Err=||AX-b||: l’erreur de la méthode. function [X,Err]=trianginf(A,b) % A est la matrice triangulaire inférieure et b le vecteur second membre. n=length(b); x(1)= for i=2:n x(i)=... end 13

Algorithme de le cas triangulaire inférieure Donner une matrice A : triangulaire inférieure et un vecteur b. La solution x est définie par :  b1   x = 1   a11   i−1

      xi

=

i X 1 h aij xj , bi − aii

Pour i = 2 jusqu’à n

j=1

     1 3 0 0 4 0 0 Application 1 : soient A1 =  9 2 0  , A2 =  7 4 0  et b =  5  , 3 8 0 3 7 −6 6 

Trouver les solutions et l’erreurs de l’équation Ai X = b pour i = 1, 2 2. Ouvrir un fichier qu’on nomme triangsup.m. Dans ce fichier, on va créer la fonction Matlab qui à partir d’une matrice supérieure A et d’un second membre b, résout le système Ax = b. La syntaxe de début du programme est la suivante : function [X,Err]=triangsup(A,b) % A est la matrice triangulaire supérieure et b le vecteur second membre. n=length(b); x(n)= for i=n-1:-1:1 % boucle incrément négatif x(i)=... end Algorithme de le cas triangulaire supérieure Donner une matrice A : triangulaire supérieure et un vecteur b. La solution x est définie par :  bn    xn =  ann       x   i

=

n i X 1 h aij xj , bi − aii j=i+1

14

Pour i = n − 1 jusqu’à 1

     3 0 −1 1 4 0 2 Application 2 : soient A1 =  0 2 4  , A2 =  0 4 2  et b =  5  , 3 0 0 6 0 0 3 

Trouver les solutions de l’équation Ai X = b, par la méthode du pivot de gauss a

7.2

Méthode de Gauss - Ouvrir un fichier qu’on nomme gauss.m. Dans ce fichier, on va créer la fonction Matlab qui à partir d’une matrice A effectue de Gauss et résout le système Ax = b. La syntaxe de début du programme est la suivante :

function [X,Err]=gauss(A,b) % A est la matrice principale et b le vecteur second membre. n=length(b); ..... ... end     2 3 3 1   1 2 3  −4 −6 3 2  ε 1  Application 3 : soient A1 =  4 5 6  , A2 =   −1 1 1 1  , A3 = 1 1 , 7 8 9 −2 −1 1 1   15    3   et b3 = 1 , avec ε = 0.01, 0.001, 0.0001. b2 =   5  2 1



 1 b1 =  2  , 3

Trouver les solutions et l’erreurs de l’équation Ai X = bi , par la méthode du pivot de Gauss.

15

Algorithme du pivot de Gauss Entrée : A = (aij )1≤i,j≤n et b = (bi )1≤i≤n . Sortie x = (xi )1≤i≤n la solution du système linéaire. % Triangulation de la matrice k=1 ; Tant que (k < n) faire  Pivot = akk ;    Si (Pivot 6= 0)      Faire i = k + 1 à n :    a   mik = ik   a kk          Faire j = k à n      aij = aij − mik akj     Fin de la boucle pour j         bi = bi − mik bk   Fin de la boucle pour i     k = k + 1;      Sinon      l = k;          Tant que (alk = 0)          l = l + 1;   

al. = ak. Pivotage Ll = Lk

% Résolution du système triangulaire supérieure bn ann Faire i = n − 1 à 1 n X
1 aij xj bi − xi = aii xn =

j=i+1

Fin de la boucle pour i

16

7.3

Méthode de décomposition LU - Le but de cette section est d’écrire une fonction ResolutionLU.m sous Matlab qui permet à partir d’une matrice A et un vecteur b, résout le système Ax = b avec de la factorisation LU de la matrice A. La syntaxe de début du programme est la suivante 1 Algorithme de LU la fonction "DecompLU.m" Algorithme de LU function [L,U] = DecompLU(A) Entrées A = (aij ), 1 ≤ i, j ≤ n (la matrice principale) Sorties : L = (lij ), 1 ≤ j ≤ i, 1 ≤ i ≤ n & U = (uij ), i ≤ j ≤ n, 1 ≤ i ≤ n l11 = 1 Faire j = 1 à n u1j = a1j % Détermination de la première ligne de U Fin de la boucle pour j Faire i = 2 à n i1 li1 = aa11 % Détermination de la première colonne de L Fin de la boucle pour i Faire i = 2 à n lii = 1 i−1 X lik uki uii = aii − k=1

Faire j = i + 1 à n i−1 X lik ukj uij = aij − k=1

i−1 i X 1 h ljk uki lji = aji − uii k=1 Fin de la boucle pour j Fin de la boucle pour i

2 Résolution d’un système linéaire avec méthode de LU (la fonction "ResolutionLU.m") La fonction ResolutionLU.m function [X, Err] = ResolutionLU(A,b) [L, U ] = DecompLU(A) ; Y=trianginf(L,b) ; [X, Err]=triangsup(U,Y) ; end



     1 2 3 1 2 3 1 Application 4 : soient A1 =  4 5 6  , A2 =  4 5 6  et b =  5  , 1 8 10 7 8 9 3 17

Trouver les solutions et l’erreurs de l’équation Ai X = b pour i = 1, 2

8

Travail à faire par les étudiants (compte rendu)

Partie I (Méthode de Cholesky ) : On rappelle la méthode de Cholesky permettant de factoriser une matrice A symétrique définie positive sous la forme A = CC T , où C est une matrice triangulaire inférieure inversible. Algorithme de Cholesky function [C] = cholesky(A) √ a11 , ai1 , i > 1. c11

c11 = ci1 =

cjj

v u j−1 u X c2jk , j > 1, = tajj − k=1

cij

=

j−1  X 1  aij − cik cjk , i = j + 1, ..., n. cjj k=1

(i) Ecrire une fonction Matlab [C]=cholesky(A) qui, étant donné la matrice A calcule la matrice C. Application : On définit   15 10 18 12  10 15 7 13   A=  18 7 27 7  12 13 7 22

calculer C. Est-ce que A est définie postive ? (ii) Ecrire une fonction Matlab [x]=resolchol(A,b)

qui, étant donnés la matrice A symétrique définie positive et le vecteur b, calcule la solution x de Ax = b en utilisant la fonction cholesky et les fonctions triangsup et trianginf. Application : Résoudre la système Ax = b avec la matrice A précédente et   53  72   b=  26  97 18

Partie II (Méthode de Richtmayer) : L’algorithme de Richtmayer permet de résoudre un système d’équations linéaires Ax = u lorsque A est une matrice tridiagonale, 

b1  a2    0   ..  A= .  ···   ..  .   0 0

c1 b2

0 c2

··· 0

··· ···

··· ···

a3 .. .

b3 .. .

c3 .. .

··· .. .

··· .. .

··· .. . ··· .. .

··· ···

··· ···

··· .. . 0 ···

0 .. . ··· .. .

0 0 .. . .. . ··· .. .

an−1 bn−1 cn−1 0 an bn

L’algorithme de Richtmayer est donnée comme suit



       .      

Algorithme de Richtmayer function [x] = Richtmayer(a,b,c,u) u1 c1 , f1 = , b1 b1 ci , pour 2 ≤ i ≤ n − 1 = − ai ei−1 + bi ui − ai fi−1 = , pour 2 ≤ i ≤ n. ai ei−1 + bi

e1 = − ei fi

xn = f n . xi

= ei xi+1 + fi , pour i = n − 1, ..., 1

Ecrire une fonction Matlab [x]=rich(a,b,c,u) qui, étant donné les vecteurs a, b, c, u, calcule la solution du système Ax = b par la méthode précédente. Application : utiliser la fonction précédente pour déterminer la solution de Ax = u dans le cas :     2 3 0 0 0 −1  1 3 1 0 0   0         A =  0 2 1 2 0  et u =   0 .  0 0 1 1 1   3  0 0 0 2 1 1

19

TP no 4 Résolution des systèmes linéairess (Méthodes Itératives)

9

Méthode de Jacobi

(a) Soit A une matrice carrée inversible, écrire une fonction Matlab [x,k]=jacobi(A,b,N,ε,x0 ) qui, étant donné la matrice A et le vecteur b, calculer la solution de Ax = b par la méthode de Jacobi où : • ε est la précision utilisée pour le test d’arrêt. On pourra tester par exemple ||x(k+1) − x(k) ||22 < ε, ||x(k+1) ||22 • N est un nombre maximal d’itérations, • x0 est un vecteur initial,

• k est le nombre d’itérations éffectuées pour obtenir la précision voulue.

Pour arrêter l’excution de la fonction en cas de convergence, vous pouvez utiliser la commande return. Lorsque le nombre maximal d’itérations est atteint avant que la précision voulue ne soit obtenue, affincher un message d’erreur signalant que la méthode de Jacobi n’a pas convergé.

9.1

Algorithme de Jacobi

pour la résolution d’un système linéaire Ax = b, où A = (aij )1≤i,j≤n est une matrice avec les aii 6= 0. On décomposer la matrice A sous cette forme : — A=D−E−F — D : diagonale de A. — E : triangulaire inférieure avec des 0 sur le diagonale. — F : triangulaire supérieure avec des 0 sur le diagonale.

En suite, l’équation :

AX = b

⇔

DX = (E + F )X + b

20

⇔

X = X − D −1 AX + D −1 b

Algorithme Méthode de Jacobi  (0) X : vecteur initial         D := diagonale de A.          ||X(k+1) − X(k) ||2  2 Tant que ≥ ε && k < N  ||X(k+1) ||22         X(k+1) = X(k) − D −1 AX(k) + D −1 b        k = k + 1;

Faire



     10 −2 1 2 −1 1 9 Application 1 : soient A1 =  −2 10 −2  , A2 =  2 2 2  et b =  12  , −2 −5 10 −1 −1 2 18 1. Trouver les solutions et les erreurs de la résolution de Ai X = b pour i = 1, 2, par la méthode de Jacobi.

10

Méthode de Gauss-Seidel

(b) Soit A une matrice carrée inversible, écrire une fonction Matlab [x,k]=gaussseidel(A,b,N,ε,x0 ) qui, étant donné la matrice A et le vecteur b, calcule la solution de Ax = b par la méthode de Gauss-Seidel où : • ε est la précision utilisée pour le test d’arrêt. On pourra tester par exemple ||x(k+1) − x(k) ||22 < ε, ||x(k+1) ||22 • N est un nombre maximal d’itérations, • x0 est un vecteur initial,

• k est le nombre d’itérations éffectuées pour obtenir la précision voulue.

Pour arrêter l’excution de la fonction en cas de convergence, vous pouvez utiliser la commande return. Lorsque le nombre maximal d’itérations est atteint avant que la précision voulue ne soit obtenue, affincher un message d’erreur signalant que la méthode de Gauss-Seidel n’a pas convergé.

21

10.1

Algorithme de Gauss Seidel

pour la résolution d’un système linéaire Ax = b, où A = (aij )1≤i,j≤n est une matrice avec les aii 6= 0. On décomposer la matrice A sous cette forme : — A=P −N — P := D − E triangulaire inférieure avec des éléments diagonaux de A sur le diagonale. — N := F triangulaire supérieure avec des 0 sur le diagonale.

En suite, l’équation :

AX = b

⇔

PX = NX + b

⇔

X = X − P −1 AX + P −1 b

Algorithme Méthode de Gauss Seidel  (0) X vecteur initial           ||X(k+1) − X(k) ||2   2  ≥ ε && k < N Tant que  ||X(k+1) ||22      X(k+1) = X(k) − P −1 AX(k) + P −1 b        k = k + 1;

Faire



     10 −2 1 2 −1 1 9 Application 2 : soient A1 =  −2 10 −2  , A2 =  2 2 2  et b =  12  , −2 −5 10 −1 −1 2 18 1. Trouver les solutions et les erreurs de la résolution de Ai X = b pour i = 1, 2, par la méthode de Gauss-Seidel. 2. Comparer les résultats.

11

Exemple d’application en physique

On considère le déplacement vertical u(x) d’un câble, représenté au repos par le segment [a, b] = [0, 1], et fixé aux extrémités, sur lequel on applique une force f (x). Ce déplacement u est la résolution de l’équation différentielle suivante :  ′′ −u (x) = f (x), pour x ∈]0, 1[, (1) u(0) = 0, u(1) = 0 Soit N ∈ N, h = N1 et xi = ih pour i = 1, ..., .N Pour approcher la solution u(x) on considère la discrétisation de l’intervalle [0, 1] en N sous-intervalles

22

[xi , xi+1 ], et on construit ui de u(xi ) par la méthode des différences finies. Cette méthode requiert de résoudre numériquement le système linéaire tridiagonale Ax = b qui suit 

2 −1 0  −1 2 −1    0 −1 2  . .. .. 1   . . . Au = b ⇔ 2  . h  ··· ··· ···   .. .. ..  . . .   0 ··· ··· 0 ··· ···

··· 0

··· ···

··· ···

0 0 .. . .. . ··· .. .

 

u1 u2 u3 .. . .. . .. .

        −1 0 · · ·     .. .. ..   . . . .   ··· ··· ···     .. .. ..   . . .   0 −1 2 −1   uN −2 ··· 0 −1 2 uN −1





f (x1 ) f (x2 ) f (x3 ) .. . .. . .. .

              =             f (xN −2 ) f (xN −1 )

              

Plus N est grand, plus l’approximation sera précise et plus la taille du système linéaire à résoudre sera élevée. 1. Résoudre analytiquement l’équation différentielle (1) pour f (x) = x(1 − x) - Indication : la solution analytique de cette équation est Z Z 1 Z x F (y)dy avec F (y) = F (y)dy + x u(x) = − 0

y

f (t)dt.

0

0

2. On suppose que la force appliquée soit f (x) = x(1−x) et on prend N = 20 intervalles. Construire la matrice A et le vecteur b correspondants, à l’aide des commandes suivantes f=’x.*(1-x)’; N=20; h=1/N; x=linspace(0,1,N-1); b=eval(f); A = (N.2 ) ∗ (diag(2 ∗ ones(n − 1, 1), 0) − diag(ones(n − 2, 1, 1)) − ...diag(ones(N − 2, 1), −1))x ;

3. Calculer la solution du système linéaire Au = b à partir de la factorisation LU de A. Utiliser les commandes tic et toc pour calculer le temps nécassaire (voir exemple ci-dessous). Ensuite comparer avec le temps de calcul pour résoudre directement : tic; u=A; toc; Refaites la même comparaision avec factorisation de Cholesky de A. ¯ 4. Programmer un programmme qui résout cette équation différentielle par les différences finies, et reprétente le déplacement u du câble aux noeuds xi de l’intervalle [0, 1]. Comparer graphiquement le déplacement approché avec la solution exacte du problème (1).

23

TP no 5 Interpolation Polynômiale (Lagrange, Phénomène de Runge)

12

Position du Problème 1. Étant donné (n + 1) points du plan (x0 , y0 ), (x1 , y1 ), . . . , (xn , yn ) trouver une fonction simple qui passe par ces (n + 1) points ; 2. Étant donné une fonction f : I −→ R connue seulement en certains points, trouver une focntion qui est une bonne approximation de f .

Un choix naturel de l’approximation est la fonction polynomiale, considérée comme étant d’autant plus simple que son degré est bas. On s’intéresse alors à présenter quelques formules numériques permettant d’approcher la valeur de la fonction à interpoler en un point donné vénient de cette approximation polynomiale. Dans ce TP, on va s’intéresser à l’implémentation des algorithmes permettant : — de calculer le polynôme d’interpolation de Lagrange d’une fonction f définit sur un intervalle [a, b]. — de constater le phénomème de Runge.

13

Interpolation de Lagrange

Soient n ∈ N et (xi , yi )0≤i≤n , (n + 1) points du plan deux à deux distinctes. Alors, il existe un unique polynôme P de degré inférieur ou égal à n qui passe par le (n + 1) points. ce polynôme présenté dans la base de Lagrange et appélé "polynôme d’interpolation de Lagrange" est donné par la formule suivante :

Pn (x) =

n X

f (xi )Li (x)

i=0

Où, pour tout i = 0, . . . , n on a

Li (x) =

n Y

j=0,j6=i

(x − xj ) , (xi − xj ) 24

0≤i≤n

Ainsi, l’algorithme du calcul d’une valeur approchée de f (t) pour t donnée, par l’interpolation de Lagrange est le suivant : Algorithme d’interpolation de Lagrange

                                    

Donner deux vecteurs X = (xi )1≤i≤n , Y = (yi = f (xi ))1≤i≤n et t P = 0; Pour i = 1 à n faire                    

L = 1; 

Pour j = 1 à i − 1    (t−x )  L = L (xi −xjj ) ;     Pour j = i + 1 à n    (t−x ) L = L (xi −xjj ) ; P = P + y(i)L;

Err =

Rb a

(f (t) − P )2 dt

1. Ouvrir un fichier qu’on nomme lagrange.m Dans ce fichier, on va créer la fonction Matlab qui à partir des points d’interpolation (xi ) et des valeurs yi = f (xi ), calcule le polynôme d’interpolation de Lagrange sur un intervalle [a, b].

13.1

Application 1 (fonction Sinus)

1. Valider votre programme par le choix d’une fonction polynomiale f (x) = x3 + x2 + 2. 2. Appliquer votre programme pour une fonction f (x) = sin(x) sur l’intervalle [a, b] = [−π, π] avec  des points d’interpolation équidistants xi = a + (b − a) Ni ), i = 0, . . . , N . 3. Tracer les deux fonctions dans une seule graphique f (x) = sin(x) et son polynôme d’interpolation de Lagrange Pn sur un intervalle [c, d] plus grand que [a, b], pour différent valeur n = 2, 10, 100, dans un Script sous matlab on nomme AppLagrangeSINUS.m.

25

13.2

Application 2 (Phénomène de Runge)

1 1. Phénomène de Runge : Refaire la même chose avec la fonction f (x) = 1+x 2 sur l’intervalle [a, b]  i aux points xi = a + N (b − a), i = 0, . . . , N . Définir un nombre de points N = 5, 10, 100, a = −5 et b = 5

2. Créer une discrétisation points.

n xi = a +

i N (b

− a),

i = 0, . . . , N

o

du segment x = [a, b] de N + 1

3. Interpoler cette fonction f à l’aide de la fonction lagrange.m. On tracera le graphe de f et de son interpolée à l’aide de la commande plot, exemple : plot(x,y) ; On prendra N de plus en plus grand. Remarque : On constate que avec l’emplacement des points équidistance (xi )0≤i≤N le polynome de Lagrange dans ce cas ne converge pas uniformenet vers f (x) qui c’est pas le cas pour la fonction f (x) = sin(x).

13.3

Application 3 (Les points de Tchebycheff)

1. Pour corriger n les erreurs dans l’exemple précédente ono définir les points de Tchebycheff sous (a−b) iπ i = 0, . . . , N cette forme x(i) = a+b 2 + 2 cos( N ),

1 2. Tracer les deux fonctions dans une seule graphique f (x) = 1+x 2 et son polynôme d’interpolation de Lagrange Pn définie par les points de Tchebycheff sur l’intervalle [a, b], pour différent valeur n = 2, 10, 100.

3. Que remarquez vous ?

13.4

Exemple d’un problème climatologie

La température moyenne de l’air au voisinage du sol dépend de la concentration K en acide carbonique (H2 CO3 ). Dans la Table 1, on donne la variation δK = θK − θK˜ de la température moyenne par rapport à une température de référence K est pour différentes latitudes et pour quatre valeurs de K. ˜ la valeur de K mesurée en 1896 et normalisée à un. On peut construire une fonction La quantité K qui, sur la base des données disponibles, permet d’approcher la température moyenne à une latitude quelconque et pour d’autres valeurs de K (voir Exemple 3.1)

26

Table 1 : Variation de la moyenne annuelle des températures sur la Terre pour différentes valeurs de la concentration K en acide carbonique à différentes latitudes (d’après Philosophical Magazine 41, 237 (1896))

Question : calculer le polynôme d’interpolation des données du Problème correspondant à K = 0.67 (première colonne de la Table 1), en utilisant seulement les valeurs de la température pour les latitudes 65, 35, 5, −25, −55.

27

TP no 6 Intégration Numérique

Objectif du TP Comme dit le proverbe, Dérive qui veut, intègre qui peut. Le calcul numérique des intégrales joue donc un rôle dans de nombreux problèmes physiques. Il s’agit ici d’écrire une série de procédures permetZ b f (x)dx, et de comparer leurs performances (vitesse de tant de calculer numériquement l’intégrale a

convergenceZ pour une précision requise). Durant ce TP, on pourra plus particulièrement tester le cas π sin(x)dx = 2, mais les procédures doivent pouvoir s’appliquer à toute fonction f définie analytique 0

(et intégrable) sur le domaine choisi.

Sujet du TP 1. Ouvrir un fichier qu’on nomme trapeze.m. Dans ce fichier, on va créer la fonction Matlab qui à partir des points d’interpolation (xi ) et des valeurs yi = f (xi ), permet d’approximer Z

b

f (x) dx

a

à l’aide de la méthode des trapèzes. 2. Tester la méthode sur f (x) = sin(x) sur [0, π]. 2

3. Tester la méthode sur f (x) = e−x sur [a = 0, b = 7] Que se passe t’il quand on augmente le nombre de points (N = 10, 100, 1000). 4. Ouvrir un fichier qu’on nomme simpson.m. Dans ce fichier, on va créer la fonction Matlab qui à partir des points d’interpolation (xi ) et des valeurs yi = f (xi ), permet d’approximer Z

b

f (x)dx,

a

à l’aide de la méthode de Simpson. 5. Tester la méthode sur f (x) = sin(x) sur [0, π]. 2

6. Tester la méthode sur f (x) = e−x sur [a = 0, b = 7] et comparer les résultats avec ceux obtenus par la méthode des trapèzes.

28