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