Sujets TP Mecanique [PDF]

Zitiervorschau

ENIM-Monastir 2019-2020

RECUEIL DES SUJETS DE TP D’ANALYSE NUMÉRIQUE

TP Analyse numérique Génie Mécanique (Meca-1)

Équipe de Mathématique de l’ENIM :2019-2020

1

TP No :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 de 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



     2 2 −4 u1 =  3  , u2 =  −1  et u3 =  −3  −1 5 6

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

Exercice 2.

Soient A1 , A2 et A3 les trois matrices suivantes       4 0 2 3 0 −1 1 3 A1 = , A2 =  0 2 4  et A3 =  2 4 2  . 5 2 0 0 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 2

- Indication : utiliser les fonctions (det, eig, rank, vrho,trace)

Exercice 3.

On note



     5/8 −1/4 1/8 1 5 0 1/4  , b =  −1  , u0 =  2  , A =  1/4 1/8 −1/4 5/8 1 −4

et on définit, pour n ≥ 0, la suite de vecteurs un+1 = Aun + b. 1. Ecrire une fonction "TP1_Exercice3.m" qui permet de calculer les termes de la suite un . La syntaxe de début du programme est la suivante : function [u] = TP1_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

0.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 3

> 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)

0.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 1.

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 (0), g(0) et g 0 (0) 2) Construire les fonctions f et g sur [0, 1]. 3) Montrer graphiquement que f et g admettent des racines dans [0, 1]

4

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.

1

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 ].

5

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. 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 function [xn ,n]=Dichotomie(f, ε, a, b) % Entrée :f, ε > 0, x0 , a et b % Sortie : xn et n.                                              

n = 0;                              

  Tant que b − a > ε Faire c=

a+b 2 ;

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

n = n + 1; xn = c; 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. 6

1.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 , 0

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

f (x) = 10 + 9e−x > 0, ∀x ∈ [0, 1], 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. Données numériques : f (x) = 10x − 9e−x , a = 0, b = 1 et ε = 10−6 . 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.

2

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]. 0

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

7

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 :

Algorithme du point fixe function [xn ,n]=PointFixe(ϕ, ε, x0 , Nmax ) % Entrée : ϕ, ε > 0, x0 et Nmax le nombre maximal d’itérations. % Sortie : xn et n.                        

n = 0; xn+1 = ϕ(xn ); 

  Tant que |xn+1 − xn | > ε et n ≤ Nmax Faire

    xn+1 = ϕ(xn );   n = n + 1; afficher la valeur approchée du point fixe de ϕ est ”xn ” afficher le nombre d’itérations pour aboutir à ce résultat est ”n”

2.1

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 (c) Données numériques : ε = 10−6 , x0 = 0.2 et Nmax = 100 8

2. Utiliser ce programme pour déterminer la racine de la fonction f2 (x) = x2 + ln(x) sur [1/2, 1] (a) Données numériques : ε = 10−7 , x0 = 0.7 et Nmax = 100

3

Méthode de Newton

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

