Introduzione a MatlabIntroduzione a Matlab

Gabriella Puppo

Che cosa è MatlabChe cosa è Matlab

Matlab èun linguaggio di programmazioneun ambiente di calcolo scientifico con

routines altamente specializzateun ambiente grafico

Argomenti trattatiArgomenti trattati

Matlab come calcolatriceInserire comandi, vettori, matriciOperazioni su vettoriCicliFile .m e functionsGraficiMetodi numerici

Matlab come calcolatricePer usare Matlab come calcolatrice, inserisco i comandi dopo il >>.Per esempio:

>> 2+1ans = 3

Oppure:>> log(4)ans = 1.3863

Matlab normalmente stampa solo le prime 4 cifre decimali,ma in realta’ne memorizza molte di piu’. Per vederle tutte:

>> format long>> log(4)ans = 1.38629436111989

Help Online

Matlab ha un ricchissimo help online, che e’ accessibile direttamente dal toolbar.In questo modo e’ facile recuperare non solo i comandi e la relativa sintassi, ma anche informazioni sugli algoritmi che sono applicati nei diversi comandi

Help onlineSpesso e’ pero’ piu’ veloce usare help direttamente dala finestra dei comandi. Per accedere alle informazioni, basta digitare help nella finestra dei comandi:

>> help

HELP topics:

matlab\general - General purpose commands.matlab\ops - Operators and special characters.matlab\lang - Programming language constructs.matlab\elmat - Elementary matrices and matrix manipulation.matlab\elfun - Elementary math functions.matlab\specfun - Specialized math functions.matlab\matfun - Matrix functions - numerical linear algebra.… etc.

Per avere informazioni su una particolare function, per esempio, eye:

>> help eye

EYE Identity matrix. EYE(N) is the N-by-N identity matrix. EYE(M,N) or EYE([M,N]) is an M-by-N matrix with 1's on the diagonal and zeros elsewhere. EYE(SIZE(A)) is the same size as A. See also ONES, ZEROS, RAND, RANDN.

Per cercare informazioni su un particolare argomento, si usa il comando lookfor (look for = cerca)

>> lookfor logarithmLOGSPACE Logarithmically spaced vector.LOG Natural logarithm.LOG10 Common (base 10) logarithm.LOG2 Base 2 logarithm and dissect floating point number.BETALN Logarithm of beta function.GAMMALN Logarithm of gamma function.LOGM Matrix logarithm.

L’output di lookfor contiene i nomi di tutte le functions che presentano la parola “logarithm”nel loro help.

Inserire comandi, vettori e matrici

Per inserire comandi, basta digitare il comando al promptper esempio:

>> pians = 3.1416

Matlab crea una variabile ans a cui assegna il valore richiesto (in questo caso pi greco). Anche qui:

>> format long>> pians = 3.14159265358979

Per inserire matrici, si usano parentesi quadre:il comando:

>> a=[2, 3; 1, 2]

produce in output:

a = 2 3 1 2

Notare che non c’è nessun bisogno di dimensionare la matrice: Matlab infatti attribuisce automaticamente la memoria richiesta. Attenzione!Matlab automaticamente stampa l’output di ogni comando: per eliminare questa risposta è necessario terminare il comando con un ;Questo comando, per esempio, non produce nessun output:

>> a=[2, 3; 1, 2];

E’ possibile costruire matrici automaticamente:

>> a=zeros(2)a = 0 0 0 0

crea una matrice 2 per 2 di zeri, mentre:

>> a=zeros(2,3)a = 0 0 0 0 0 0

crea una matrice 2 per 3.N.B. Le functions di Matlab (come zeros) possono accettare un numero variabile di elementi in input.Analogamente funzionano le functions ones (che genera matrici di 1), rand (che genera matrici di numeri casuali), eye (che genera le matrici identità).

Column notationIl carattere : indica un ciclo implicito, che si usa per creare vettori:

>> x=1:5x = 1 2 3 4 5Si può introdurre anche un incremento non intero:

>> x=1:.1:2x = Columns 1 through 8 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 Columns 9 through 11 1.8000 1.9000 2.0000

Operazioni su vettoriMatlab esegue automaticamente le operazioni algebriche sulle matrici:

>> a=ones(2,3);>> b=ones(2,3);>> a+bans = 2 2 2 2 2 2

