COURTE INTRODUCTION A MATLAB

Contenu

Mise en route
Elements de l'interface graphique
Philosophie
Variables, vecteurs, matrices
Aide
Principales commandes et syntaxe sur les vecteurs et matrices
Fonctions et scripts
Graphiques
Resoudre un systeme d'equations differentielles ordinaires
Codes

Mise en route

Pour creer un raccourci de l'application Matlab, ajouter dans votre fichier .profile les lignes suivantes

# Set PATH variable

export PATH=/Applications/MATLAB_R2012b.app/bin:$PATH

Si le fichier .profile n'existe pas vous pouvez le creer. Relancez ensuite la fenetre de terminal pour que le raccourci fonctionne.

A l'invite de commande dans votre fenetre de terminal, tapez matlab pour lancez Matlab en mode graphique. Pour lancez Matlab en ligne de commande seulement (sans fenetre), tapez matlab -nodesktop. Ce tutoriel suppose que vous avez lance Matlab en mode graphique.

Elements de l'interface graphique

La fenetre Matlab est divisee en plusieurs espaces qui peuvent etre configures. L'espace le plus important est la Command Windows. Cette fenetre permet d'executer des commandes Matlab. L'autre espace important est l'editeur de texte, qui vous permet d'editer des fichier qui pourront etre executes par Matlab.

Philosophie

Matlab, qui est l'abbreviation de matrix laboratory, est optimise pour manipuler des matrices. Pour en tirer partie, il faut vectorialiser les expressions. En C, pour initialiser un vecteur x qui prendra les valeurs de l a 1000000 (un million), on procede de la facon suivante:

double *x;
int i, sizex = 1000000;
x = malloc(size*sizeof(double));
for ( i = 0; i < sizex; i++)
{
    x[i] = i+1.0;
}

Sous Matlab, la meme procedure se fait de la facon suivante:

x = 1:1000000;

L'initialisation est vectorialisee: le membre de droite (1:1000000) est le vecteur 1 a 1000000 est l'assignation se fait de vecteur a vecteur.

Notez qu'il n'est pas necessaire de declarer la variable x, elle est cree lors de l'initialization. Le type sera automatiquement un array de double:

class(x)

ans =

double

Contrairement a C, les expressions n'ont pas a se terminer par un ;. Si le ; est omis, le resultat de l'expression sera affiche dans la Command Windows. L'expression

1+1;

ne produit aucin effet (mais le resultat est stocke dans une variable appelee ans). L'expression

1+1

ans =

     2

renvoie le resultat escompte: 2.

Les valeurs scalaires sont des arrays 1x1, et ne sont pas different des arrays representant des vecteurs ou des matrices.

On aurait pu initialiser x comme en C, mais cela aurait ete beucoup plus long. Les commandes tic et toc permettent de calculer un intervalle de temps ecoule:

clear x; % on efface la variable x
tic
x = 1:1000000;
toc

Elapsed time is 0.003285 seconds.

Maintenant, avec la boucle for

clear x;
tic
for i = 1:1000000
    x(i) = i;
end
toc

Elapsed time is 0.573735 seconds.

La boucle for est beaucoup plus lente que l'assignation vectorielle! Conclusion, si on peut, on evite les boucles for en Matlab.

Variables, vecteurs, matrices

La type de variable de base de Matlab est l'array de double. Il existe d'autres types (class en Matlab) comme les char, cell, struct, int8, etc, mais les double sont utilisees la plupart du temps. Comme les variables creees sont de classe double par defaut, on suppose qu'on travaillera en double dans ce tutoriel. Un array de double A peut comporter plusieurs dimensions. Pour un array de dimension 2, c.-a-d. un matrice (le defaut), on accede a l'entree ligne i, colonne j par A(i,j). (Comparer a C: A[i-1][j-1].)

Voici un liste de commandes pour creer des vecteurs (mx1 pour un vecteur colonne, et 1xn pour un vecteur ligne).

m = 10; n = 10;
xcol  = ones(m,1);      % cree un vecteur colonne de dimension 10 avec UNs partout
xline = ones(1,n);      % cree un vecteur ligne de dimension 10 avec UNs partout
ycol = (1:m)';          % cree un vecteur colonne des entiers de 1 a m. Le ' est l'adjoint
yline = 1:n;            % cree un vecteur ligne des entiers de 1 a n.
z1 = linspace(0,1,n);   % cree le vecteur ligne des n doubles espaces lineairement entre 0 et 1.
h = 1/(n-1);            % definit la distance entre deux elements de z1
z2 = 0:h:1;             % autre facon de definir le meme vecteur, mais attention aux erreurs d'arrondis!
isequal(z1,z2)          % on teste pour savoir si |z1 == z2|, le resultat est FAUX.

ans =

     0