f (x) f 0 (x) si n) xn − ff0(x (xn )

ϕ(x) = x − xn+1 =

f 0 (x) 6= 0

Algorithme de Newton function [xn ,n]=Newton(f, ε, x0 , Nmax ) % Entrée : f, ε > 0, x0 et Nmax le nombre maximal d’itérations. % Sortie : xn et n.                         

n = 0; xn+1 = xn − 

f (xn ) f 0 (xn ) ;

  Tant que |xn+1 − xn | > ε et n ≤ Nmax Faire

     xn+1 = xn −   n = n + 1;

f (xn ) f 0 (xn ) ;

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

3.1

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. 9

3. Ecrire une programme qui permet de comparer entre les trois méthodes 4. Données numériques : f (x) = 10x − 9e−x , ε = 10−6 , x0 = 0.3 et Nmax = 1000

4

Quelques problème physique à résoudre

4.1

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 :



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[.

4.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 α. 10

Figure 1 – Le système de quatre barres du problème statique

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 a2 + a22 − a23 + a24 a1 cos(β) − cos(x) − cos(β − x) + 1 a2 a4 2a2 a4

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

11

TP No : 3 Résolution des systèmes linéaires (Méthodes Directes)

5

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.

6 6.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. 12

n=length(b); x(1)= for i=2:n x(i)=... end Algorithme dans le cas triangulaire inférieure function [x,Err]=trianginf(A,b) % 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.  b1   ; x1 =   a11    Faire i=2 à n     i−1   X 1  xi = bi − aij xj ; aii   j=1    Fin boucle i        Err = norm(Ax − b);



     4 0 0 3 0 0 1 Application 1 : soient A1 =  9 2 0  , A2 =  7 4 0  et b =  5  7 −6 6 8 0 3 3 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

13

Algorithme dans le cas triangulaire supérieure function [x,Err]=triangsup(A,b) % 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.

 bn   ; xn =   ann     Faire i=n-1 à 1   n   X  1  bi − xi = aij xj ; aii  j=i+1     Fin boucle i        Err = norm(Ax − b);

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

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

6.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 =  , A3 = , −1 1 1 1  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. 14

Algorithme du pivot de Gauss function [x,Err]=Gauss(A,b) % 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   akk          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 xn =

bn ; ann

Faire i = n − 1 à 1 n X
1 xi = bi − aij xj ; aii j=i+1

Fin de la boucle pour i

15

6.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 ≤ i, j ≤ n, & U = (uij ), 1 ≤ i, j ≤ 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 li1 =

ai1 a11 ;

% Détermination de la première colonne de L

Fin de la boucle pour i Faire i = 2 à n lii = 1; uii = aii −

i−1 X

lik uki ;

k=1

Faire j = i + 1 à n uij = aij −

i−1 X

lik ukj ;

k=1 i−1

lji =

 X 1  aji − ljk uki ; 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")

16

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      4 5 6 5  Application 4 : soient A1 = 4 5 6 , A2 = et b = 1 8 10 7 8 9 3 Trouver les solutions et l’erreurs de l’équation Ai X = b pour i = 1, 2

7

Travail à faire par les étudiants (compte rendu)

Les comptes-rendus sont à rendre au plus tard lors de la prochaine séance.

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.

17

Algorithme de Cholesky function [C] = Cholesky(A) % Entrées : A = (aij ), 1 ≤ i, j ≤ n (la matrice principale). % Sorties : C = (Cij ), 1 ≤ i ≤ n, 1 ≤ j ≤ n. C11 =

√

a11 ;

Faire i = 2 à n Ci1 =

ai1 C11 ;

% Détermination de la première colonne de C

Fin de la boucle pour i Faire j = 2 à n v u j−1 u X 2 ; Cjk Cjj = tajj − k=1

Faire i = j + 1 à n j−1

Cij =

 X 1  aij − Cik Cjk ; Cjj k=1

Fin de la boucle pour i Fin de la boucle pour j

(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 = fn . 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éaires (Méthodes Itératives)

8

Méthode de Jacobi 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é.

8.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 function [x,k]=jacobi(A,b,N,ε,x0 ) % Entrées : A = (aij ), 1 ≤ i, j ≤ n, b = (bi )1≤i≤n , N ,ε et x0 . % Sorties : x = (xi ), 1 ≤ i ≤ n, et k le nombre d’iteration.  D := diagonale de A;         k = 0;          X(k+1) = X(k) − D−1 AX(k) + D−1 b;      ||X(k+1) − X(k) ||2  2  ≥ ε && k < N Tant que    ||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 2 2  et b =  12  , Application 1 : soient A1 =  −2 10 −2  , A2 =  2 −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.   1 −10  0  2. Données Numériques : N = 1000, ε = 10 , x0 = 0

9

Méthode de Gauss-Seidel 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, 21

• 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é.

9.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 function [x,k]= GaussSeidel(A,b,N,ε,x0 ) % Entrées : A = (aij ), 1 ≤ i, j ≤ n, b = (bi )1≤i≤n , N ,ε et x0 . % Sorties : x = (xi ), 1 ≤ i ≤ n, et k le nombre d’iteration.  P := Triangulaire inférieure de A;         k = 0;          X(k+1) = X(k) − P −1 AX(k) + P −1 b;      ||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 2 2  et b =  12  , Application 2 : soient A1 =  −2 10 −2  , A2 =  2 −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.   1 2. Données Numériques : N = 1000, ε = 10−10 , x0 =  0  0 3. Comparer les résultats. 22

10

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 :  00 −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 [xi , xi+1 ], et on construit Xi 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 1 Z Z x u(x) = − F (y)dy + x F (y)dy avec F (y) = 0

0

y

f (t)dt.

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); b=eval(f)’; A = (N.∧ 2) ∗ (diag(2 ∗ ones(N, 1), 0) − diag(ones(N − 1, 1), 1) − diag(ones(N − 1, 1), −1)) ; 3. Calculer la solution du système linéaire AX = 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

Statistiques descriptives Réalisation de tableaux, graphiques, droite de Régression, corrélation

11

Forme des données

En statistique on veillera toujours à ce que les données aient la forme suivante : les variables en colonne et les individus en ligne. La première ligne devra comporter le nom des variables. Si les données n’ont pas initialement sous cette forme on veillera à la leur donner. Ceci simplifiera grandement les traitements ultérieurs dans Excel ou tout autre logiciel statistique.

11.1

Les fonctions sous Excel :

Moyenne, variance, écart-type, médiane, quartiles, centiles

- La fonction MOYENNE permet de calculer la moyenne. - La fonction VAR permet de calculer la variance (attention les autres variantes produisent des valeurs légèrement différentes). - La fonction ECARTYPE permet de calculer l’écart-type. - La fonction MEDIANE permet de calculer la médiane. - La fonction QUARTILE.INCLURE permet des calculer les quartiles, elle réclame en second argument le numéro du quartile à calculer 0,1,2,3, ou 4. (exemple Q1 = QUARTILE.INCLURE(matrice,1) et Q3 =QUARTILE.INCLURE(matrice,3) ) - La fonction COEFFICIENT.CORRELATION(matrice1 ;matrice2) permet de calculer le coefficient de corrélation. Tous ces indicateurs peuvent aussi être calculés simplement à partir de formules. Il n’existe en général pas dans Excel de fonction bien adaptée au traitement des données par classe. Dans ce cas, il faut repartir du cours et adapter les formules. Ces indicateurs peuvent aussi être facilement obtenus à l’aide de l’utilitaire d ?analyse statistique (Données, puis Utilitaire d’analyse, puis Statistiques descriptives), attention de l’avoir installé au préalable.

12

Partie I : Statistiques descriptives 1D

Exercice 1 "Caractères Quantitatifs discrets (Diagrammes)"

24

Dans une classe, les notes obtenues à un devoir sont : Notes : 8; 8; 8; 11; 11; 13; 13; 14; 14; 16; 16; 16; 17; 18; 18 Notes Effectif

8

11

13

14

16

17

18

1) Saisir le tableau dans un fichier "Exercice1" puis représenter le diagramme en bâtons correspondant a ce tableau Exercice 2 "Caractères Quantitatifs continus(Histogrammes)"

On a demandé la taille des élèves dans une classe de 33 élèves. On obtient les résultats suivants : Taille(en cm) Effectif

150-160 3

160-170 12

170-175 9

175-180 6

180-190 2

190-200 1

1) Saisir le tableau dans un nouveau fichier "Exercice2", puis tracer l’Histogramme correspondant à ce tableau. Exercice 3 "Caractères Qualitatifs (Diagrammes Circulaires)"