o anche:>> a=2*eye(2)a = 2 0 0 2

Naturalmente, le operazioni richieste devono essere ben definite:

>> a*b??? Error using ==> *Inner matrix dimensions must agree.

Perché il prodotto fra matrici è definito solo quando il numero di colonne della prima matrice e il numero di righe della seconda coincidono.Posso invece moltiplicare a per la trasposta di b. Per calcolare la trasposta:

>> b'ans = 1 1 1 1 1 1

ora il prodotto è ben definito:

>> a*b'ans = 3 3 3 3

Si possono calcolare funzioni di matrici:

>> a=zeros(1,2)a = 0 0>> b=cos(a)b = 1 1

Con questo sistema è possibile calcolare in modo vettoriale i valori di una funzione:

>> x=1:0.1:2;>> fx=cos(3*x)+2;

Per calcolare un prodotto, una potenza o un quoziente, Matlab distingue due operatori diversi. Nel caso del prodotto per esempio:

* denota il prodotto fra due matrici .* denota il prodotto fra le singole componenti

La stessa distinzione vale per / (quoziente) e ^ (potenza)

Esempio

>> x=1:0.1:2;>> fx=cos(3*x)*exp(x);??? Error using ==> *Inner matrix dimensions must agree.

Il comando corretto è:

>> fx=cos(3*x).*exp(x);

Un altro esempio:

>> x=ones(2,2);>> x^2ans = 2 2 2 2>> x.^2ans = 1 1 1 1

Infatti X^2 indica il prodotto della matrice X con sé stessa, che è definito solo per matrice quadrate, cioè X^2 = X*X,mentre A=X.^2 indica la matrice con elementi A i,j = ( X i,j ) 2 .

Operatori relazionali

Gli operatori relazionali più comuni sono:== uguale~= diverso da< minore di<= minore o uguale etc.Esempi:

>> x=2;>> x==0ans = 0>> x==2 ans = 1

(questa relazione e’falsa:)=> ans=0

(questa relazione è vera:)=> ans=1

Gli operatori relazionali possono essere applicati anche alle matrici:

>> a=[1 2; 0 -1];>> a>0 (qui i primi due elementi sono veri)ans = 1 1 0 0

>> a>=0 (qui i primi tre elementi sono veri)ans = 1 1 1 0

Operatori logiciGli operatori logici più comuni sono:& and logico| or logico~ not logicoEsempi:

>> x=1; y= -1;>> x>0 & y>0 (questa relazione è falsa)ans = 0>> x>0 | y>0 (questa relazione è vera)ans = 1

Un esempio piu’ complicato:

if x=='domenica' | x=='sabato 'display('Evviva!')elseif x=='venerdi 'display('Torno a casa')elsedisplay('Vado al Poli')end

Se imposto:

ans =Evviva!

>> x='sabato '

ottengo:

In questo esempio, x deve sempre essere impostato come variabile di 8 caratteri

Ciclo if … elseif …endIl ciclo basato su if ha la struttura:if espressione istruzioniendEsempio:

>> a=[1,4];>> if a>0 sqrt(a) endans = 1 2

Le istruzioni vengono eseguite solo se espressione è vera, cioè se espressione è diversa da zero.

>> if cos(2)display('ciao')end

Esempio:

Ciclo for … endIl ciclo for ha la struttura:for variabile = espressione istruzioniendDi solito espressione è un vettore:

>> s=0;>> for i=1:10s=s+i;end>> ss = 55

calcola la somma dei primi 10 numeri interi

I cicli for possono essere uno dentro l’altro:

>> n=4;>> for i=1:n for j=1:i a(i,j) = 1; end end

Crea una matrice triangolare inferiore:

>> aa = 1 0 0 0 1 1 0 0 1 1 1 0 1 1 1 1

Non sempre i cicli hanno indici interi.Per esempio:

>> for x=[pi, 51, -72.1]display(x)end

in output produce questo risultato:

x = 3.14159265358979x = 51x = -72.1

Ciclo while … endIl ciclo while ha la seguente strutturawhile espressioneistruzioniendEsempio

>> i=1;>> while i<5 i=i+1; end>> ii = 5

Calcolare la precisione di macchina

Devo calcolare un numero x=2^p tale che x=x+1

