Matlab 1
Matlab 1
Matlab 1
Matlab est un logiciel de calcul numérique produit par MathWorks (voir le site web http://
www.mathworks.com/). Il est disponible sur plusieurs plateformes.
Matlab est un langage simple et très efficace, optimisé pour le traitement des matrices, d’où son
nom. Pour le calcul numérique, Matlab est beaucoup plus concis que les “vieux” langages (C,
Pascal, Fortran, Basic). Un exemple: plus besoin de programmer des boucles modifier pour un
à un les éléments d’une matrice. On peut traiter la matrice comme une simple variable. Matlab
contient également une interface graphique puissante, ainsi qu’une grande variété d’algorithmes
scientifiques.
On peut enrichir Matlab en ajoutant des “boîtes à outils” (toolbox) qui sont des ensembles de
fonctions supplémentaires, profilées pour des applications particulières (traitement de signaux,
analyses statistiques, optimisation, etc.).
L’Université de Genève dispose d’une centaine de licences Matlab qui sont à la disposition de
la communauté universitaire. Il suffit d’être connecté au réseau de l’université pour pouvoir
l’utiliser. Plusieurs toolbox sont aussi disponibles (voir Appendice 1), leur nombre dépend de
l’environnement, il est parfois plus restreint.
Ces notes ne constituent pas une référence exhaustive sur MatLab, mais soulignent les no-
tions principales de manière succincte et servent de guide pour le travail en laboratoire. Elles
devraient inciter l’étudiant à chercher lui-même les compléments d’information qui lui sont
nécessaires, soit avec les outils d’aide en ligne, soit dans les ouvrages suivants:
• Introduction à Matlab
J.-T. Lapresté (Ellipses, 1999)
• Mastering Matlab 6
D. Hanselman B. Littlefield (Prentice Hall, 2001)
• Apprendre et maîtriser Matlab
M. Mokhtari A. Mesbah, (Springer, 1997)
• Solving problems in scientific computing using Maple and Matlab
W. Gander, J. Hrebicek (Springer, 1995, second edition)
• Numerical Methods Using Matlab
G. Lindfield J. Penny (Prentice Hall, 2nd edition : 2000)
Matlab est relativement coûteux car il contient un nombre impressionnant de fonctions. Il existe
une version étudiant à un prix abordable et un clone (Octave), disponible en freeware, dont la
compatibilité avec Matlab est assez bonne: (http://www.octave.org/ ou dans la distribution
SuSE de Linux).
Les exemples donnés dans le cours sont disponibles sur le serveur servtp3-1 des Travaux Pra-
tiques Avancés, dossier public/infophys/matlab/cours/. Les solutions d’exercices sont
dans le dossier public/infophys/matlab/exercices/.
Nous nous bornons ici à décrire le language Matlab qui est indépendant de la plateforme uti-
lisée. L’environnement de travail offre plusieurs fonctionalités assez conviviales qui sont utiles
au développement. L’utilisateur les mettra à profit sans difficulté.
MatLab1.fm (15October2004) 2
1. Aspects élémentaires
1.1 Aides
help -> donne de l’aide sur une fonction ou un toolkit (help help)
helpdesk -> documentation en hypertexte (requiert Netscape ou autre)
helpwin -> aide en ligne dans une fenêtre séparée
lookfor -> recherche d’un mot clé (lent)
which -> localise fonctions et fichiers
what -> liste des fichiers matlab dans le répertoire courant
exist -> check si une fonction ou une variable existe dans le workspace
who, whos -> liste des variables dans le workspace
>> autre=3;
Opérations élémentaires:
+ - * / or \ ^
>> 4/2
ans =
2
>> 2\4
ans =
2
MatLab1.fm (15October2004) 3
1.4 Variables spéciales
pi inf i or j realmin realmax eps ans flops (# d’opérations effectuées)
>> eps
ans =
2.2204e-16
>> realmax
ans =
1.7977e+308
>> realmin
ans =
2.2251e-308
c2 = 3*(2-sqrt(-1)*3)
c2 =
6.0000 - 9.0000i
real(c1) imag(c1) abs(c1)
>> angle(c1)
ans =
-1.1071
autres fonctions:
conj isreal
besselj besselh .... beta ... erf ... gamma ... legendre
MatLab1.fm (15October2004) 4
cross (produit vectoriel)
dot (produit scalaire)
>> v1=[1 3 5]
v1 =
1 3 5
>> v2=[2 4 6]
v2 =
2 4 6
>> cross(v1,v2)
ans =
-2 4 -2
>> whos
Name Size Bytes Class
ans 1x1 8 double array
c1 1x1 16 double array (complex)
c2 1x1 16 double array (complex)
v1 1x3 24 double array
v2 1x3 24 double array
1.7 Affichage
FORMAT Set output format.
All computations in MATLAB are done in double precision.
FORMAT may be used to switch between different output
display formats as follows:
FORMAT SHORT (default) Scaled fixed point format with 5 digits.
FORMAT LONG Scaled fixed point format with 15 digits.
FORMAT SHORT G Best of fixed or floating point format with 5 digits.
FORMAT LONG G Best of fixed or floating point format with 15 digits.
Autres format : cf help format
Spacing:
FORMAT COMPACT Suppress extra line-feeds.
FORMAT LOOSE Puts the extra line-feeds back in.
examples:
>> pi
ans =
3.1416
>> format long g
>> pi
ans =
3.14159265358979
MatLab1.fm (15October2004) 5
1.8 Entrées-sorties
Deux commandes utiles pour gérer le workspace, dont la taille dépend de votre espace de swap:
>> save % écrit toutes les variables du workspace dans le fichier matlab.mat
>> load % charge dans le workspace toutes les variables du fichier matlab.mat
MatLab1.fm (15October2004) 6
2. Vecteurs
2.1 Création de vecteurs
Par défaut, le vecteur est une ligne à plusieurs colonnes
a) vecteur ligne par énumération des composantes:
>> v = [1 3.4 5 -6]
v =
1.0000 3.4000 5.0000 -6.0000
c) vecteur colonne:
>> xcol = x’
xcol =
0
0.2856
0.5712
0.8568
1.1424
1.4280
1.7136
1.9992
2.2848
2.5704
2.8560
3.1416
MatLab1.fm (15October2004) 7
2.2 Adressages et indexages
>> x(3) % 3ème élément du vecteur x
ans =
0.5712
>> x([8 3 9 1]) % une sélection de composantes (on les désigne avec un autre
vecteur!)
ans =
1.9992 0.5712 2.2848 0
>> c = [a b]
c =
1 2 3 10 20 30
Tableau 1: Notations
[] énumération d’éléments
: descripteur d’éléments de vecteur/matrice
() ensemble d’arguments
, séparateur d’arguments
; séparateur des lignes dans les matrices
supression du résultat de l’évaluation d’une instruction
’ transposition de matrice
. force l’opérateur à s’appliquer sur chaque élément du vecteur/matrice
% délimitateur de commentaire
... continuation de l’instruction sur la ligne suivante
MatLab1.fm (15October2004) 8
3. Matrices
3.1 Création de matrices
Un matrice est un ensemble de lignes comportant toutes le même nombre de colonnes
Matlab, depuis la version 5, supporte les matrices à n dimensions (n>2)
On peut étendre aux matrices les autres manières de définir des vecteurs.
Par exemple:
>> m2 = [1:1:3 ; 11:1:13]
m2 =
1 2 3
11 12 13
3.2 Transposition
l’opérateur apostrophe utilisé pour créer un vecteur colonne est en fait l’opérateur transposition:
>> m2’
ans =
1 11
2 12
3 13
Une exception:
>> m2^2
??? Error using ==> ^
Matrix must be square.
Dans ce cas, Matlab veut calculer le produit matriciel m2 * m2
La solution est l’usage du point qui force l’opération sur chaque élément:
>> m2 .^ 2
ans =
MatLab1.fm (15October2004) 9
1 4 9
121 144 169
a) Multiplications
>> m1 % rappelons la définition de m1
m1 =
1 2 3
4 5 6
7 8 9
>> m1 * m2’ % le produit matriciel n’est possible que lorsque les dimensions
sont cohérentes
ans =
14 74
32 182
50 290
>> m1 * m2
??? Error using ==> *
Inner matrix dimensions must agree.
b) Divisions
>> m2/m3 % division matricielle à droite
ans =
1.0000 -0.0000
9.5406 -1.5960
>> m2\m3 % division matricielle à gauche (cf algèbre linéaire!)
ans =
-0.5000 -0.8257 -0.4500
0 0 0
0.5000 0.9419 1.1500
MatLab1.fm (15October2004) 10
>> m3./m2 % chaque élément de m3 est divisé par l’élément équivalent m2
ans =
1.0000 1.0000 1.0000
0.0909 0.2635 0.7692
>> zeros(2,5)
ans =
0 0 0 0 0
0 0 0 0 0
MatLab1.fm (15October2004) 11
- Changez la valeur de l’élément dans la 2ème ligne, 6ème colonne, que se passe-t-il?
- Mettez tous les éléments de la 4ème colonne à 4
- Créez B en prenant les lignes de A en sens inverse
- Créer C en accolant toutes les lignes de la première et troisième colonne de B à la droite de A
- Créer D sous-matrice de A faite des deux premières lignes et les deux dernières colonnes de
A. Trouvez aussi une manière de faire qui ne dépende pas de la taille de A.
Note: chacun de ces exercices se fait en une seule instruction, sans boucles ittératives.
>> A = [1 2 3 ; 4 5 6 ; 7 8 9 ]
A =
1 2 3
4 5 6
7 8 9
MatLab1.fm (15October2004) 12
>> tril(A) % triangle inférieur
ans =
1 0 0
4 5 0
7 8 9
Il y a encore bien d’autres fonctions pour travailler les matrices, voir la liste dans l’annexe.
c) Exercice (avancé):
Sans utiliser de boucles d’ittération, ajouter aux éléments d’une matrice l’indice de leur col-
onne.
sparse peut aussi être utilisé pour créer directement une matrice clairsemée:
>> S = sparse([2 1 3 4], [1 2 3 1], [4 5 6 7], 4, 3)
S =
(2,1) 4
(4,1) 7
(1,2) 5
(3,3) 6
Le gain de place dans le workspace est d’autant plus significatif que la matrice est grande (on
utilise bucky, une matrice clairsemée prédéfinie):
>> clear
>> B_sparse = bucky;% matrice clairsemée
>> B_full = full(B_sparse); % matrice complète
>> whos
Name Size Bytes Class
B_sparse 60x60 28800 double array
B_full 60x60 2404 sparse array
MatLab1.fm (15October2004) 13
4. Programmer en Matlab
4.1 Opérateurs logiques et de relation
< plus petit
> plus grand
<= plus petit ou égal
>= plus grand ou égal
== égal
~= pas égal
& et
| ou
~ not
xor(x,y) ou exclusif
any(x) retourne 1 si un des éléments de x est non nul
all(x) retourne 1 si tous les éléments de x sont nuls
isequal(A,B), ischar etc...
b) While
while expression
(commands)
end
c) If-then-else
if expression1
(commandes à exécuter si expression1 est“vrai”)
elseif expression2
(commandes à exécuter si expression2 est“vrai”)
else
(commandes à exécuter si aucune expression est “vrai”)
end
MatLab1.fm (15October2004) 14
4.3 M-Files ou scripts
Un script (ou M-file) est un fichier (message.m par exemple) contenant des instructions Matlab.
Voici un exemple de script:
Matlab vous offre un éditeur pour écrire et mettre au point vos M-files:
>> edit % lance l’éditeur de MatLab. Voir aussi la barre d’outils
Tout autre éditeur de texte convient aussi.
Les M-files sont exécutés séquentiellement dans le “workspace”, c’est à dire qu’ils peuvent ac-
céder aux variables qui s’y trouvent déjà, les modifier, en créer d’autres etc.
Exercice:
Dans un script “pluspetitnombre.m”, utiliser la structure “while” pour évaluer le plus petit
nombre tel que, ajouté à 1, on obtient un nombre supérieur à 1.
4.4 Fonctions
On peut écrire des fonctions MatLab que l’on peut ensuite appeler depuis un script.
Voici une fonction” temps” définie dans le fichier temps.m:
function y=temps(x)
% TEMPS(X) affiche un message suivant le temps qu’il fait
% et retourne le paramètre d’entrée X changé de signe
if length(x)>1 error(’X doit être un scalaire’); end
if x~=0
disp(’Hello, il fait beau’)
else
disp(’Espérons que demain sera meilleur!’)
end
y=-x;
return
MatLab1.fm (15October2004) 15
voir les utiliser depuis l’extérieur de la fonction. Par exemple:
global x y
Matlab offre plusieurs moyens de vérifier les arguments d’entrées et de sorties d’une fonction:
nargin retourne le nombre d’arguments d’entrée
nargout retourne le nombre d’arguments de sortie
nargchk vérifie le nombre d’arguments d’entrée
error Affiche un message d’erreur
inputname retourne le nom d’un argument d’entrée
Exercice:
Ecrivez une fonction “index_of_max” qui retourne l’index de l’élément le plus grand d’un vec-
teur A donné comme argument d’entrée.
Examples:
>> pwd
E:\matlab
>> cd E:\infophys\matlab\cours
>> ls
ReO3_GX.m diffsplines.m rex.m test_filter.m
ballistic.m message.m temps.m
ballistic_ODE.m myfun1.m test.dat
>> cd ..
>> cd exercices
>> pwd
E:\infophys\matlab\exercices
MatLab1.fm (15October2004) 16
5. Graphisme
Nous donnons ici les indications minimum. Utilisez help et les autres commandes d’aide pour
affiner vos connaissances et vos graphiques.
5.1 Graphiques à 2D
>> help graph2d % intro. au graphisme 2D et tableau des fonctions disponibles
a) courbes: plot
>> x=linspace(0,2*pi,30);
>> y=sin(x);
% un plot rudimentaire:
>> plot(x,y)
>> title(’Trigonométrie’) % MatLab1_fig1
MatLab1_fig1
% quelques améliorations:
>> grid on
>> axis([ 0 2*pi -1 1])
% add a second curve
>> hold on
>> z=cos(x)
>> plot(x,z,’c+’) % MatLab1_fig2
>> clf % efface la figure
MatLab1_fig2
semilogx, semilogy et loglog sont semblables à plot mais permettent de faire des plots log.
MatLab1.fm (15October2004) 17
Exercice :
Ecrivez un script pour montrer l’évolution graphique de vos dépenses pour la santé en fonction
du montant de vos frais médicaux annuels. Evaluez 2 situations: 1) Prime mensuelle: 300.- ; la
caisse rembourse 90% des frais avec une franchise de 150.- 2) rabais de 37% sur vos primes si
la franchise est de 1’500.- Quand êtes-vous gagnant?
b) histogrammes: hist
>> x=-2.9 : 0.2 : 2.9; %défini l’intervalle et la largeurs des canaux de l’his-
togramme
>> y=randn(50000,1); % génère des nombres aléatoires répartis selon une distr.
normale
>> hist(y,x) % dessine l’histogramme MatLab1_fig3
MatLab1_fig3
5.2 Graphiques à 3D
>> help graph3d % introduction au graphisme 3D et tableau des fonctions dis-
ponibles
MatLab1.fm (15October2004) 18
MatLab1_fig4
Résultat de testmeshgrid.m :
>> x =
1 10 100
MatLab1.fm (15October2004) 19
y =
0 25 50 75 100
X =
1 10 100
1 10 100
1 10 100
1 10 100
1 10 100
Y =
0 0 0
25 25 25
50 50 50
75 75 75
100 100 100
P_4_2 =
10
75
MatLab1_fig5
d) Surface avec illumination: surfl
>> clf
>> [X,Y,Z] = peaks(30);
>> surfl(X,Y,Z) % graphique avec illumination
>> shading interp % meilleure interpolation
>> colormap pink %choix d’une palette de couleurs prédéfinie
MatLab1.fm (15October2004) 20
>> view(-37.5+90, 30) % changement point de vue: view(azimut, elevation)
MatLab1_fig6
MatLab1_fig6
MatLab1_fig7
% pour transformer ce graphe en échelle de couleurs
>> colormap(’gray’)
>> pcolor(X,Y,Z) % plot pseudo-couleur de Z sur la grille X,Y
>> shading interp
>> axis(’square’) % MatLab1_fig8
MatLab1.fm (15October2004) 21
>> colormap(’default’)
MatLab1_fig8
5.3 Animations
moviein: initialisation d’un buffer pour les images
getframe: capture une figure dans une image de l’animation
movie: affiche l’animation
Exemple:
% RotatePeaks.m - Exemple d’animation
%
% On crée plusieurs frames dans lesquels on met un même graphique 3D
% vu sous des angles différents de manière à créer une impression de
% rotation
%
clear
%
% préparation du graphique qui va servir de base à l’animation
[X,Y,Z] = peaks(30);
surfl(X,Y,Z)
axis([-3 3 -3 3 -10 10])
axis off
shading interp
MatLab1.fm (15October2004) 22
colormap (hot)
%
% construction de l’animation
m=moviein(15); %initialise l’animation (15 images)
for i=1:15
view(-37.5+24*(i-1),30) % change le point de vue
m(:,i)=getframe; % capture une image; Notez que m(i)= ... va ausi
end
%
% présentation du résultat
movie(m,5) % affiche l’animation 5+1 fois
Exercice:
Construisez une animation du jeu de la vie.
On part avec un espace 2D composé de (m*n) cellules qui sont aléatoirement vivantes ou
mortes. A chaque génération (itération), cette population évolue en suivant les deux règles suiv-
antes:
1) une cellule morte devient vivante si elle est entourée exactement par 3 cellules vivantes parmi
les 8 cellules qui l’entourent.
2) une cellule vivante survit seulement si elle est entourée par deux ou trois cellules vivantes.
Conseil: utiliser la fonction graphique “image” pour visualiser l’espace dans lequel évolue la
population de cellules.
MatLab1.fm (15October2004) 23
6. Quelques applications de Matlab
6.1 Fits et interpolations
a) Fits de polynômes: (interpolation.m)
>> x = [0 : 0.1 : 1];
>> y = [ -.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2];
>> p = polyfit(x,y,2) % fit un polynôme du deuxième degré
p =
-9.8108 20.1293 -0.0317
La solution est y = -9.8108x2 + 20.1293x - 0.0317.
Regardons le résultat:
>> xi = linspace(0, 1, 100);
>> yi = polyval(p, xi);
>> plot(x, y, ’o’, xi, yi, ’--’) % MatLab1_fig9 (plus bas)
MatLab1_fig9
MatLab1.fm (15October2004) 24
6.2 Intégrations
Note : “humps” est une fonction numérique prédéfinie dans MatLab.
a) Intégrer un vecteur:
>> x = linspace(-1, 2, 1000);
>> y = humps(x);
>> area=trapz(x,y)
area =
26.3450
6.3 Différentiations
Contrairement à l’intégration, la différentiation numérique est difficile, puisqu’elle est très sen-
sible à la forme locale de la fonction. La meilleure manière de faire de la différentiation est de
fitter un polynôme ou une fonction spline. Dans les deux cas, on connaît l’expression analytique
de la dérivée.
a) Dérivée approximative
La fonction DY = diff(Y) retourne en DY la dérivée de Y calculée simplement en prenant la
différence entre les éléments adjacants de Y. “diff” retourne un vecteur moins long que Y, ce
qui peut être gênant. Dans ce cas, vous pouvez utiliser “gradient”.
MatLab1.fm (15October2004) 25
% DIFFSPLINES.M Différentiation avec des cubic-splines
% Fitte un ensemble de valeurs avec des splines et calcule la dérivée
clear
clf
%Il faut commencer par fitter des cubic splines (voir plus haut)
xi = linspace(0, 1, 100); % mesh fin entre 0 et 1
yi = spline(x, y, xi); % fit des cubic splines.
plot(x, y, ’o’, xi, yi, ’-’);
hold on
% différentiation:
pp = spline(x,y); % obtient les polynômes partiels de la spline
if pp.pieces(1)~=length(y)-1
error(’pas de pp-form pour cette spline cubique’)
end
[br, coefs, l, k] = unmkpp(pp); % obtient les coefficients des polynômes
sf = k-1:-1:1; % scaling factors
dco = sf(ones(l,1),:) .* coefs(:,1:k-1); % évalue les coeffs de dériva-
tion
ppd = mkpp(br, dco); %construit la pp-form pour la différentiation
yid = ppval(ppd, xi); % calcul de la dérivée
plot( xi, yid/4, ’--’); MatLab1_fig10
grid on
title(’cubic spline (ligne continue); dérivée (pointillé)’)
MatLab1_fig10
MatLab1.fm (15October2004) 26
Exercice 6.3.1 :
Ecrivez une fonction derpline(x, y, xi) qui calcule la dérivée par la méthode des splines,
d’une courbe définie par un ensemble de paires (x,y) et retourne la dérivée aux différentes va-
leurs spécifiées dans un vecteur xi. Cette fonction doit pouvoir donner un message d’erreur si
les paramètres d’entrées ne sont pas ceux qu’elle attend. (Indication: s’inspirer du script diff-
splines.m ci-dessus)
d) Gradient, Laplacien
>> gradient (F) % dérivée de F si F est un vecteur
>> gradient(H) % gradient de H, si H est une matrice
>> del2(U) % approximation du Laplacien de U
---r1------r4------r7------r10---
| | | | |
v0 r2 r5 r8 r11
| | | | |
---r3------r6------r9------r12---
Note:
Selon le système d’équations linéaires, on doit faire attention à la méthode utilisée. Si inv ne
marche pas, voici un petit guide pour tenter de trouver une solution :
Si A est une matrice triangulaire régulière, essayer \
sinon si A est définie positive ou carrée symétrique ou hermitienne, essayer la factorisation
de Cholesky : chol
sinon si A est non carrée, essayer la décomposition QR (Q= othogonal unitary matrix,
R=upper triangular matrix, ) : qr
MatLab1.fm (15October2004) 27
sinon pour A clairsemée, il faut essayer la famille d’outils que matlab met à disposition dans
ce cas (voir la discussion dans Numerical Methods using Matlab, 2nd edition,
G. Lindfield and J.Penny, Prentice Hall, 2000)
Ψ( r ) = ∑ bi Φi ( r )
i=1
L’équation de Schrödinger se transforme en
∑ ( Hi, j – Eδ i, j )bi = 0
i
qui n’a de solutions que si le déterminant de (Hi,j-Eδi,j)=0 (voir Mécanique Quantique).
Pour un Hi,j donné, on veut trouver les valeurs propres E de l’énergie et les vecteurs propres bi
qui donnent les fonctions d’onde.
La fonctions eig permet d’obtenir vecteurs et valeurs propres.
Illustration:
Voici comment se décrivent les états électroniques dans un cristal de ReO3 (approximation des
liaisons fortes, voir cours de Physique du Solide):
MatLab1.fm (15October2004) 28
0 0 0 d1 sy sz 0 0 0 0 0 0 0 0;
0 0 0 0 b1 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 b1 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 d1 sz sx 0 0 0 0 0;
0 0 0 0 0 0 0 b1 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 b1 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 d4 0 Sx Sy Sz;
0 0 0 0 0 0 0 0 0 0 d4 -sqrt(3)*Sx sqrt(3)*Sy 0;
0 0 0 0 0 0 0 0 0 0 0 a1 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 a1 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 a1];
E(:,j)=sort(eig(A+A’));%A+A’= Hamiltonien,eig=valeurs propres,
%sort=classement
end
[nb nk]=size(E); %nb: nombre de valeurs propres
figure(1)
hold; grid;
for n=1:nb
plot(kx,E(n,:)) % graphe des bandes (valeurs propres de l’Hamiltonien H)
end
axis([min(kx) max(kx) -0.6 0.8])
xlabel(’kx’)
ylabel(’Energy’)
title(’Structure de bandes du ReO3 : G-X’);
Cet exemple illustre admirablement le language synthétique et les outils qu’offre MatLab!
Exercice 6.4.2:
Déterminez les fréquences propres (naturelles) du système de ressorts couplés ci-dessous.
% __ __
% //| | | k4 | |
% //| k1 | |--------/\⁄ \⁄ \⁄ \------- | |
% //|--/\⁄ \⁄ \--| | __ |m3|
% //| |m1| k2 | | k3 | |
% //| | |--/\⁄ \⁄ \--|m2|--/\⁄ \⁄ \--| |
% //| |__| |__| |__|
%
% |--> q1 |--> q2 |--> q3
%
%
MatLab1.fm (15October2004) 29
ode113 :Adams-Bashforth-Moulton.
MatLab offre aussi une démo (pas très explicite) sur la solution des ODE:
>> odedemo
La stratégie pour résoudre les ODE est d’écrire une fonction (par exemple la fonction “myode”,
contenue dans le M-file myode.m) qui décrive le membre de droite de l’ équation à résoudre
(voir help odefile).
On peut alors obtenir la solution y ( t ) , dans un intervalle t=[tinitial tfinal] et pour des conditions
initiales définies dans y0 en utilisant ode45.
ode45 (et les autres) permet de résoudre des systèmes de plusieurs ODE couplés. On va utiliser
cette propriété dans l’exemple ci-dessous.
MatLab1.fm (15October2004) 30
0
-g ];
clear; clf;
global g
g = 9.81; % const. de gravitation
% conditions initiales du tir:
h = 1; % altitude de départ
v0 = 25; % vitesse initiale
theta = pi/180*15; % angle de tir
% conditions initiales:
xin = [0, h, v0*cos(theta), v0*sin(theta)];
% temps de vol jusqu’à atteindre l’altitude 0
tmax = (xin(4) + sqrt(xin(4)^2 + 2*g*xin(2)))/g;
% solution numérique des ODE :
[t x] = ode23(’ballisticODE’, [0 tmax], xin);
%plot de la solution avec une interpolation faite avec des splines cubiques
N = max(x(:, 1)); xi = 0:N/100:N;
axis([0,max(x(:,1)), 0, max(x(:,2))])
grid on; hold on;
plot(xi, spline(x(:,1), x(:, 2), xi), ’:r’);
hold off;
Exercice 6.5.1 :
Reprenez l’exemple du tir ballistique et introduisez une force de frottement proportionnelle à
la vitesse (voir MatLab1_fig11 ).
MatLab1_fig11
Exercice 6.5.2 :
Résolvez les équations du mouvement pour les ressorts couplés de l’exercice 6.4.2
MatLab1.fm (15October2004) 31
6.6 Zéros et minima de fonctions; fits
a) Recherche des zéros d’une fonction
hump est une fonction prédéfinie dans MatLab (voir MatLab1_fig12 )
MatLab1_fig12
b) Minimalisations
Recherche des extrema d’une fonction.
Très utile quand vous voulez fitter une fonction paramétrique à une série de mesures afin de dé-
terminer les paramètres de votre modèle.
MatLab1.fm (15October2004) 32
3.14159268915185
>> fminbnd(’cos’, 0, 15)
ans =
9.42477753165182
Exercice 6.6.1:
Ecrivez un script fit.m qui
1) génère une mesure m(x) décrite par une gaussienne avec du bruit (utiliser randn).
2) fitte un modèle gaussien g(x) en utilisant fminsearch. Les paramètres du modèle seront:
la position x0 du centre de la gaussienne,
la largeur de la gaussienne σ
l’amplitude de la gaussienne A
Pour utiliser fminsearch, il faudra écrire une fonction fitgaus qui retourne la somme des résid-
us entre la mesure et le modèle. Le résidu au point x est: (m(x) - g(x))2 . Utiliser “sum”.
3) Construire une figure en deux parties l’une montrant la mesure avec le modèle de départ et
l’autre la mesure avec le fit obtenu. (MatLab1_fig13 )
Rappel : gaussienne :
2
( x – x0 )
g ( x ) = A exp – --------------------
-
2σ
2
MatLab1_fig13
MatLab1.fm (15October2004) 33
6.7 Filtrage de signaux
Le filtrage des signaux est un très vaste domaine. On trouvera une série de filtres prêts à l’em-
ploi dans divers toolbox. Citons les fonctions:
conv - convolution
filter - 1D Digital filter
filter2 - 2D Digital filter
latcfilt - Lattice filter
fir1, fir2 ...-Window based FIR filter design (Finite Impulse Response)
medfilt1 - 1D median filtering
besself - Bessel analog filter
clear
figure(1)
clf
subplot(2,1,1)
N = 1001;
bruit = 0.1
t = [0 : N-1]’/(N-1);
Smooth = exp( -100*(t-1/5).^2) + exp( -500*(t-2/5).^2) + ...
exp(-2500*(t-3/5).^2) + exp(-12500*(t-4/5).^2) ;
Noisy = Smooth + bruit* rand(size(t)); % avec bruit
plot(t, Noisy)
hold on
title(’Filtre fenêtre constante’)
subplot(2,1,2)
plot(t, Smooth-Filtered)
title(’Smooth-Filtered’)
MatLab1.fm (15October2004) 34
Cet exemple montre clairement l’effet de retard introduit par filter.
MatLab1_fig14
MatLab1.fm (15October2004) 35
Exemple 2 : Utilisation plus subtile de filter:
Filtre de Savitzky-Golay : ce filtre revient à fitter localement un polynôme.
clear
figure(2)
clf
subplot(2,1,1)
N = 1001;
bruit = 0.1
t = [0 : N-1]’/(N-1);
Smooth = exp( -100*(t-1/5).^2) + exp( -500*(t-2/5).^2) + ...
exp(-2500*(t-3/5).^2) + exp(-12500*(t-4/5).^2) ;
Noisy = Smooth + bruit* rand(size(t)); % avec bruit
plot(t, Noisy)
hold on
title(’Filtre de Savitzky-Golay’)
M = 3 % 4 degré du polynome
nl= 8 % 5 nombres de points à gauche
nr= 8 % 5 nombres de points à droite
A = ones(nl+nr+1, M+1);
for j = M:-1:1
A(:,j) = [-nl:nr]’ .* A(:,j+1);
end
[Q, R] = qr(A) ; % Orthogonal-triangular decomposition
c = Q(:, M+1) / R(M+1, M+1);
Filtered = filter(c(nl+nr+1:-1:1), 1, Noisy);
% compensation du retard introduit par le filtre :
Filtered(1:nl) = Noisy(1:nl);
Filtered(nl+1:N-nr) = Filtered(nl+nr+1:N);
Filtered(N-nr+1:N) = Noisy(N-nr+1:N);
plot(t, Filtered+0.2, ’-r’)
subplot(2,1,2)
plot(t, Smooth-Filtered)
title(’Smooth-Filtered’)
MatLab1.fm (15October2004) 36
MatLab1_fig15
MatLab1.fm (15October2004) 37
Exercice 6.7.1:
a) Ecrivez un script TEST_CONV.M qui filtre le signal de TEST_FILTER.M par convolution
avec une fonction gaussienne. Utiliser conv.
b) Déconvoluez le résultat pour tenter de retrouver le signal de départ. Utiliser deconv.
A titre d’exemple, nous appliquons un filtre non linéaire simple à une mesure de comptage
simulée. Le but du filtre est de mettre en évidence un petit signal gaussien superposé à un grand
bruit de fond et noyé dans le bruit statistique.
Le principe de ce filtre est de prendre la dérivée d’ordre n du signal et de l’élever à une puissance
p. On peut choisir interactivement m et n, ainsi que la statistique du signal (en donnant le taux
de comptage accumulé à son pic). Le bruit est réparti selon la racine carée du signal, comme
c’est le cas pour toute mesure de comptage.
MatLab1.fm (15October2004) 38
maxfin=max(fin);
minfin=min(fin);
fin=fin./(maxfin-minfin);
fin=fin.*an;
plot(x,y,’-’,xdi,fin,’-’)
grid
Exercice 6.7.2:
Observez comment se comporte le filtre rex.m
- a) sans bruit (dans ce cas, le taux de comptage dans le pic ne joue pas de rôle). Essayez no-
tamment n=2 p=1; n=3 p=1 et n=3 p=3. (voir MatLab1_fig16 )
- b) avec du bruit, en fonction du taux comptage accumulé dans votre signal
MatLab1_fig16
Exercice 6.7.4:
Construisez un filtre plus performant que celui implémenté dans rex.m et rex2.m
MatLab1.fm (15October2004) 39
6.8 Optimisation
Au paragraphe 6.6 nous avons abordé la minimalisation. On appelle généralement optimisation
une minimalisation sous contraintes. Pour traiter ce type de problème avec Matlab, il faut avoir
à disposition la “Optimization Toolbox” qui ne fait pas partie de Matlab lui-même. Cette toolbox
est à disposition sur le réseau de l’Université, mais le nombre de licences est assez restreint.
help optim donne les grandes lignes des outils de cette toolbox
helpdesk fournit une description plus détaillée.
Au travers d’un exemple et d’un exercice, nous allons étudier deux des fonctions d’optimisation
disponibles.
Exemple d’optimisation :
clear
disp(’Problème: minimiser cx, sous les contraintes Ax <= b’)
disp(’avec 3 constraintes sur x’)
disp(’ ’)
disp(’Coefficients c de la fonction objectif f(x)’)
c = [-5 -4 -6]
disp(’Constraintes : membres de gauche’)
A = [1 -1 1; 3 2 4; 3 2 0]
disp(’Contraintes: membres de droite’)
b = [20 42 30]
%
% Limites inférieures et supérieures de x :
n = max(size(c));
L = zeros(n,1);
U = 10^10*ones(n,1);
%
% Optimisation :
[x,lambda,exitflag] = linprog(c,A,b,[],[],L,U);
if exitflag<= 0
disp(’ERROR, exitflag =’)
exitflag
else
disp(’Solution optimale pour x’)
x
disp(’Valeur optimale de f(x)’)
z = c*x
disp(’Vérification des contraintes :’)
for i=1:n
A(i,:)*x
end
end
MatLab1.fm (15October2004) 40
Exercice 6.8.1:
Déterminer un acier aux propriétés données, à partir de relations empiriques.
Le problème est exposé ci-dessous:
MatLab1.fm (15October2004) 41
% demandé s’il était possible d’utiliser une stratégie plus adéquate.
%
% La méthode présentée ici est basée sur la théorie de l’optimisation
% (minimalisation sous contraintes) qui consiste, à partir d’un acier de
% départ, à descendre dans l’espace de minimalisation jusqu’au minimum qui
% donnera l’acier aux propriétés optimum. Un bémol : la solution peut toujours
% n’être que un minimum local. Il faut donc la vérifier avec soin.
% Une difficulté résolue par cette méthode : le respect de contraintes, qui
% peuvent être de plusieurs types : inégalités, égalités, limites inférieures
% et/ou supérieures.
%
% En Matlab, les problèmes d’optimisation peuvent se résoudre aisément avec
% une des fonctions du "toolbox" d’optimisation. Nous utiliserons fminimax.
%
% Pour fminimax, il faut fournir une fonction de "résidus", ce que l’on peut
% faire en évaluant:
%
% residu d’une propriété de l’acier x =
% (propriété désirée - propriété obtenue pour l’acier x)^2
%
% La fonction (ResiduProprietes dans l’exemple ci-dessous) doit fournir un
% vecteur contenant le résidu pour chaque propriété. Comme les paramètres
% de la fonction sont définis, on devra passer les autres quantités
% nécessaires par l’intermédiaire de variables globales.
MatLab1.fm (15October2004) 42
Appendice 1 :
Les principales “toolbox” disponibles sur le serveur de MatLab à l’Université
MatLab1.fm (15October2004) 43
Appendice 2 : Tableaux de fonctions
Extrait de : Mastering Matlab 5, D. Hanselman, B. Littlefield, Prentice Hall, 1998.
Operations
Syntaxe d’adressage
MatLab1.fm (15October2004) 44
Matrices spéciales
Fonctions matricielles
MatLab1.fm (15October2004) 45
Fonctions matriceilles (suite)
MatLab1.fm (15October2004) 46
Tratiement de données
MatLab1.fm (15October2004) 47
Transformées de Fourier
MatLab1.fm (15October2004) 48
Table des Matières
Introduction 2
1. Aspects élémentaires
1.1 Aides 3
1.2 Variables scalaires, workspace, opérations élémentaires 3
1.3 Commentaires, ponctuation 3
1.4 Variables spéciales 4
1.5 Nombres complexes 4
1.6 Fonctions mathématiques 4
1.7 Affichage 5
1.8 Entrées-sorties 6
1.9 Terminer Matlab 6
1.10 Personnaliser Matlab 6
2. Vecteurs
2.1 Création de vecteurs 7
2.2 Adressages et indexages 8
2.3 Combinaison de vecteurs 8
3. Matrices
3.1 Création de matrices 9
3.2 Transposition 9
3.3 Opérations scalaires-matrices 9
3.4 Opérations entre matrices 10
3.5 Matrices particulières 11
3.6 Caractéristiques des matrices 11
3.7 Manipulations de matrices et sous-matrices 11
3.8 Matrices clairsemées 13
4. Programmer en Matlab
4.1 Opérateurs logiques et de relation 14
4.2 Contrôler l’exécution 14
4.3 M-Files ou scripts 15
4.4 Fonctions 15
4.5 Gestion du système de fichiers 16
5. Graphisme
5.1 Graphiques à 2D 17
5.2 Graphiques à 3D 18
5.3 Animations 22
MatLab1.fm (15October2004) 49
6. Quelques applications de Matlab
6.1 Fits et interpolations 24
6.2 Intégrations 25
6.3 Différentiations 25
6.4 Algèbre linéaire 27
6.5 Equations différentielles ordinaires (ODE) 29
6.6 Zéros et minima de fonctions; fits 32
6.7 Filtrage de signaux 34
6.8 Optimisation 40
Appendices
A.1 Les toolbox disponibles à l’Université 43
A.2 Tableaux de fonctions 44
MatLab1.fm (15October2004) 50