Dans une entreprise, on a demandé aux employés leur moyen de transport pour venir au travail. Les résultats sont les suivants : Moyen utilisé Effectif

à pied 50

en voiture 110

en métro 200

1) Saisir le tableau dans un nouveau fichier "Exercice3", puis tracer le Diagramme circulaire correspondant à ce tableau. Rappel : — Fréquence : on appelle fréquence le rapport entre l’effectif d’une valeur et l’effectif total. Si ni N = n1 + n2 + · · · + nt , alors fi = . N — Fréquence cumulée : on note gk la fréquence cumulée du caractère xk , c’est-à-dire gk = f1 + ni · · · + fk où fi = est la fréquence du caractère xi . La courbe des fréquences cumulées est celle N obtenue en joignant les points (xi , gi ). Exercice 4 "Courbe de Fréquence Cumulée"

25

Notes Effectif Fréquence Fréquence cumulée

1 2

2 1

3 6

4 5

5 2

6 9

7 7

8 1

9 0

10 1

1) Saisir le tableau dans un nouveau fichier "Exercice4" puis compléter. 4) Construire la courbe de fréquences cumulées.

13

Partie II : Statistique à deux dimensions

Dans la première partie, les séries statistiques étudiées étaient des séries simples : on étudiait un seul caractère dans une population. Il peut être utilise de considérer en même temps plusieurs caractères d’une même population. Par exemple, la température et la pression d’un milieu à différentes heures d’une journée.

