Calcolo Numerico I - a.a. 2010-2011
Laboratorio 3 - Introduzione a MATLAB
Grafica 2D in Matlab
Per tracciare in Matlab il grafico di una funzione y = f (x) per
un intervallo assegnato di valori della variabile x, utilizziamo il
comando plot. Ad esempio
>> x=linspace(0,2*pi,20); % vettore delle ascisse
>> y=x.*sin(x); % vettore delle ordinate
>> plot(x,y)
Se vogliamo assegnare un colore e uno stile della linea (ad esempio
rosso (red) ’r’, e linea tratteggiata, ’--’)
>> plot(x,f(x),’r--’)
>> grid % tracciamo le linee di griglia
Il comando doc LineSpec mostra una lista di tutte le opzioni
disponibili per specificare lo stile della linea.
Per rappresentare sul medesimo grafico più di una funzione
>>
>>
>>
>>
x=linspace(0,2*pi,20);
y1=x.*sin(x);
y2=x.*cos(x);
plot(x,y1,x,y2), grid
Alternativamente
>>
>>
>>
>>
>>
>>
x=linspace(0,2*pi,20);
y1=x.*sin(x);
plot(x,y1,’r’), grid
hold on
y2=x.*cos(x);
plot(x,y2,’g’)
Il comando hold on funziona come un interruttore acceso/spento;
fino a quando non si digita hold off tutti i grafici vengono disegnati nella medesima finestra.
Per disegnare più grafici, ciascuno su una specifica finestra
>>
>>
>>
>>
>>
>>
x=linspace(0,2*pi,20);
y1=x.*sin(x);
plot(x,y1,’r’), grid
figure(2)
y2=x.*cos(x);
plot(x,y2,’g’)
Disegnamo ora il grafico di sin(x) per x ∈ [−π, π] e completiamolo
con un titolo e delle etichette sugli assi
>>
>>
>>
>>
>>
x=linspace(-pi,pi,200);
y=sin(x);
plot(x,y)
title(’Grafico di sin(x)’)
xlabel(’Asse x’), ylabel(’Asse y’)
2
Per ridimensionare gli assi
>> xmin=-2; xmax=2; ymin=-1; ymax=1;
>> axis([xmin xmax ymin ymax]);
Il comando doc axis visualizza tutti i possibili formati di scalatura e di modifica degli assi del grafico.
La barra di comandi al di sopra della finestra del grafico permette
numerose operazioni, fra cui lo zoom (icona “lente”) e la modifica
delle proprietà degli oggetti - spessore delle linee, colore delle linee,
carattere delle etichette...- del grafico stesso (icona “freccia”)
Funzioni simboliche (inline)
Esiste in Matlab una sintassi che permette di definire una funzione direttamente nello spazio di lavoro (ovvero “in linea”) senza
ricorrere ad un file esterno.
L’istruzione relativa è inline, seguita dall’espressione della funzione e dall’elenco dei parametri da cui essa dipende.
Per esempio, definiamo la funzione f (x) = (sin(x) + x)2:
>> f=inline(’(sin(x)+x).^2’,’x’)
dove abbiamo indicato esplicitamente che f è funzione di x. Attenzione alla sintassi con gli apici e i punti e attenzione all’operazione
di elevamento a potenza componente per componente!
3
Ad una funzione cosı̀ definita non sono associati dei valori numerici
(verificare con whos f), ma solo un’espressione matematica. Se ora
vogliamo associare dei valori numerici, scriviamo
>> x=linspace(0,2*pi,20);
>> y=f(x); %assegna a f dei valori numerici in
%corrispodenza del vettore x
Verificare con whos y che y è un vettore numerico!
Riprendiamo l’esempio del disegno del grafico della funzione y =
x sin(x) e usiamo il comando inline
>> x=linspace(0,2*pi,20);
>> f=x.*sin(x);
>> plot(x,f(x)), grid
Attenzione: perché il comando
>> plot(x,f)
non funziona?
Elementi di Grafica 3D in Matlab
Per disegnare una curva nello spazio 3D, ad esempio un’elica circolare definita in forma parametrica
4
Figura 1: Grafico dell’elica circolare ottenuto con il comando plot3.
>> t=linspace(0,10*pi,1000);
>> plot3(sin(t), cos(t),t) % sintassi generale: plot3(x,y,z)
Il risultato è rappresentato in Fig. 1
Per tracciare grafici di superfici, ad esempio della funzione z =
sin(x) sin(y), dapprima dobbiamo creare le matrici X e Y delle coordinate x e y. A questo scopo, usiamo il comando meshgrid che
converte una coppia di vettori x,y in una coppia di matrici X,Y,
come segue:
>> x=linspace(-pi,pi,40);
y=linspace(-pi,pi,40);
5
>> [X,Y]=meshgrid(x,y);
quindi calcoliamo in corrispondenza delle coordinate contenute in
tali matrici i valori della funzione z (attenzione alla sintassi .!)
>> Z=sin(X).*sin(Y);
e infine rappresentiamo la superficie
>> surf(X,Y,Z) %superficie interpolata
>> mesh(X,Y,Z) % superficie "wireframe"
I risultati ottenuti con i due comandi sono riportati in Fig. 2 La
modalità con cui viene colorata la superficie generata da surf può
essere modificata con il comando shading.
Il comando caxis([min,max]) imposta i valori che corrispondono
al minimo e al massimo della scala di colori e il comando colormap
modifica la scala di colori.
Per disegnare le linee di livello della funzione, possiamo utilizzare
infine il comando contour(X,Y,Z,N), eventualmente specificando nel parametro N il numero di linee di livello da disegnare. Il
risultato ottenuto è riportato in Fig. 3
6
Figura 2: Grafico 3D della funzione z = sin(x) sin(y) ottenuto con i comandi surf (figura
in alto) e mesh (figura in basso).
7
Figura 3: Linee di livello della funzione z = sin(x) sin(y) ottenute con il comando contour
.
Grafici in scala logaritmica
Spesso è utile rappresentare il grafico di una certa funzione in scala
logaritmica, su entrambi gli assi cartesiani o su uno solo dei due.
In Matlab esistono i seguenti comandi:
• loglog(x,y): grafico in scala logaritmica su entrambi gli assi
• semilogx(x,y): grafico in scala logaritmica sulle ascisse
• semilogy(x,y): grafico in scala logaritmica sulle ordinate
Osservazione importante: funzioni del tipo
y = kxn,
8
dove k e n sono due parametri reali, sono rappresentate da una
retta in scala logaritmica. Infatti, se prendiamo il logaritmo in
base 10 di entrambi i membri, otteniamo
Log10y = Log10k + nLog10x.
Ponendo Y = Log10y, X = Log10x, P = Log10k, tale espressione
si può scrivere come
Y = nX + P
Abbiamo ottenuto una retta, il cui coefficiente angolare della retta
è l’esponente n!
Ad esempio, disegnamo nella medesima finestra il grafico delle funzioni y = x, y = x2, y = x3, y = x4, prima in scala lineare e poi
in scala logaritmica su entrambi gli assi.
Per creare un vettore di ascisse x che comprenda diversi ordini
di grandezza conviene usare il comando logspace(p1,p2,N) che
genera un vettore di N elementi “logaritmicamente equispaziati”,
ovvero ottenuti come se usassimo il comando 10.ˆ(linspace(p1,p2,N))
>> x=logspace(-1,5,10);
>> plot(x,x,x,x.^2,x,x.^3,x,x.^4)
>> figure(2), loglog(x,x,x,x.^2,x,x.^3,x,x.^4)
Nel grafico in scala lineare riusciamo a distinguere solo la potenza
maggiore, mentre le altre risultano “schiacciate” vicino all’asse x,
mentre in scala logaritmica sono molto più chiari gli ordini di
9
grandezza. Inoltre, in scala logaritmica l’esponente di ciascuna funzione è immediatamente riconoscibile dalla pendenza della retta in
scala logaritmica (si veda la Fig. 4).
Figura 4: Grafico delle funzioni y = x, y = x2 , y = x3 , y = x4 in scala lineare (sinistra) e
in scala logaritmica su entrambi gli assi (destra).
Inoltre, funzioni del tipo
y = kanx,
sono rappresentate da una retta in un piano con asse y logaritmico. Se prendiamo il logaritmo in base 10 di entrambi i membri,
otteniamo infatti
Log10y = Log10k + nxLog10a.
Ponendo Y = Log10y, A = Log10a, P = Log10k, tale espressione
si può scrivere come
Y = nAx + P.
10
Quando usare la scala logaritmica ?
• per rappresentare un insieme di dati che si estende in un
intervallo molto vasto di valori;
• per mettere in evidenza particolari andamenti nella variazione
dei dati.
Ricordiamo che:
• in un grafico in scala logaritmica non si possono rappresentare
quantità negative o zero (scegliere opportunamente l’intervallo
delle ascisse!)
• Attenzione: Matlab assegna alle etichette dei segni di graduazione degli assi i valori rappresentati e non i loro logaritmi
11
Polinomi in Matlab
Matlab non prevede un oggetto particolare di tipo polinomio, ma
rappresenta i polinomi attraverso vettori che contengono i coefficienti del polinomio stesso, ordinati partendo dal coefficiente del
termine di grado maggiore (detto anche coefficiente direttivo) fino
ad arrivare al coefficiente del termine di grado minore (detto anche
termine noto)
Per esempio, i polinomi p(x) = x3 −x2 +4x−1 e q(x) = x4 −3x−1
di grado, rispettivamente, n = 3 e n = 4, diventano i vettori (di
lunghezza n + 1)
>> p=[1 -1 4 -1];
>> q=[1 0 0 -3 -1];
% zero per i termini con coeff. nullo
Valutiamo un polinomio con coefficienti contenuti in p in un punto
x e assegniamo alla variabile w tale valore
>> x=2;
>> w=polyval(p,x);
Possiamo anche valutare il polinomio in un vettore di punti contenuto in x
>> x=linspace(-1,1);
>> y=polyval(p,x);
>> plot(x,y) % disegnamo il grafico di p(x) in [-1,1]
12
Per sommare algebricamente i due polinomi p1(x) = 3x3 − 2x2 +
x + 15 e p2(x) = x − 2
>> p1=[3 -2 1 15];
>> p2=[0 0 1 -2];
>> s=p1+p2;
NB: i due vettori dei coefficienti degli addendi devono avere la
stessa lunghezza, quindi abbiamo dovuto “allungare” il vettore dei
coefficienti di p2 con degli zeri!
Per calcolare il prodotto di due polinomi, ad esempio p(x) di grado
m e q(x) di grado n
>> p=[3 -2 1 15]; %m=3
>> q=[1 -2];
%n=1
>> r=conv(p,q); % coeffs del polinomio di grado m+n=4
dove r è il vettore dei coefficienti del polinomio r(x) = p(x)q(x).
Per calcolare la derivata di un polinomio p(x) rispetto alla sua
variabile indipendente, p0(x), usiamo il comando polyder. Ad
esempio, calcoliamo p0(x), con p(x) = 3x3 − 2x2 + x + 15
>> p=[3 -2 1 15]; %m=3
>> polyder(p)
9
-4
1
dove il risultato è ancora un vettore di coefficienti con il significato
usuale.
13
Per calcolare la primitiva di un polinomio, si usa il comando polyint(p,k)
dove k rappresenta il valore a cui si vuole porre la costante dell’integrale indefinito. Essa può essere omessa nel caso in cui si voglia
k = 0. Ad esempio
Z
(x2 − 2x − 1) dx
>> p=[1 -2 -1];
>> P=polyint(p);
dove il risultato è ottenuto a meno di una costante (fissata qui a
zero). Se vogliamo invece calcolare l’integrale definito
Z 3
(x2 − 2x − 1) dx
0
possiamo combinare il comando polyint con il comando polyval
>> p=[1 -2 -1];
>> P=polyint(p);
>> I=polyval(I,3)-polyval(I,0)
14
Rappresentazione dei numeri reali:
È impossibile rappresentare su una macchina l’infinità dei numeri
reali. Si lavora con un sottoinsieme di dimensione finita: l’insieme
dei numeri floating-point.
Ad esempio:
2/3 = 0.6
>> 2/3
ans =
0.6667
1/3 = 0.3
>> 1/3
ans =
0.3333
I numeri reali vengono arrotondati.
Usare solo 4 cifre decimali per rappresentare i numeri sarebbe tuttavia molto grossolano, ma questo è solo un possibile formato del
sistema per visualizzare un numero reale che non coincide con la
sua rappresentazione interna.
Altri possibili formati sono:
>> format long
>> 2/3
15
ans =
0.66666666666667
>>1/3
ans=
0.3333333333333
>> 1/100
ans =
0.01000000000000
>> format long e
>> 1/100
ans =
1.000000000000000e-02
Altre variabili predefinite:
• eps (epsilon macchina),
• realmax, realmin (più grande e più piccolo numero floating
point),
• Inf (+∞ è il risultato di operazioni come 1/0 o di un overflow,
ad es. exp(1000))
• Nan (Not-a-Number, è il risultato di operazioni come 0/0 o
inf-inf).
16