MatLab le nouveau logiciel de programmation pour les mathematiciens

I. Le logiciel MatLab

MatLab, contraction de Matrix et de Laboratory,
est un logiciel optimisé pour la calcul matriciel. Cela a pour
conséquence qu'un programme, pouvant être fait soit à l'aide d'une ou
plusieurs boucles, soit à l'aide d'un calcul matriciel, s'exécutera
généralement plus vite en utilisant un calcul matriciel puisque, dans
ce cas, MatLab accède aux données dans un ordre adapté à leur stockage
en mémoire.
La gestion des boucles sous MatLab s'est améliorée avec les
dernières versions ; cependant, il est souvent plus avantageux de
vectorialiser ses calculs de manière à obtenir dans un code une
majorité de calculs matriciels à la place des boucles.

II. Configuration matérielle pour les tests

Processeur: 1,6 GHz

RAM: 2 Go

Swap: 4 Go


III. Comment vectorialiser un calcul

Tout d'abord, la plupart des fonctions MatLab peuvent agir sur des
vecteurs et non seulement sur des simples nombres ; par exemple, pour
calculer les valeurs du logarithme népérien pour les entiers de 1 à 10^6, stockées dans un vecteur ligne v, il suffit de faire:

v=log((1:1e6));

Vous pouvez alors vérifier, à l'aide des quelques lignes suivantes, que cette opération est plus rapide qu'une boucle.

tic;
v=zeros(1,1e6);
for p=1:1e6
v(p)=log(p);
end

IV. Expérience sur un exemple concret

Dans le but d'illustrer l'avantage de la vectorialisation, travaillons
sur un exemple concret. Savoir ce que fait le programme que nous allons
utiliser n'est pas une nécessité absolue, il a été choisi de manière à
montrer que vectorialiser un calcul n'est pas un artifice destiné à
améliorer la vitesse d'exécution de simples "cas d'école" mais peut
servir dans des cas réels de programmation. Il est donc possible de
passer directement aux sous-sections [Vous devez être inscrit et connecté pour voir ce lien] et [Vous devez être inscrit et connecté pour voir ce lien] sans que la compréhension de l'astuce en soit altérée.


IV.1. Exemple concret choisi

Le programme utilisé pour les tests calcule le champ magnétique généré dans l'espace par une boucle de courant.



Le champ magnétique généré par une telle boucle est donné par la loi de Biot et Savart:



IV.2. Programme non vectorialisé

function B = compute_B_loop(R,I,x,n)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calcule le champ magnétique généré par une boucle de courant
%% R est le rayon de la boucle
%% I est l'intensité du courant
%% x est la matrice des points où on veut calculer le champ B
%% La taille de x est 3 x nombre_de_points
%% ( x1 ... ...)
%% x = ( y1 ... ...)
%% ( z1 ... ...)
%% n est le nombre de points de discrétisation pour la boucle de courant
%% B est une matrice 3 x nombre_de_points
%% ( Bx1 ... ...)
%% B = ( By1 ... ...)
%% ( Bz1 ... ...)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu0 = 4*pi*10^(-7);

nx=size(x,2);

B=zeros(3,nx);

for p=1:nx
for q=0:n-1
normxmy3=((x(1,p)-R*cos(2*pi/n*q))^2+(x(2,p)-R*sin(2*pi/n*q))^2+x(3,p)^2)^(3/2);
B(1,p)=B(1,p)+cos(2*pi/n*q)*x(3,p)/normxmy3;
B(2,p)=B(2,p)+sin(2*pi/n*q)*x(3,p)/normxmy3;
B(3,p)=B(3,p)...
+(-sin(2*pi/n*q)*(x(2,p)-R*sin(2*pi/n*q))-...
cos(2*pi/n*q)*(x(1,p)-R*cos(2*pi/n*q)))/normxmy3;
end
end


B = mu0*R*I/(2*n)*B;

IV.3. Programme vectorialisé


function B = compute_B(R,I,x,n)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calcule le champ magnétique généré par une boucle de courant
%% R est le rayon de la boucle
%% I est l'intensité du courant
%% x est la matrice des points où on veut calculer le champ B
%% La taille de x est 3 x nombre_de_points
%% ( x1 ... ...)
%% x = ( y1 ... ...)
%% ( z1 ... ...)
%% n est le nombre de points de discrétisation pour la boucle de courant
%% B est une matrice 3 x nombre_de_points
%% ( Bx1 ... ...)
%% B = ( By1 ... ...)
%% ( Bz1 ... ...)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu0 = 4*pi*10^(-7);