L’échantillon est noté Ω : l’effectif de l’échantillon est n = card(Ω). Dans le cadre de la corrélation, nous  travaillons sur un échantillon de n observation, constituées de couples (xi , yi ) c-a-d Ω = (xi , yi ), i = 1...,n . n

1X xi n i=1 n 1X (xi − x)2 — Variance : V (x) = n i=1 n 1X — Moyenne : y = yi n i=1 n 1X — Variance : V (y) = (yi − y)2 n — Moyenne : x =

i=1

n 1X

— Covariance : Cov(x, y) = E [X − E(X)][Y − E(Y )] = x i − x yi − y . n



i=1

Propriétés : Cov(X, Y ) = E(XY ) − E(X)E(Y ) Cov(X, Y ) = Cov(Y, X) Cov(aX + b, cY + d) = acCov(X, Y ) X et Y indépendants −→ Cov(X, Y ) = 0, la réciproque étant fausse. Cov(x, y) — Coefficient de corrélation : r = p V (x)V (y) — — — —

Exercice 5 (Corrélation) 26

La corrélation entre deux variables X et Y mesure le lien linéaire entre deux variables. La fonction correl calcule la corrélation entre deux variables . On considère deux vecteurs X = (xi ) et Y = (yi ) avec 1 ≤ i ≤ 50 - Tracer le nuage de points xi , yi 1) 2) 3) 4)




1 ≤ i ≤ 50

Vérifier que si xi = i et yi = 2xi alors le coefficient de corrélation entre X et Y vaut 1. Vérifier que si xi = i et yi = −2xi alors le coefficient de corrélation entre X et Y vaut −1. √ Calculer le coefficient de corrélation entre X et Y si xi = i et yi = xi . Calculer le coefficient de corrélation entre X et Y si xi = i et yi = x2i .

Exercice 6 (Coefficient de Corrélation)

1) Simuler deux échantillons X = (xi ) et Y = (yi ) avec 1 ≤ i ≤ 20 indépendants et de même loi normale de moyenne 0 et de variance V = 1, N (0, 1).
2) Tracer le nuage de points xi , yi i ∈ 1, . . . , 20 . 3) Calculer le coefficient de corrélation entre les deux échantillons aléatoires X et Y . 4) Commenter.  5) Construire la série Z = (zi ), i ∈ 1, . . . , 20 définie par z1 = x1 et pour j > 1 zj = azj−1 + xj où a = 0, 85  6) Tracer le nuage de points (z
i , zi+1 ), i ∈ 1, . . . , 19 et calculer le coefficient de corrélation entre
z1 , . . . , z19 et z2 , . . . , z20 . 7) Commenter. Exercice 7

Dans un circuit électrique, on a relevé l’évolution temporelle de la tension daux bornes d’un appareil inconnu, repris dans le tableau ci dessous T (en ms) U (en Volt)

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 5.7 5.5 5.3 5.1 4.9 4.8 4.5 4.2 3.8 3.4 3.0 2.5 2 1.4 0.7

1) Représenter ce nuage de points. 2) On ajuste cette série par le modèle linéaire U = aT + b. Déterminer les coefficients de ce modèle. 3) Déterminer , au niveau de 95% l’intervalle de prévision de la tension pour un temps de 3.3ms Exercice 8 (Régression Linéaire)

27

En 1973, F. J. Anscombe a publié dans le numéro 27 de American Statistician un jeu de données très intéressantes pour montrer les pièges du calcul "aveugle" du coefficient de corrélation linéaire. Il s’agit de 11 individus sur lesquels sont mesurées 8 variables. On notera que les variables x1 , x2 , x3 sont identiques. x1 10 8 13 9 11 14 6 4 12 7 5

1 2 3 4 5 6 7 8 9 10 11

1) 2) 3) 4)

x2 10 8 13 9 11 14 6 4 12 7 5

x3 10 8 13 9 11 14 6 4 12 7 5