%Calcola la precisione di macchinap = 0;epsilon=1;while 1~=1+epsilon epsilon = epsilon/2; p = p+1;endepsilon=epsilon*2, p=p-1

N.B. Il ciclo viene eseguito una volta di troppo, per questo nell’ultimariga il valore di epsilon viene corretto

Eseguendo il programma precedente, troviamo:

epsilon = 2.2204e-016p = 52

Questo è lo stesso valore contenuto nella variabile intrinsecaeps, che contiene appunto la precisione di macchina:

>> epsans = 2.2204e-016

Calcolare il piu’ piccolo numero floating point della forma x=2^p

Devo trovare un numero della forma x=2^p tale che x sia considerato 0.

% Calcola il piu' piccolo numero floating point della % forma xmin=2^px=1;while x>0 xmin = x; x=x/2;endxmin

Risultati:

Il piu’ piccolo numero floating point è>> xminxmin = 4.9407e-324Notare che se dimezzo xmin, trovo:>> xmin/2ans = 0

Il più grande numero floating point

Il più grande numero floating point della forma 2^p è>> xmaxxmax = 8.9885e+307

Se raddoppio xmax, trovo>> xmax*2ans = Inf

File .m e functions

Un file .m (M-file) è un programma riconoscibile da Matlab. La scrittura di files .m permette di:

Sperimentare con un algoritmo, senza dover reintrodurre una lunga lista di comandi

Ottenere una documentazione permanente per un lavoro

Ottenere programmi che possono essere riutilizzati, per esempio cambiando solo i dati

Scambiare programmi con altri utenti

Struttura di un file .m

Script M-files: sono files di comandi. Non hanno variabili in entrata e in uscita e operano sulle variabili del workspace

function M-files: sono files di comandi, che hanno argomenti in entrata e in uscita. Le variabili interne a questi programmi non influenzano le variabili del workspace

I files .m sono di due tipi:

Commenti

Sia gli scripts che le functions devono contenere righe di commento.

I commenti sono segnalati da %: Matlab ignora tutti i caratteri di una riga dopo il %

Le prime righe di commento di uno script o di una function diventano parte dello help online

Esempio: file radice.m% Questo file calcola la radice degli elementi di % una matrice a, se a>0, altrimenti da' un messaggio di erroreif a>=0 a=sqrt(a)else display('errore')end

Attenzione: nel workspace deve essere stata definita una variabile a. Inoltre l’esecuzione di questo script modifica il contenuto della variabile a

Function M-files

Esempiofunction a=radfunz(x)% RADFUNZ(X) calcola la radice degli elementi di X% se X>=0, altrimenti stampa un messaggio di errore%if x>=0 a=sqrt(x)else display('errore')end

Questo file deve essere salvato come radfunz.m

Struttura di una function

La function inizia con una riga che ne specifica il nome (nell’esempio radfunz), le variabili di input e le variabili di output.

La function deve essere salvata in un file con lo stesso nome (nell’esempio radfunz.m)

I commenti dopo la prima riga faranno parte dello help on-line

Seguono le istruzioni con eventuali altri commenti

Un altro esempio:function [xmin,xmax]=minmax(a)%MINMAX(A,M,N) calcola l'elemento minimo, XMIN, e l’elemento % massimo, XMAX della matrice A.xmin=Inf; xmax=-Inf;% ricava le dimensioni della matrice A:[m,n] = size(a);for i=1:m for j=1:n if a(i,j) > xmax xmax = a(i,j); end if a(i,j) < xmin xmin = a(i,j); end endend

function [out1,out2,…]=funz(in1,in2,….)

Gli argomenti in output vanno a sinistra dell’ =, fra [ ]

Gli argomenti in input vanno a destra dell’ = , fra ( )

Posso usare un numero di argomenti minore di quello indicato nella definizione della function, sia in entrata che in uscita.

Per esempio: a=funz(b), assegna a “in1” il valore “b”, e ad “a” il valore “out1”

La function precedente ha la seguente struttura

Grafici

Preparare un vettore di ascissePreparare un vettore di ordinateFare il graficoEsempio: grafico di cos(4x)*exp(x), su [0,2]

Per ottenere il grafico di una funzione, devo:

>> x=0:0.01:2;>> f=cos(4*x).*exp(x);>> plot(x,f)

