Universite Aix-Marseille 1

Licence de mathematiques
Analyse numerique (2LMAT6B)
TP Matlab
Ce TP est une adaptation (avec pas mal d'ajouts) pour le public de licence de l'université de Marseille écrite par
B. Andreianov, R. Herbin , F. Hubert et N. Jullian, des pages de TP de l'université du New Hampshire écrites
par Kelly Black, Pour voir la version originale : http://www.math.unh.edu/~mathadm/tutorial/software/matlab/

La commande pour lancer matlab est matlab (!).

Pour sortir de matlab, il faut faire quit ou exit.

Attention il faut ABSOLUMENT sortir de matlab avant de fermer la fenètre unix ou vous l'avez appelé e
avant de faire le logout de votre terminal.

En effet, sinon le processus matlab n'est pas interrompu et vous bloquez un des 30 accès a matlab... Si vous
n'êtes pas sorti "proprement" de matlab, il faut vous relogger, récuperer le numéro du processus matlab par la
commande top, par exemple, et faire kill -9 n0 du processus.

On va utiliser matlab comme logiciel de programmation pour se simplifier la vie, puisque de nombreuses
procédures existent déja dans matlab. Le langage de base est proche des notations utilisées en algèbre linéaire et
analyse numérique, mais il faut connaître quelques extensions de syntaxe que nous allons introduire dans les
feuilles qui suivent. On commence par les notions les plus simples sur les vecteurs et les matrices pour arriver
jusqu'à la programmation...

Warning : il est possible (probable, certain...) qu'il comporte des fautes de frappe (inéxactitudes, erreurs,
bourdes colossales ...) que les auteurs vous seraient reconnaissants de leur signaler... [email protected]

Table des matieres:

1. Les vecteurs
2. Les matrices
3. Les opérations sur les vecteurs
4. Les boucles
5. Le graphisme
6. Les procédures
7. Les fonctions
8. Les matrices creuses

Les vecteurs
Commencons par la manipulation des vecteurs dans Matlab. Matlab est un logiciel de programmation qui rend
en particulier plus facile la manipulation des matrices et vecteurs. L'interface de Matlab est écrit dans un langag
qui ressemble énormement aux notations employées en algèbre linéaire et en analyse numérique, On commence
ici par donner les bases pour la manipulation des vecteurs.

Pour démarrer matlab, ouvrir une fenètre unix (ou linux) et taper la commande : matlab

Cette commande démarre le logiciel, qui attend ensuite que vous lui entriez une commande matlab. Dans le texte
qui suit, toute ligne démarrant par deux signes "strictement superieur a" ( >> ) est une ligne de commande
matlab ou sera écrite votre commande.

Presque toutes les commandes de base de Matlab s'articulent autour de l'utilisation des vecteurs. Pour simplifier
la création des vecteurs, on peut spécifier la première composante, l'increment et la dernière composante. Matlab
retrouvera automatiquement le nombre de composantes et leurs valeurs. Par exemple, pour créer un vecteur
dont les composantes sont 0, 2, 4, 6, and 8, on peut taper la ligne suivante :
>> [0,2,4,6,8]

ans =

0 2 4 6 8

ou bien de manière plus compacte (on va de 0 à 8 avec un pas ègal à 2) :

>> 0:2:8
ans =

0 2 4 6 8

Matlab garde en mémoire le dernier résultat. Dans l'exemple précédent, une variable "ans" a été créee
(pour "answer" = reponse en anglais...) . Pour obtenir la transposée du précédent vecteur, taper la commande
suivante (dans les notations anglo-saxonnes, le ' signifie la transposée):
>> ans'

ans =

0
2
4
6
8
Pour pouvoir garder en mémoire les vecteurs crées, il faut leur donner un nom. Par exemple on peut créer le
vecteur ligne v en tapant la commande suivante :
>> v = [0:2:8]

v =

0 2 4 6 8

>> v

v =

0 2 4 6 8

>> v;
>> v'
ans =

0
2
4
6
8
Remarquer que dans l'exemple précédent, lorsqu'on termine la ligne par un point virgule, le résultat n'est pas
affiché.
Ceci sera utile lorsqu'on travaillera sur des très gros systèmes d'équations.

Matlab vous permet de travailler sur certaines parties d'un vecteur. Par exemple si vous voulez seulement les tro
premières composantes du vecteur, vous utilisez la meme notation que pour la création du vecteur.

Pour afficher les trois premières composantes de v :

>> v(1:3)

ans =

0 2 4

Pour afficher les composantes 1 et 3 de v :

>> v(1:2:3)
ans =

0 4

ce qui donne le même résultat que :

>> v([1,3])
ans =

0 4

ou encore :
>> v(1:2:4)
ans =

0 4

Pour afficher le vecteur transposé du précèdent :

>> v(1:2:4)'

ans =

0
4

Vous pouvez aussi utiliser un incrément négatif pour définir un vecteur en partant de la dernière composante.
Ainsi au lieu de définir
>> v=[0,-2,-4,-6,-8]
v =

0 -2 -4 -6 -8

vous pouvez définir:

>> v=[0:-2:-8]

v =

0 -2 -4 -6 -8

Une fois que vous maitrisez les notations,vous pouvez effectuer des opérations, par exemple:
>> v(1:3)-v(2:4)

ans =

2 2 2

La commande clear permet d'effacer toutes les variables (attention, matlab ne demande pas de confirmation, et
une fois effacées, les variables sont irrécuperables...) :

>> clear
>> v
??? Undefined function or variable 'v'.

retour au sommaire

Les matrices
On donne ici les outils de base pour la définition et manipulation des matrices. On suppose connues les première
notions sur les vecteurs données au paragraphe précèdent.

On définit une matrice de manière semblable à un vecteur. Pour définir une matrice, on peut en effet la voir
comme une colonne de vecteurs lignes (attention à laisser les espaces...)
>> A = [ 1 2 3; 3 4 5; 6 7 8]

A =

1 2 3
3 4 5
6 7 8

On peut aussi la voir comme une ligne de vecteurs colonnes