x4 8 8 8 8 8 8 8 19 8 8 8

y1 8.04 6.95 7.58 8.81 8.33 9.96 7.24 4.26 10.84 4.82 5.68

y2 9.14 8.14 8.74 8.77 9.26 8.10 6.13 3.10 9.13 7.26 4.74

y3 7.46 6.77 12.74 7.11 7.81 8.84 6.08 5.39 8.15 6.42 5.73

y4 6.58 5.76 7.71 8.84 8.47 7.04 5.25 12.50 5.56 7.91 6.89

Calculer la moyenne et la variance des 8 variables x1 , . . . , x4 et y1 , . . . , y4 . Calculer les coefficients de corrélation des couples (x1 , y1 ), . . . , (x4 , y4 ). Que constate t-on ? Tracer la représentation des couples (x1 , y1 ), . . . , (x4 , y4 ). Commenter Tracer la droite de régression décrite par la méthode des MOINDRES CARRÉS.

Exercice 9 (Droite de régression et points atypiques)

Douze personnes sont inscrites à une formation. Au début de la formation, ces stagiaires subissent une épreuve A notée sur 20. A la fin de la formation, elles subissent une épreuve B de niveau identique. Les résultats sont donnés dans le tableau suivant : Epreuve A Epreuve B

3 8

4 9

6 10

7 13

9 15

10 14

9 13

11 16

12 13

13 19

15 6

4 19

1) Représenter le nuage de points. Déterminer la droite de régression. Calculer le coefficient de détermination. Commenter. 2) Deux stagiaires semblent se distinguer des autres. Les supprimer et déterminer la droite de régression sur les dix points restants. 3) Calculer le coefficient de détermination. Commenter

28

TP No :6

Recherche Opérationnelle Méthode graphique & Algorithme de Simplexe

14

Introduction

Dans ce TP, on s’intéresse à la résolution du système linéaire z =< X, C > sous des contraintes AX ≤ b où A est une matrice de taille n × m, z la fonction objectif et X le vecteur inconnu. Pour ce faire, on utilisera les deux méthodes suivante : 1 Méthode graphique : en hachurant la zone correspondant aux contraintes dite domaine réalisable, et en traçant les lignes de niveaux (des lignes parallèles dont chacune correspond à une valeur constante de z) de la fonction à maximiser. On obtient ainsi la solution optimale elle correspond à un "sommet" du domaine réalisable. 2 Méthode du simplexe : c’est un algorithme fondamental pour la résolution des programmes linéaires, dont notamment les problèmes économiques. Il a été mis au point par le mathématicien américain George Dantzig en 1947. Ainsi, l’algorithme permet de maximiser ou minimiser une relation économique dite fonction objective de plusieurs variables sous certaines contraintes.

Exemple : Soit P1 le programme suivant

max [Z] = x1 + 3x2  L1 : x1 + x2 ≤ 14    L2 : −2x1 + 3x2 ≤ 12 Sous les contraintes : L : 2x1 − x2 ≤ 12    3 L4 : x1 ≥ 0, x2 ≥ 0 Avec - Z La fonction objectif. - L1 à L3 Les contraintes. - L4 Les contraintes de non négativité. Question : Trouver la solution du programme linéaire (P1 ) ? — On commencer par résoudre les contraintes L1 à L4 la zone des solutions réalisable. — Ensuite, Traçant les lignes de niveaux qui correspondant à l’équation Z = cste. — Enfin, résoudre le problème. 29

15

Méthode graphique

On trace les droites d’équations : L1 : x + y = 14, L2 : −2x + 3y = 12 et L3 : 2x − y = 12, ensuite en hachurant la zone correspondante aux contraintes, et en traçant les lignes de niveaux (des lignes parallèles) de la fonction à maximiser. On obtient ainsi la solution optimale qui correspond à un "sommet" du domaine réalisable. En utilisant la fonction drawfr.m fournie aux étudiants pendant la séance de TP pour tracer l ?ensemble des contraintes L1 à L4 .

Exemple pour le programme P1 : Le programme dans l’éditeur de Matlab est : clear c =[1 3] ; A =[1 1 ; -2 3 ; 2 -1] ; rel =’< < > >’ Le second membre des contraintes

Dans ce cas on trouve le resultat suivant : Feasible region and a level line with the objective value = 8.1667 8 7 6