EsercizioScrivere una function che calcoli la funzione esponenziale, utilizzando i primi N termini della serie di Taylor.

Valutare l’ errore per x fissato, utilizzando come confronto la funzione exp di Matlab, e disegnare un grafico dell’ errore in funzione di N.

La function che calcola l’esponenziale può essere scritta come:

function ex=esponenziale(n,x)% EX=ESPONENZIALE(N;X) Calcola l'esponenziale di e^x% utilizzando i primi N termini dello sviluppo in serie% dell'esponenzialeex=1;for k=1:n den=factorial(k); ex=ex+x.^k/den;end

Mentre la function che calcola l’ errore può essere scritta così:

function err=errore_exp(n,x)% ERR=ERRORE_EXP(N,X) Disegna un grafico dell'errore fra% la formula di Taylor calcolata per i primi K termini,% per k=1:N, per l' esponenziale e il% risultato fornito dalla function EXP% ERR(K) contiene il l'errore commesso con i primi K % termini K=1,...,N.

for k=1:n err(k)=abs(exp(x)-esponenziale(k,x));endsemilogy(err)

Si nota che l’ errore diminuisce rapidamente, se aumento N:

Tuttavia, per valori di N più elevati, l’ errore non diminuisce più,perché si è a livello dell’ errore di macchina.

Esempio. Grafico di una circonferenza

>> t=0:0.01:2*pi;>> x=cos(t);>> y=sin(t);>> plot(x,y,'g+')>> axis equal

Calcolo della norma 2 di una matrice 2 X 2

function c=norma(a)%NORMA(A) fornisce una stima della norma 2 di una matrice 2 per 2

%costruisce il cerchio unitario e ne stampa il grafico:t=0:0.05:2*pi;x=cos(t);y=sin(t);plot(x,y)axis equalhold on

%calcola a*[x;y] e stampa il grafico di ogni puntoc=0;for i=1:length(t) b=a*[x(i); y(i)]; plot(b(1),b(2),'g+') nb = sqrt( b(1)^2 +b(2)^2 ); if nb>c c=nb; endend

…il listato continua ...

Esempio di cancellazione numerica

Calcolare (1-x)^6 con le due formule:

y1 = (1-x)^6y2 = 1-6x +15x^2 -20x^3 +15x^4 -6x^5 +x^6

e confrontare y1 e y2 in un intorno di 1

Esempio di cancellazione numerica

%Calcola (1-x)^6 con le due formule:%y1 = (1-x)^6%y2 = 1-6x +15x^2 -20x^3 +15x^4 -6x^5 +x^6%e confronta y1 e y2 in un intorno di uno%k = 0for delta = [0.1, 0.01, 0.005, 0.0025] h = delta/100; x = 1-delta:h:1+delta; y1 = (1-x).^6; y2 = 1 -6*x +15*x.^2 -20*x.^3 +15*x.^4 -6*x.^5 + x.^6; k = k+1; subplot(2,2,k) plot(x,y1) hold plot(x,y2,'g') axis([1-delta 1+delta -max(abs(y2)) max(abs(y2)) ])end

L’output dello script precedente è:

Istruzione subplot

L’istruzione subplot(M,N,K) crea una figura contenente M*Ngrafici, distribuiti su M righe ed N colonne.

L’indice K indica che le istruzioni seguenti si riferiscono al K-esimo grafico, con K che varia fra 1 e M*N

Esercizi

1) Scrivere una function cha calcoli il valor medio di una sequenza di numeri assegnati.

2) Scrivere una function che produca una matrice quadrata di zeri, con degli 1 sulla diagonale che scende da destra verso sinistra.3) Scrivere una function che disegni un rettangolo, centrato sull’origine, con i lati assegnati in input.

4) Scrivere una function che disegni l’effetto di una matrice 2 X 2 assegnata, sul segmento S: 0;21 yxS

Risoluzione di sistemi lineari con Matlab

Matlab risolve il sistema lineare Ax=b con il comando:x=A\b.Se A è quadrata, questo comando implica i seguenti passi:- Calcola la fattorizzazione LU con pivoting sulle colonne- Risolvi i due sistemi triangolari

Se A è rettangolare (sistemi sovradeterminati), Matlab calcola lasoluzione nel senso dei minimi quadrati, cioè:- Calcola la fattorizzazione QR con pivoting sulle colonne- Risolve il sistema triangolare superiore