pour connaitre les dimensions des variables xcol et xline, on peut utiliser les commandes length et size. La commande size renvoie la taille de chaque dimension

size(xcol)
size(xline)

La commande length renvoie la plus grande taille seulement, ce qui est pratique pour des expressions du type

mcols = length(xcol);

L'autre facon de connaitre l'etat des variables est la commande whos

whos

  Name        Size                   Bytes  Class     Attributes

  ans         1x2                       16  double              
  h           1x1                        8  double              
  i           1x1                        8  double              
  m           1x1                        8  double              
  n           1x1                        8  double              
  x           1x1000000            8000000  double              
  xcol       10x1                       80  double              
  xline       1x10                      80  double              
  ycol       10x1                       80  double              
  yline       1x10                      80  double              
  z1          1x10                      80  double              
  z2          1x10                      80  double

Cette commande affiche toutes les variables du workspace courant. Le workspace contient les variables definies par l'utilisateur.

De la meme facon que les vecteurs, on peut creer des matrices mxn.

A  = ones(m,n);         % cree une matrice 10x10 avec UNs partout
Ar = rand(m,n);         % cree une matrice 10x10 avec nombres aleatoires entre 0 et 1
A  = vander(z1);        % cree une matrice carree de Vandermonde de taille length(z1)

Aide

Pour de l'aide sur une commande, tapez

help nom_de_la_commande

Par exemple, pour afficher l'aide sur la commande ones

help ones

 ONES   Ones array.
    ONES(N) is an N-by-N matrix of ones.
 
    ONES(M,N) or ONES([M,N]) is an M-by-N matrix of ones.
 
    ONES(M,N,P,...) or ONES([M N P ...]) is an M-by-N-by-P-by-... array of
    ones.
 
    ONES(SIZE(A)) is the same size as A and all ones.
 
    ONES with no arguments is the scalar 1.
 
    ONES(..., CLASSNAME) is an array of ones of class specified by CLASSNAME.
 
    Note: The size inputs M, N, and P... should be nonnegative integers. 
    Negative integers are treated as 0.
 
    Example:
       x = ones(2,3,'int8');
 
    See also EYE, ZEROS.

    Overloaded methods:
       distributed/ones
       codistributor2dbc/ones
       codistributor1d/ones
       codistributed/ones
       gpuArray/ones

    Reference page in Help browser
       doc ones

La ligne See also est utile pour trouver des commandes similaire. Ici on trouve les commandes zeros et eye qui permettent d'initialiser des matrices (ou vecteurs) de zeros et la matrice identite.

Principales commandes et syntaxe sur les vecteurs et matrices

a = 2.0         % assigne la valeur 2 a la varialble a
(a == 2)        % teste l'egalite entre 2 et a
find(x == 200)  % trouve entrees i verifiant l'expression logique x(i) == 200
find(x < 10)    % trouve entrees i verifiant l'expression logique x(i) < 10
A*xcol          % multiplication matrice*vecteur
xline*A         % multiplication vecteur*matrice
A\x             % division a gauche (equivalent a A^(-1)*x)
A(:,1:2:end)    % colonnes impaires de A
A([1 2 5 6],:)  % lignes 1, 2, 5, et 6

Fonctions et scripts

Les fonctions Matlab sont le plus souvent ecrites dans un fichier texte separe. Le nom du fichier est le nom de la fonction. La fonction du fichier mafonction.m sera donc appelee par la commande mafonction. Construisons un fonction qui prend en argument deux vecteurs de meme taille, x et y, et qui renvoie les coefficients du polynome de degre n qui interpole au sens des moindres carres y aux points x: y ~ p(x).

edit monpolyfit.m

La premiere ligne doit etre

function p = monpolyfit(x,y,n)

Le mot cle function indique que le fichier est une fonction. p est la la variable retournee par la fonction et (x,y,n) sont les arguments de la fonction. La notation function [p1, p2] = mafonction(x,y) est utilisee s'il y a plusieurs variables de sortie et function mafunction(x,y) est utilisee si il n'y a pas de variable de sortie (pour une fonction qui genererait un graphique ou un fichier, par exemple).

La deuxieme ligne devrait etre un ligne descriptive de la fonction, et est en commentaire:

% MONPOLYFIT interpolation polynomiale par moindres carres

Cette ligne est utilisee pour afficher les descriptions courtes des fonction lorque, par exemple la commande lookfor est utilisee. Les lignes suivantes, en commentaire, decrivent l'usage de la fonction. C'est le moment de mettre le See also, pour faire un lien vers d'autres fonctions Matlab. Ici on pourrait faire un lien vers la fonction polyfit de Matlab.

% See also POLYFIT

(Les majuscules ne se verront pas avec la commande help).

Ensuite viens le corps de la fonction. Tout d'abord on verifie que les deux vecteurs x et y sont des vecteurs de meme taille, et que n+1 est inferieur ou egal a la taille de y.