— La ligne bleue est une ligne de niveau (Z = cst).

5

x2

— La zone rouge est l’ensemble des solution réalisables du programme P1 .

4 3 2 1 0 0

1

2

3

4

x1

5

6

7

8

Figure 2 – zone des solutions réalisables du P1

Donc la solution est situé dans le sommet en haut de coordonnée (x1 , x2 ) = (6, 8) ce qui donne la valeur de la fonction objectif z = 1 × 6 + 3 × 8 = 6 + 24 = 30

Exercice 1.

30

On considérer les deux programmes linéaires suivants :

Le programme linéaire P2 max [Z] = 7x + 5y  L1 : x ≤ 300      L2 : y ≤ 400 L3 : x + y ≤ 500 Les contraintes :   L : 2x + y ≤ 700    4 L5 : x ≥ 0, y ≥ 0

Le programme linéaire P3 max [Z] = x + 2y  L1 : x + 3y ≥ 12    L2 : 2x + y ≥ 14 Les contraintes : L : 4x + y ≥ 8    3 L4 : x ≥ 0, y ≥ 0

1) À l’aide de la fonction drawfr.m tracer l’ensemble des solutions réalisables et la ligne de niveau. 2) Déterminer les solutions des programmes P2 et P3 par la méthode graphique.

Exercice 2.

On considère les deux programmes linéaires P4 et P5 suivants : Le programme linéaire P5

Le programme linéaire P4 max [Z] = 9x + 3y   L1 : 5x + y ≥ 180 L2 : 5x + y ≥ 110 Les contraintes :  L5 : x ≥ 0, y ≥ 0

min [Z] = 12x − 20y    L1 : −3x + 2y ≤ 8  L2 : x + 2y ≤ 3 Les contraintes : L : y ≤ 10    3 L4 : x ≥ 0, y ≥ 0

1) Résoudre les programmes linéaires P4 et P5 à l’aide de la fonction drawfr.m sous Matlab. 2) Interpréter les résultats dans les deux cas ? Remarque : 1. Avec la méthode graphique on distingue les différents cas qui peuvent se produire lors de la résolution d’un programme linéaire. En effet, si le domaine réalisable est l’ensemble vide alors il n’y a pas de solution, le problème est non réalisable. Si la zone est infinie le problème est alors réalisable, mais non borné. 2. la résolution graphique est assez simple ici car la fonction n’est composée que de deux variables. Si nous avions 3 variables il aurait fallu faire une résolution avec un graphique en trois dimension ce qui est plus complexe pour déterminer les sommets du polyèdre. 3. Un problème de dimension 4 ou plus est impossible à résoudre graphiquement. Or notre algorithme doit être en mesure de résoudre des problèmes de dimension n ≥ 4 donc il nous faut une méthode de résolution plus adaptée.

31

16

Méthode de simplexe

Dans la plupart des problèmes réels, on a plus que deux variables à déterminer. Une procédure algébrique pour résoudre les programmes linéaires avec plus que deux variables fera l’objet de cette partie. C’est la méthode de simplexe. L’idée de la méthode du simplexe est de passer d’une solution de base admissible à une autre solution de base admissible adjacente en améliorant la fonction objective. On démarre d’une solution de base admissible quelconque (si possible l’origine). Le programme linéaire est sous la forme suivante :

  min ou max Z =< C, X > A.X ≤ b  X≥0

En introduisant les variables d’écart le problème linéaire deviendra de la forme suivante :



x1  ..  .   xp Avec X =   sp+1   ..  .

 min Z =< C, X >       a11 x1 + · · · + a1p xp + sp+1 = b1 .. .. .. . . .     an1 x1 + · · · + a1p xp + sn = bn   xi ≥ 0 ∀1 ≤ i ≤ n

     .    

sn

Rappel : Algorithme de Simplexe (exemple) Soit P1 le programme suivant max [Z] = x1 + 3x2  L1 : x1 + x2 ≤ 14    L2 : −2x1 + 3x2 ≤ 12 Sous les contraintes : L : 2x1 − x2 ≤ 12    3 L4 : x1 ≥ 0, x2 ≥ 0 1 Solution de base réalisable : Une solution de base telle que toutes les variables prennent des valeurs non-négatives est appelée solution de base réalisable. 2 Méthode du simplexe : partir d’une solution de base admissible et passer à une solution de base voisine qui améliore la valeur de l’objectif. - Solution voisine : changement d ?une variable en base. — Détermination de la variable entrante. — Détermination de la variable sortante. — Pivotage. 32