>> B = [ [1 2 3]' [2 4 7]' [3 5 8]']

B =
 1 2 3
2 4 5
3 7 8

(Attention encore aux espaces.)

Si vous avez utilisé un paquet de variables pendant votre session matlab et que vous voulez en avoir la liste et la
taille, vous pouvez utiliser la commande whos pour en avoir la liste.

>> whos

Name Size Bytes Class

A 3x3 72 double array
B 3x3 72 double array

Grand total is 18 elements using 144 bytes

Comme on l'a déjà mentionné, les notations utilisées par Matlab sont proches de celles d'algèbre linéaire. La
multiplication matrice-vecteur est facile à faire, attention toutefois de respecter les tailles des matrices et des
vecteurs à multiplier, sinon, gare aux messages d'erreur !!. Il y a en particulier une contradiction entre la
définition habituelle des vecteurs en algèbre linéaire, qui sont des vecteurs colonnes, et la définition des vecteurs
par matlab, qui sont des vecteurs lignes. Définissons par exemple le vecteur v ainsi
>> v = [0:2:8]

v =

0 2 4 6 8

Pour matlab, v est un vecteur ligne (ce qui est raisonnable puisqu'on l'a écrit sous forme d'une ligne), i.e. une
matrice (1,3); donc si vous faites le produit A v, vous courez a la catastrophe:
>> A*v(1:3)
??? Error using ==> *
Inner matrix dimensions must agree.

Ce message d'erreur signifie que vous pouvez multiplier une matrice (m,p) par une matrice (n,q) à condition que
p = n (auquel cas vous obtenez une matrice (m,q). Il faut ici multiplier la matrice A par le vecteur colonne v', qui
peut être vu comme une matrice (3,1).
>> A*v(1:3)'

ans =

16
28
46

Habituez-vous a ce message d'erreur, une fois qu'on se met a multiplier des matrices et des vecteurs dans tous le
sens, il est facile d'oublier la taille des objets qu'on a définis, en particulier en raison du problème signalé sur la
definition des vecteurs...

On peut aussi travailler avec certaines parties des matrices, exactement comme avec les vecteurs. Une fois
encore, faites attention à ce que les tailles des objets manipulés soient compatibles.
>> A(1:2,3:4)
??? Index exceeds matrix dimensions.

>> A(1:2,2:3)

ans =

2 3
4 5

>> A(1:2,2:3)'

ans =

2 4
3 5

Une fois que vous savez créer et manipuler une matrice, vous pouvez effectuer des opérations standard dessus,
par exemple trouver son inverse. Exemple où ca marche bien, la matrice B définie plus haut :
>> inv(B)

ans =

-3.0000 5.0000 -2.0000

-1.0000 -1.0000 1.0000

2.0000 -1.0000 0

Attention toutefois quand vous calculez l'inverse d'une matrice : matlab ne sait identifier un 0 qu'à la precision de
la machine près. Il va donc vous inverser des matrices carrées qui ne sont pas forcément inversibles, mais il vous
avertit que son calcul risque de ne pas être terrible en vous donnant en particulier le conditionnement de la
matrice (calcule avec la precision machine).

Exemple ou matlab sait que votre matrice n'est pas inversible (et donc il vous le dit) :
>> C = [1 1 1; 2 1 2; 0 0 0]

C =

1 1 1

2 1 2

0 0 0

>> inv(C)

Warning: Matrix is singular to working precision.

ans =

Inf Inf Inf

Inf Inf Inf

Inf Inf Inf

(Inf = Infinity). Voyons maintenant un exemple où il ne sait pas que votre matrice n'est pas inversible, prenons la
matrice A définie plus haut et cherchons son inverse par matlab.
>> inv(A)

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 4.565062e-18

ans =

1.0e+15 *

-2.7022 4.5036 -1.8014

5.4043 -9.0072 3.6029
-2.7022 4.5036 -1.8014

Or si on calcule le déterminant de A :
>> det(A)

ans =

Vous pouvez vérifier ce calcul a la main : la matrice A est effectivement singulière.

Attention aussi au fait que matlab fait la distinction entre majuscules et minuscules, ce qui peut poser des
problèmes lorsqu'on commence à programmer des algorithmes un peu plus compliqués.
>> inv(a)
??? Undefined function or variable a.

Dans l'exemple précèdent, on avait défini auparavant la matrice A, mais on n'a jamais défini la variable a.

Matlab comporte des procédures qui permettent de trouver la solution d'un système d'équations linéaires. Par
exemple, si vous voulez trouver le vecteur x solution de Ax=b, une possibilité (certainement pas la plus
efficace...) est de calculer l'inverse de A et de multiplier le résultat par b. En fait cette procédure est lente et peu
précise : elle est à proscrire (d'ailleurs quand vous calculez 21/7, vous viendrait-il a l'idée de commencer par
calculer 1/7 et de le multiplier ensuite par 21 ?) Les méthodes de résolution de systèmes linéaires utilisées par
Matlab sont du type de celles que vous avez vues en cours (décomposition LU avec pivot ou décomposition
LDLt dans le cas des matrices symétriques)

On définit un vecteur colonne v :

>> v = [1 3 5]'

v =

1
3
5

et on cherche a calculer la solution du système B x = v où la matrice B est définie plus haut (on sait qu'elle est
inversible) :
>> x=B\v

x =

-1

Aucun problème... Par contre si on lui demande la même chose avec A :

>> y = A\v

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 4.565062e-18

y =

1.0e+15 *

1.8014
-3.6029
1.8014

Voyons ce que donne le produit B par x :

>> B*x

ans =

1
3
5

Ouf, on a bien Bx = v (ne pas oublier le * pour le produit dans matlab). On peut aussi résoudre des systèmes de
la forme x A = b mais attention alors aux dimensions de x, A et b, (il faut qu'elles soient compatibles...)
>> x1 = v'/B

x1 =

4.0000 -3.0000 1.0000

>> x1*B

ans =

1.0000 3.0000 5.0000

Attention : la commande A/B effectue l'opération A*inv(B) alors que la commande A\B effectue la commande
inv(A)*B.

Exercice 1 : Calculez (A/B)*B pour vérifier.

On peut aussi vouloir calculer les valeurs propres et les vecteurs propres d'une matrice. Il existe deux procédure
dans matlab, l'une calcule seulement les valeurs propres, l'autre calcule les valeurs propres et les vecteurs propre
Si vous ne vous souvenez plus quelle est la syntaxe d'utilisation de cette procédure, vous pouvez taper help eig
après le prompt de matlab (eig pour eigenvalue ou eigenvector (valeur propre ou vecteur propre en anglais) ou
consulter l'aide en ligne de matlab accessible en haut à droite de votre fenêtre.

>> eig(A)

ans =

14.0664
-1.0664
0.0000

>> [v,e] = eig(A)

v =

-0.2656 0.7444 -0.4082

-0.4912 0.1907 0.8165
-0.8295 -0.6399 -0.4082

e =

14.0664 0 0
0 -1.0664 0
0 0 0.0000

>> diag(e)

ans =

14.0664
-1.0664
0.0000

retour au sommaire

Les opérations sur les matrices et les vecteurs

Matlab permet de créer facilement les vecteurs et les matrices. La vraie puissance de Matlab réside dans la facilit
avec laquelle on peut manipuler les vecteurs et les matrices. Nous supposons ici connues les définitions et les
opérations de base de l'algèbre linéaire, ainsi que les deux premiers paragraphes de ce TP.

Nous illustrons d'abord les manipulations simples telle que addition, soustraction et multiplication matricielle.
Ensuite, nous verrons les opérations "composante par composante". Puis nous combinons ces opérations afin
d'illustrer qu'un effort très modique permet de définir des opérations relativement compliquées.

Notons qu'un vecteur n'est rien d'autre qu'une matrice et que tout réel est considéré comme une matrice de taille
1*1. La commande size vous permet d'obtenir la taille d'une matrice, la commande length la longueur d'un
vecteur.
>> v = [1 2 3]';size(v)
ans =
3 1

>> v = [1 2 3];size(v)
ans =
1 3

>> length(v)
ans =
3

>> a=1; size(a)

ans =
1 1

ADDITION, SOUSTRACTION

Les notations sont classiques pour l'algèbre linéaire (notez que les vecteurs doivent être définies comme
colonnes!). Définissons deux vecteurs et faisons leur différence :
>> v = [1 2 3]'

v =

1
2
3

ou bien

>> v = [1; 2; 3]

v =

1
2
3

ce qui est plus long à taper.

>> b = [2 4 6]'

b =

2
4
6

>> v+b

ans =

3
6
9

>> v-b

ans =
 -1
-2
-3

MULTIPLICATION

La multiplication des vecteurs et des matrices est possible si et seulement si les dimensions sont compatibles ;
dans le cas contraire un message d'erreur apparaît.
>> v*b
Error using ==> *
Inner matrix dimensions must agree.

>> v*b'

ans =

2 4 6
4 8 12
6 12 18

>> v'*b

ans =

28
Avec la dernière commande, vous avez effectué le produit scalaire des vecteurs v et b. Notez encore une fois qu
les vecteurs sont des COLONNES.

Souvent, on veut faire une opération avec chaque élément d'une matrice. Matlab permet d'effectuer ses
opérations très rapidement (beaucoup plus rapidement que si l'on utilise une boucle!). Par exemple, supposons
que l'on veut multiplier chaque élément du vecteur v par l'élément correspondant du vecteur b (c'est-à-dire
trouver le vecteur de composantes v(1)*b(1), v(2)*b(2), v(3)*b(3) ). Le symbole de l'opération en question est
noté ".*" dans Matlab.
>> v.*b

ans =

2
8
18

>> v./b

ans =

0.5000
0.5000
0.5000
En fait, on peut mettre le point devant toute opération pour dire à Matlab que vous voulez faire cette opération
composante par composante. Bien sur, cela ne change rien pour l'addition et la soustraction. Mais :
>> A=[1 2; 3 4], B=[1 -1; 2 2]

A =
 1 2
3 4

B =

1 -1
2 2

>> A*B, A.*B

ans =

5 3
11 5

ans =

1 -2
6 8

>> A/B, A./B

ans =

-0.5000 0.7500
-0.5000 1.7500

ans =

1.0000 -2.0000
1.5000 2.0000

CONCATENATION

On peut aussi facilement concaténer des vecteurs ou des matrices là encore, attention à la compatibilité des
dimensions :
>> [A,B]

ans =

1 2 1 -1
3 4 2 2

>> [A B]

ans =

1 2 1 -1
3 4 2 2

>> [A ;B]

ans =

1 2
 3 4
1 -1
2 2

Matlab sait appliquer des FONCTIONS non-linéaires à un vecteur ou une matrice composante par composante.
En fait, c'est ce qu'il fait par défaut.
>> sin(v)

ans =

0.8415
0.9093
0.1411

>> log(v)

ans =

0
0.6931
1.0986
>> x=[0:pi/2:2*pi]

x =

0 1.5708 3.1416 4.7124 6.2832

>> sin(x)

ans =

0 1.0000 0.0000 -1.0000 -0.0000

Remarquez que les FONCTIONS DES MATRICES sont définies différemment :
>> A=[0.5 0.2; -0.1 0.3]

A =

0.5000 0.2000
-0.1000 0.3000

>> EA=exp(A), EMA=expm(A)

EA =

1.6487 1.2214
0.9048 1.3499

EMA =

1.6333 0.2979
-0.1489 1.3354

La matrice EA est l'exponentielle de A "élément par élement", tandis que la matrice EMA est l'exponentielle
matricielle de A. En guise de vérification :
>> EMA-eye(2)-A-1/2*A^2-1/6*A^3-1/24*A^4
ans =

1.0e-03 *

0.1349 0.2012
-0.1006 -0.0663

où la commande eye(n) produit une matrice identité de taille n.

Certaines fonctions standard sont définies de la même façon :

>> log(EA), logm(EMA)

ans =

0.5000 0.2000
-0.1000 0.3000

ans =

0.5000 0.2000
-0.1000 0.3000

Matlab permet de travailler avec des fonctions "composante par composante". Grâce à cela, nous pouvons défin
rapidement et facilement de nouvelles opérations. Dans l'exemple suivant, nous définissons et manipulons un
vecteur "très long". (Remarquez que ";" à la fin d'une commande permet de ne pas afficher le résultat de son
execution.)
>> x = [0:0.1:100]

x =

Columns 1 through 7

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000

(...encore un paquet de lignes...)

Columns 995 through 1001

99.4000 99.5000 99.6000 99.7000 99.8000 99.9000 100.0000

>> y = sin(x).*x./(1+cos(x));
Dans ce cas, il est préférable d'afficher le résultat GRAPHIQUEMENT. Pour cela, on utilise la commande plot
de Matlab.
Dans cet exemple, on rencontre également la commande la plus utile de Matlab... c'est le "HELP".
>> plot(x,y)
>> plot(x,y,'rx')
>> help plot

>> help plot

PLOT Linear plot.

PLOT(X,Y) plots vector Y versus vector X. If X or Y is a matrix,
then the vector is plotted versus the rows or columns of the matrix,
whichever line up. If X is a scalar and Y is a vector, length(Y)
 disconnected points are plotted.

PLOT(Y) plots the columns of Y versus their index.

If Y is complex, PLOT(Y) is equivalent to PLOT(real(Y),imag(Y)).
In all other uses of PLOT, the imaginary part is ignored.

Various line types, plot symbols and colors may be obtained with
PLOT(X,Y,S) where S is a character string made from one element
from any or all the following 3 colunms:

y yellow . point - solid

m magenta o circle : dotted
c cyan x x-mark -. dashdot
r red + plus -- dashed
g green * star
b blue s square
w white d diamond
k black v triangle (down)
^ triangle (up)
< triangle (left)
> triangle (right)
p pentagram
h hexagram

For example, PLOT(X,Y,'c+:') plots a cyan dotted line with a plus

at each data point; PLOT(X,Y,'bd') plots blue diamond at each data
point but does not draw any line.

PLOT(X1,Y1,S1,X2,Y2,S2,X3,Y3,S3,...) combines the plots defined by

the (X,Y,S) triples, where the X's and Y's are vectors or matrices
and the S's are strings.

For example, PLOT(X,Y,'y-',X,Y,'go') plots the data twice, with a

solid yellow line interpolating green circles at the data points.

The PLOT command, if no color is specified, makes automatic use of

the colors specified by the axes ColorOrder property. The default
ColorOrder is listed in the table above for color systems where the
default is yellow for one line, and for multiple lines, to cycle
through the first six colors in the table. For monochrome systems,
PLOT cycles over the axes LineStyleOrder property.

PLOT returns a column vector of handles to LINE objects, one

handle per line.

The X,Y pairs, or X,Y,S triples, can be followed by

parameter/value pairs to specify additional properties
of the lines.

See also SEMILOGX, SEMILOGY, LOGLOG, GRID, CLF, CLC, TITLE,

XLABEL, YLABEL, AXIS, AXES, HOLD, COLORDEF, LEGEND, and SUBPLOT.

>> plot(x,y,'y',x,y,'go')
>> plot(x,y,'y',x,y,'go',x,exp(x+1),'m--')

Exercice 2 : Tracer la courbe de la fonction x -> sin(x) entre 0 et 100 avec le pas 0.1.
EXEMPLE
Supposons que vous voulez approcher la dérivée d'une fonction donnée par les différences finies. (Notez que, à
la différence de Maple, Matlab ne calcule pas la dérivée symboliquement; il est utilisé pour des calculs
numériques matriciels, éventuellement très volumineux). Dès que vous avez les points de discretisation et les
valeurs de la fonction en ces points, il est facile de récupérer le vecteur de différences finies (par exemple,
décentré à droite) :
>> diff = (y(2:1001)-y(1:1000))./(x(2:1001)-x(1:1000));
>> whos
Name Size Bytes Class

ans 1x1 8 double array

diff 1x1000 8000 double array
x 1x1001 8008 double array
y 1x1001 8008 double array

Grand total is 3003 elements using 24024 bytes

Comparons avec la dérivée exacte (regardez préalablement que doit faire la commande max,
en tapant help max) :
>> diff-cos(x(1:1000));
>> max(ans)

ans =

0.0500

Les différences centrées sont plus précises :

>> diff1= (y(3:1001)-y(1:999))./(x(3:1001)-x(1:999));
>> diff1-cos(x(2:1000));
>> max(ans)

ans =

0.0017

Voici les commandes qui permettent de discrétiser la dérivée première et la dérivée seconde de la fonction
sin(x)^4 sur la partition de [0,100] avec le pas 0.1, et mettre en evidence l'erreur maximale commise.
>> z=sin(x).^4;
>> diff1= (z(3:1001)-z(1:999))./0.2;
>> diff1 - 4*(sin(x(2:1000)).^3).*cos(x(2:1000));
>> max(ans)

ans =

0.0181

>> diffdoubles= (z(3:1001)-2*z(2:1000)+z(1:999))./(0.1)^2;

>> diffdoubles - 4*( 3*(sin(x(2:1000)).^2).*(cos(x(2:1000)).^2) - sin(x(2:1000)).^4);
>> max(ans)

ans =

0.0332

Exercice 3:

1. Tracer la fonction x*atan(x) sur l'intervalle [-5,5] avec la précision "visuelle".

2. Tracer sa dérivée sur le même intervalle.

3. Définir le vecteur de discrétisation de [0,1] avec le pas alterné 0.01/0.03 : les noeuds seront en 0, 0.01,
0.04, 0.05, 0.08, 0.09, 0.12 etc...
Réponse possible:
>> f=[0:0.04:1];
>> g=zeros(1,51);
>> g(1:2:51)=f(1:26);
>> g(2:2:50)=f(1:25)+0.01;
>> g

4. Proposer une formule pour approcher la dérivée seconde d'une fonction à partir des valeurs de la
fonction sur cette partition (utiliser le développement de Taylor pour ajuster les coefficients).
Réponse possible :
>> diffdoubles=2*( (g(2:50)-g(1:49)).*y(3:51) - (g(3:51)-g(1:49)).*y(2:50) + (g(3:51)-
g(2:50)).*y(1:49) )./( 0.01*0.03*(g(3:51)-g(1:49)) );

5. Comparer graphiquement la dérivée seconde numérique de la fonction sin(2x) avec le graphique exac
.
Réponse :
>> f=[0:0.04:1];
>> g(1:2:51)=f(1:26);
>> g(2:2:50)=f(1:25)+0.01;
>> y=sin(2*g);
>> diffdoubles=2*( (g(2:50)-g(1:49)).*y(3:51) - (g(3:51)-g(1:49)).*y(2:50) + (g(3:51)-
g(2:50)).*y(1:49) )./( 0.01*0.03*(g(3:51)-g(1:49)) );
>> plot(g(2:50),diffdoubles,g(2:50),-4*sin(2*g(2:50)))

retour au sommaire

Les boucles
Ici nous montrons comment utiliser les boucles for et while. Nous commençons par les boucles for avec
l'application à la méthode d'Euler pour approximation des équations différentielles ordinaires. Nous donnons
ensuite des exemples avec while.

Les boucles forpermettent de répéter une sequence de commandes. Toutes les structures "boucle" de Matlab
commencent par un mot clé (for, while) et finissent avec le mot end. Pour avoir le syntaxe complet, tapez help
for.

Le boucle suivante tourne 4 fois et imprime ce qu'elle fait (pour la rendre muette, il faut mettre ; après j^2) :
>> for j=1:4, j^2, end

ans =

ans =

4
ans =

ans =

16

On peut aussi définir un vecteur, puis modifier chaque composante à l'aide d'une boucle :
>> v=[0:0.5:pi], M=length(v)

v =

0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000

M =

>> for i=1:M, v(i)=v(i)*i; end, v

v =

0 1 3 6 10 15 21

C'est un simple exemple d'utilisation de boucles en Matlab. Aussi simple qu'inutile...

DO NOT DO THIS IN PRACTICE !!!!

L'execution de boucles est BEAUCOUP plus lente que la modification de vecteurs à l'aide des opérations
"composante par composante".

Illustrons cela avec la montre automatique (faites help tic pour voir comment elle marche) :
>> v=[0:0.01:100*pi]; M=length(v)

M =

31416

>> tic, w=v.*[1:M]; TimeDirect=toc

TimeDirect =

0.0159

>> tic, for i=1:M, v(i)=v(i)*i; end, TimeBoucle=toc

TimeBoucle =

0.7330

>> TimeBoucle/TimeDirect
ans =

46.2328

>> max(v-w)

ans =

Le temps d'execution d'une boucle est plus que 40 fois plus long !!! Et pourtant on calcule bien la même chose,
car la dernière commande montre que toutes les composantes de (v-w) valent 0. En utilisant, bien entendu, une
opération "composante par composante" en non pas un boucle !

Toutefois, les boucles sont plus faciles à programmer lorsque il s'agit d'une opération récursive. Dans l'exemple
suivant, nous commencons par la seconde ligne de la matrice et retranchons la ligne précèdente, puis descendons
une ligne plus bas et exécutons la même opération, et ainsi de suite.
>> A=rand(4)

A =

0.0196 0.5028 0.1897 0.5417

0.6813 0.7095 0.1934 0.1509
0.3795 0.4289 0.6822 0.6979
0.8318 0.3046 0.3028 0.3784

>> for i=2:4, A(i,:)=A(i,:)-A(i-1,:), end

A =

0.0196 0.5028 0.1897 0.5417

0.6616 0.2067 0.0038 -0.3908
0.3795 0.4289 0.6822 0.6979
0.8318 0.3046 0.3028 0.3784

A =

0.0196 0.5028 0.1897 0.5417

0.6616 0.2067 0.0038 -0.3908
-0.2822 0.2222 0.6784 1.0887
0.8318 0.3046 0.3028 0.3784

A =

0.0196 0.5028 0.1897 0.5417

0.6616 0.2067 0.0038 -0.3908
-0.2822 0.2222 0.6784 1.0887
1.1140 0.0824 -0.3757 -0.7103

La commande rand(n) crée une matrice dont les coefficients sont choisit de façon aléatoir
entre 0 et 1.
Exemple 1 : le pivot de Gauss.

L'algorithme de la méthode de Gauss (sans permutation) peut être programmé en deux boucles et une seule
opération :
>> B=rand(5), A=B;

B =

0.6932 0.1988 0.4199 0.3678 0.5692

0.6501 0.6252 0.7537 0.6208 0.6318
0.9830 0.7334 0.7939 0.7313 0.2344
0.5527 0.3759 0.9200 0.1939 0.5488
0.4001 0.0099 0.8447 0.9048 0.9316

>> for j=2:5, for i=j:5, A(i,:) = A(i,:) - A(j-1,:)*A(i,j-1)/A(j-1,j-1); end, end, A

A =

0.9830 0.7334 0.7939 0.7313 0.2344

0 -0.3184 -0.1400 -0.1479 0.4039
0 0 0.6485 0.7413 0.4700
0 0 -0.0000 -0.7600 0.0159
0 0 0.0000 0 0.5311
Est-ce le bon résultat ? Souvenez-vous que la matrice triangulaire obtenue par la méthode de Gauss est
exactement la matrice U dans le l'algorithme LU. Matlab sait faire la décomposition LU (voir help lu pour son
utilisation) :
>> [L,U]=lu(B); U

U =

0.9830 0.7334 0.7939 0.7313 0.2344

0 -0.3184 -0.1400 -0.1479 0.4039
0 0 0.6485 0.7413 0.4700
0 0 0 -0.7600 0.0159
0 0 0 0 0.5311

En fait, Matlab optimise le choix de pivot en réarrangant les lignes et les colonnes de la matrice à l'aide de la
matrice P mentionnée dans help lu; nous avons triché pour avoir une matrice avec P=eye(5), en modifiant la
"vraie" matrice que random nous donne au départ.

Exercice 4 : Trouvez comment on a fait.

Exemple 2 : approximation des EDO

Construisons une solution approchée du problème de Cauchy y'=x^2-y^2, y(0)=1, en utilisant la méthode d'Eule
explicite. On définit le pas de discrétisation h, puis on trouve les points de la partition {x1,x2,x3,...} de
l'intervalle [0,2] avec le pas uniforme égal à h, puis on trouve les valeurs yi qui approchent y(xi). Les valeurs
xi et yi sont stockées en vecteurs x et y .
>> h = 0.1;
>> x = [0:h:2];
>> y = 0*x;
>> y(1) = 1;
>> for i=2:length(x),
y(i) = y(i-1) + h*(x(i-1)^2 - y(i-1)^2);
end
>> plot(x,y)
>> plot(x,y,'go')
>> plot(x,y,'go',x,y)
Si vous n'aimez pas les boucles for, vous pouvez utiliser while à la place. Le boucle while répète une séquence
de commandes tant qu'une certaine condition est satisfaite. Dans cet exemple le problème de Cauchy y'=x-|y|,
y(0)=1 est approché par la méthode d'Euler améliorée :

>> h = 0.005;
>> x = [0:h:4];
>> y = 0*x;
>> y(1) = 1;
>> i = 1;
>> while(i<length(x))
z=y(i) + h* (x(i)-(abs(y(i))));
y(i+1) = y(i) + h/2*( (x(i)-(abs(y(i))))+ (x(i+1)-(abs(z))) );
i = i + 1;
end
>> plot(x,y,'go')
>> plot(x,y)

Exercice 5 :

1. Définir en une ligne le vecteur à composantes sin(1), sin(4), sin(7), sin(10), ..., sin(10000).
Faire la même chose en utilisant un boucle. Comparer les temps d'execution.

2. Considérons les deux transformations d'un vecteur v :

>> for i=2:length(v), v(i)=v(i)-v(i-1); end
>> for i=1:length(v)-1, v(i)=v(i)-v(i+1); end
Une des deux peut être programmée avec une grande économie de temps d'execution.
Laquelle? Comment?
Réécrire l'autre transformation en utilisant un boucle "while".

3. Considérons la transformation suivante :

>> x=[1:1000].^2;
>> for i=2:1000, x(i)=x(i)+x(i-1); end ?
C'est une vraie récurrence, donc l'usage d'une boucle semble justifié.

(a) Expliquer alors le résultat suivant :

>> x=[1:1000].^2;
>> tic, y=[1:1000].*([1:1000]+1).*(2*[1:1000]+1)/6; Time1=toc;
>> tic, for i=2:1000, x(i)=x(i)+x(i-1); end, Time2=toc;
>> x(1:7),y(1:7),Time2/Time1

ans =

1 5 14 30 55 91 140

ans =

1 5 14 30 55 91 140

ans =

12.3685

(b) Améliorer l'algorithme suivant :

 >> x=[1:1000].;
>> for i=2:1000, x(i)=x(i)+x(i-1); end.

4. Montrer que toute boucle "for" peut être réécrite comme un boucle "while".
Proposer un exemple de boucle "while" qui ne peut pas être réécrite comme
un boucle "for". (Suggestion : la méthode de Newton vue en Maple...)

retour au sommaire

Le graphisme
Dans cette partie, on introduit les opérations élémentaires nécessaires à la création de graphiques. Afin de
montrer l'utilisation de la commande plot, une approximation par la methode d'Euler de l'équation differentielle
y'= 1/y avec y(0)=1 est calculée, puis tracée. On donne le pas de discrétisation h = 1/16, puis on exécute la
méthode d'Euler. Enfin, on saisit la solution exacte afin d'effectuer une comparaison avec la solution approchée.
(Cet exemple est tiré de la partie qui précède sur les boucles.)
>> h = 1/16;
>> x = 0:h:1;
>> y = 0*x;
>> size(y)

ans =

1 17

>> max(size(y))

ans =

17

>> y(1) = 1;
>> for i=2:max(size(y)),
y(i) = y(i-1) + h/y(i-1);
end
>> true = sqrt(2*x+1);
Maintenant, on a l'approximation et la solution exacte. Pour les comparer, l'approximation est tracée en vert (g
pour "green") avec des 'o', alors que la solution exacte est affichée avec les paramètres par défaut). C'est la
commande plot qui produit les graphiques sous matlab. Elle interprète une suite de triplets (ou couple, la
troisième composante étant facultative) formées de deux vecteurs contenant respectivement les abscisses et
ordonnées des points, et enfin des options d'affichage. Pour obtenir la liste de ces options (si on desire modifier
l'affichage par défaut), il suffit de les demander par l'intermédiaire de la commande help plot.
>>plot(x,y,'go',x,true)
Bien, mais il serait aussi intéressant de tracer l'erreur :
>> plot(x,abs(true-y),'mx')
Pour afficher plusieurs graphiques dans une même figure, on utilisera la commande subplot. Matlab peut
interpréter une fenètre comme un tableau de graphiques. Ici on désire avoir une ligne et deux colonnes, pour
placer dans la première case (première au sens naturel de la lecture) la courbe représentant la fonction exacte et
approchée, alors que le second emplacement recevra la représentation de l'erreur.
>> subplot(1,2,1);
>> plot(x,y,'go',x,true)
>> subplot(1,2,2);
>> plot(x,abs(true-y),'mx')
Figure 1. Les deux graphiques de la première approximation

Recommencons au début. On calcule un nouvelle approximation en ayant divisé en deux le pas de discrétisation.
Mais avant, on efface et on réinitialise la figure en utilisant la commande clf. (Remarquez qu'on utilise de
nouveaux vecteurs x1 et y1.)
>> clf
>> h = h/2;
>> x1 = 0:h:1;
>> y1 = 0*x1;
>> y1(1) = 1;
>> for i=2:max(size(y1)),
y1(i) = y1(i-1) + h/y1(i-1);
end
>> true1 = sqrt(2*x1+1);
La nouvelle approximation est tracée, mais attention! Comme matlab nous le fait remarquer, les vecteurs passés
la commande plot doivent avoir la même taille.
>> plot(x,y1,'go',x,true1)
??? Error using ==> plot
Vectors must be the same lengths.
On peut donner des titres aux axes et aux graphiques. L'exemple suivant met en valeur l'utilisation de la
commande subplot pour l'activation d'un graphique particulier au sein d'un tableau de graphiques.
>> plot(x1,y1,'go',x1,true1)
>> plot(x1,abs(true1-y1),'mx')
>> subplot(1,2,1);
>> plot(x,abs(true-y),'mx')
>> subplot(1,2,2);
>> plot(x1,abs(true1-y1),'mx')
>> title('Error for h=1/32')
>> xlabel('x');
>> ylabel('|Error|');
>> subplot(1,2,1);
>> xlabel('x');
>> ylabel('|Error|');
>> title('Error for h=1/16')
Figure 2. Les erreurs pour deux approximations

Enfin, si on désire faire une impression ou une sauvegarde dans un fichier (sans utiliser les menus déroulants des
figures), on peut utiliser la commande print. L'exemple suivant crée un fichier postscript erreur.ps, il ne restera
plus qu'a l'imprimer depuis un prompt UNIX en utilisant la commande lp.
>> print -dps erreur.ps
retour au sommaire

Les procédures
Dans cette partie, on va introduire les opérations de base nécessaires à la creation d'une procédure. C'est un
fichier dans lequel une succession de commandes (interprétables par matlab) va être inscrite. On pourra ainsi
exécuter cette liste de commandes sans avoir à les retaper, juste en les appelant par le nom du fichier. Par
exemple, on va écrire une procédure effectuant une approximation d'une équation différentielle par la méthode
d'Euler.

On peut utiliser l'éditeur de texte proposé par matlab en cliquant en haut à droite de la fenêtre de commandes
matlab sur open file ou encore utiliser n'importe quel éditeur proposé par unix. Pour ouvrir un éditeur de texte
sous unix, si vous n'en avez encore jamais utilisé, vous pouvez suivre le guide :

Bien, maintenant il reste à écrire la procedure et à la sauvegarder dans le répertoire courant sous un nom suivi d
l'extension.m par exemple simpleEuler.m. Il est extrèmement important de commenter vos programmes, de
manière à pouvoir les relire rapidement pour savoir ce qu'ils font.

% fichier: simpleEuler.m
% ce fichier matlab calcule l'approximation de
%
% dy/dx = 1/y
% y(0) = starty
%
%
% Pour executer cette procedure il faut preciser :
% h : la taille du pas
% starty : la valeur initiale
%
% La procedure cree trois vecteurs.
%
% x contient les points de discretisation,
% il commence a xo et son pas est h
%
% y contient l'approximation de l'equation differentielle
% de condition initiale starty aux points x
%
%
% true contient la solution exacte de l'equation differentielle
%
% si vous voulez completer l'aide, il suffit de rajouter des
% signes de pourcentage et de completer les commentaires.
%
x = [0:h:1];

y = 0*x;
y(1) = starty;

for i=2:max(size(y)),
y(i) = y(i-1) + h/y(i-1);
end

true = sqrt(2*x+1);

Une fois la sauvegarde effectuée il suffit de revenir a la fenètre matlab et taper le nom (sans l'extention .m) de la
procédure (c'est-à-dire : simpleEuler).
>> clear
>> simpleEuler
??? Undefined function or variable h.

Error in ==> /home/utilisateur/simpleEuler.m

On line 29 ==> x = [0:h:1];
Ici, on avait pris soin d'effacer toutes les variables alors connues de matlab (à l'aide de la commande clear) afin
d'observer son comportement. Le message est clair, matlab ne sait pas à quoi correspond "h" et il lui a été
demandé de travailler avec ce fameux "h" à la ligne 29 du programme simpleEuler.m, situé dans le répertoire
/home/utilisateur/ (on suppose que c'est notre répertoire de travail actuel). Recommençons maintenant, en ayant
précisé les variables nécessaires à notre procédure.
>> h = 1/16;
>> starty = 1;
>> simpleEuler
>> whos
Name Size Bytes Class

ans 1x140 280 char array

h 1x1 8 double array
i 1x1 8 double array
starty 1x1 8 double array
true 1x17 136 double array
x 1x17 136 double array
y 1x17 136 double array

Grand total is 194 elements using 712 bytes

>> plot(x,y,'rx',x,true)
Une fois les variables en mémoire, notre procédure peut fonctionner. L'utilisation de la commande whos nous
permet de constater que les variables x,y et starty ont bien été créees. Mais on remarque aussi la création de la
variable i, utilisée pendant la procédure. Ainsi, une procédure agit sur les variables locales (celles qui ont été
créees en ligne) et en crée d'autres si elles n'existent pas encore. C'est un inconvénient, en effet, si on avait
sauvegardé un calcul dans i avant l'execution de la procédure, ce calcul aurait été perdu. En fait une procédure
execute la liste des commandes comme si elles avaient été tapées en ligne.

Si on désire de nouveau executer les commandes contenues dans simpleEuler mais avec un nouveau pas de
discrétisation, il suffit de changer celui-ci. Mais attention, la procédure va sauvegarder ses nouveaux calculs dan
les mêmes variables que tout à l'heure. Si on désire les conserver, il faut prendre soin de les sauvegarder sous de
noms qui ne seront pas utilisés durant celle-ci.
>> x0 = x;
>> y0 = y;
>> true0 = true;
>> h = h/2;
>> simpleEuler
>> whos
Name Size Bytes Class

ans 1x140 280 char array

h 1x1 8 double array
i 1x1 8 double array
starty 1x1 8 double array
true 1x33 264 double array
true0 1x17 136 double array
x 1x33 264 double array
x0 1x17 136 double array
y 1x33 264 double array
y0 1x17 136 double array

Grand total is 293 elements using 1504 bytes

>> plot(x0,abs(true0-y0),'gx',x,abs(true-y),'yo');

On dispose maintenant des deux approximations, sauvegardées dans x0 et y0 pour le pas de 1/16, et x et y pour
le pas de 1/32.

retour au sommaire

Les fonctions
Lorsqu'on désire répeter une suite de commandes, mais qu'on veut également pouvoir les appliquer à des
variables dont le nom changera, sans risquer d'effacer des données sauvegardées dans une variable dont le nom
est utilisé au cours de cette liste de commandes, on utilise une fonction.

C'est comme une procédure, sauf qu'on va preciser les variables entrantes et sortantes. Ces variables sont muette
; le choix de leurs noms n'importe qu'au sein de la fonction. De plus, toutes les autres variables utilisées ne sont
définies que localement, leur utilisation ne modifiera pas des variables existantes, elle ne sont gardées en mémoir
que durant l'execution de la fonction.

On va prendre pour exemple une fonction effectuant une élimination de Gauss. On désire pouvoir passer en
argument la matrice A et le vecteur b, pour obtenir le vecteur x résolvant le système Ax=b. La première ligne de
la fonction doit préciser quelles sont les variables entrantes et sortantes, ainsi que leurs noms. Dans notre cas,
c'est ce que fait la ligne :
function [x] = gaussElim(A,b)
Pour avoir plus de variables en entrée ou sortie, il suffit de les rajouter en les séparant par une virgule.

Il ne reste plus qu'à copier le texte suivant dans un fichier que l'on sauvegardera sous gaussElim.m.

function [x] = gaussElim(A,b)

% Fichier gaussElim.m
% Cette fonction effectue l'elimination de Gauss
% sur le systeme donne
% i.e., soient A et b elle peut trouver x,
% Ax = b
%
% Pour executer cette fonction, il faut specifier :
% A - la matrice du membre de droite
% b - le vecteur du membre de gauche
%
% La fonction retourne x par elimination de Gauss.
% ex: [x] = gaussElim(A,b)
%

N = max(size(A));

% Elimination de Gauss

for j=2:N,
for i=j:N,
m = A(i,j-1)/A(j-1,j-1);
A(i,:) = A(i,:) - A(j-1,:)*m;
b(i) = b(i) - m*b(j-1);
end
end

% algorithme de remontee
x = zeros(N,1);
x(N) = b(N)/A(N,N);

for j=N-1:-1:1,
x(j) = (b(j)-A(j,j+1:N)*x(j+1:N))/A(j,j);
end

Maintenant que notre fonction existe, créons une matrice et un vecteur (de tailles cohérentes) afin de résoudre le
système. Et pour bien vérifier que les variables utilisées dans la fonction et celles déjà existantes, sont bien
indépendantes, choisissons d'autres noms.
>> M = [1 2 3 6; 4 3 2 3; 9 9 1 -2; 4 2 2 1]

M =

1 2 3 6
4 3 2 3
9 9 1 -2
4 2 2 1

>> v = [1 2 1 4]'

v =

1
2
1
4

>> [res] = gaussElim(M,v)

res =
 0.6809
-0.8936
1.8085
-0.5532

>> whos
Name Size Bytes Class

M 4x4 128 double array

v 4x1 32 double array
res 4x1 32 double array

Grand total is 24 elements using 192 bytes

Notre système a été résolu, et on remarque bien que les variables utilisées au cours de la fonction n'ont pas laissé
de "trace" (pas de i,j,N, etc ...).

Si on désire, on peut aussi demander à une fonction d'appeler une autre fonction. Ca pourrait être le cas de la
sélection de la méthode de résolution parmi vos algorithmes préferés lors de la résolution d'un problème. Ici, on
se sert de cette fonctionnalité pour la sélection de la fonction f pour la résolution de l'équation différentielle
y'=f(x,y) par la méthode d'Euler.

Comme d'habitude, il reste à copier le programme et à le sauvergarder. Le nom de l'enregistrement sera bien sur
eulerApprox.m.

function [x,y] = eulerApprox(startx,h,endx,starty,func)

% Fichier: eulerApprox.m
% Cette fonction calcule l'approximation de l'equation
% differentielle donnee par
% y' = func(x,y)
% y(startx) = starty
%
% Pour executer cette fonction vous devez specifier
% startx : la valeur initiale pour x
% h : le pas de discretisation
% endx : la valeur finale pour x
% starty : la condition initiale
% func : nom de fonction pour le calcul du second
% membre de l'equation differentielle.
% Doit etre une chaine de caracteres
%
% ex: [x,y] = eulerApprox(0,1,1/16,1,'f');
% Retourne l'approximation de l'equation differentielle
% ou x va de 0 a 1 avec un pas de 1/16.
% La condition initiale est 1, et le membre de gauche
% est calcule par la fonction donnee par f.m.
%
% La fonction cree deux vecteurs :
% - les points de discretisation : x
% - l'approximation de l'EDO aux points de discretisation : y
%

x = [startx:h:endx];

y = 0*x;
y(1) = starty;

for i=2:max(size(y)),
y(i) = y(i-1) + h*feval(func,x(i-1),y(i-1));
end

Dans cet exemple, on approche l'équation différentielle y'=1/y. Pour ce faire, on crée la fonction, qu'on enregistr
sous f.m.

function [f] = f(x,y)

% Evaluation du membre de gauche de l'equation differentielle

f = 1/y;

Voila, qui est fait. Maintenant on passe a l'utilisation. Au fait, quand on met des commentaires dans nos
fonctions, ils peuvent être visionnés comme aide en ligne avec la commande help. On lit alors que la fonction
doit être définie comme une chaine de caractères. C'est ce qu'on va faire en passant comme argument 'f' pour
indiquer l'utilisation de la fonction f.m (la commande interne feval va permettre d'appeler f.m).
>> help eulerApprox

Fichier: eulerApprox.m
Cette fonction calcule l'approximation de l'equation
differentielle donnee par
y' = func(x,y)
y(startx) = starty

Pour executer cette fonction vous devez specifier

startx : la valeur initiale pour x
h : le pas de discretisation
endx : la valeur finale pour x
starty : la condition initiale
func : nom de fonction pour le calcul du second
membre de l'equation differentielle.
Doit etre une chaine de caracteres

ex: [x,y] = eulerApprox(0,1,1/16,1,'f');

Retourne l'approximation de l'equation differentielle
ou x va de 0 a 1 avec un pas de 1/16.
La condition initiale est 1, et le membre de gauche
est calcule par la fonction donnee par f.m.

La fonction cree deux vecteurs :

- les points de discretisation : x
- l'approximation de l'EDO aux points de discretisation : y

>> [x,y] = eulerApprox(0,1/16,1,1,'f');

>> plot(x,y)
L'execution de la fonction produit deux vecteurs et les enregistre sous les noms specifiés (ici : x et y).

retour au sommaire
Les matrices creuses
La modélisation de problèmes physiques conduit à résoudre de très grands systèmes. Mais la matrice associée à
ce système possède généralement un grand nombre de zéros. Afin d'économiser la place mémoire nécessaire au
stockage de ces matrices, dites creuses (sparse en anglais), une série de commandes leur a été dediée. En fait ell
consiste en quelque sorte à ne stocker que les termes non nuls, l'emplacement de ces termes et la taille de la
matrice.

Le but de ces commandes est aussi d'automatiser la création de la matrice. Il n'est pas rare que le nombre de
termes se compte en centaines. Pour commencer voyons quelques matrices pleines qui pourront nous servir par
la suite.
>> ones(3,4)

ans =

1 1 1 1
1 1 1 1
1 1 1 1

>> zeros(2,5)

ans =

0 0 0 0 0
0 0 0 0 0

>> eye(3,5)

ans =

1 0 0 0 0
0 1 0 0 0
0 0 1 0 0

Ces exemples parlent d'eux-mêmes, passons maintenant aux matrices creuses. La commande est sparse. Elle
permet entre autres de transformer le mode de stockage d'une matrice pleine en celui d'une matrice creuse.
>> clear
>> A=eye(3,5), B=sparse(A)

A =

1 0 0 0 0
0 1 0 0 0
0 0 1 0 0

B =

(1,1) 1
(2,2) 1
(3,3) 1

>> whos
Name Size Bytes Class

A 3x5 120 double array

 B 3x5 60 sparse array

Grand total is 18 elements using 180 bytes

On constate bien que la matrice B prend moins de place en mémoire. Mais elle est moins lisible, sans l'utilisation
de la commande whos. L'affichage ne nous permet pas de déterminer sa taille. La matrice B est en fait composée
de coordonnées et de valeurs. C'est un autre mode d'utilisation de la commande sparse. On donne trois vecteurs
contenant respectivement des numéros de lignes, de colonnes et enfin des valeurs. Bien sur, ces vecteurs doiven
avoir la même taille.
>> C=sparse([1 2 3],[1 2 3], [1 1 1])

C =

(1,1) 1
(2,2) 1
(3,3) 1

>> whos B C
Name Size Bytes Class

B 3x5 60 sparse array

C 3x3 52 sparse array

Grand total is 6 elements using 112 bytes

Les matrices B et C retournent les mêmes lignes, mais elles n'ont pas la même taille. Lorsque les matrices sont d
taille "raisonnable", on peut utiliser la commande full pour les visualiser.
>> full(B) , full(C)

ans =

1 0 0 0 0
0 1 0 0 0
0 0 1 0 0

ans =

1 0 0
0 1 0
0 0 1

Bon, c'est bien mais il y a encore plus simple pour la création de B et C.

>>B2=speye(3,5) , C2=speye(3,3) , D=speye(3)

B2 =

(1,1) 1
(2,2) 1
(3,3) 1

C2 =

(1,1) 1
(2,2) 1
(3,3) 1
D =

(1,1) 1
(2,2) 1
(3,3) 1

En fait, speye est la version creuse de eye. On remarque aussi au passage que la donnée d'un seul argument
produit une matrice carrée. Il faut faire attention, pour matlab, certains arguments sont facultatifs mais peuvent
exister par défaut. Par exemple, la commande sparse demande le nombre de lignes et de colonnes en quatrième e
cinquième arguments.
>> A=sparse([1:3],[1:2:6],[1:3:9]) , B=sparse([1:3],[1:2:6],[1:3:9],4,10)

A =

(1,1) 1
(2,3) 4
(3,5) 7

B =

(1,1) 1
(2,3) 4
(3,5) 7

>> full(A) , full(B)

ans =

1 0 0 0 0
0 0 4 0 0
0 0 0 0 7

ans =

1 0 0 0 0 0 0 0 0 0
0 0 4 0 0 0 0 0 0 0
0 0 0 0 7 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0

Si on ne précise pas la taille, matlab prend la plus petite permettant à la matrice de contenir tous ses coefficients
non nuls. On va profiter de la forme de saisie ci-dessus pour augmenter la taille de la matrice.
>> n=3, B=sparse([1:n],[1:2:2*n],[1:3:3*n],n+1,3*n+1)

n =

B =

(1,1) 1
(2,3) 4
(3,5) 7

>> n=30 , B=sparse([1:n],[1:2:2*n],[1:3:3*n],n+1,3*n+1);

n =
 30

Cette fois-ci, pas question d'afficher le resultat, mais on peut observer le placement des termes non nuls de la
matrice grace a la commande spy.
>> spy(B)

Cette commande, qui fonctionne aussi avec les matrices pleines, peut être utile lors de la recherche d'erreurs au
sein d'un programme, en revelant un problème de remplissage.

Les matrices creuses se manipulent commes les pleines. Pour s'en convaincre, il suffit de les manipuler.

Une autre commande utile à la saisie de matrice est spdiags. Elle permet de créer des matrices à partir de leurs
"diagonales". Les arguments de la commande spdiags sont : une matrice dont les colonnes correspondent aux
diagonales (aux élements de bords près, voir plus loin), un vecteur ligne (ayant autant de colonnes que la matrice
précèdente) dont les valeurs déterminent le placement des diagonales, et enfin la taille de la matrice. Voici un
exemple :
>> A=spdiags([ 1; 2; 3], [ 0 ], 3, 3);full(A)

ans =

1 0 0
0 2 0
0 0 3
>> A=spdiags([ 1; 2; 3], [ 0 ], 4, 3);full(A)

ans =

1 0 0
0 2 0
0 0 3
0 0 0

Dans le vecteur de placement, l'indice 0 correspond a un placement sur la diagonale principale. Le terme est ici
abusif, il faut comprendre comme diagonale principale , l'ensemble des coefficients ai,i. On remarque que la
commande est acceptée pour des matrices non carrées. En fait, pour que la commande spdiags(A,v,n,p) soit
acceptee, il suffit que le nombre de lignes de la matrice A soit p et qu'elle ait autant de colonnes que le vecteur v.
>> d1=(11:14)';d2=(21:24)';
>> A=spdiags([d1 d2], [-1 1], 4, 4);full(A)

ans =

0 22 0 0
11 0 23 0
0 12 0 24
0 0 13 0

Le coefficient -1 dans le vecteur de position indique que le vecteur colonne d1 correspond aux coefficients ai,i-1.
On observe que le dernier élément de d1 n'est pas pris en compte dans la construction de la matrice, de même
que le premier élément de d2. Il faut donc faire attention, lors de la construction, aux termes qui risquent de ne
pas apparaître. Voici un exemple de création de matrice en fonction de sa taille.
>> n=5, e=ones(n,1); A=spdiags([e 2*e e], [-1 0 1], n, n); full(A)

n =

ans =

2 1 0 0 0
1 2 1 0 0
0 1 2 1 0
0 0 1 2 1
0 0 0 1 2

Ainsi qu'un exemple qui, en plus, combine plusieurs outils.

>> n=8, A=spdiags([[1:n]' [ones(n,1)*[1:n-1]]], [0:n-1], n+1, n); full(A)

n =

ans =

1 1 2 3 4 5 6 7
0 2 1 2 3 4 5 6
0 0 3 1 2 3 4 5
0 0 0 4 1 2 3 4
0 0 0 0 5 1 2 3
0 0 0 0 0 6 1 2
 0 0 0 0 0 0 7 1
0 0 0 0 0 0 0 8
0 0 0 0 0 0 0 0

Pour manipuler le contenu des matrices (creuses ou pleines), on utilise la commande find. Elle donne la liste des
indices qui vérifient un test donné. Ainsi, la liste des indices dont le coefficient est 4 dans la matrice A est :
>> [i,j]=find(A==4)

i =

4
1
2
3
4

j =

4
5
6
7
8

On profite de cet exemple pour observer qu'une matrice peut etre considerée par matlab comme un vecteur
(construit à partir des colonnes de la matrice). Ici, l'indice 31 correspond au couple ligne-colonne (4,4). Enfin on
remplace tous ces "4" par des "44".
>> ind=find(A==4)

ind =

31
37
47
57
67

Enfin on remplace tous ces "4" par des "44".

>> A(ind)=44*ones(size(ind)); full(A)

ans =

1 1 2 3 44 5 6 7
0 2 1 2 3 44 5 6
0 0 3 1 2 3 44 5
0 0 0 44 1 2 3 44
0 0 0 0 5 1 2 3
0 0 0 0 0 6 1 2
0 0 0 0 0 0 7 1
0 0 0 0 0 0 0 8
0 0 0 0 0 0 0 0

Ici, la commande size permet d'obtenir la taille du vecteur ind, afin que l'affectation soit cohérente (A(ind) et
44*ones(size(ind)) ont bien la même taille). On peut aussi utiliser cette commande pour créer des fonctions
capables d'agir sur des vecteurs. Créons la fonction F, telle que F(x)=exp(-1/x) si x>0 et F(x)=0 sinon. Pour cela
on crée le fichier F.m qui contient les lignes suivantes.
function y=F(x)

% F agit terme a terme sur les vecteurs

%
% F(x)=exp(-1/x) si x>0
% F(x)=0 sinon
%

y=zeros(size(x));
ind=find(x>0);
y(ind)=exp(-1./x(ind));

Il ne reste plus qu'à effectuer un tracé de la fonction F.

>> x=-1:0.0001:1; plot(x,F(x))

Cette fonction s'utilise bien sur des vecteurs. L'utilisation de la commande find rend la fonction plus rapide que s
elle avait ete programmée sans. Ecrivons la même fonction dans un style plus classique de programmation, puis
enregistrons-la sous F2.m.

function y=F2(x)

% F2 agit terme a terme sur les vecteurs

%
% F2(x)=exp(-1/x) si x>0
% F2(x)=0 sinon
%

for i=1:length(x)
if x(i)>0, y(i)=exp(-1./x(i));
else y(i)=0;
end
end

On peut maintenant comparer les vitesses d'execution des deux fonctions.

>> x=-1:0.0001:1;
>> tic, F2(x); toc

elapsed_time =

30.8921

>> tic, F(x); toc

elapsed_time =

0.0473

La manipulation matricielle est bien plus rapide (la vitesse d'execution dépend aussi de la machine, ici vos
résultats seront probablement différents).

Voila vous avez en votre possession une partie des outils de matlab, des compléments seront donnés pour les ca
particuliers mais vous pouvez aussi découvrir d'autres commandes à partir de l'aide en ligne (helpwin).
retour au sommaire