sx = size(x);
sy = size(y);
if min(sx)>1 % alors pas un vecteur
    error('monpolyfit:tailledonnees','X n''est pas un vecteur');
end
if min(sy)>1 % alors pas un vecteur
    error('monpolyfit:tailledonnees','Y n''est pas un vecteur');
end
if sx(2) ~= 1 % alors pas un vecteur colonne (~= : pas egal)
    x = x.'; % on transpose (pas conjugue)
    sx = size(x); % on recalcule la taille de x
end
if sy(2) ~= 1 % alors pas un vecteur colonne (~= : pas egal)
    y = y.'; % on transpose (pas conjugue)
    sy = size(y); % on recalcule la taille de y
end
if (sx(1)~=sy(1)) % verifions les tailles de X et Y
    error('monpolyfit:tailledonnees','Les tailles de X et Y ne correspondent pas');
end
if (n+1)>length(x)
    error('monpolyfit:tailledonnees','Le degre n du polynome a ajuster est trop eleve, il doit etre inferieur a la longueur de X');
end

Maintenant on etablit notre equation normale pour les moindres carres. On genere un matrice de Vandermonde A a l'aide de x:

A = vander(x);

A est une matrice de taille length(x). On ne retient de la matrice que les n+1 dernieres lignes

m = sx(1);
A = A(:,(m-n):m);

On cherche les coefficients du polynome p qui approximent le mieux A*p=y. L'equation normale (en notation Matlab) est A'*A*p = A'*y. Si on fait une decomposition QR de la matrice A, on obtient l'equation R*p = Q'*y. On resoud ces equations a l'aide des fonctions householder, Qadj et backsubs.

[W,R]=householder(A);
Qy = Qadjb(W,y);
p = backsubs(R(1:(n+1),1:(n+1)),Qy(1:(n+1)));

Le code est termine, pas besoin de terminer le fichier avec un return et un end. Le vecteur p est a ete defini et sera retourne par la fonction. Toutes les variables creees dans la fonction sont locales a cette fonction et ne peuvent etre vues a l'exterieur de celle-ci. Les variables sont detruites lors du retour de la fonction, a moins de declarer ces variables persistent. Cela pourrait etre utile, par exemple si on voulait faire plusieurs interpolations de meme degre avec le meme vecteur x. Il serait alors plus efficace de ne pas recalculer A est sa decomposition a chaque appel de la fonction.

Pour resoudre un probleme, il est utile d'ecrire dans un fichier la liste des commandes a faire. Ce fichier n'est qu'une liste d'expressions a evaluer par Matlab. Tout fichier avec l'extension .m dans le chemin de Matlab peut etre appele de l'invite de commande. Un fichier .m est un fichier script si il ne commence pas par function. Les fichiers script sont executes ligne par ligne et les variables sont dans le workspace.

On cree un fichier script appele fit_polynome.m

edit fit_polynome.m

La premiere ligne devrait etre une description du fichier en commentaire:

%% fit_polynome: script pour ajuster un polynome a des donnees

Ensuite on definit les vecteur x et y.

x = 10*rand(34,1);      % 34 points aleatoires dans l'intervalle [0, 10];
y = 10*rand(34,1) + x;  % 34 points aleatoires correles avec x
n = 3;                  % interpolation avec un polynome de degre 3
p = monpolyfit(x,y,n);  % on trouve le polynome qui s'ajuste le mieux

En executant le script fit_polynome, on obtient le polynome p demande. Toutes les variables definies dans le script se retrouvent dans le workspace.

Graphiques

Matlab possede plusieurs commandes pour visualiser des donnees. Pour tracer des donnes en 2D, on utilise la commande plot. Rajoutons dans le fichier script fit_polynome les commandes pour visualiser les nuages de points (x,y) et le graphe du polynome p.

%% Visualisation 2D des donnees
figure(1);              % cree ou active la figure 1
clf;                    % clf efface la figure
plot(x,y,'s')           % trace les points avec des carres ('s' pour square)
hold on                 % les prochains appels a plot n'effaceront pas les autres
r = max(x)-min(x);      % l'etendue de x
a = min(x) - 0.1*r;     % limite inferieure sur x pour tracer le polyome
b = max(x) + 0.1*r;     % limite superieure sur x pour tracer le polyome
xp = linspace(a,b,100); % 100 points equirepartis sur [a,b]
plot(xp,polyval(p,xp))  % trace le graphe de p aux valeurs xp
xlabel('x')             % etiquette abscisse
ylabel('y')             % etiquette ordonnees
legend('donnees','polynome',2) % legende, les etiquettes sont dans
                               % l'ordre d'aparition des plots

Pour voir le resultat, on execute le script,

fit_polynome