v.base x3 x4 x5 ∆j

A1 1 -2 2 -1

A2 1 3 -1 -3

A3 1 0 0 0

A4 0 1 0 0

A5 0 0 1 0

B 14 12 12 0

L1 L2 L3 L4

Le premier tableau du programme P1 : Élection de la variable entrante et sortante de la base : 1) La variable qui rentre en base : — On choisit la colonne dont sa valeur dans la ligne ∆j est le plus grand entre tous les négatifs. Dans ce cas, cela serait la variable x1 de coefficient −1. — La colonne de la variable qui rentre en base qu’on appelle colonne-pivot. Une fois obtenue la variable entrante en base, on détermine quelle sera la variable sortante. 2) La variable qui sort en base : — On prend la décision sur la base d’un simple calcul : diviser chaque terme indépendant (colonne B) entre l’élément correspondant de la colonne-pivot, à condition que les deux éléments soient strictement positifs (supérieures à zéro). On choisit la variable dont le résultat est minimal. b2 b3 12 12 1 — Dans cet exemple on a : xb11 = 14 1 = 14, x21 = −2 < 0, Non (on prend pas) et x31 = 2 = 6. =⇒ le minimum dans la ligne 3. On choisit la variable sortante x5 . 3) Condition d’arrêt : — Si l’objectif est la maximisation, quand dans la dernière ligne (ligne indicatrice) n’existe aucune valeur negative parmi les coûts réduits (colonnes B en avance), on obtient la condition d’arrêt. Dans ce cas, on attend la fin de l’algorithme puisqu’on ne peut pas l’améliorer. Deuxième tableau du simplexe Pivotage de Gauss 2 est le pivot (intersection de la colonne qui rentre et de la ligne qui sort)

L01 = L1 − 1L03 L03 = L23

L02 = L2 − (−2)L03 L04 = L4 + 1L03

De même pour les autres tableaux — Dans cet exemple on a ∆02 = −3.5 =⇒ la colonne de pivot A2 y x2 variable qui rentre. b3 8 6 1 2 — Ensuite on a : xb12 = 1.5 = 24 , xb22 2 = 12, et x31 = −0.5 < 0 on prend pas =⇒ le minimum dans 33

v.base x3 x4 x1 ∆j

A1 0 0 1 0

A2 1.5 2 -0.5 -3.5

A3 1 0 0 0

A4 0 1 0 0

la ligne 1. y x3 variable qui sort. — La colonne de pivot est 2 et la ligne de pivot 1

34

A5 -0.5 1 0.5 0.5

B 8 24 6 6

L1 L2 L3 L4

Troisième tableau v.base x2 x4 x1 ∆j

A1 0 0 1 0

A2 1 0 0 0

A3 2/3 -4/3 1/3 7/3

A4 0 1 0 0

A5 -1/3 5/3 1/3 -2/3

B 16/3 40/3 26/3 74/3

L1 L2 L3 L4

A3 2/5 -4/5 3/5 9/5

A4 1/5 3/5 -1/5 2/5

A5 0 1 0 0

B 8 8 6 30

L1 L2 L3 L4

La colonne de pivot est 5 et la ligne de pivot 2 Dernier tableau v.base x2 x5 x1 ∆00j

A1 0 0 1 0

A2 1 0 0 0

