Poly TP Maths 2023 2024

Année 2023-2024
MATHS–EI1
Françoise Foucher
Mazen Saad
Marie Billaud-Friess
Mathématiques
Énoncés TP
Table des matières
1 Initiation Matlab – Equation de la chaleur 1D 5
2 Résolution d’un système linéaire tridiagonal par plusieurs méthodes (LU, Ja-
cobi, Gauss-Seidel) 19
3 Optimisation sans contrainte 25
4 Optimisation avec contraintes 31
5 Lois usuelles – Convergence de suites de variables aléatoires 39
6 Evaluation d’une proportion – Moindres carrés avec contraintes linéaires 45
3
Mathématiques – Énoncés TP
4/49
Chapitre 1
Initiation Matlab – Equation de

la chaleur 1D
Objectif
— Initiation au calcul scientifique avec Matlab. Matlab est un logiciel de calcul numérique
produit par MathWorks (https :// www.mathworks.com/) pour le traitement des matrices
et pour le calcul numérique. Matlab est l’abréviation en anglais de "matrix laboratory",
tous les objets en jeu sont des matrices, y compris les scalaires (matrices 1×1). Tout ce qui
est fait dans ce TP peut être réalisé avec le logiciel libre de calcul numérique
OCTAVE (https ://www.gnu.org/software/octave/), du projet informatique
GNU.
— Ecriture et test d’un programme qui résout, comme vu en TD par la méthode des dif-
férences finies, l’équation de la chaleur en une dimension −u00 = f sur ]0, 1[ avec les
conditions au bord u(0) = ug et u(1) = ud .
Travail à rendre
— La fiche de résultats complétée sur papier au fur et à mesure des exercices.
— Le code des fonctions et des scripts, dans une archive TP1_GROUPE_Nom1_Nom2.zip.
— Déposer l’archive sur hippocampus, cours MATHS_S5 ou MATHS_S6 suivant le semestre.
Il faudra auparavant s’inscrire dans votre groupe via l’activité "Inscription dans son groupe
pour déposer les TP".
1.1 Premières commandes, variable ans, point-virgule ;

Une fois Matlab lancé, se placer dans son dossier de travail Z:\MATLAB, puis utiliser la
fenêtre de commandes pour taper les instructions. Chaque ligne de commande commence par
le sigle >>. Par exemple :
>>1+1
ans =
5
2
>>t=1+1
t =
2
>>t=1+1;
>>t
t=
2
>>u=sin(t)
u =
0.9093
>>v=exp(u)
v =
2.4826
>>format long
>>v
v =
2.48257772801500
>>format short
>>v
v =
2.4826
>>who
Your variables are:
ans t u v
>>whos
Name Size Bytes Class
ans 1x1 8 double array
t 1x1 8 double array
u 1x1 8 double array
v 1x1 8 double array
>> clear
>> help format
format Set output format.
format with no inputs sets the output format to the default appropriate
for the class of the variable. For float variables, the default is
format SHORT.
...
On remarque que :
1. par défaut, tout calcul est affecté dans la variable ans,
2. la commande format permet de modifier l’affichage du format des différentes variables :
short pour 5 chiffres affichés, long pour 16 chiffres affichés,
3. les commandes who, whos permettent de lister l’ensemble des variables utilisées,
4. la commande clear efface le contenu de tous les variables utilisées,
5. on peut rappeler la commande précédente à l’aide de la flèche ↑ du pavé numérique,
6. on obtient de l’aide avec la commande help ou doc,
7. avec l’instruction t=1+1;, le calcul de t a été effectué mais il n’est pas affiché à l’écran.
C’est l’instruction t de la ligne suivante qui permet l’affichage.
6/49
1.2 Variables prédéfinies et codage double précision des

réels
Certaines variables sont déjà définies, par exemple pi, i :
>> pi
ans =
3.1416
>> i
ans =
0 + 1.0000i
et aussi des variables relatives à la représentation des réels par Matlab : realmin, realmax, eps,
inf, NaN, . . .Matlab utilise un codage double précision des réels : codage binaire sur 64 bits (1
pour le signe, 11 pour l’exposant, 52 pour la mantisse), ce qui donne la représentation d’un réel
x:
e 1 1 1
x = (±1) 2 a0 + a1 + a2 2 + · · · + a52 52
2 2 2
avec a0 = 1, a1 ,. . .,a52 égaux à 0 ou 1, et l’exposant e ∈ [−210 , 210 − 1] = [−1024, 1023].
Le plus petit réel (en valeur absolue) codable est noté realmin :
realmin ' 2−1024 ' 10−308
>> realmin
ans =
2.2251e-308
Pour les grandes valeurs réelles, on a le plus grand réel codable noté realmax :
realmax ' 21023 ' 10308
>> realmax
ans =
1.7977e+308
au delà duquel Matlab renvoie la variable notée Inf, abréviation de Infinity. On parle aussi dans
ce cas de overflow :
>> v=1.e+400
v =
Inf
Quand le nombre ne peut pas être évalué, Matlab renvoie la variable NaN qui signifie "Not a
Number" :
>> v=O/0
v =
NaN
7/49
Il ne faut pas confondre le plus petit réel codable realmin avec l’erreur relative entre deux réels
qui ont la même représentation, notée eps :
eps = 2−52 ' 10−16
>> eps
ans =
2.2204e-16
On l’appelle aussi "zéro machine". C’est le plus grand réel qui vérifie 1+eps=1 au sens où 1+eps
et 1 sont codés identiquement :
>> 1. + eps
1.0000
Comme eps ' 10−16 , on peut considérer que les réels sont ainsi codés avec 16 chiffres signifi-
catifs.
Conséquence : dans une addition x+y, si le rapport |y/x| est inférieur à eps alors l’addition ne
retiendra que la valeur de x.
Exercice 1
Définir les variables x = 1020 et y = 103 et faire calculer les expressions suivantes,
qu’on comparera aux valeurs exactes. Compléter la fiche de résultats.
>> x = 1e20; y = 1e3;

>> y/x
>> x+y
>> x+y-x
>> x-x+y
>> x*(1+y/x)-x
>> y/(x+y-x)
1.3 Nombres complexes

>> c1 = 1-2i
c1 =
1.0000 - 2.0000i
>> c2 = 3*(2-sqrt(-1)*3)
c2 =
6.0000 - 9.0000i
>> c3 = conj(c2);% complexe conjugué
>> a1 = real(c1);% partie réelle
>> b2 = imag(c2);% partie imaginaire
>> r = abs(c2);% module
>> t = angle(c1);% argument
1.4 Les commentaires, les affichages de message

Il est important de pouvoir écrire des commentaires lors de l’élaboration d’un programme. Pour
cela, sur une ligne, tout ce qui est après le symbole % n’est pas pris en compte.
8/49
>> pi % pour connaitre la valeur de pi
De même, il est important d’afficher des informations au cours de l’exécution d’un programme.
Pour cela, on peut utiliser la fonction disp :
>> disp(’la valeur de pi avec 15 décimales vaut : ’); disp(pi);

la valeur de pi avec 15 décimales vaut :
3.141592653589793
ou la fonction fprintf, avec %f pour un format décimal et %e pour un format scientifique :
>> fprintf(’la valeur de pi avec 10 décimales vaut pi = %15.10f\n’,pi);

la valeur de pi avec 10 décimales vaut pi = 3.1415926536
>> fprintf(’la valeur de pi avec 10 décimales vaut pi = %15.10e\n’,pi);
la valeur de pi avec 10 décimales vaut pi = 3.1415926536e+00
1.5 Vecteurs - Matrices

Les tableaux (vecteurs, matrices, ...) s’utilisent dans la fenêtre de commande Matlab ou dans des
programmes (fonctions ou scripts).
Les indices des vecteurs et matrices commencent toujours à 1. L’indice 0 n’existe
pas en Matlab.
La ième composante d’un vecteur x est notée x(i) ; le coefficient de la ième ligne, jème colonne
d’une matrice A est noté A(i,j). Il n’y a pas à déclarer les tableaux, mais il est cependant
recommandé de réserver une place mémoire pour les vecteurs et matrices en les initialisant à 0
(voir plus loin).
1.5.1 Création d’un vecteur ligne
Un vecteur ligne peut être défini de plusieurs façons :

(i) soit par une instruction de la forme
x = [x1 x2 ... xn] ou x = [x1, x2, ..., xn]
où les éléments xi sont séparés par des espaces ou des virgules et mis entre [ ].
>> x1 = [1.1 2.3 3.5 -4. -8.] % (vecteur de 5 composantes)

x1 =
1.1000 2.3000 3.5000 -4.0000 -8.0000
>> x2 = [3 -5.2 2*sqrt(3)] % (vecteurs de 3 composantes).
x2 =
3.0000 -5.2000 3.4641
>>x2(2)
ans =
-5.2000
(ii) soit par une subdivision uniforme d’un intervalle donné [a,b] en sous-intervalles de
longueur donnée h :
>> x = a:h:b; % vecteur ligne, x(i)=a+(i-1)h
>> x = a:b; % par défaut h = 1
9/49
>> x = 5 :-1 : 1
x = 5 4 3 2 1
>> x = 1: 1.1 : 5
x = 1.0000 2.1000 3.2000 4.3000
>> x = 5:1
x =
Empty matrix: 1-by-0
>> x = 0 : pi/2 : 2 * pi
x =
0 1.5708 3.1416 4.7124 6.2832
(iii) soit par une subdivision uniforme d’un intervalle donné [a,b] comprenant un nombre
de points donné n :
>> x = linspace(a,b,n); % vecteur ligne, x(i)=a+(i-1)h, h=(b-a)/(n-1), i = 1,...,n
>> x = linspace(a,b); % par défaut n = 100
>> y = linspace(0,pi,9)
y =
0 0.3927 0.7854 1.1781 1.5708 1.9635 2.3562 2.7489 3.1416
(iv) soit par l’utilisation d’une fonction appliquée à un vecteur ligne, par exemple, sin(x)
ou cos(x) avec x un vecteur ligne.
>> x = (0 : 0.2 : 1); % x est un vecteur ligne
>> y = exp(x); % y est un vecteur ligne, y(i)=exp(x(i))
(v) soit par initialisation à un vecteur ligne avec des 0 ou des 1, de longueur donnée n :
>> x = zeros(1,n); % vecteur ligne, x(i)=0, i = 1,...,n
>> x = ones(1,n); % vecteur ligne, x(i)=1, i = 1,...,n
1.5.2 Vecteur transposé, création d’un vecteur colonne
Un vecteur colonne peut être défini de plusieurs façons :
(i) comme le transposé d’un vecteur ligne : x’ définit dans Matlab le vecteur transposé du
vecteur x.
>> y = (1 : 4)’; % vecteur colonne de composantes 1,2,3,4
>> y = [3 -4.5 2.1]’; % vecteur colonne de composantes 3,-4.5,2.1
>> y = linspace(0,1,5)’; % vecteur colonne de composantes 0,0.25,0.5,0.75,1
(ii) ou directement, en utilisant des points virgules ; au lieu de virgules ou espaces :
x = [x1; x2; ...; xn]
>> z = [3.1452 ; -3 ; 4. ; 5.256];
(iii) ou par l’utilisation d’une fonction appliquée à un vecteur colonne x, par exemple, sin(x)
ou cos(x).
>> x = (0 : 0.2 : 1)’; % x est un vecteur colonne
>> y = exp(x); % y est un vecteur colonne, y(i)=exp(x(i))
(iv) ou par initialisation à un vecteur colonne avec des 0 ou des 1, de longueur donnée n :
>> x = zeros(n,1); % vecteur colonne, x(i)=0, i = 1,...,n
>> x = ones(n,1); % vecteur colonne, x(i)=1, i = 1,...,n
1.5.3 Opérations sur les vecteurs
Quelques opérations sur les vecteurs :
10/49
u+v : addition de deux vecteurs u et v de même longueur,

dot(u, v) : produit scalaire de deux vecteurs u et v,
u’*v : produit scalaire de deux vecteurs colonnes u et v,
u.*v : multiplie deux vecteurs u et v composantes par composantes,
u./u : divise deux à deux les composantes de deux vecteurs u et v,
u.ˆ v : élève les composantes de u à la puissance des composantes de v,
abs(u) : vecteur donnant les valeurs absolue des composantes d’un vecteur u,
sum(u) : somme des composantes d’un vecteur u,
mean(u) : moyenne des composantes d’un vecteur u,
length(u) : longueur d’un vecteur u,
min(u) : plus petite composante d’un vecteur u,
max(u) : plus grande composante d’un vecteur u,
norm(u) : norme euclidienne d’un vecteur u,
u(n1:n2) : extrait un sous-vecteur de composantes u(n1) à u(n2).
>> u = (1 : 4) , v = (2 : 5)
u =
1 2 3 4
v =
2 3 4 5
>> z = u .* v % z(i) = u(i)*v(i)
z =
2 6 12 20
>> w = v./u % w(i) = v(i)/u(i)
w =
2.0000 1.5000 1.3333 1.2500
>> t = v .^u % t(i) = v(i) ^ u(i)
t =
2 9 64 625
>> x = (0 : 0.2 : 1)
x =
0 0.2000 0.4000 0.6000 0.8000 1.0000
>> z = exp(x).* cos(x)
z =
1.0000 1.1971 1.3741 1.5039 1.5505 1.4687
\%z est un vecteur de composantes z(i) = exp(x(i)) cos(x(i)) pour i = 1,...,6
>> y = z(3:5)
y =
1.3741 1.5039 1.5505
1.5.4 Création de matrices

Une matrice de n lignes et p colonnes peut être définie :
(i) par une relation de la forme :
A = [a11 a12 a1p ; a21 a22 a2p ; · · · an1 an2 anp ]

>> A = [1 2 3 4 ; 5 6 7 8 ; 9 10 11 12] % A a 3 lignes et 4 colonnes
A =
1 2 3 4
5 6 7 8
9 10 11 12
11/49
(ii) ou à partir de vecteurs ou de matrices plus petites par concaténation :

>> A = [1 2 3; 4 5 6]
A =
1 2 3
4 5 6
>> B = [A; 7 8 9]
B =
1 2 3
4 5 6
7 8 9
(iii) ou par initialisation à des matrices particulières :
>> A = zeros(n,p); % matrice n lignes et p colonnes, A(i,j) = 0
>> A = eye(n); % matrice identité d’ordre n
>> A = ones(n,p); % matrice n lignes, p colonnes, A(i,j) = 1
>> A = rand(n,p) % matrice de nombres aléatoires, 0 < A(i,j) < 1
>> A = diag(u); % matrice diagonale, de diagonale égale au vecteur u
>> A = diag(u,1); % matrice de sur-diagonale égale à u, nulle ailleurs
1.5.5 Opérations sur les matrices

Quelques opérations possibles
A’ : transposée de A si A est réelle, adjointe si A complexe,
+ , * : addition, multiplication de deux matrices compatibles,
x = A(1, :) : première ligne de A,
y = A(:, 2 : 3) : deuxième et troisième colonnes de A,
w = A(3, i : j) : du i-ième au j-ième elements de la ligne 3,
u = A(:,2) : deuxième colonne de A,
z = A(:) : mise sous forme d’une colonne,
c*A : multiplie tous les éléments de A par le scalaire c,
A.ˆ m : élévation à la puissance m de chaque élément de la matrice A,
Aˆ m : élévation à la puissance m de la matrice A,
[n,p] = size(A) : donne le nombre n de lignes et le nombre p de colonnes de la matrice A,
n = size(A,1) : nombre de lignes de la matrice A,
p = size(A,2) : nombre de colonnes de la matrice A,
eig(A) : vecteur donnant les valeurs propres de la matrice A,
det(A) : déterminant de la matrice A,
rank(A) : rang de la matrice A,
trace(A) : trace de la matrice A,
cond(A) : conditionnement de la matrice A,
spy(A) : représentation graphique de la matrice A,
inv(A) : inverse de la matrice A si elle existe,
x = A\b : résolution du système linéaire Ax = b,
x = linsolve(A,b) : résolution du système linéaire Ax = b (équivalent à x = A\b).
Exemples avec les fonctions spy, size et diag :
>>B = rand(20,30);
>>B(3:6,10:18) = 0.;
>>spy(B)
>>[n,m] = size(B)
12/49
>>NbLigneB = size(B,1)
>>NbColonneB = size(B,2)
>>A = diag([1 2 3],-1)+diag([4 5 6],1)
Exercice 2
A l’aide des fonctions eye, diag et ones, créer une matrice tridiagonale d’ordre 10
ci-dessous. Reporter la formule utilisée sur la fiche de résultats.
 
2 −1
−1 2 −1 0 
 

 −1 2 −1 

A=  .. .. .. 
 . . . 

 . . . .

 0 . . −1
−1 2
1.6 M-Files ou scripts

Un script (ou M-file) est un fichier contenant des instructions Matlab et dont le nom se termine
par le suffixe .m. Matlab contient un éditeur pour écrire et mettre au point les scripts. Pour
exécuter le script il suffit de cliquer sur run dans le menu de l’éditeur, ou bien, dans la fenêtre
de commande, d’entrer le nom du script :
>> premier
par exemple pour le script premier de l’exercice suivant. Les scripts sont exécutés séquentiellement
dans le "workspace", et ils peuvent accéder aux variables qui s’y trouvent déjà, les modifier, en
créer d’autres, etc.
Pour calculer le temps d’exécution d’une partie d’un script, mettre l’instruction tic au début et
l’instruction toc à la fin.
Exercice 3
Télécharger le fichier premier.m dans le dossier de travail (dossier dans lequel on
exécute les commandes Matlab : Z:\MATLAB\TP1 recommandé). Attention : ce script
contient des erreurs. Le déboguer. S’aider du carré en haut à droite, rouge quand il
reste au moins une erreur, orange quand il reste des alertes (par exemple quand il
manque ; car on souhaite un affichage de contrôle), vert quand la syntaxe est correcte.
Il peut encore rester des erreurs à l’exécution. S’aider aussi des variables, visibles dans
la fenêtre de droite, ou prévoir des affichages supplémentaires dans le code pour suivre
ce que fait le programme. Comprendre ce que fait le programme et compléter la fiche
de résultats.
Exercice 4
Télécharger le fichier TP1.m, dans le dossier de travail. Il contient un script qui va
permettre de calculer une solution approchée de l’équation de la chaleur en dimension
1 par la méthode des différences finies. Commencer à compléter ce script où cela est
indiqué pour l’exercice 4 : calcul du pas h, des vecteurs X et Xint et de la matrice A.
Compléter la fiche de résultats avec les valeurs obtenues pour N s = 3. Vérifier qu’on
obtient les résultats attendus.
13/49
1.7 Fonctions
Une fonction est un script admettant des variables d’entrées et des variables de sorties, par
exemple f (x, y) = x2 + y 2 admet x et y comme variables d’entrées et le résultat z = x2 + y 2
comme argument de sortie. Le nom de la fonction sera le nom du fichier. Voici une fonction "fonc"
définie dans un fichier fonc.m :
function [z,t] = fonc(x,y,m)

z = x^2+y^2;
t=x-y+m;
end
Les variables x, y et m sont les variables d’entrées, elles sont utilisées par la fonction sans être
modifiées. Les variables z et t sont les arguments de sortie, la fonction doit affecter une valeur
pour z et t.
Utilisation de cette fonction : Une fonction est appelée depuis un script (un programme
principal ou une autre fonction) ou dans la fenêtre de commande.
>> [z,t] = fonc(1., 0., -4.)

>> [y1,y2]= fonc(-1., sin(1.), sqrt(2.) )
Définition locale d’une fonction simple : On peut définir une fonction localement à l’in-
térieur du script qui l’utilise, à condition qu’elle ne renvoie qu’un seul résultat et que son ex-
pression soit simple (sans boucles itératives ou de choix). Par exemple, pour définir la fonction
f (x, y) = x2 + y 2 :
>> f = @(x,y) x^2+y^2;

>> z = f(2,1)
z =
5
Cela est utile en particulier pour redéfinir localement une fonction d’une seule variable à partir
d’une fonction plus complexe définie dans un fichier .m, par exemple pour tracer son graphe, voir
plus loin.
1.8 Boucles et choix

1.8.1 Opérateurs logiques de comparaison
< plus petit
> plus grand
<= plus petit ou égal
>= plus grand ou égal
== égal ( compare si deux nombres sont égaux ou non)
∼= différent
&& et logique
|| ou logique
∼ non (pour le contraire)
1 vrai
0 faux
14/49
1.8.2 Choix simple ou alternatif (if)

Faire >> help if pour une description complète. Exemple :
a= log(2.); b= sqrt(2.)/2.;
if (a > b)
disp(’ a est plus grand que b’)
else
disp(’a est plus petit que b’)
end
1.8.3 Choix multiple selon un paramètre (switch)

Faire >> help switch pour une description complète. Exemple :
switch m % faire selon la valeur de m

case 0 % si m=0
x=0
case {1, 2, 9} % si m=1 ou m=2 ou m=9
x= m^2
case {-3,-1, 99} % si m=-3 ou m=-1 ou m= 99
x= m+ m^2
otherwise % sinon
x=-9999999
end
Exercice 5
Ecrire dans un fichier f.m (qu’on crée avec l’éditeur de Matlab) la fonction f(x,m)
définie par : 
0.

 si m = 0
−2. si m = 1



f (x, m) = x2 si m = 2

4π 2 sin(2π x) si m = 3




10. e−x

si m = 4
de façon à ce que, si x est un vecteur de coordonnées xi , alors le résultat soit le vecteur
des valeurs f (xi , m). S’aider de la trame ci-dessous pour écrire la fonction. Tester la
fonction dans la fenêtre de commande, pour les différentes valeurs de m, avec l’instruc-
tion f([0;0.5;1],m) qui doit renvoyer le vecteur colonne des valeurs de f en 0, 0.5 et
1. Compléter la fiche de résultats.
function y = f(x,m)
% Terme source
switch m
case 0
y = zeros(size(x)); %vecteur de zéros, de même dimension que x
case 1
y = ....
case 2
y = .....
case 3
15/49
y = .....
case 4
y = .....
otherwise
error(’Argument m indéfini dans f’);
end
end
Exercice 6
Ecrire dans un fichier solex.m (qu’on crée avec l’éditeur de Matlab) la fonction
solex(x,ug,ud,m) qui donne la solution exacte de l’équation −u00 = f associée à
la fonction f(x,m) et aux valeurs ug et ud. S’aider de la trame ci-dessous, qui donne la
solution pour le cas m = 0 (f = 0). Indication : on calcule ces solutions "à la main"
en intégrant deux fois la fonction f . Ce calcul peut être fait plus tard, il suffira alors de
compléter cette fonction. Exécuter dans la fenêtre de commande, pour les différentes
valeurs de m, l’instruction solex([0;0.5;1],1,2,m) qui doit renvoyer le vecteur co-
lonne des solutions exactes en 0, 0.5 et 1 dans le cas ug=1, ud=2. Compléter la fiche de
résultats.
function y = solex(x, ug, ud, m)

% Solutions exactes
switch m
case 0
y = (ud - ug) * x + ug;
case 1
y = ....
case 2
y = .....
case 3
y = .....
case 4
y = .....
otherwise
error(’Argument m indéfini dans f’);
end
end
Exercice 7
Compléter le script TP1.m avec le calcul du second membre et le calcul de la solution
exacte à l’aide des fonctions f et solex. Rendre actives les instructions de calcul des
erreurs. Vérifier que tout fonctionne.
1.8.4 Boucle itérative pour (for)

Faire >> help for pour une description complète. Exemple :
A = zeros(4,5)
for i = 1: 4
for j=5:-1:1
A(i,j) = i+j^2 ;
end
16/49
end
A
1.8.5 Boucle itérative tant que (while)

Faire >> help while pour une description complète. Exemple :
% recherche l’indice du 1er coefficient non nul (numériquement) d’un vecteur u

n = length(u);
test = 0;
i = 0;
while (~ test) && (i < n)
i = i + 1;
test = (abs(u(i)) > eps);
end
i
1.9 Graphismes
Faire >> help graph2d pour l’aide sur les graphes en 2D.
Dans le script courbes.m ci-dessous sont données différentes instructions qui pourront être testées
et utilisées plus tard.
% courbes.m
% Exemples de graphismes en 2D
% graphe 1
figure
x=0:0.05:5;
y=sin(x.^2);
plot(x,y);
xlabel(’Time’)
ylabel(’Amplitude’)
title(’ma première courbe’)
legend(’toto’)
%graphe 2
figure
z = -2.9:0.2:2.9;
bar(z,exp(-z.*z));
%graphe 3
figure
subplot(2,1, 1)
plot(x,y);
subplot(2,1,2)
bar(z,exp(-z.*z));
% graphe 4
figure
Y=[0:0.05: 1];Z1=sin(2*pi*Y);Z2=cos(2*pi*Y);
17/49
plot(Y,Z1,’:b’,Y,Z2,’+k’);
title(’Exemple de courbes’);
xlabel(’Y’);ylabel(’Z’);
legend(’sin’,’cos’);
Exercice 8
Finaliser le script TP1.m en rendant actives les instructions qui permettent de tracer
les graphes après le calcul des erreurs. Vérifier que tout fonctionne.
1.10 Quelques fonctions mathématiques

Quelques fonctions usuelles :
— sin, cos, tan, sinh, cosh, tanh,...
— asin, acos, atan, asinh, acosh, atanh, ...
— exp, log, log10, sqrt, ...
et aussi :
fix(x) : plus proche entier inférieur ou égal à x en valeur absolue
floor(x) : partie entière de x
ceil(x) : plus proche entier plus grand ou égal à x
round(x) : entier le plus proche de x
mod(x,y) : le reste de la division euclidienne de x par y
sign : signe de x (vaut -1, 0 ou 1)
Exercice 9
1. A l’aide du script TP1, vérifier que la résolution fonctionne bien : choisir m = 0 ou
1 (car l’erreur de la méthode est nulle pour ces exemples où f est un polynôme de
degré inférieur ou égal à 1), h = 0.1 (c’est-à-dire N s = 9). Les graphes des solu-
tions exactes et approchées doivent se superposer, et l’erreur doit être de l’ordre
du zéro machine. Sinon, vérifier chaque étape de calcul. Reporter les valeurs des
erreurs sur la feuille de résultats.
2. Etudier la convergence de la méthode : pour chacun des exemples m = 2, 3 et
4, vérifier sur le graphe pour h = 0.1 que les solutions sont proches, et reporter
sur la feuille de résultats les erreurs obtenues avec h = 0.1, 0.01, 0.001. Vérifier
la convergence numérique d’ordre 2, c’est-à-dire que l’erreur diminue à peu près
proportionnellement au carré de h.
18/49
Chapitre 2
Résolution d’un système linéaire

tridiagonal par plusieurs
méthodes (LU, Jacobi,
Gauss-Seidel)
Objectif
— Programmer, tester et comparer plusieurs méthodes pour résoudre un système linéaire
tridiagonal : la méthode directe LU et les méthodes itératives de Jacobi et Gauss-Seidel.
On comparera aussi aux résultats donnés par la fonction linsolve de Matlab, qui utilise
une factorisation LU (voir doc linsolve).
— Savoir programmer une méthode itérative.
— Observer que la précision des résultats peut varier suivant la matrice du système linéaire
et suivant la méthode de résolution utilisée. Connaitre les paramètres qui donnent des
informations sur la qualité du résultat.
— Comprendre les résultats par une étude théorique.
Tout ce qui est fait dans ce TP peut être réalisé avec le logiciel libre de calcul
numérique OCTAVE (https ://www.gnu.org/software/octave/).
Travail à rendre
— La fiche de résultats complétée sur papier.
— Le code des fonctions et du script dans une archive TP2_GROUPE_Nom1_Nom2.zip
— Déposer l’archive sur hippocampus, cours MATHS_S5 ou MATHS_S6 suivant le semestre.
Il faudra auparavant s’inscrire dans votre groupe via l’activité "Inscription dans son groupe
pour déposer les TP" si cela n’a pas été déjà fait.
19
2.1 Fonctions de résolution d’un système linéaire tridiago-

nal
1. Ecrire les fonctions suivantes de résolution d’un système linéaire tridiagonal par les mé-
thodes LU, Jacobi et Gauss-Seidel, suivant les algorithmes étudiés en TD :
(a) function x = lu_tridiag(n,A,b) :

Cette fonction calcule la sous-diagonale l de la matrice L, la diagonale d et la sur-
diagonale u de la matrice U (cf. TD) puis calcule la solution x en résolvant succes-
sivement les systèmes Ly = b puis U x = y (cf. TD). Attention : ne pas calculer les
matrices L et U entièrement (pour ne pas encombrer inutilement la mémoire avec des
zéros), et ne pas utiliser la fonction linsolve de Matlab.
(b) function [x,k,converge] = jacobi_tridiag(n,A,b,itermax,tolerance) :

Suivre la trame de la fonction donnée ci-dessous. Appliquer le schéma écrit en TD
pour le calcul de xsuiv en fonction de x. Attention : ne pas calculer les matrices D,
E, F (cela encombrerait inutilement la mémoire et risquerait d’augmenter le nombre
d’opérations).
(c) function [x,k,converge] = gauss_seidel_tridiag(n,A,b,itermax,tolerance) :

Reprendre la fonction jacobi_tridiag dans laquelle il suffit de modifier le calcul de
xsuiv en fonction de x, en suivant le schéma écrit en TD. Attention : de même que
pour la méthode de Jacobi, ne pas calculer les matrices D, E, F .
Indication : trame de la fonction jacobi_tridiag :

%__________________________________________________________________________
function [x,k,converge] = jacobi_tridiag(n,A,b,itermax,tolerance)
%--------------------------------------------------------------------------
% resout le systeme tridiagonal Ax=b par la methode de Jacobi
% itermax = nombre d’iterations maximal
% tolerance = ecart maximal entre les deux derniers iteres
%
% k = nombre d’iterations effectuees
% la methode converge si ecart entre deux derniers iteres < tolerance
% converge = 1 si la methode converge, converge = 0 sinon
% ------------------------------------------------------------------------
x = zeros(n,1); xsuiv = x;
k = 0;
converge = 0 ;
while ((~converge)&&(k<itermax))
% calcul de xsuiv en fonction de x suivant le schema de Jacobi
.....
ecart = norm (xsuiv-x);
converge = (ecart < tolerance);
k = k+1;
x = xsuiv;
end
20/49
2. Vérifier le bon fonctionnement de ces fonctions par exemple pour le système Ax = b avec :
     
2 −1 0 0 1 1
−1 2 −1 0  0 1
A=  0 −1 2 −1 et b = 0 , de solution x = 1 .
    
0 0 −1 2 1 1
2.2 Tests avec le script TP2.m

Les tests seront faits sur deux systèmes Ax = b, de dimension N qu’on fera varier. Le second
membre b aura toutes ses composantes égales à 1. La matrice A sera au choix A1 ou A2 :
   
2 −1 4 −1
−1 2 −1 0  −1 4 −1 0 
   

 −1 2 −1 


 −1 4 −1 

A1 =  .. .. ..  A2 =  .. .. .. 

 . . . 


 . . . 

 .. ..   .. .. 
 0 . . −1  0 . . −1
−1 2 −1 4
1. Télécharger le fichier TP2.m depuis hippocampus.

2. Compléter l’initialisation de la matrice A avec au choix A1 ou A2 , et du second membre b,
vecteur colonne avec toutes ses coordonnées égales à 1.
3. Compléter la résolution par les méthodes LU, Jacobi et Gauss-Seidel à l’aide des fonctions
précédemment écrites. Pour chaque méthode, utiliser tic et toc pour obtenir le temps de
résolution et calculer le résidu relatif kAx − bk2 /kbk2 , comme cela est fait pour la résolution
par linsolve.
4. Pour la méthode de Jacobi, on rappelle que la matrice d’itération est J = D−1 (E + F ) =
D−1 (D − A). On obtient respectivement pour A1 et A2 (cf. calcul analogue en TD) :
   
0 1 0 1
1 0 1 0  1 0 1 0 
   
 1 0 1   1 0 1 
1  1 
J1 =  .. .. ..  J2 =  .. .. .. 
2
 . . . 
 4
 . . . 

 .. ..   .. .. 
 0 . . 1  0 . . 1
1 0 1 0
Définir ces deux matrices dans le script puis faire le calcul du rayon spectral en utilisant
que ρ(J) = kJk2 car J est une matrice réelle symétrique. La norme matricielle indice 2
d’une matrice peut s’obtenir avec la fonction Matlab norm(J,2).
5. Pour la méthode de Gauss-Seidel, ajouter le calcul du rayon spectral de la matrice d’ité-
ration G en utilisant que ρ(G) = ρ(J)2 , résultat vu (admis) en TD pour une matrice A
tridiagonale.
6. Tests et tableaux de résultats. Pour chacune des deux matrices A1 et A2 , compléter
les tableaux de la fiche de résultats pour les différentes méthodes (linsolve, lu_tridiag,
jacobi_tridiag, gauss_seidel_tridiag), successivement pour N = 10, 100, 1000 : temps de
21/49
calcul, résidu relatif, et nombre d’itérations (pour les méthodes itératives). Pour les mé-
thodes de Jacobi et Gauss-Seidel, choisir tolerance assez petite pour que le résidu obtenu
soit du même ordre de grandeur que celui obtenu par la méthode LU (pour comparer des
résolutions de précisions similaires), et itermax assez grand pour que la précision souhaitée
soit atteinte (le calcul peut prendre plusieurs minutes dans certains cas).
2.3 Analyse des résultats

1. Valeur du conditionnement et son influence sur les résultats : à l’aide des tableaux de ré-
sultats pour chaque matrice A1 et A2 , observer si le conditionnement de la matrice aug-
mente et si cela a une influence sur le résidu obtenu avec la résolution LU. Vérifier que ces
résultats sont en accord avec les résultats théoriques.
Indications sur les résultats théoriques : en reprenant les résultats de l’exercice 1.4 des TD
concernant la matrice A = h12 A1 = (N + 1)2 A1 , on obtient :

k π
Spectre(A1 ) = {4 sin2 ; k = 1, . . . , N }
N +1 2
et
4
cond2 (A1 ) = cond2 (A) ∼ (N + 1)2 quand N → +∞
π2
Puis, en remarquant que A2 = A1 + 2I, on obtient :

2 k π
Spectre(A2 ) = {4 sin + 2 ; k = 1, . . . , N }
N +1 2
et
λmax (A2 ) 4+2
cond2 (A2 ) = ∼ π2
∼ 3 quand N → +∞
λmin (A2 ) +2
(N +1)2
2. Convergence des méthodes itératives Jacobi et Gauss-Seidel : à l’aide des tableaux de ré-
sultats pour chaque matrice A1 et A2 , observer pour les méthodes Jacobi et Gauss-Seidel, le
lien entre la valeur du rayon spectral de la matrice d’itération et la vitesse de convergence.
Vérifier que ces résultats sont en accord avec les résultats théoriques.
Indications sur les résultats théoriques : en reprenant ce qui a été fait dans l’exercice 2.3
des TD pour la matrice A = h12 A1 = (N + 1)2 A1 , on obtient pour A1 :
1
J = D−1 (E + F ) = D−1 (D − A1 ) = I − A1
2
de valeurs propres :

2 k π k
µk = 1 − 2 sin = cos π , k = 1, . . . , N
N +1 2 N +1
d’où
π
ρ(J) = cos → 1 quand N → +∞
N +1
et
2 2 π
ρ(G) = ρ(J) = cos → 1 quand N → +∞
N +1
22/49
Puis, pour A2 , on obtient :

1
J = D−1 (E + F ) = D−1 (D − A2 ) = I − A2
4
de valeurs propres :

1 2 k π 1 k
µk = 1 − 4 sin +2 = cos π , k = 1, . . . , N
4 N +1 2 2 N +1
d’où
1 π 1
ρ(J) = cos → quand N → +∞
2 N +1 2
et
1 π 1
ρ(G) = ρ(J) = cos2 2
→ quand N → +∞
4 N +1 4
3. Facultatif : Etudier l’influence du point initial pour les deux méthodes itératives.
4. Facultatif : pour pouvoir augmenter davantage la taille N de la matrice, transformer la
fonction lu_tridiag en remplaçant la matrice pleine A par trois vecteurs contenant res-
pectivement sa sous-diagonale, sa diagonale et sa sur-diagonale. Refaire les tests avec cette
nouvelle fonction pour observer l’évolution du résidu quand N augmente. On doit observer
que le résidu augmente beaucoup pour A1 tandis qu’il reste petit pour A2 , ce qui est en
accord avec les résultats théoriques sur la valeur du conditionnement et son influence sur
la résolution du système linéaire.
23/49
24/49
Chapitre 3
Optimisation sans contrainte
Objectif
— Programmer, tester et comparer plusieurs méthodes pour minimiser une fonction sans
contrainte : la méthode de gradient à pas constant et la méthode de Newton (méthodes à
programmer). On comparera aussi aux résultats données par la méthode de Nelder-Mead
(fonction fminsearch de Matlab).
— Connaître les propriétés des différentes méthodes (hypothèses requises sur la fonction à
minimiser, vitesse de convergence, ...).
— Observer les résultats sur plusieurs fonctions ayant des propriétés différentes (α-convexe
ou non, coercive ou non, admettant un ou plusieurs points de minimum local, ...).
— Être capable de choisir une méthode d’après les propriétés de la fonction à minimiser.
Fichiers à utiliser
Télécharger dans un dossier de travail (Z:\MATLAB\TP3 par exemple), à partir d’hippocampus

(cours MATHS, CHAPITRE 3), les fichiers f.m, gradf.m, hessienf.m, test_fmin.m.
numérique OCTAVE (https ://www.gnu.org/software/octave/).
Travail à rendre
— Un compte-rendu dans un fichier TP3_GROUPE_Nom1_Nom2.pdf, comprenant
les résultats des tests (tableaux et graphes) et les réponses aux questions.
— Le code des fonctions fgradient.m et fnewton.m et du script test_fmin.m, dans une
archive TP3_GROUPE_Nom1_Nom2.zip.
— Déposer le compte-rendu et l’archive du code sur hippocampus, cours MATHS_S5 ou
MATHS_S6 suivant le semestre.
25
3.1 Résultats théoriques sur les exemples de fonctions à

minimiser
Fonction no 0 : f : R2 → R, f (x) = (x1 − 1)2 + x22
Solution évidente : point de minimum (global) en x = (1, 0), minimum qui vaut 0.
Fonction no 1 : f : R2 → R, f (x) = (x1 − 12 )2 (x1 + 1)2 + (x2 + 1)2 (x2 − 1)2

On calcule le gradient de f en x :
4(x1 − 21 )(x1 + 1)(x1 + 41 )

∇f (x) =
4(x2 + 1)(x2 − 1)x2
On en déduit les points stationnaires :
1 1 1 1 1 1
(−1, −1), (−1, 0), (−1, 1), (− , −1), (− , 0), (− , 1), ( , −1), ( , 0), ( , 1).
4 4 4 2 2 2
On calcule la matrice hessienne en x :

λ1 (x) 0
∇2 f (x) =
0 λ2 (x)
de valeurs propres λ1 (x) = 12x21 + 6x1 − 23 et λ2 (x) = 12x22 − 4. On calcule :
(
λ1 (−1, x2 ) = 29 , λ1 (− 14 , x2 ) = − 49 , λ1 ( 12 , x2 ) = 92
λ2 (x1 , −1) = 8, λ2 (x1 , 0) = −4, λ2 (x1 , 1) = 8
Conclusion : f admet un minimum local aux points (−1, −1), (−1, 1), ( 21 , −1), ( 21 , 1), qui vaut 0
(c’est en fait un minimum global car f (x) ≥ 0, ∀x), un maximum local au point (− 41 , 0), et un
point selle aux quatre autres points stationnaires.
Fonction no 2 : f : R2 → R, f (x) = exp(x1 ) − x1 + x22 + x2

On calcule le gradient de f en x :

exp(x1 ) − 1
∇f (x) =
2x2 + 1
On en déduit l’unique point stationnaire (0, − 21 ).
On calcule la matrice hessienne de f en x :

2 exp(x1 ) 0
∇ f (x) =
0 2
de valeurs propres λ1 (x) = exp(x1 ) et λ2 (x) = 2, strictement positives. Donc f strictement
convexe. De plus, on vérifie que f est coercive (cf. TD).
Conclusion : f admet un unique point de minimum, le point (0, − 21 ), et le minimum vaut
f (0, − 12 ) = 0.75.
Fonction no 3 : f : R3 → R, f (x) = x21 + x1 x2 + x22 − x2 x3 + 23 x23 − 9x1 − x2 − 7x3

On vérifie que c’est une fonction quadratique, de la forme f (x) = 21 (Ax, x) − (b, x) avec :
   
2 1 0 9
A = 1 2 −1 et b = 1
0 −1 3 7
26/49
On peut vérifier que la matrice symétrique A est définie-positive, par exemple en calculant ses
valeurs propres avec la fonction eig de Matlab. On en déduit que f admet un point de mini-
mum global strict en x solution du système linéaire Ax = b. La résolution de ce système donne
x = (5, −1, 2).
Conclusion : f admet un point de minimum global strict en (5, −1, 2), minimum qui vaut
f (5, −1, 2) = −29.
Les fichiers f.m, gradf.m et hessienf.m contiennent respectivement les fonctions :
— function y = f(x,choixf) qui calcule f (x) pour la fonction no choixf,

— function g = gradf(x,choixf) qui calcule le vecteur gradient g(x) pour la fonction
no choixf,
— function H = hessienf(x,choixf) qui calcule la matrice hessienne H(x) pour la fonc-
tion no choixf.
3.2 Minimisation par la fonction fminsearch de Matlab

function [x,fval,exitflag,output] = fminsearch(f,x0,options)
Voir la documentation avec la commande doc fminsearch.

Cette fonction utilise la méthode itérative de Nelder-Mead dont le principe est de partir d’un
simplexe de l’espace (c’est-à-dire un triangle dans R2 , une pyramide dans R3 , ...), et à chaque
itération de remplacer un des sommets par un point voisin où la fonction est plus petite. Quand
l’algorithme converge, on obtient un simplexe de taille très petite autour du point de minimum.
Compléter le fichier test_fmin.m pour appliquer cette méthode aux fonctions no 0, 1, 2, et 3 en

suivant les indications données en commentaire.
La fonction fminsearch est appelée avec en entrée un point initial x0 et une précision souhai-
tée tolerance (champ ’TolX’ du paramètre options, initialisé par la fonction optimset) : le
calcul des itérés s’arrête quand la distance entre les deux derniers itérés est inférieure à cette va-
leur. En sortie, on indique l’algorithme utilisé (champ algorithm du résultat output, donné par
output.algorithm), le point de minimum obtenu x, la valeur de ce minimum fval et le nombre
d’itérations effectuées (champ iterations du résultat output, donné par output.iterations).
Pour les fonctions no 0, 1 et 2 qui sont définies sur R2 , le script test_fmin.m propose le tracé de
graphes qui peuvent aider à localiser les minima recherchés :
— surfaces y = f (x) avec les fonctions meshgrid et surfl,
— courbes de niveaux de f avec la fonction contour,
— visualisation avec la fonction quiver des vecteurs gradients approchés à l’aide de la fonc-
tion gradient,
— possibilité de mettre plusieurs graphes dans la même fenêtre avec la fonction subplot.
Exécuter le script test_fmin.m pour faire les tests suivants :
27/49
Fonction 0 : vérifier sur cet exemple simple que le script fonctionne bien.
Fonction 1 : tester avec une tolérance de 10−4 en choisissant successivement les points initiaux
x0 = (0, 0), (0, 0.5), (0.4, 0.5), qui sont dans la zone du point de minimum (0.5, 1). Si besoin,
augmenter le nombre maximum d’itérations (champ ’MaxIter’ du paramètre options) et le
nombre maximum d’appels à la fonction f (champ ’MaxFunEvals’ du paramètre options), à
l’aide de la fonction optimset. S’aider de la documentation avec la commande doc fminsearch.
Questions :
1. Faire un tableau de résultats dans lequel on indiquera pour chaque point initial : le point
de minimum obtenu, la valeur en ce point, le nombre d’itérations, la nature du point
(minimum local, maximum local ou point selle).
2. Est-ce que la méthode converge vers le point de minimum (0.5, 1) avec chacun des points
x0 choisis ? Est-ce que le nombre d’itérations varie beaucoup ?
3. Facultatif : compléter l’étude en choisissant un autre point x0 .
Fonction 2 : tester avec une tolérance de 10−4 , et le point initial x0 = (−1.5, −1.5) puis un autre
point initial de son choix.
Questions :
1. Faire comme pour la fonction no 1 un tableau regroupant les résultats obtenus.
2. L’unique point de minimum (0, −0.5) est-il correctement obtenu quel que soit le choix
de x0 ? Est-ce que le nombre d’itérations varie beaucoup ?
Fonction 3 : tester avec une tolérance de 10−4 , et le point initial x0 = (0, 0, 0) puis un autre point
initial de son choix.
Questions :
1. Faire comme pour les fonctions no 1 et 2 un tableau regroupant les résultats obtenus.
2. L’unique point de minimum (5, −1, 2) est-il correctement obtenu quel que soit le choix
de x0 ? Est-ce que le nombre d’itérations varie beaucoup ?
3.3 Minimisation par la méthode du gradient à pas constant

En partant d’un point initial x0 ∈ Rn , on calcule les itérés xk par la formule de récurrence :
xk+1 = xk − ρ∇f (xk )
où ρ > 0 est un réel fixé (indépendant de k).
Si f est de classe C 1 , α-convexe, et que ∇f est lipschitzienne, alors on sait qu’il existe ρ > 0
suffisamment petit tel que la suite (xk ) converge vers l’unique point de minimum x∗ de f , quel
que soit le point de départ x0 (cf. TD et EXOS_SUP).
Si ces conditions suffisantes de convergence ne sont pas vérifiées, on peut quand même appliquer
la méthode mais la suite (xk ) ne convergera pas forcément vers un point de minimum, et le
résultat pourra dépendre du point initial.
Ecrire une fonction Matlab :
28/49
function [x,fval,cvge,iter] = fgradient(choixf,x0,rho,tolerance,itermax)
qui réalise le calcul les itérés xk de la méthode du gradient, avec le critère de convergence
kxk+1 − xk k ≤ tolerance, en un nombre maximum d’itérations itermax, c’est-à-dire qu’on
calcule les itérés tant que kxk+1 − xk k > tolerance et k < itermax. Les résultats seront x le
dernier itéré calculé, fval la valeur de f en x, cvge un indicateur de convergence, et iter le
nombre d’itérations effectuées. Indication : suivre le même algorithme que celui utilisé pour la
fonction jacobi_tridiag du TP2.
Compléter cette fonction fgradient pour qu’elle affiche le chemin suivi par les itérés successifs sur
le graphe des courbes de niveaux, pour les exemples dans R2 (fonctions no 0, 1 et 2). Indication :
plot([x(1),xsuiv(1)],[x(2),xsuiv(2)],’k-*’) tracera un segment noir entre les points x
et xsuiv et placera des * sur ces deux points. Adapter la couleur et le symbole (* ou + ou autre)
pour que le tracé se voit bien.
Compléter le script Matlab test_fmin.m pour appliquer cette fonction aux différents exemples
étudiés.
Questions :
1. Faire les mêmes tests qu’avec la fonction fminsearch. On pourra choisir ρ égal à 10−1 ,
et le diminuer s’il ne convient pas pour obtenir la convergence.
2. Regrouper les résultats dans un tableau, sans oublier de préciser la valeur de ρ.
3. Donner les graphes pour les fonctions no 1 et 2.
4. Est-ce que les points de minimum recherchés sont correctement obtenus avec chacun des
points x0 choisis ?
5. Dans le cas où la méthode converge vers un point qui n’est pas un point de minimum,
expliquer ce qui se passe.
6. Est-ce que le nombre d’itérations varie beaucoup en fonction de x0 ? En fonction de ρ ?
7. Déduire plusieurs conseils pour utiliser la méthode du gradient.
3.4 Minimisation par la méthode de Newton

Soit f : Rn −→ R deux fois dérivable. Un extremum local de f vérifie la condition nécessaire
∇f (x) = 0. La méthode de Newton recherche les solutions de ∇f (x) = 0, autrement
dit les points stationnaires de f . En partant d’un point initial x0 ∈ Rn , on calcule les itérés
xk ∈ Rn par le schéma itératif :
(
H(xk )qk = ∇f (xk ),
xk+1 = xk − qk
où H(xk ) est la matrice hessienne formée des dérivées partielles secondes de f au point xk .
Indication : le calcul de qk se fait en résolvant le système linéaire H(xk )qk = ∇f (xk ) (pour éviter
le calcul, plus coûteux, de l’inversion de la matrice). Attention : le vecteur ∇f (xk ) doit être en
colonne.
Il restera à caractériser la solution trouvée : minimum local, maximum local, point selle.
Ecrire une fonction Matlab :
29/49
function [x,fval,cvge,iter] = fnewton(choixf,x0,tolerance, itermax)
qui réalise le calcul des itérés, avec le critère de convergence kxk+1 − xk k ≤ tolerance, en un
nombre maximum d’itérations itermax. Les résultats seront x le dernier itéré calculé, fval la
valeur de f en x, cvge un indicateur de convergence, et iter le nombre d’itérations effectuées.
Prévoir aussi que la fonction fnewton affiche le chemin suivi par les itérés successifs sur le graphe
des courbes de niveaux (pour les exemples dans R2 ).
Compléter à nouveau le script Matlab test_fmin.m pour appliquer cette méthode sur les diffé-
rents exemples.
Questions :
1. Faire les mêmes tests qu’avec les fonctions fminsearch et fgradient, des tableaux de
résultats analogues et donner les graphes pour les fonctions no 1 et 2.
2. Pour lesquelles des fonctions la méthode a-t-elle permis d’obtenir un point de mini-
mum pour chacun des points x0 choisis ? En déduire les propriétés d’une fonction qui
conviennent pour rechercher un point de minimum par la méthode de Newton.
3. Dans ces situations, que remarque-t-on sur le nombre d’itérations par rapport aux mé-
thodes précédentes ?
4. Est-ce que le nombre d’itérations varie beaucoup en fonction de x0 ? Dans le cas de la
fonction no 3, le nombre d’itérations est très faible et ne varie pas. Expliquer pourquoi.
5. Déduire des questions précédentes les situations où il est intéressant d’utiliser la méthode
de Newton.
6. Dans les autres situations, faut-il conseiller la méthode de Nelder-Mead (fonction fmin-
search) ou la méthode du gradient ? Donner des arguments pour guider dans ce choix.
30/49
Chapitre 4
Optimisation avec contraintes
Objectif
— Savoir choisir une fonction de minimisation sous contraintes, suivant le logiciel de calcul
utilisé (Matlab ou Octave) et suivant le type problème (linéaire ou non).
— Savoir écrire les contraintes sous la forme attendue par la fonction de minimisation utilisée.
— Comprendre ce que signifie qu’une contrainte soit active ou inactive.
— Adapter la méthode du gradient, programmée dans le TP précédent pour un problème
sans contrainte, aux cas où il est possible à chaque itération de projeter l’itéré sur le
domaine des contraintes.

(cours MATHS, CHAPITRE 4), les fichiers Matlab test_linprog.m, fobj.m, nonlincon.m,
test_fmincon.m.
Tout ce qui est fait dans ce TP peut être réalisé avec le logiciel libre de calcul numé-
rique OCTAVE (https ://www.gnu.org/software/octave/) mais les fonctions de mi-
nimisation sont un peu différentes. Dans ce cas, télécharger les fichiers Octave test_glpk.m,
fobj.m, nonlincon.m, nonlinconeq.m, test_sqp.m.
Travail à rendre
— Un compte-rendu dans un fichier TP4_GROUPE_Nom1_Nom2.pdf comprenant
les résultats des tests (valeurs, tableaux, graphes) et les réponses aux questions.
— Le code des fonctions et des scripts dans une archive TP4_GROUPE_Nom1_Nom2.zip
31
4.1 Optimisation linéaire

On minimise une fonction linéaire sous contraintes linéaires, en distinguant les contraintes en
inégalité, les contraintes en égalité et les contraintes de borne, comme ceci :

Ax ≤ b

min J(x) = (f, x) sous les contraintes Aeq x = beq
x 
lb ≤ x ≤ ub

4.1.1 Résolution d’un premier problème de logistique

Une ferme souhaite utiliser un capital de 18 000 euros pour semer du colza, du tournesol et des
pois sur une surface de 12 ha, tout en limitant la quantité de fertilisant utilisé à 1600 kg au total.
Les caractéristiques de ces cultures sont données dans le tableau suivant. On cherche comment
répartir les 12 hectares entre colza, tournesol et pois pour maximiser le revenu qui reviendra à
la ferme.
coût en euros/ha fertilisant en kg/ha revenu pour la ferme en euros/ha

colza 1600 270 155
tournesol 1300 150 160
pois 1700 110 150
Table 4.1 – une ferme
Modélisation
Ce problème peut s’écrire sous la forme du problème d’optimisation linéaire suivant :

16 x1 + 13 x2 + 17 x3 ≤ 180


270 x + 150 x + 110 x ≤ 1600
1 2 3
max J(x) = 155 x1 + 160 x2 + 150 x3 sous les contraintes
x x1 + x2 + x3 = 12


x1 ≥ 0, x2 ≥ 0, x3 ≥ 0

(4.1)
Résolution par la fonction linprog de Matlab
La fonction linprog fait partie de la boîte à outils optimization de Matlab.
function [x,fval,exitflag,output,lambda] =
linprog(f,A,b,Aeq,beq,lb,ub,x0,options)
Pour l’utiliser, s’aider de la documentation, qui s’obtient avec doc linprog, et compléter le script
test_linprog.m :
— Compléter l’initialisation de f,A,b,Aeq,beq,lb,ub pour le problème (4.1).
— Si l’indicateur exitflag vaut 1, ce qui indique le succès du calcul, faire afficher les résultats
suivants :
— le nombre d’itérations effectuées output.iterations,
— la solution x et la valeur du min obtenu fval,
— le vecteur lambda.ineqlin, relatif aux contraintes Ax ≤ b,
32/49
— le vecteur lambda.lower, relatif aux contraintes lb ≤ x.
Résolution par la fonction glpk de Octave (si impossibilité d’utiliser Matlab)
function [x,fval,errnum,extra] =
glpk(f,A,b,lb,ub,ctype,vartype,s)
Pour l’utiliser, s’aider de la documentation, qui s’obtient avec doc glpk, et compléter le script
test_glpk.m :
— Compléter l’initialisation de f,A,b,lb,ub,ctype,vartype pour le problème (4.1).
— Si l’indicateur errnum vaut 0, ce qui indique le succès du calcul, faire afficher les résultats
suivants :
— la solution x et la valeur du min obtenu fval,
— le vecteur extra.lambda, relatif aux contraintes écrites dans A et b,
Questions :
1. Donner les résultats obtenus. En déduire la répartition optimale des 12 hectares en colza,
tournesol et pois pour maximiser le revenu qui reviendra à la ferme.
2. Les coordonnées des vecteurs lambda.ineqlin et lambda.lower de la fonction linprog
de Matlab, ou du vecteur extra.lambda de la fonction glpk de Octave indiquent si les
contraintes sont actives ou non :
— coordonnée nulle pour une contrainte inactive,
— coordonnée non nulle pour une contrainte active : dans ce cas la solution est sur la
frontière de cette contrainte, la contrainte est vérifiée avec égalité.
Est-ce que les coordonnées non nulles de ces vecteurs correspondent à des contraintes
pour lesquelles l’égalité est vérifiée avec la solution x ?
3. Sélectionner une contrainte inactive et résoudre à nouveau le problème sans cette
contrainte. Donner les résultats obtenus. Est-ce qu’on obtient la même solution ?
4.2 Résolution d’un deuxième problème de logistique

Le but est d’optimiser la livraison d’une marchandise . On dispose de M entrepôts, indicés par
i = 1, · · · M , disposant chacun d’un stock si . Il faut livrer N clients, indicés par j = 1, N , qui
ont commandé chacun une quantité rj . Le coût de transport entre l’entrepôt i et le client j est
donnée par ci,j . Les inconnues sont les M N quantités vi,j de marchandise partant de l’entrepôt
i vers le client j. On veut minimiser le coût de transport tout en satisfaisant les commandes des
clients et sous la contrainte de limite des stocks. Ce problème se formule ainsi :
PN
X M X N Pj=1 vi,j ≤ si , i = 1, ..., M

M
min J(v) = ci,j vi,j sous les contraintes i=1 vi,j = rj , j = 1, ..., N (4.2)
v∈RM N 
i=1 j=1 
vij ≥ 0, i = 1, ..., M, j = 1, ..., N
On va considérer ce problème pour les données suivantes :

— M = 2 entrepôts, de stocks s1 = 200 et s2 = 500,
— N = 5 clients, ayant commandé les quantités r1 = 22, r2 = 38, r3 = 8, r4 = 44, r5 = 28,
— coûts de transport
33/49
. de l’entrepôt 1 vers les clients 1 à 5 : c1,1 = 10, c1,2 = 5, c1,3 = 8, c1,4 = 1, c1,5 = 15,
. de l’entrepôt 2 vers les clients 1 à 5 : c2,1 = 8, c2,2 = 12, c2,3 = 11, c2,4 = 20, c2,5 = 9.
Transformation du problème
On utilise un vecteur x ∈ R10 pour représenter les quantités recherchées vi,j :
x1 = v1,1 , . . . , x5 = v1,5 , x6 = v2,1 , . . . , x10 = v2,5 ,
10
et un vecteur f ∈ R pour représenter les coûts ci,j :
f1 = c1,1 , . . . , f5 = c1,5 , f6 = c2,1 , . . . , f10 = c2,5 .
Écrire les contraintes sous forme matricielle :
PN
1. les contraintes j=1 vi,j ≤ si , i = 1, ..., M sous forme matricielle Ax ≤ b,
PM
2. les contraintes i=1 vi,j = rj , j = 1, ..., N sous forme matricielle Aeq x = beq ,
3. et enfin les contraintes vij ≥ 0, i = 1, ..., M, j = 1, ..., N sous la forme lb ≤ x.
Résolution par la fonction linprog de Matlab ou glpk de Octave

Compléter le script test_linprog.m de Matlab, ou le script test_glpk.m de Octove pour pouvoir
résoudre ce problème.
Questions :
1. Quelles sont pour chaque client les quantités de marchandises venant de chacun des
deux entrepôts ? Regrouper les résultats dans un tableau par exemple. Est-ce que les
contraintes en inégalité sont actives ou non ? Dire ce que cela signifie pour le stock des
entrepôts.
2. Résoudre à nouveau en multipliant par 5 les quantités demandées par les clients. Ana-
lyser la situation et les résultats.
3. Mêmes questions en multipliant les quantités demandées par 6.
4.3 Optimisation non linéaire

On minimise une fonction quelconque sous contraintes quelconques, en distinguant les contraintes
linéaires, les contraintes de borne et les contraintes non linéaires, comme ceci :



 Ax ≤ b
Aeq x = beq



min J(x) sous les contraintes lb ≤ x ≤ ub
x 
c(x) ≤ 0





c (x) = 0
eq
4.3.1 Résultats théoriques sur les exemples à étudier

Problème no 1

1 2 x 1 ≥ 0

min J(x) = −3x1 + x2 sous les contraintes x2 ≥ 0
x 2 
 2
x1 + x22 ≤ 1
34/49
L’ensemble des contraintes est un quart de disque, fermé borné convexe non vide. La fonction J
est continue et convexe, non strictement. J étant continue sur K fermé borné non vide de R2 ,
on sait que la fonction J admet au moins un point de minimum (global) sur K. De plus, J et
K étant convexes, les points de minimum de J sur K sont points de minimum global et ce sont
les solutions des conditions (KKT). On obtient un et un seul point de minimum (global) x∗ de
J sur K :
x∗ = (1, 0), J(x∗ ) = −3
Problème no 2

2 4 −3
min J(v) = (Av, v) sous la contrainte kvk = 1 avec A =
v −3 10
2
√ 1 de R . La matrice A est symétrique, définie-
L’ensemble des contraintes est le cercle de rayon
positive, de plus petite valeur propre 7 − 3 2. D’après la résolution faite en TD, J admet
exactement deux points de minimum (global) x∗ et −x∗ sous ces contraintes :
 √
x∗ = ( √ 1 √ , √ 2−1√ ) ' (0, 9239; 0, 3827)
4−2 2 4−2 2
J(x∗ ) = J(−x∗ ) = 7 − 3√2 ' 2, 7574
Problème no 3
On veut installer une nouvelle antenne pour améliorer un réseau de communication. Cette an-
tenne doit se trouver le plus près possible de 4 clients particuliers, en donnant la priorité à
ceux qui ont utilisé le plus de temps de communication. D’autre part, cette antenne doit se
trouver à au moins 10 km des deux autres antennes déjà existantes du réseau. On représente la
localisation des antennes et des clients dans un repère orthonormée du plan (unité = km). Les
coordonnées des 2 antennes déjà existantes sont (−5; 10) et (5; 0). Les coordonnées des clients
sont {(5; 10), (10; 5), (0; 12), (12; 0)} et leurs nombres d’heures de communications respectives
{200; 150; 200; 300}.
La variable du problème est la position x de la nouvelle antenne. On choisit comme fonction
objectif la somme des carrés des distances entre x et les quatre clients, pondérée par les nombres
d’heures de communications :
J(x) = 2((x1 − 5)2 + (x2 − 10)2 ) + 1, 5((x1 − 10)2 + (x2 − 5)2 )
+2(x21 + (x2 − 12)2 ) + 3((x1 − 12)2 + x22 )
Enfin, les contraintes s’écrivent :
(
(x1 + 5)2 + (x2 − 10)2 ≥ 100
(x1 − 5)2 + x22 ≥ 100
L’ensemble K des contraintes est le plan R2 privé de la réunion de deux disques ouverts. C’est
un fermé non vide, non borné, et non convexe. On peut vérifier que J est coercive sur K. J
étant continue et coercive sur K fermé non vide de R2 , on en déduit que J admet au moins un
minimum (global) sur K.
La résolution des conditions (KKT) nécessairement vérifiées par tous les points de minimum
(local ou global) conduit aux résultats suivants :
min global en x∗ ' (8, 3307; 9, 4112) qui vaut ' 515,




min local en (−5; 0) qui vaut 1980,


 min local en ' (1, 6193; −9, 4112) qui vaut ' 2704,
min local en ' (−14, 5140; 13, 0794) qui vaut ' 4718

35/49
4.3.2 Résolution par la fonction fmincon de Matlab
La fonction fmincon fait partie de la boîte à outils optimization de Matlab.
function [x,fval,exitflag,output,lambda] =
fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
Pour l’utiliser, s’aider de la documentation donnée avec doc fmincon et des commentaires au
début du script test_fmincon.m.
1. Utiliser la fonction fobj.m :

function [J,gJ,HJ] = fobj(x,choixpb)
qui, suivant le problème choisi choixpb, calcule en x la valeur J de J, le vecteur gJ gradient
de J et la matrice HJ matrice hessienne de J.
2. Compléter la fonction nonlincon.m :

function [c,ceq,gc,gceq] = nonlincon(x,choixpb)
qui, suivant le problème choisi choixpb, calcule en x les contraintes en inégalité c et
les contraintes en égalité ceq et leurs gradients respectifs gc et gceq. S’il y a plusieurs
contraintes, par exemple en inégalité, c est alors un vecteur de contraintes, et gc une ma-
trice contenant les gradients de ces contraintes, rangés par colonne.
3. Compléter le script test_fmincon.m pour résoudre les trois problèmes par la fonction
fmincon. En sortie, on affiche l’algorithme utilisé, le point de minimum obtenu, la valeur
de ce minimum et le nombre d’itérations effectuées. Ajouter sur le graphe des courbes de
niveau l’affichage du point initial choisi x0 et du point de minimum obtenu.
Indication : comme nos fonctions fobj et nonlincon dépendent de x et de choixpb, utili-
ser @(x)fobj(x,choixpb) à la place de fun et @(x)nonlincon(x,choixpb) à la place de
nonlcon.
4. Résoudre les trois problèmes, avec la tolérance 10−4 (déjà fixée au début de test_fmincom.m
par le champ ’TolX’ du paramètre options, initialisé par la fonction optimoptions) et
avec différents choix du point initial x0 :
— pour le problème no 1, x0 = (0; 0, 5), (5; 5),
— pour le problème no 2, x0 = (10; 10), (0; −2),
— pour le problème no 3, x0 = (−15; −15), (−5; −15), (0; −15), (−15; 15), (15; 15).
Rassembler les résultats dans des tableaux.
Questions :
1. Regrouper les résultats obtenus dans des tableaux. Donner un graphe correspondant
pour chaque problème.
2. Pour le problème no 1, est-ce que le point initial x0 a une forte influence sur le nombre
d’itérations ?
3. Pour le problème no 2, pour lequel il existe plusieurs points de minimum, quel conseil
donner ?
4. Pour le problème no 3, expliquer les difficultés pour trouver le point de min global.
36/49
4.3.3 Résolution par la fonction sqp de Octave, si impossibilité d’utiliser Matlab

function [x,fval,info,iter,nf,lambda] =
sqp(x0,fobj,nonlinconeq,nonlincon,lb,ub)
Pour l’utiliser, s’aider de la documentation donnée avec doc sqp.

1. La fonction objectif et ses dérivées seront calculées par la fonction fobj.m, comme avec
Matlab.
2. Les contraintes non linéaires sont définies dans deux fonctions disctinctes à compléter :
function [ceq,gceq] = nonlinconeq(x,choixpb)
function [c,gc] = nonlincon(x,choixpb)
qui, suivant le problème choisi choixpb, calculent en x respectivement les contraintes en
égalité ceq et leurs gradients gceq et les contraintes en inégalité c et leurs gradients gc.
S’il y a plusieurs contraintes, ceq et c sont des vecteurs, et gceq et gc sont des matrices
contenant les gradients des contraintes, rangés par colonne. Attention : la fonction sqp
considère la contrainte sous la forme c(x) ≥ 0, tandis qu’en cours on a utilisé la convention
c(x) ≤ 0.
3. Compléter le script test_sqp.m pour résoudre les trois problèmes par la fonction sqp, en
suivant les mêmes instructions que pour le script test_fmincon avec Matlab. Indication :
comme nos fonctions fobj, nonlincon et nonlinconeq dépendent de x et de choixpb, uti-
liser @(x)fobj(x,choixpb) à la place de fobj et analogue pour les deux autres fonctions.
4. Faire les mêmes tests que ceux demandés dans la résolution par fmincon de Matlab.
4.3.4 Résolution par la méthode du gradient projeté à pas constant

En partant d’un point initial x0 ∈ Rn , on calcule les itérés xk par la formule de récurrence :
xk+1 = p(xk − ρ∇f (xk ))
où ρ > 0 est un réel fixé (indépendant de k) et p la projection orthogonale (pour la norme
euclidienne) sur l’ensemble K des contraintes.
Si K est un convexe fermé non vide, alors la projection est bien définie. Si de plus f est de classe
C 1 , α-convexe, et que ∇f est lipschitzienne, alors on sait qu’il existe ρ > 0 suffisamment petit
tel que la suite (xk ) converge vers l’unique point de minimum x∗ de f sur K, quel que soit le
point de départ x0 (cf. EXOS_SUP).
Aucun des trois problèmes ne vérifie l’ensemble de ces conditions suffisantes de convergence.
Cependant, on peut facilement définir le projeté sur K pour les deux premiers problèmes. On
pourra donc calculer les itérés dans ces deux cas mais la suite (xk ) ne convergera pas forcément
vers un point de minimum, et le résultat pourra dépendre du point initial.
1. Écrire la fonction :
function y = projete(x,choixpb)
qui calcule le projeté orthogonal y sur K d’un point x suivant que choixpb vaut 1 ou 2.
Rappel : le projeté orthogonal de x sur K est le point y de K tel que kx − yk est minimum,
avec k.k norme euclidienne. Indication : résoudre graphiquement à la main pour déduire
un algorithme convenable.
37/49
2. Écrire la fonction (en reprenant la fonction fgradient du TP3) :

function [x,fval,cvge,iter] =
fgradientprojete(choipb,x0,rho,tolerance,itermax)
qui réalise le calcul les itérés xk de la méthode du gradient projeté, avec le critère de
convergence kxk+1 − xk k ≤ tolerance, en un nombre maximum d’itérations itermax,
c’est-à-dire qu’on calcule les itérés tant que kxk+1 − xk k > tolerance et k < itermax.
Les résultats seront x le dernier itéré calculé, fval la valeur de f en x, cvge un indicateur
de convergence, et iter le nombre d’itérations effectuées. De plus, la fonction fgradient
affichera le chemin suivi par les itérés successifs sur le graphe des courbes de niveaux.
3. Compléter le script Matlab test_fmincon.m (ou le script Octave test_sqp.m) pour appli-
quer cette méthode aux problèmes 1 et 2. Faire les mêmes tests qu’avec la fonction fmincon
(ou sqp), avec tolerance = 1e-4 et en choisissant rho=1e-1, qu’on pourra diminuer si
besoin pour essayer d’obtenir la convergence de la méthode. Rassembler les résultats dans
des tableaux, sans oublier d’indiquer la valeur du paramètre rho.
Questions :
1. Regrouper les résultats obtenus dans des tableaux. Donner les graphes correspondants.
2. Est-ce que la méthode de gradient projeté donne des résultats similaires à ceux de la
fonction fmincon pour résoudre les problèmes no 1 et 2 ?
3. Plus généralement pour un problème avec un domaine des contraintes sur lequel il facile
de projeter les itérés, quel peut être l’avantage d’utiliser la méthode de gradient projeté ?
38/49
Chapitre 5
Lois usuelles – Convergence de

suites de variables aléatoires
Objectif
— Comprendre ce qui est représenté par un histogramme.
— Savoir obtenir des simulations de variables aléatoires suivant des lois usuelles.
— Etudier différents types de convergences d’une suite de variables aléatoires.
(cours MATHS, CHAPITRE 5), le script Matlab TP5.m.
numérique OCTAVE (https ://www.gnu.org/software/octave/). Toutes les fonctions
utilisées sont identiques sauf pour la loi de Poisson (randp au lieu de poissrnd).
Travail à rendre
les résultats des tests et les réponses aux questions.
5.1 Simulations et histogrammes de lois usuelles

Supposons qu’on dispose d’un vecteur x d’un grand nombre de valeurs xi , i = 1..N d’une variable
alétoire X donnée. Pour obtenir une représentation approchée de la loi de distribution suivie par
39
X, on va représenter un histogramme de 10 classes de même amplitude associé aux valeurs xi .

Pour cela, on réalise les opérations suivantes :
— on utilise la fonction Matlab [nelements,centers] = hist(x), qui pour un vecteur
donné x, détermine les vecteurs nelements et centers donnant respectivement les nombres
d’éléments et les centres des classes obtenues en divisant l’intervalle [mini xi , maxi xi ] en
10 intervalles de même longueur h.
— on souhaite que les surfaces des rectangles de l’histogramme soient égales aux fréquences
des valeurs xi obtenues dans les classes correspondantes : h × hauteur = nelements/N .
Pour cela on définit un vecteur hauteurs égal au vecteur nelements divisé par N h.
— on trace l’histogramme avec la fonction Matlab bar(centers,hauteurs)
Ainsi, les surfaces des rectangles de l’histogramme donnent des approximations de P (X ∈ I)
où I est une des classes. En particulier les valeurs des hauteurs de l’histogramme donnent des
approximations de la densité de probabilité quand X est une variable aléatoire continue.
Questions : Pour les lois usuelles qui suivent, faire les graphes et histogrammes demandés.
Observer si les histogrammes obtenus sont cohérents avec les valeurs des espérances et va-
riances. Dans le cas des lois continues, vérifier que la hauteur des histogrammes correspond
approximativement à la densité de probabilité.
5.1.1 Loi continue uniforme sur [0,1]

Utiliser la fonction rand de Matlab pour simuler la réalisation de N = 1000 valeurs d’une variable
aléatoire X qui suit la loi continue uniforme U(0, 1). Faire un graphe donnant le nuage de points
obtenus (i, xi ), i = 1, . . . , N , et un graphe donnant l’histogramme associé à ces valeurs xi (avec
la fonction subplot de Matlab les deux graphes sont dans la même figure). Refaire ce test avec
N = 5000 ou N = 10000.
5.1.2 Loi binomiale

Soit X une variable aléatoire discrète qui suit une loi binomiale B(n, p). Une réalisation de X est
obtenue en additionnant n réalisations yj de variables aléatoires indépendantes Yj qui suivent
la loi de Bernoulli B(p). Pour obtenir une réalisation yj de Yj , on utilise une réalisation u d’une
variable aléatoire U de loi U(0, 1) obtenue avec la fonction rand. Comme on a :
P (Yj = 1) = p = P (U ∈ [0, p])
on fixera yj = 1 si u ≤ p et yj = 0 sinon. Ecrire une fonction :

function x = function binomiale(N,n,p)
qui calcule un vecteur x de N réalisations de X qui suit la loi B(n, p) pour n et p donnés. Puis, à
l’aide de cette fonction, dessiner le nuage de points et l’histogramme obtenus avec N (N = 1000
par exemple) réalisations de X pour n = 20 et p = 0, 3. Refaire ce test avec un autre couple
(n, p).
5.1.3 Loi exponentielle

1
Si U est une variable aléatoire qui suit la loi uniforme U(0, 1), alors pour λ > 0, X = − ln U est
λ
une variable aléatoire qui suit la loi exponentielle E(λ) (cf. TD). En utilisant ce résultat, écrire
une fonction :
40/49
function x = exponentielle(N,lambda)
qui calcule un vecteur x de N réalisations de X qui suit la loi E(λ) pour λ donné. Puis, à l’aide
de cette fonction, dessiner le nuage de points et l’histogramme obtenus avec N (N = 1000 par
exemple) réalisations de X pour λ = 3. Faire de même pour une autre valeur de λ.
5.1.4 Loi de Poisson
Utiliser la fonction poissrnd de Matlab (toolbox Statistics) pour obtenir N (N = 1000 par
exemple) réalisations de X qui suit la loi P(λ) pour λ = 3. Dessiner le nuage de points et
l’histogramme correspondants. Faire de même pour une autre valeur de λ.
5.1.5 Loi normale
Utiliser la fonction randn de Matlab, qui simule des réalisations d’une variable aléatoire de loi
normale centrée réduite N (0, 1), pour écrire une fonction :
function x = normale(N,m,v)
qui calcule un vecteur x de N réalisations de X de loi N (m, v) pour une moyenne m et une
variance v données. Puis, à l’aide de cette fonction, dessiner le nuage de points et l’histogramme
obtenus avec N (N = 1000 par exemple) réalisations de X pour m = 3 et v = 10.
5.2 Loi des grands nombres, convergence presque sûre et

convergence en probabilité
5.2.1 Loi forte des grands nombres
On considère une suite de variables aléatoires (Xn ) indépendantes, suivant une même loi de pro-
babilité et telles que E(|Xn |) < +∞. On note m = E(Xn ) (même espérance pour tout n car
n
1X
même loi pour tous les Xn ), et on note X̄n = Xi . Alors la suite (X̄n ) converge presque
n i=1
sûrement vers m, ce qui signifie que l’événement {ω ∈ Ω / X̄n (ω) → m} est de probabilité 1.
Etant donné qu’on ne peut pas simuler un événement de probabilité nulle, ce résultat signifie
que toute réalisation de la suite (X̄n ) donne une suite de réels qui converge vers m.
Questions : Vérifier qu’on observe ce résultat en réalisant les tests suivants :

1. Pour nmax assez grand (au moins 100), calculer des réalisations xn , n = 1, . . . , nmax,
d’une variable aléatoire X qui suit U(0, 1) et en déduire les valeurs correspondantes x̄n
de X̄n , n = 1, . . . , nmax. On pourra utiliser la fonction cumsum de Matlab. Dessiner les
points (n, x̄n ), n = 1, . . . , nmax et relever la valeur de x̄nmax . Faire plusieurs fois ce test.
2. Faire la même chose avec une loi de Poisson P(λ).
41/49
5.2.2 Convergence presque sûre et convergence en probabilité
Soit (Xn )n≥1 la suite de variables aléatoires, où chaque Xn prend les valeurs n ou 0, avec les
probabilités :
1 1
P (Xn = n) = √ et P (Xn = 0) = 1 − √
n n
Cette suite converge en probabilité vers 0. En effet, pour tout ε > 0, on a (à partir du rang
n > ε) :
1
P (|Xn | > ε) = P (Xn = n) = √ → 0 quand n → +∞
n
En revanche, cette suite ne converge pas presque sûrement vers 0. On peut le démontrer en
utilisant le résultat suivant (qu’on admet) :
Une suite de variables aléatoires indépendantes
X (Xn )n≥1 converge presque surement vers 0 si et
seulement si pour tout ε > 0, la série P (|Xn | > ε) est convergente.
n≥1
Questions :
1. Calculer des réalisations xn de Xn , n = 1, . . . , nmax, pour nmax assez grand (au moins
100). Tracer les points (n, xn ) obtenus. Comment observe-t-on que (Xn ) ne converge
pas sûrement vers 0 ?
2. Facultatif : Proposer un test pour observer que la probabilité pour que Xn = n diminue
quand n augmente.
3. On remplace la probabilité √1n par n12 . Montrer que Xn converge presque sûrement vers
0. Refaire le test de la quaetion 1. Est-ce qu’on observe la convergence de presque sûre
de (xn ) vers 0 ?
5.3 Théorème central limite, convergence en loi, approxi-

mation de lois
5.3.1 Théorème central limite
Si on considère une suite de variables aléatoires indépendantes (Xn ) suivant une même loi de
n
1X
probabilité, de moyenne m et de variance σ 2 finie, et qu’on note X̄n = Xi , alors la suite
n i=1
X̄n − m
Zn = converge en loi vers une variable aléatoire Z qui suit la loi normale N (0, 1).
√σ
n
En pratique, on approche la loi de Zn par la loi normale N (0, 1) quand n > 30.
Questions :
1. Pour une suite de variables aléatoires indépendantes (Xn ) qui suivent la loi uniforme
U(0, 1), représenter l’histogramme obtenu avec N (N = 1000 par exemple) réalisations
de Znmax , pour nmax assez grand (au moins 100). Comparer à un histogramme obtenu
avec la loi N (0, 1), avec le même nombre N de réalisations.
2. Faire la même chose avec une loi de Poisson P(λ).
42/49
5.3.2 Approximation d’une loi binomiale par une loi de Poisson

Si on considère une suite de variables aléatoires indépendantes (Xn ) suivant la loi binomiale
B(n, p), alors quand n → +∞, de telle sorte que np → λ, λ > 0, la suite (Xn ) converge en loi
vers une variable aléatoire X qui suit la loi de Poisson P(λ).
En pratique, on approche la loi B(n, p) par la loi P(np) quand n > 50 et np < 18.
Question : Mettre en évidence ce résultat à l’aide d’histogrammes d’une loi B(n, p) et de la
loi P(np), pour un exemple de valeurs n et p convenables.
5.3.3 Approximation d’une loi binomiale par une loi normale

Si on considère une suite de variables aléatoires indépendantes (Xn ) suivant la loi binomiale
Xn − m
B(n, p), alors quand n → +∞, la suite Zn = √ , où on note q = 1 − p, converge en loi vers
npq
une variable aléatoire Z qui suit la loi normale N (0, 1).
En pratique, on approche la loi B(n, p) par la loi N (np, npq) quand n > 50 et np > 18.
Question : Mettre en évidence ce résultat à l’aide d’histogrammes d’une loi B(n, p) et de la
loi N (np, npq), pour un exemple de valeurs n et p convenables.
5.3.4 Approximation d’une loi de Poisson par une loi normale

Si on considère une suite de variables aléatoires indépendantes (Xn ) suivant la loi de Poisson
Xn − n
P(n), alors quand n → +∞, la suite Zn = √ converge en loi vers une variable aléatoire Z
n
qui suit la loi normale N (0, 1).
En pratique, on approche la loi P(λ) par la loi N (λ, λ) quand λ > 18.
Question : Mettre en évidence ce résultat à l’aide d’histogrammes d’une loi P(λ) et de la loi
N (λ, λ), pour un exemple de valeur de λ convenable.
43/49
44/49
Chapitre 6
Evaluation d’une proportion –

Moindres carrés avec contraintes
linéaires
Objectif
— Déterminer un intervalle de confiance et comprendre l’information qu’il donne.
— Savoir résoudre un problème de moindres carrés avec contraintes linéaires, et étudier un
exemple d’application.
Il n’y a pas de fichier à télécharger.
numérique OCTAVE (https ://www.gnu.org/software/octave/). Toutes les fonctions
utilisées sont identiques sauf pour la fonction (qp au lieu de lsqlin).
Travail à rendre
les résultats des tests et les réponses aux questions.
6.1 Évaluation d’une proportion

On s’intéresse à une élection entre deux candidats C0 et C1 . On réalise un sondage auprès d’un
grand nombre d’électeurs, pour obtenir la réalisation (x1 , . . . , xn ) d’un échantillon (X1 , . . . , Xn )
45
de taille n, où Xi vaut 1 ou 0 suivant que le vote déclaré est pour C1 ou C0 (les votes blancs ou
indécis ne sont pas retenus pour l’échantillon).
Le problème est d’évaluer, à partir de l’échantillon (X1 , . . . , Xn ), la proportion p des suffrages

pour C1, parmi tous les suffrages exprimés. Les variables aléatoires Xi , i = 1 . . . , n, sont indé-
pendantes et suivent la loi de Bernoulli B(p). On va utiliser un estimateur de p pour prédire le
résultat de l’élection, et construire un intervalle de confiance In pour p, de niveau asymptotique
(1 − α), pour α donné. On rappelle que pour X de loi B(p), on a E(X) = p et V (X) = p(1 − p).
6.1.1 Construction de l’estimateur

On note p̄n la proportion
Pn de l’échantillon qui possède la modalité 1, c’est-à-dire la proportion
Pn des
xi égaux à 1 : p̄n = n1 i=1 xi . C’est une réalisation de la moyenne d’échantillon X̄n = n1 i=1 Xi ,
qui est un estimateur de la proportion recherchée p, et qui a les propriétés suivantes :
— sans biais : E(X̄n ) = p,
— convergent p.s., d’où la convergence de la suite de réels (p̄n ) vers p, pour toute réalisation
de l’échantillon, quand n → +∞,
— convergent L2 avec V (X̄n ) = p(1−p) n → 0 quand n → +∞, donc la dispersion de X̄n
autour de son espérance E(X̄n ) = p tend vers 0 quand n → +∞.
6.1.2 Construction d’un intervalle de confiance

On cherche un intervalle In tel que P (p ∈ In ) ' 1 − α. On dispose du résultat suivant (obtenu
par le théorème Central Limite et la convergence p.s. de X̄n vers p, cf TD) :
X̄n − p
La variable aléatoire Zn = q converge en loi vers une variable aléatoire Z qui suit la
X̄n (1−X̄n )
n
loi normale N (0, 1).
On considère le réel z α2 tel que P (|Z| ≤ z α2 ) = 1 − α pour Z qui suit la loi N (0, 1). Indication :
pour α = 5%, z α2 ' 1.96. Avec la convergence en loi de Zn vers Z, on obtient :
 s 
X̄n (1 − X̄n ) 
P |X̄n − p| ≤ z α2 → P (|Z| ≤ z α2 ) = 1 − α quand n → +∞
n
On en déduit l’intervalle de confiance pour p, de niveau asymptotique (1 − α) :

 s s 
X̄n (1 − X̄n ) X̄n (1 − X̄n )
In = X̄n − z α2 , X̄n + z α2 
n n
Avec une réalisation p̄n de X̄n , on calcule un intervalle In . On pourra considérer que n est assez
grand pour utiliser l’approximation P (p ∈ In ) ' 1 − α, si np̄n ≥ 5 et n(1 − p̄n ) ≥ 5.
6.1.3 Mise en oeuvre sur des échantillons d’intentions de vote

Écrire un script proportion.m pour faire les tests suivants dans le cas où l’intention
réelle de vote de la population entière pour le candidat C1 est p = 51% :
46/49
1. Pour N = 10000 simuler la réalisation (x1 , . . . , xN ) d’un échantillon (X1 , . . . , XN ) des

intentions de vote. Sur un graphe avec n en abscisse et les probabilités en ordonnée, re-
présenter les estimations p̄n et les intervalles de confiance In de niveau asymptotique 95%,
obtenus pour n = 1000, 2000, 3000, ..., 10000 (en prenant les 1000 premiers éléments de
l’échantillon de taille N , puis les 2000 premiers, etc, comme si on agrandissait la taille de
l’échantillon au fur et à mesure en sondant davantage de personnes). Utiliser la fonction
errorbar pour tracer les intervalles de confiance.
Question : Donner le graphe obtenu. Quelles observations peut-on faire quand la taille n
de l’échantillon augmente ?
2. Refaire ce test et le graphe correspondant plusieurs fois (à chaque fois la réalisation
(x1 , . . . , xN ), N = 10000, est nouvelle).
Question : Observer qu’avec certains échantillons, on peut se tromper dans la prédiction,
c’est-à-dire prédire que le candidat C1 perd l’élection. Est-ce que ce problème semble ap-
paraitre pour n’importe quelle valeur de la taille n de l’échantillon ? Illustrer par plusieurs
exemples de graphes obtenus.
3. Déterminer l’intervalle de confiance IN , pour p de niveau asymptotique 95%, uniquement
pour N = 10000, pour 100 réalisations (x1 , . . . , xN ) de l’échantillon (X1 , . . . , XN ). Comp-
ter le nombre de fois où la valeur exacte de p (p = 51%) appartient à l’intervalle IN ?
Question : Donner les résultats obtenus. Que signifie que l’intervalle de confiance a un
niveau asymptotique de 95% ? Est-ce que ce test permet d’observer cela ?
6.2 Moindres carrés avec contraintes linéaires
6.2.1 Présentation du problème

1 Ax ≤ b

min J(x) = kCx − dk22 sous les contraintes Aeq x = beq
x 2 
lb ≤ x ≤ ub

On va s’intéresser ici au cas avec contraintes égalités uniquement (Aeq x = beq ), c’est-à-dire aux
problèmes de la forme :
1
min J(x) = kCx − dk22 sous la contrainte Aeq x = beq (6.1)
x 2
Rappel du résultat du cours (Chapitre 4) :

On suppose que 
C ∈ Mm,n ,

d ∈ Rm ,

m≥n

et pour les contraintes : 

Aeq ∈ Mp,n ,

beq ∈ Rp ,

p≤n

47/49
Le système Aeq x = beq représente p contraintes affines sur les n coordonnées de x.
Si C est de rang n et Aeq de rang p, alors J admet un et un seul point de minimum x∗ sur
l’ensemble des contraintes, et il existe un vecteur de multiplicateurs de Lagrange λ∗ ∈ Rp tel que
(x∗ , λ∗ ) est solution du système linéaire par blocs suivant :
CT C ATeq
T
x C d
= (6.2)
Aeq 0 λ beq
6.2.2 Résolution du problème

1. Ecrire une fonction :
function [x,Jval,residu,lambda,exit] = lsqlineq(C,d,Aeq,beq)
qui écrit et résout le système (6.2). En sortie, la fonction donnera le point x où le min est
atteint, la valeur Jval=J(x), le vecteur residu=Cx-d, le vecteur lambda des multiplicateurs
de Lagrange et un indicateur exit qui vaut 1 si le calcul du min a été possible, 0 sinon.
2. Écrire un script mcarres.m pour tester cette fonction sur le problème suivant :
On cherche x ∈ K = {x ∈ R3 ; Aeq x = beq } qui vérifie au mieux (au sens des moindres
carrés) le système sur-déterminé Cx = d, avec :
   
1 0 0 65, 5
1 1 1 −142
   
1 2 4  82, 5  4 0 1 7
C= , d=
 , Aeq = , beq = .
 −12 
1 3 9  3 1 0 7
 
1 4 16  25 
1 5 25 −45
On doit obtenir la solution x∗ = (2; 1; −1) et le vecteur de multiplicateurs de Lagrange

λ∗ = (−1; 2).
3. Facultatif : comparer avec les résultats obtenus par la fonction lsqlin de Matlab (qui fait
partie de la boîte à outils optimization de Matlab) ou par la fonction qp de Octave.
6.2.3 Approximation-interpolation d’un nuage de points
On cherche un polynôme de degré inférieur ou égal à trois qui interpole les deux points (−3, 2) et
(3, 5), et passe au mieux, au sens des moindres carrés discrets, par les points (xi , yi ), i = 1, ..., 5
rassemblés dans le tableau suivant :
xi -2 -1 0 1 2
yi 6 4 5 3 1
48/49
Questions :
1. Écrire ce problème sous la forme (6.1) d’un problème de moindres carrés avec contraintes
linéaires.
2. Compléter le fichier mcarres.m pour résoudre ce problème. Quel polynôme obtient-on ?
3. Faire une figure avec les points d’interpolation, les points d’approximation et le graphe
du polynôme sur l’intervalle [-3,3].
49/49

Poly TP Maths 2023 2024

Transféré par

Droits d'auteur :

Formats disponibles

Poly TP Maths 2023 2024

Transféré par

Informations du document

Description originale:

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Poly TP Maths 2023 2024

Transféré par

Droits d'auteur :

Formats disponibles

Année 2023-2024

1 Initiation Matlab – Equation de la chaleur 1D 5

3 Optimisation sans contrainte 25

4 Optimisation avec contraintes 31

5 Lois usuelles – Convergence de suites de variables aléatoires 39

6 Evaluation d’une proportion – Moindres carrés avec contraintes linéaires 45

Initiation Matlab – Equation de

1.1 Premières commandes, variable ans, point-virgule ;

1.2 Variables prédéfinies et codage double précision des

realmin ' 2−1024 ' 10−308

realmax ' 21023 ' 10308

eps = 2−52 ' 10−16

>> x = 1e20; y = 1e3;

1.3 Nombres complexes

1.4 Les commentaires, les affichages de message

>> pi % pour connaitre la valeur de pi

>> disp(’la valeur de pi avec 15 décimales vaut : ’); disp(pi);

ou la fonction fprintf, avec %f pour un format décimal et %e pour un format scientifique :

>> fprintf(’la valeur de pi avec 10 décimales vaut pi = %15.10f\n’,pi);

1.5 Vecteurs - Matrices

1.5.1 Création d’un vecteur ligne

Un vecteur ligne peut être défini de plusieurs façons :

>> x1 = [1.1 2.3 3.5 -4. -8.] % (vecteur de 5 composantes)

1.5.2 Vecteur transposé, création d’un vecteur colonne

Un vecteur colonne peut être défini de plusieurs façons :

1.5.3 Opérations sur les vecteurs

Quelques opérations sur les vecteurs :

u+v : addition de deux vecteurs u et v de même longueur,

1.5.4 Création de matrices

A = [a11 a12 a1p ; a21 a22 a2p ; · · · an1 an2 anp ]

(ii) ou à partir de vecteurs ou de matrices plus petites par concaténation :

1.5.5 Opérations sur les matrices

Exemples avec les fonctions spy, size et diag :

1.6 M-Files ou scripts

function [z,t] = fonc(x,y,m)

>> [z,t] = fonc(1., 0., -4.)

>> f = @(x,y) x^2+y^2;

1.8 Boucles et choix

1.8.2 Choix simple ou alternatif (if)

1.8.3 Choix multiple selon un paramètre (switch)

switch m % faire selon la valeur de m

function y = solex(x, ug, ud, m)

1.8.4 Boucle itérative pour (for)

1.8.5 Boucle itérative tant que (while)

% recherche l’indice du 1er coefficient non nul (numériquement) d’un vecteur u

1.10 Quelques fonctions mathématiques

Résolution d’un système linéaire

2.1 Fonctions de résolution d’un système linéaire tridiago-

(a) function x = lu_tridiag(n,A,b) :

(b) function [x,k,converge] = jacobi_tridiag(n,A,b,itermax,tolerance) :

(c) function [x,k,converge] = gauss_seidel_tridiag(n,A,b,itermax,tolerance) :

Indication : trame de la fonction jacobi_tridiag :

2.2 Tests avec le script TP2.m

1. Télécharger le fichier TP2.m depuis hippocampus.

2.3 Analyse des résultats

Puis, pour A2 , on obtient :

Optimisation sans contrainte

Télécharger dans un dossier de travail (Z:\MATLAB\TP3 par exemple), à partir d’hippocampus