nx=size(x,2);

%% x est la matrice des points de l'espace, on la répète
x=repmat(reshape(x.',1,nx,3),[n 1 1]);

%% y est la matrice des points de la boucle
y=R*[cos(2*pi/n*(0:n-1));sin(2*pi/n*(0:n-1));zeros(1,n)];

%% On répète la matrice y
y=repmat(reshape(y.',n,1,3),[1 nx 1]);

%% On peut alors calculer matriciellement tous les vecteurs SM (i.e. yx)
xmy=x-y; clear x y;

%% On calcule tous les vecteurs dl de la formule (on multipliera par R à la fin)
dy=[-sin(2*pi/n*(0:n-1));cos(2*pi/n*(0:n-1));zeros(1,n)];

%% On répète cette matrice
dy=repmat(reshape(dy.',n,1,3),[1 nx 1]);

%% On peut alors calculer matriciellement tous les vecteurs dl^SM de la formule (i.e. dy^yx)
cdyxmy=cross(dy,xmy,3); clear dy;

%% On calcule au passage toutes les normes de SM (yx)
norm_xmy_3=repmat(dot(xmy,xmy,3).^(3/2),[1 1 3]); clear xmy;

%% Il ne reste plus qu'à assembler B
B=cdyxmy./norm_xmy_3; clear cdyxmy norm_xmy_3;
B=sum(B);
B = mu0*R*I/(2*n)*B;
B = squeeze(B);
B = B.';
IV.4. Temps d'exécution pour des calculs importants

Exécutons maintenant ces deux programmes pour les points de [-5,5]x[-5,5]x[-5,5] sur une grille régulière 50x50x50 avec 50 points de discrétisation sur la boucle de courant.


L'exécution du programme non vectorialisé

[X Y Z]=meshgrid(linspace(-5,5,50));
x=[reshape(X,1,50^3);reshape(Y,1,50^3);reshape(Z,1,50^3)];
clear X Y Z;
tic;B_loop=compute_B_loop(1,1,x,50);toc;
clear B_loop; nous donne un temps d'exécution de 20,8s.


L'exécution du programme vectorialisé

[X Y Z]=meshgrid(linspace(-5,5,50));
x=[reshape(X,1,50^3);reshape(Y,1,50^3);reshape(Z,1,50^3)];
clear X Y Z;
tic;B=compute_B(1,1,x,50);toc;
clear B; nous donne un temps d'exécution de 14,4s.


Le programme non vectorialisé prend donc près de 1,5 fois plus de temps que le calcul vectorialisé.


Imaginez donc le résultat sur un programme qui tourne en quelques
heures. C'est toujours plus plaisant d'avoir les résultats après
"seulement" 3h que d'attendre 4h30. D'autant plus si les résultats sont faux et qu'il faut relancer le code après correction... ;-)



%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





toc; Le temps d'exécution pour cette boucle est de 3,7s, alors que le temps d'exécution du même calcul écrit de cette manière

tic;v=log((1:1e6));toc; est de 0,3s.


Soit un facteur 10 entre les deux temps d'exécution.


Ensuite, beaucoup de fonctions MatLab permettent de s'affranchir de
calculs à l'aide de boucles. En voici une liste non exhaustive:

• Les opérateurs arithmétiques *, /, ^
  peuvent agir sur des tableaux de même taille en effectuant l'opération
  terme à terme ; il suffit de rajouter un point devant l'opérateur pour
  préciser qu'on veut réaliser cette opération, c'est-à-dire .*, ./, .^. Par exemple, (1:10).^2 donnera les 10 premiers carrés parfaits.
  
• repmat - très simple d'utilisation, tapez help repmat ou doc repmat
  pour avoir une aide quant à cette fonction - permet de répéter un
  tableau (ce tableau peut être un vecteur ou une matrice) pour en
  construire un plus grand.
  
• reshape permet de redimensionner un tableau ;
  l'utilisation de cette fonction implique de connaître quelque peu
  comment MatLab "replace" les éléments lors d'un redimensionnement. De
  nouveau, l'aide MatLab donnée par help reshape ou doc reshape vous sera utile, vous pouvez consulter l'astuce sur les [Vous devez être inscrit et connecté pour voir ce lien] pour obtenir des explications différemment formulées.
  
• Beaucoup d'autres !
  

:Oo: Euh c'est difficile a déchiffré mais j'arrive avec difficulté de comprendre

Un peut difficile mais efficase