— Tous les ∆00j ≥ 0, ∀j = 1, ...5 =⇒ Condition d’arrêt. — La solution optimale se lit dans le tableau : les variables qui Pne sont pas dans le tableau sont hors-base donc à zéro. La solution optimale est donc : Z = 5i=1 ci xi la valeur correspondante de B soit Z = 6 × 1 + 3 × 8 + 8 × 0 = 6 + 24 = 30 atteinte en S = (6, 8, 0, 0, 8). Rappel : Algorithme de Simplexe (cas général) On choisit parmi les variables hors base, la variable xj , telle que : ∆j = max(∆i < 0; i = 1; . . . ; n) avec ∆i = −ci On détermine parmi les variables de base la variable sortante xi telle que : abiji = min{ abiji ; i = 1; . . . ; n} avec bi > 0 et aij > 0 On élimine xj par pivotage (élimination de Gauss) : 1. Diviser la contrainte i par le pivot xij 2. Modifier la contrainte k, pour tout k 6= i, en multipliant la contrainte i par akj et en la soustrayant de la contrainte k. 3. Modifier la fonction économique en multipliant la contrainte i par cj et en la soustrayant de l’expression de z. La solution optimale est atteinte lorsque tous les ∆i sont positifs. On utilisant la fonction simplex2p.m fournie aux étudiants pendant la séance de TP pour résoudre le programme linéaire suivant : Exemple pour le programme P1 : Le programme dans l’éditeur de Matlab est :

35

clear c =[1 3] ; A =[1 1 ; -2 3 ; 2 -1] ; rel =’< < > >’ Le second membre des contraintes Maximisation (max) ou minimisation (min)

Dans ce cas on trouver le resultat avec tous les tableaux de simplexe

Exercice 3.

On considére les deux programmes linéaires P4 et P5 suivantes. Le programme linéaire P5

Le programme linéaire P4 max [Z] = 9x + 3y   L1 : 5x + y ≥ 180 L2 : 5x + y ≥ 110 Les contraintes :  L5 : x ≥ 0, y ≥ 0

min [Z] = 12x − 20y    L1 : −3x + 2y ≤ 8  L2 : x + 2y ≤ 3 Les contraintes : L : y ≤ 10    3 L4 : x ≥ 0, y ≥ 0

1) Résoudre les programmes linéaires P4 et P5 à l’aide de la fonction simplex2p.m sous Matlab. 2) Comparer les résultat avec la méthode graphique dans les deux cas ?

Exercice 4.

On considère les deux programmes linéaires P6 et P7 suivantes. Le programme linéaire P6

Le programme linéaire P7

max [Z] = x + 3y + z min [Z] = 14x + 12y + 12z   L : 5x + y ≤ 100  1    L1 : x − 2y + 2z ≥ 1 L2 : y + 3z ≤ 300 L2 : x + 3y − z ≥ 3 Les contraintes : Les contraintes : L : 2x − y + z ≤ 900    L3 : x ≥ 0, y ≥ 0, z ≥ 0  3 L4 : x ≥ 0, y ≥ 0, z ≥ 0 1) Résoudre les programmes linéaires P6 et P7 à l’aide de la fonction simplex2p.m sous Matlab. 2) Comparer les résultat dans les deux cas ?

36

17

Fonction "linprog"

Les problèmes de cette partie sont à traiter à l’aide de la fonction linprog définir sous MATLAB. La commande help linprog donne accès à l’aide sur cette a fonction. Il est important de remarquer que linprog résout un problème de minimisation. Donc, le maximum d’une fonction z sera cherché comme le minimum de la fonction −z.

Exercice 5.

On utilisant la commande linprog résoudre les programmes linéaires P1 à P7 . Voici l’exemple du programme P1 : clear all clc c =[-1 -3] ; ←− A =[1 1 ; -2 3 ; 2 -1] ; ←− b=[14 12 12] ; ←− [X, Zmax]=linprog(c,A,b,[ ],[

Les paramètres de la fonction objectif Z. Matrice des contraintes. Le second membre des contraintes ],[0 0 0 0 0],[ ],[ ],1)

Exercice 6.

Une usine fabrique deux produits P 1 et P 2 à l’aide de trois matières premières M 1, M 2 et M 3 dont on dispose en quantité limitée. On se pose le problème de l’utilisation optimale de ce stock de matières premières c’est-à-dire la détermination d’un schéma, d’un programme de fabrication tel que : — les contraintes de ressources en matières premières soient respectées, — le bénéfice réalisé par la vente de la production soit maximum. Modèle mathématique : — Données numériques des contraintes. La disponibilité en matières premières est de 18 unités de M 1, 8 unités de M 2 et 14 unités de M 3. — Caractéristiques de fabrication. Elles sont données dans le tableau ci-dessous :

P1 P2

M1 1 3

M2 1 1

M3 2 1

1) Écrire le programme linéaire correspondant au problème 2) Résoudre le programme par les différents méthodes

37