Se il condizionamento di A è elevato, stampa un messaggio diwarning, ma calcola ugualmente la soluzione.

Esempio>> a=[2 -1 0; -1 2 -1; 0 -1 2];>> b=ones(3,1);>> x=a\bx = 1.5000 2.0000 1.5000

Attenzione! Matlab calcola una soluzione anche quando il sistema non ammette soluzione.Nell’esempio seguente, a è singolare:

>> a=[1 2 3; 4 5 6; 7 8 9];>> b=[1; 1; 0];>> x=a\b;Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.541976e-018.>> norm(a*x-b)ans = 5.0990

Non è chiaro il modo in cui Matlab calcola la soluzione in questocaso.

Passare una function come argomento

Molte funzioni “lavorano” su altre funzioni, predefinite dall’utente.Supponiamo di voler scrivere una function che stampa il grafico di una funzione sull’intervallo [a,b].La function richiesta deve avere in input il nome della function di cui vogliamo il grafico e l’intervallo [a,b].La sua intestazione sarà quindi:

function plot_funz(funz,a,b)

La difficoltà nasce dal fatto che per Matlab il parametro “funz” deve essere una stringa di caratteri, oppure uno “handle” contenente il nome della function.

Supponiamo quindi di aver creato un M-file f.m che contiene la function di cui vogliamo il grafico:

function f=f(x)% Calcola la funzione exp(x)*cos(4*x)f=exp(x).*cos(4*x);

Poiché il nome della function deve essere passato come stringa di caratteri, la chiamata a plot_funz per disegnare il grafico di f(x) sull’intervallo [0,5] sarà:

>> plot_funz('f',0,5)

Oppure, in Matlab 6, si puo’ usare anche uno handle (@):

>> plot_funz(@f,0,5)

Function feval

All’interno della function plot_funz, è necessario valutare la funzione passata come argomento sui valori x delle ascisse.Per far questo devo usare la function feval:

f=feval(funz,x);

Questa istruzione trasforma la stringa funz , associando a funzil file funz.m . Poi, esegue il comando funz(x).Possiamo ora dare il listato della function plot_funz

Listato di plot_funz

function plot_funz(funz,a,b)% PLOT_FUNZ(FUNZ,A,B): disegna il grafico della funzione% FUNZ sull'intervallo [a,b], usando 201 punti.% FUNZ deve essere una stringa di caratterih=(b-a)/200;x=a:h:b;f=feval(funz,x); %Valuta funz nei valori xplot(x,f)

Inline functions

Per costruire funzioni semplici, posso usare l’istruzione inline, che genera funzioni di una riga. Consideriamo il file f.m:

function f=f(x)% Calcola la funzione exp(x)*cos(4*x)f=exp(x).*cos(4*x);

Un modo più semplice di generare la stessa funzione è:

>> f=inline('exp(x).*cos(4*x)')f = Inline function: f(x) = exp(x).*cos(4*x)

Attenzione!

L’oggetto f creato dall’istruzione inline è una stringa di caratteri.Quindi se voglio disegnare il grafico di f usando la function plot_funz, ho due possibilità:1) -creo un file f.m, che contiene la funzione desiderata. - chiamo plot_funz con il comando:

>> plot_funz('f',0,5)

2) -creo f come inline function e poi chiamo plot_funz:

>> f=inline('exp(x).*cos(4*x)');>> plot_funz(f,0,5)

Esempio

Risolvere l’equazione:

x^3 - 1 – exp(x) = 0

Definisco la funzione e traccio un grafico:

>> f=inline('x^3-1-exp(x)')f = Inline function: f(x) = x^3-1-exp(x)>> fplot(f,[-2,5])>> grid on

Ottengo la seguente figura:

Ci sono due zeri: uno nell’intervallo [1,3], l’altro in [4,5]

Quindi, per trovare gli zeri di f, posso usare i seguenti comandi:

>> [x1,res] = fzero(f,[1,3])x1 = 2.0811res = 0>> [x2,res] = fzero(f,[4,5])x2 = 4.5037res = 1.4211e-14

Se f fosse stata memorizzata nel file f.m, il comando sarebbe stato:

>> [x1,res] = fzero('f',[1,3])