Resoudre un systeme d'equations differentielles ordinaires

Matlab possede plusieurs solveurs EDO, dont une suite edo*. Le solveur multi-usage a utiliser par default est ode45 (prononce ode quatre cinq) un solveur Runge-Kutta a pas de temps adaptatif d'ordre 4 et 5. Les autres solveurs utiles sont ode23, un solveur d'ordre ode15s. Tous les solveurs ont la meme syntaxe:

[t,x] = ode45(@monedo,[t0 tfinal],CI);
sol   = ode45(@monedo,[t0 tfinal],CI);

La premiere syntaxe renvoie la solution x au temps t. La deuxieme syntaxe renvoie une structure sol qui contient la solution ainsi que des informations sur le solveur utilise. La fonction monedo est une fonction Matlab qui definit le systeme d'EDO. Le systeme est integre sur l'intervalle [t0 tfinal], avec conditions initiales CI. L'equation differentielle est de la forme

$\frac{dx}{dt}(t) = F(t,x).$

avec $x(t)$ un vecteur de taille $n$ . La fonction EDO doit suivre une structure bien definie, avec une premiere ligne

function dxdt = monedo(t,x)

La fonction prend en argument le temps t et l'etat du systeme x. La variable x est un vecteur colonne de taille n. La fonction retourne la derivee de x au temps t dans dxdt un vecteur colonne de taille n. Les conditions initiales sont un vecteur de la meme dimension que x. Comme premier exemple, on prend le systeme de Lotka-Volterra. La fonction lotkavolterra.m specifie le systeme suivant:

$\frac{dx}{dt} = x (a - b y),$$ $$\frac{dy}{dt} = -y (c - d x).$

Le systeme est resolu de la facon suivante:

CI = [30; 2];
tspan = [0 20];
sol = ode23(@lotkavolterra,tspan,CI);

On trace la solution a l'aide de deval, qui permet d'evaluer la solution sol a n'importe quel temps.

xint = linspace(0,20,200);
yint = deval(sol,xint);
figure(2); clf;
plot(xint,yint)
xlabel('temps'); ylabel('population');

Dans cet exemple, il n'est pas possible de facilement changer les parametres du systeme, qui sont code dans le fichier lotkavolterra.m. De plus, pour chaque resolution du systeme, il faut definir les conditions initiales et l'intervalle d'integration. Un fichier script peut aider, mais ne permet pas d'acceder aux parametres, qui sont variables locales a lotkavolterra. Un moyen de s'en sortir est d'utiliser une fonction imbriquee. Une fonction imbriquee est une sous-fonction, mais qui partage certaines variables avec la fonction principale dans laquelle elle se trouve.

A lieu de definir un script pour resoudre le systeme d'EDO, on definit une fonction run_lotka.m

function sol = run_lotka()
% RUN_LOTKA simulation du systeme de Lotka-Volterra

% parametres d'integration
CI = [30; 2];
tspan = [0 20];

% parametres du systeme
pars = [0.1 1 2 ];      % Les variables a, b, c, d sont partagees
b = 0.01;   % avec la fonction imbriquee lotka
c = 1;
d = 0.02;

% resolution du systeme pour 3 jeux de parametres differents
for i = 1:3
    a = pars(i); % on met a jour la valeur de a
    sol{i} = ode23(@lotka,tspan,CI);
end

% plot des 3 solutions
nx = 200;
xint = linspace(tspan(1),tspan(2),nx);
figure(2); clf;
newplot;            % met en place un nouveau plot
hold on             % qui ne s'effacera pas a chaque appel de plot
cols = 'kbr';       % couleurs k: noir, b: bleu, r:rouge
for i = 1:3
    yint = deval(sol{i},xint);
    plot(xint,yint,'Color',cols(i))
end
xlabel('temps'); ylabel('population');

     function dxdt = lotka(~,xx)  % fonction imbriquee
     % LOTKA systeme d'EDO

         x = xx(1); y = xx(2);

         dxdt = [ x * (a - b * y)
             -y * (c - d * x)];
     end                         % on termine d'abord la fonction imbriquee

end                             % on termine ensuite la fonction principale

De cette facon, on peut resoudre le systeme pour n'importe quelle condition initiale, intervalle de temps et parametres du systeme. Le fichier run_lotka.m peut etre lance:

run_lotka;

Codes

Les fichiers matlab peuvent etre telecharges ici: householder.m, Qadjb.m, backsubs.m, lotkavolterra.m, run_lotka.m

La démo enti�re peut �tre téléchargée ici: tutorielmatlab.m. Ce fichier est un fichier script, qui peut exécuté par la commande tutorielmatlab, mais vous pouvez aussi l'exécuter en mode démo avec la commande echodemo tutorielmatlab. Assurez-vous d'avoir tous les fichiers nécessaires dans votre répertoire courant.