Tecnologie digitali
Scheda di laboratorio n. 1/2008
MatLab : che fare?
Avendo a disposizione gli elementi di base di un linguaggio di programmazione che permette di fare
calcoli con quantità numeriche, sorge il problema di quali applicazioni siano possibili. È bene tenere
presente che qualunque settore della matematica, in particolare se applicata alla fisica, usa il calcolo
numerico per eseguire conti altrimenti proibitivi1 . Il calcolo numerico, grazie anche alla efficienza
raggiunta nei moderni elaboratori, è dunque una tecnica tra le più importanti nel lavoro scientifico
moderno. Illustriamo una applicazione, in una situazione elementare, ma di notevole importanza storica:
l’uso di formule approssimate per il calcolo dell’integrale definito di una funzione di variabile reale.
Inoltre impareremo a rappresentare graficamente, anche in forma animata, leggi orarie di punti materiali
in movimento.
Vettori
Il tipo di dato fondamentale di MatLab è il vettore. In senso informatico esso è una raccolta di elementi
dello stesso tipo (numeri per il momento) ordinati in una successione finita. L’elemento generico del
vettore si indica con il suo nome-radice, ad esempio a, tramite la scrittura a(i) che indica l’elemento
che si trova all’i-simo posto: a(3) indica l’elemento che si trova al terzo posto a partire dall’elemento
iniziale a(1).
Esistono casi particolari in cui l’assegnazione di un vettore è facilitata dalla sintassi di MatLab.
Innanzitutto vi è il caso di vettori di numeri reali che rappresentino la equipartizione di un intervallo
chiuso I = [a, b] in un numero prefissato di N sottointervalli di ampiezza dx = (b − a)/N . L’istruzione
è, ad esempio:
x = linspace (a, b, N+1)
N+1 rappresenta il numero di punti equispaziati, compresi i due estremi in cui è suddiviso l’intervallo.
Questa istruzione è molto utile insieme con la proprietà delle funzioni di MatLab di agire attraverso i
vettori, ossia di restituire, quando hanno come argomento un vettore x, un vettore che ha per elementi
i valori della funzione f (x) calcolati per ciascun elemento di x.
√
Es.1 Sia N = 10, I = [0, 10] f (x) = x. Si assegni a=0, b=1, N=10 e si costruisca il vettore equispaziato x =
linspace(a, b, N+1), e da questo si costruisca il vettore
p f= sqrt(x). Si verifichi, scrivendo sul quaderno la tabella,
che i due vettori x ed f sono le righe di una tabella x, (x). Mediante l’istruzione
plot(x,f)
si tracci, usando i vettori cosı̀ costruiti, il grafico cartesiano della funzione radice quadrata.
1
Uno dei luoghi capitali della moderna scienza computazionale è il Courant Institute of Mathematical Sciences della
università di New York, intitolato al matematico Richard Courant, uno degli ideatori del metodo degli elementi finiti, tra i
maggiori progressi tecnici nel calcolo numerico. Nei programmi del Courant Institute, la teoria matematica si intreccia ai
risvolti applicativi nelle più diverse discipline per il tramite delle più avanzate tecniche computazionali. Notizie sull’istituto
sono reperibili in rete, all’indirizzo http://math.nyu.edu/.
1
2
Leggi orarie con MatLab
Descrivere il moto di un punto materiale significa individuare la posizione del punto nello spazio ad
ogni istante di tempo. Questa descrizione è usualmente chiamata legge oraria. Il moto generale di un
punto materiale nello spazio si ottiene assegnando una dipendenza dal tempo del suo vettore posizione2
−
→
r :
−
legge oraria → →
r (t)
Rappresentazione computazionale di una legge oraria: caso 1D
Rappresentare (ai fini del calcolo numerico) una legge oraria x(t) in una dimensione (ovvero di un
moto rettilineo) significa assegnare due vettori:
1. t = (t1 , t2 , ..., ti , ..., tn ) che elenca i tempi a cui si conosce la posizione del punto mobile.
2. x = (x1 , x2 , ..., xi , ..., xn ) che elenca la posizione (espressa come coordinata rispetto ad un’origine arbitraria) occupata dal corpo agli istanti elencati in t.
Si osservi che questa rappresentazione è anche quella che si otterrebbe se sottoponessimo la posizione
del punto mobile ad una operazione di misura. Per quanto rapida sia una singola misurazione
essa richiede un certo tempo Tc (detto tempo di campionamento) per cui la misura successiva
corrisponderà ad un tempo separato dal precedente all’incirca da Tc . L’insieme di dati sperimentali ( o
campioni) sarà quindi ordinabile nelle colonne di una tabella, ossia in vettori.
Es.2 Si consideri la legge oraria del moto armonico:
x = x0 cos(ωt) +
v0
sin(ωt)
ω
Si scriva una function che, avendo come input la frequenza del moto ω, la posizione x0 e la velocità v0 iniziali restituisca
i vettori x e t che rappresentano la legge oraria in un numero di punti N , anche esso dato come input, per un intervallo
di tempo che comprende circa 5 oscillazioni.
Es.3 Con i vettori x e t calcolati nell’esercizio precedente si tracci il grafico orario del moto (x vs. t)
Velocità
→
Come è noto la definizione di velocità si fa a partire dalla nozione di posizione −
r (t) tramite l’operazione
di derivata:
−
d→
r
−
→
v (t) =
dt
Se è nota la espressione analitica della legge oraria si possono applicare le regole di derivazione e
calcolare l’espressione analitica della velocità. Se, invece, la legge oraria è nota in forma computazionale,
ossia per valori discretizzati del tempo, tutto quello che si può fare è utilizzare i valori conosciuti per
approssimare il valore della derivata ossia della velocità. Nel caso unidimensionale note le coppie
(ti , xi ) del tempo e della posizione, la velocità al tempo corrispondente all’indice i si può approssimare
mediante:
xi+1 − xi
(1.1)
vi =
ti+1 − ti
Riducendo l’intervallo tra un istante ed il successivo, ossia aumentando il numero n dei campioni della
legge oraria, l’approssimazione diviene sempre migliore. Per eseguire l’operazione richiesta risulta utile
la funzione di MatLab
2
Si intende per vettore posizione del punto P la freccia che ha la coda in un punto O scelto come origine e la punta in
−−→
→
P. Si indica anche con −
r = OP .
A. Di Lieto, F. Maccarrone - Tecnologie Digitali I - 2008/09
3
diff(X)
la quale restituisce un vettore che ha come componente i-sima la differenza X(i+1)-X(i). Con questa
funzione, se x e t sono i vettori che rappresentano la legge oraria il vettore velocità calcolato applicando
la formula (1.1), si calcola mediante:
v = diff(x)./diff(t)
(1.2)
Es.4 Se x e t hanno N elementi, quanti elementi ha il vettore v calcolato mediante la (1.2)?
Es.5 Si modifichi la funzione dell’es. 2 in modo che oltre ad x e t restituisca anche i vettori v e v approx; il primo di
questi sia calcolato in corrispondenza degli stessi valori di t usati per la x(t) eseguendo direttamente la derivata analitica
della legge oraria:
vx = −ωx0 sin(ωt) + v0 cos(ωt)
e l’altro applicando la formula (1.2). Con i vettori t v e v approx si tracci il grafico cartesiano che confronta il valore
esatto della velocità con quello approssimato. Cosa si osserva al variare del numero N di elementi? In quali istanti la
approssimazione è meno accurata? Quale spiegazione daresti?
Integrali
Consideriamo una funzione f (x) senza particolari patologie, ossia che possa essere calcolata in tutti i
punti di un intervallo I = [a, b]. Esistono vari modi di definire l’integrale definito di f(x) tra a e b
che si indica con la scrittura:
Z
b
I=
f (x)dx
a
Una delle interpretazioni di questa operazione consiste nel calcolo dell’area compresa tra le rette x = a,
x = b, l’asse delle ascisse ed il grafico della curva (vedi figura). Una approssimazione dell’area si può
ottenere con il metodo dei rettangoli, ovvero suddividendo l’intervallo I in sottointervalli uguali di
ampiezza dx mediante i punti x1 = a, x2 = a + dx, ..., xN +1 = b e calcolando l’area dei rettangoli che
hanno altezza uguale al valore della funzione nell’i-simo intervallino (ad esempio nell’estremo sinistro
xi dove vale fi = f (xi )) e base uguale alla spaziatura dx. Evidentemente la somma approssima l’area a
meno dei triangoli curvilinei che sono sopra o sotto la curva a seconda che nell’intervallino la funzione
sia crescente o decrescente3 . In simboli l’operazione è:
I'
N
X
fi ∗ dx
(1.3)
i=1
3
La funzione può anche essere in parte sopra e in parte sotto. Si faccia un disegno in cui ciò avviene. Discutere
graficamente l’approssimazione cercando di illustrare qualitativamente in quali casi sia migliore o peggiore.
A. Di Lieto, F. Maccarrone - Tecnologie Digitali I - 2008/09
4
Es.6 Si scriva una function che, avendo come input il numero N di sottointervalli e gli estremi di integrazione a e b,
definisca due vettori: x equispaziato tra a e b con N + 1 elementi e f che contiene i valori della funzione f = x3 calcolati
nei punti di x. Mediante una opportuna istruzione si eseguano le operazioni espresse dalla eq. (1.3) e si assegni il risultato
alla variabile I approx che deve essere restituita dalla function. Si ricordi che in MatLab l’istruzione S = sum (x)
restituisce la somma delle componenti del vettore x. Il risultato è una approssimazione dell’integrale
Z b
b4 − a4
x3 dx =
4
a
Si consideri il calcolo fatto con a = 0, b = 1 e diversi valori di N . Descrivere qualitativamente come varia il valore approssimato al crescere di N . Quali sono i valori di N affinché l’approssimazione differisca dal valore esatto rispettivamente
alla terza o alla sesta cifra decimale?
Es.7 Cerchiamo di rendere quantitativa la dipendenza dello scarto tra l’approssimazione I approx ed il risultato esatto
I ex = 1/4. Per valori crescenti di N si annoti in una tabella a tre colonne il valore di N, della ampiezza dx degli elementi
della partizione di I e dello scarto delta I = abs(I approx - I ex). Dopo avere raccolto un numero sufficiente di dati
si faccia il grafico di delta I in funzione di dx. Cosa si osserva?
Il metodo dei rettangoli è estremamente semplice, ma non è particolarmente efficiente e preciso. Sono
possibili molti altri metodi che richiedono il calcolo della funzione in un insieme di punti (anche non
necessariamente equispaziati) dell’intervallo di integrazione I. Un esempio è il metodo dei trapezi.
Esso consiste nel approssimare l’area tracciando il segmento che unisce gli estremi di ogni sottointervallo.
In questo modo per ogni sottointervallo c’è da calcolare l’area di un trapezio rettangolo (vedi figura),
da cui il nome. Usando la formula per il calcolo dell’area di un trapezio le cui basi sono i valori della
funzione f agli estremi dell’intervallo e l’altezza è l’ampiezza del sottointervallo, il metodo dei trapezi
si traduce nella relazione:
N
X
1
I'
(fi + fi+1 ) ∗ dx
(1.4)
2
i=1
Es.8 Si ripeta l’esercizio 6 usando la regola dei trapezi. Si assegni il risultato alla variabile I2 approx. Quale è il valore
di N affinché l’approssimazione differisca dal valore esatto alla terza cifra decimale? E alla sesta? Trarre una conclusione
sensata, tra le molte che si potrebbero trovare, dal confronto dei due metodi.
Es.9 Si ripeta l’esercizio 7 usando la regola dei trapezi.
Es.10 Si tracci sullo stesso grafico l’andamento dello scarto della approssimazione dal risultato esatto in funzione di N
per i due metodi. Quale conclusione si trae?
Fast food
Come per gran parte dei problemi numerici MatLab si preoccupa di fornire funzioni che eseguano il
lavoro richiesto. Ciò che si deve fare è imparare la sintassi di tali funzioni. Naturalmente è altamente
sconsigliato per un fisico utilizzare le funzioni a scatola chiusa in particolare se dal risultato del calcolo
dipende una conclusione interessante. Nei casi in cui si conosce il risultato esatto può essere tuttavia
utile fare esercizio con le funzioni predefinite valutando l’attendibilità delle risposte. Per l’integrazione
numerica MatLab fornisce varie funzioni. Esaminiamone una la cui sintassi minimale è:
quad(FUN,a,b)
dove a e b sono gli estremi di integrazione e FUN è la funzione da integrare. La funzione FUN può essere
definita in forma simbolica tramite l’istruzione
FUN = inline(’f(x)’,’x’)
A. Di Lieto, F. Maccarrone - Tecnologie Digitali I - 2008/09
5
essendo f(x) l’espressione di MatLab che calcola la funzione da integrare, considerando l’argomento
x come vettoriale e quindi usando il . per indicare le operazioni componente per componente quando
necessario. Per illustrare il funzionamento si esegua il seguente esercizio:
Es.11
Si debba calcolare l’integrale
Z
1
I=
0
2x
dx
1 + x2
il cui risultato esatto4 è I = ln(2)..Si scriva uno script composto di quattro istruzioni. Le prime due devono definire gli
estremi di integrazione a = 0 e b = 1, la terza definisca la funzione F tramite F = inline(’2*x./(1+x.∧2)’,’x’);
e l’ultima l’integrale
I=quad(F,a,b); . Si confronti il risultato dello script con quello esatto.
Grafici orari in 2D
→
Rappresentare ai fini di calcolo numerico una legge oraria −
r (t) in due dimensioni significa assegnare:
1. Il vettore t = (t1 , t2 , ..., ti , ..., tn ) che elenca i tempi a cui si conosce la posizione del punto mobile.
2. I due vettori x = (x1 , x2 , ..., xi , ..., xn ) e y = (y1 , y2 , ..., yi , ..., yn ) che individuano la posizione
(espressa, ad esempio, come coordinate cartesiane rispetto ad un’origine arbitraria di una coppia
di assi ortogonali) occupata dal corpo agli istanti elencati in t.
Es.12
Si consideri la legge oraria del moto di caduta libera (uniformemente accelerato in direzione −y) :
1
x = x0 + v0x t; y = y0 + v0y t − gt2
2
Si scriva una function che, avendo come input il valore g del campo di gravità, le coordinate x0 e y0 e le componenti
della velocità v0x e v0y iniziali, tracci i grafici orari del moto (x vs. t e y vs. t) uno sotto l’altro (servirsi della funzione
subplot) nella stessa finestra grafica.
Animazione grafica
L’istruzione O = PATCH(X,Y,C) ha come argomenti due vettori di lunghezza uguale k e disegna in 2D
un poligono avente per vertici consecutivi i punti (X[i], Y[i]). Il parametro C può essere una lettera
che indica il colore della faccia del poligono, ad esempio ’r’ per il rosso, . . . 5 .
Es.13 Il centro di un poligono regolare di k vertici ha coordinate x0 , y0 ed il raggio del cerchio circoscritto C sia r.
Scegliamo una coppia di assi coordinati cartesiani con l’origine in C e si orienti l’asse x in modo che uno dei vertici sia
nel punto P1 = (r, 0). Muovendosi in senso antiorario su C si incontra il vertice adiacente a P1 quando si è percorso un
angolo:
2π
θ2 =
k
. Le coordinate cartesiane di P2 sono:
x2 = r cos θ2
y2 = r sin θ2
Analogamente si vede che lo i-simo vertice ha coordinate.
xi = r cos θi
con θi =
.
yi = r sin θi
i−1
2π. Il vettore degli angoli dei vertici si potrebbe assegnare in MatLab con teta = [0:2*pi/k:(k-1)/h*2*pi
k
Si scriva una function che avendo come input le coordinate del centro del poligono, il numero di vertici ed il suo raggio
restituisca le coordinate dei suoi vertici.
4
5
La primitiva della funzione integranda è, a meno di una costante F (x) = ln(1 + x2 )
In generale, C può essere una terna RGB e questo permette una maggiore scelta per il colore della faccia del poligono.
A. Di Lieto, F. Maccarrone - Tecnologie Digitali I - 2008/09
6
Es.14 Si prepari una function in cui si usano gli assi predefiniti (disponibili con l’istruzione axis) e si disegni un
poligono regolare di raggio r = 1/100 della lunghezza degli assi, con k = 50 vertici intorno al punto (x, y) = (0.5, 0.5),
facendosi restituire l’handle O all’oggetto grafico.
Un oggetto grafico è dotato di molte proprietà che ne rendono possibile la restituzione, ad esempio
sullo schermo corrente. In MatLab l’accesso a queste proprietà avviene attraverso la funzione set che
ha la seguente sintassi minimale:
SET(H,’PropertyName’,PropertyValue)
Qui H è l’handle all’oggetto, ’PropertyName’ è una parola, indicata tra apici, che denota la proprietà
sulla quale si vuole intervenire e PropertyValue il nuovo valore6 che deve avere la proprietà.
Es.15 Usando, dalla finestra dei comandi la funzione set applicata all’oggetto generato nell’esercizio precedente si
modifichino le proprietà XData ed YData, che contengono i vettori dei vertici del poligono traslando i punti secondo il
−
→
vettore T = (0.2, 0.2)
Es.16 Mettere insieme gli esercizi 2, 14 e 15 già svolti in un unica funzione. Si disegni il poligono nell’origine e poi,
rendendo generici i parametri Tx e Ty dell’esercizio precedente, si inserisca l’istruzione set in un ciclo che in successione
renda la coordinata x del centro del poligono quelle del punto mobile calcolate nell’esercizio 2. L’istruzione set va fatta
seguire dall’istruzione drawnow che forza il programma a applicare le modifiche al grafico immediatamente7 . Si dovrebbe
vedere una animazione del moto oscillatorio.
Traiettoria in 2D
Si può disegnare un nuovo oggetto grafico sovrapposto ad uno già esistente impiegando l’istruzione
hold. Questa modifica una proprietà degli assi correnti che riguarda appunto la persistenza degli oggetti
presenti in caso di un nuovo comando grafico applicato agli stessi assi. Normalmente la persistenza è
off e un nuovo comando grafico sostituisce quello precedente, sopprimendolo. Se si vuole mantenere
il vecchio oggetto grafico occorre modificare la proprietà facendola diventare on: questa istruzione ha
proprio questo funzione8 .
Si dice traiettoria del moto la curva che viene tracciata dal corpo in movimento. La sua rappresentazione in 2D sono i due vettori x e y, che rappresentano la funzione y(x) (o x(y)).
Es.17 Si scriva una function che, utilizzando quella scritta nell’esercizio precedente e modificata in modo che restituisca
i vettori x ed y, tracci con questi la traiettoria del moto. Mantenendo la curva tracciata (istruzione: hold) si ripeta
quanto fatto con l’esercizio 15, realizzando sul grafico della traiettoria una animazione del moto del proiettile.
6
Il tipo di valore, ad es. numerico, letterale, vettoriale etc., deve essere tra quelli legittimi della proprietà.
Senza questa istruzione le modifiche vengono applicate tutte alla fine della esecuzione e rimane solo l’ultima modifica.
8
In effetti la istruzione è di tipo toggle. Ciò significa che agisce come un pulsante che passa allo stato 1 se lo stato
attuale è 0 e allo stato 0 se quello attuale è 1. In altre parole se la luce è accesa il pulsante toggle la spegne, mentre se è
spenta, lo stesso pulsante la accende.
7
A. Di Lieto, F. Maccarrone - Tecnologie Digitali I - 2008/09