Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
Initiation à MATLAB
TP4
Exercice 1
Traduire en langage Matlab l’algorithme PGCD suivant :
a,b entiers positifs
tant-que a �b faire
si a > b alors a �a - b
sinon b � b - a
fin tant-que
PGCD � a
Corrigé de l’exercice 1
Pour traduire en langage Matlab l’algorithme PGCD considéré dans cet exercice, nous avons besoin
des structures de contrôle, en l’occurrence la boucle while et le test if-else-end. Nous proposons le
M-file suivant pour le programme traduisant cet algorithme:
clear all;
close all;
clc;
a=input('Entrer l''entier positif a ');
b=input('Entrer l''entier positif b ');
a0=a;
b0=b;
while a~=b
if a>b
a=a-b;
else
b=b-a;
1
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
end
end
str=['Le plus petit diviseur commun de ' num2str(a0) ' et '
num2str(b0) ' est ' num2str(a)];
disp(str);
Remarques:
*Le programme M-file constituant la réponse à cet exercice n’est pas unique. Chacun a sa façon
personnelle de programmer. Le programme proposé est donc un exemple seulement de solution. Il
existe cependant des critères à satisfaire: le programme doit être sans redondances, facile à lire et
optimisé.
*On a fait appel aux variables auxiliaires a0 et b0 pour stocker les valeurs initiales de a et b qui sont
modifiées dans la boucle while et ne gardent pas donc leurs valeurs d’entrées à la sortie de cette
boucle. Elles sont en fait écrasées dans les lignes 10 et 12.
*On peut tester le programme en faisant appel à la commande gcd de Matlab qui calcule le plus
grand diviseur commun de deux nombre entiers positifs. Pour cela, il suffit de taper dans notre cas
(a0=18 ; b0=5).
>> gcd(a0,b0)
ans =
* Il faut faire la distinction entre la commande disp qui affiche la valeur de la variable sans écrire le
nom de la variable et la commande display qui écrit le nom de la variable avant de l’afficher. Taper
pour cela les deux commandes suivantes:
>> display(str)
str =
Le plus petit diviseur commun de 18 et 25 est 1
>> disp(str)
Le plus petit diviseur commun de 18 et 25 est 1
2
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
Exercice 2
Ecrire la fonction vec2col.m qui transforme tout vecteur (ligne ou colonne) passé en argument
en vecteur colonne. Un traitement d’erreur testera que la variable passée à la fonction est bien
un vecteur et non une matrice (utiliser la commande size).
Corrigé de l’exercice 2
On propose le M-file de type fonction suivant :
function [y]=vect2col(x)
n=size(x);
n1=n(1);
n2=n(2);
if n1==1
y=x.';
elseif n2==1
y=x;
else
erreur=['la variable entrée n''est pas un vecteur ligne ou
colonne'];
display(erreur)
end
On peut le tester avec les commandes
>> vect2col([1 1 1])
ans =
>> vect2col([1; 1; 1])
ans =
3
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
>> vect2col([1 1; 1 2])
erreur =
la variable entrée n'est pas un vecteur ligne ou colonne
Exercice 3
Ecrire un script qui permet selon le choix de transformer les cordonnées cartésiennes (x,y,z)
en coordonnés cylindriques (r, q,z) ou bien en coordonnées sphériques (r, q, j) .
Corrigé de l’exercice 3
Un exemple de script répondant à l’énoncé est donné dans la suite. Remarquer que l’on sort les
coordonnées cylindriques selon l’ordre : ( r, q,z ) et les coordonnées sphériques selon l’ordre
( r, q, j) .
function [a,b,c]=transformation(x,y,z,choix)
switch choix
case 'cylindrique'
a=sqrt(x^2+y^2);
b=atan(y/x);
c=z;
case 'sphérique'
a=sqrt(x^2+y^2+z^2);
b=atan(y/x);
c=atan(sqrt(x^2+y^2)/z);
otherwise
display('Vous n''avez pas préciser la nature de la
tansformation')
end
4
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
end
Exemples d’utilisation :
>> [a,b,c]=transformation(1,1,1,'cylindrique')
a=
1.4142
b=
0.7854
c=
>> [a,b,c]=transformation(2,3,4,'sphérique')
a=
5.3852
b=
0.9828
c=
0.7336
Avec Matlab vous avez quatre commandes qui vous permettent de faire la transformation des
coordonnées cartésiennes en coordonnées cylindriques ou bien des coordonnées cartésiennes en
coordonnées sphériques, ainsi que leurs inverses.
[theta,rho,z]=cart2pol(x,y,z) % transforamtion catésiennes en cylindriques
[x,y,z]=pol2cart(theta,rho,z) % transforamtion cylindriques en cartésiennes
[theta,phi,r]=cart2pol(x,y,z) % transforamtion catésiennes en sphériques
[x,y,z]=sph2cart(theta,phi,r) % transforamtion sphériques en catésiennes
5
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
Attention:
* Ici tous les angles sont calculés en radians (Matlab utilise cette unité
par défaut).
* Matlab sort les coordonnées cylindriques selon l’ordre ( q,r,z ) et les
coordonnées sphériques selon l’ordre ( q, j,r ) .
* Dans Matlab l’angle phi est compté à partir du plan horizontal et non pas
à partir de l’axe des z comme le montre la figure suivante :
Notre
choix
Définition des coordonnées sphériques selon Matlab :
phi est compté à partir du plan: phi=90-
Exercice 4
On considère l’algorithme suivant pour calculer π: on génère n couples { ( x k ,y k ) } de nombres
aléatoires dans l’intervalle [0, 1], puis on calcule le nombre m de ceux qui se trouvent dans le
premier quart du cercle unité, c'est-à-dire vérifiant la condition: x 2k + y 2k �.
1
4m
Naturellement, π est la limite de la suite u n = . Ecrire un programme Matlab pour calculer
n
cette suite et observer comment évolue l’erreur quand n augmente.
Corrigé de l’exercice 4
On peut proposer le script suivant :
6
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
clear all;
close all;
clc;
format long;
N=input('Entrer N? ')
x=rand(N,1);
y=rand(N,1);
z=x.^2+y.^2;
m=sum(z<=1);
uN=4*m/N;
display(uN);
erreur=(pi-uN)/pi;
display(erreur);
A l’exécution, on obtient :
Entrer N? 2.e7
N=
20000000
uN =
3.141851800000000
erreur =
-8.248886433788238e-005
La convergence est très lente avec le processus de Monte Carlo. Sur ma machine, on ne peut pas
dépasser N = 2.e7 car la mémoire se sature.
Ceux qui préfèrent les boucles et les tests peuvent être intéressés par la variante suivante:
7
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
clear all;
close all;
clc;
format long;
N=input('Entrer N? ')
m=0;
tic;
for k=1:N
xk=rand(1);
yk=rand(1);
if xk^2+yk^2<=1
m=m+1;
else
end
end
uN=4*m/N;
display(uN);
erreur=(pi-uN)/pi;
display(erreur);
toc;
En exécutant ce script, on obtient :
Entrer N? 2.e7
N=
20000000
uN =
3.141854000000000
erreur =
8
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
-8.318914608747151e-005
Elapsed time is 106.232164 seconds.
On a placé tic et toc pour compter la durée d’exécution, c'est-à-dire celle que l’on mesurerait en
consultant une horloge étalonnée déclenchée avant de lancer le programme et arrêtée à la fin de
l’exécution du programme (A ne pas confondre avec le temps CPU que l’on peut obtenir avec la
commande cputime). On a attendu dans ce cas à peu près 106 s.
Avec la première version, on obtient
Entrer N? 2.e7
N=
20000000
uN =
3.141189800000000
erreur =
1.282322803158388e-004
Elapsed time is 3.452282 seconds.
La durée d’exécution est meilleure avec la première version du programme. Elle est seulement de
3.3% de celle qui est exigée avec la deuxième version. Cependant la deuxième version ne souffre
pas de limitation au niveau de l’occupation de la mémoire car on peut dépasser facilement un
nombre de tirages supérieur à N = 2.e7 , du fait que l’on ne stocke pas en mémoire les variables: on
les écrase.
Un programme n’est jamais unique, mais il existe des critères qui font qu’une version de script est
meilleure que l’autre: la clarté, la vitesse d’exécution et l’occupation de la mémoire sont des
éléments importants.
L’étude de la convergence est résumée dans le tableau suivant:
N 1.e3 1.e4 1.e5 1.e6 1.e7 2.e7
Erreur 1.4513e-2 4.3267e-3 1.0973e-3 2.5612e-4 -7.0494e-6 1.5446e-4
Il faut remarquer que ces erreurs varient d’un tirage à l’autre et que le tableau représente le
résultat d’un tirage particulier. Par ailleurs la convergence n’est pas monotone comme le montre les
deux dernières colonnes.
9
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
Exercice 5
Ecrire un script qui permet de calculer par la méthode de dichotomie la racine de l’équation
nonlinéaire suivante: 5sin(x) - exp(x / 2) = 0 sur l’intervalle [0,1] .
Le script devra permettre de calculer par l’algorithme de dichotomie une approximation à h
près de la solution de l’équation f(x) = 0 dans l’intervalle [a,b], a, b et h étant saisis par
l’utilisateur et la fonction f étant définie par la commande inline.
Principe de la méthode de dichotomie :
Soit f une fonction continue sur [a, b] telle que f(a)f(b) < 0. Nécessairement f a au moins un
zéro dans ]a, b[ (ce résultat est un cas particulier du théorème des valeurs intermédiaires).
Supposons pour simplifier qu’il est unique et notons le α (dans le cas où il y a plusieurs zéros,
on peut localiser graphiquement, à l’aide de la commande plot, un intervalle qui n’en contient
qu’un).
La méthode de dichotomie consiste à diviser en deux un intervalle donné, et à choisir
le sous-intervalle où f change de signe. Plus précisément, si on note I(0) =]a, b[ et I(k) le
sous-intervalle retenu à l’étape k, on choisit le sous-intervalle I(k+1) de I(k) pour lequel f a un
signe différent à ses deux extrémités. En répétant cette procédure, on est assuré que chaque
I(k) ainsi construit contiendra α. La suite {m(k)} des milieux des intervalles I(k) convergera
vers α puisque la longueur de ces intervalles tend vers zéro quand k tend vers l’infini.
La méthode est initialisée en posant
a(0) = a, b(0) = b, I(0) =]a(0), b(0)[, m(0) = (a(0) + b(0))/2.
A chaque étape k ≥ 1, on choisit le sous-intervalle I(k)=]a(k), b(k)[ de
I(k−1)=]a(k−1), b(k−1)[ comme suit
étant donné m(k−1) = (a(k−1) + b(k−1))/2,
si f(m(k−1)) = 0,
alors α = m(k−1)
et on s’arrête ;
sinon, si f(a(k−1))f(m(k−1)) < 0
poser a(k) = a(k−1), b(k) = m(k−1);
si f(m(k−1))f(b(k−1)) < 0
poser a(k) = m(k−1), b(k) = b(k−1).
On définit alors m(k) = (a(k) + b(k))/2 et on incrémente k de 1.
A l’itération N la racine de f(x) = 0 est approximée par m(N) à h près si b(N) - a(N) �h .
�b - a �
On démontre que le nombre minimum des itérations est N min �log 2 � � -1 .
�h �
Corrigé de l’exercice 5
�b - a �
En fait, on démontre que le nombre minimum des itérations est N min �log 2 � � +1 et non
�h �
�b - a �
pas N min �log 2 � � -1 . Une solution possible pour cet exercice est donnée par le script
�h �
suivant :
10
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
clear all;
close all;
clc;
format long;
f=inline('5*sin(x)-exp(x/2)');
eta=1.e-12;
a0=0;
b0=1;
N=log2((b0-a0)/eta)+1;
for k=1:N
m0=(a0+b0)/2;
if f(m0)==0
elseif f(a0)*f(m0)<0
a1=a0; b1=m0;
else
a1=m0; b1=b0;
end
a0=a1;
b0=b1;
end
erreur=b0-a0;
display(erreur);
solution=m0;
display(solution);
%On peut vérifier que f(solution) est à peu près égal à zéro
valeurdef=f(m0);
display(valeurdef);
Noter que la fonction logarithme en base 2 s’écrit log2.
11
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
On n’était pas obligé de fixer le nombre des itérations a priori. En utilisant la boucle while sur
l’erreur avec un mouchard qui indique les itérations effectuées, le script devient
clear all;
close all;
clc;
format long;
f=inline('5*sin(x)-exp(x/2)');
eta=1.e-12;
a0=0;
b0=1;
N=log2((b0-a0)/eta)+1;
erreur=b0-a0;
n=0;
while erreur>eta
n=n+1;
m0=(a0+b0)/2;
if f(m0)==0
elseif f(a0)*f(m0)<0
a1=a0; b1=m0;
else
a1=m0; b1=b0;
end
a0=a1;
b0=b1;
erreur=b0-a0;
end
display(n);
display(erreur);
12
Ecole Nationale Des Sciences
Appliquées de Tétouan- ENSA 2AP-1
– 2015/2016
solution=m0;
display(solution);
%On peut vérifier que f(solution) est à peu près égal à zéro
valeurdef=f(m0);
display(valeurdef);
13
