Scuola di Specializzazione in Fisica Medica Introduzione a ... · A. Andreazza – Introduzione a...

Introduzione a MATLAB (A. Andreazza)

Scuola di Specializzazione in Fisica Medica

A. Andreazza – Introduzione a MATLAB

MATLAB

•  Linguaggio di programmazione ad alto livello per applicazioni scientifiche •  Ambiente di sviluppo interattivo •  Disponibilità di funzione matematich per un’ampia varietà di applicazioni •  Strumenti per grafica 2D e 3D •  Possibilità di costruire interfacce grafiche ed interfacciarsi con librerie in altri

linguaggi di programmazione. –  https://it.mathworks.com/products/matlab.html

•  Disponibile su Windows, Linux e Mac •  L’Università di Milano ha una licenza campus, per cui è accessibile a tutti gli

studenti –  Possibile installarlo sul proprio PC –  Istruzioni su: http://www.unimi.it/ateneo/80207.htm

SCUOLA DI SPECIALIZZAZIONE IN FISICA MEDICA 2


Avvio MATLAB

•  Aprire una finestra terminale •  (consigliato) Creare una directory per MATLAB

–  mkdir MATLAB •  Entrare nella directory

–  cd MATLAB •  Far partire MATLAB

–  matlab2015 & –  La versione usata in laboratorio è la 2015b


& fa partire il processo in background, permettendo di continuare ad usare il terminale

A. Andreazza – Introduzione a MATLAB SCUOLA DI SPECIALIZZAZIONE IN FISICA MEDICA 4

Pannello dei comandi

Pannello delle variabili Directory di lavoro


Primi comandi

•  Provare ad eseguire dei calcoli: –  costanti: 5*3 2^10 pi/2 –  funzioni: sin(pi/2) sind(90) exp(-1) –  e cosa succede se facciamo sqrt(-3)?

•  Definire della variabili (scalari) –  r = sqrt(2) –  C = 2*pi*r –  year = 2018

•  Cambiare il formato di stampa –  format long

•  Documentazione, cancellazione variabili e finestra –  help format –  clear –  clc


Costanti predefinite: π: pi unità immaginaria: i j 1i 1j

Variabile ans definita automaticamente

Per documentazione: in linea: help comando dettagliata: doc comando

Gioco 1 format HEX mostra la rappresentazione interna dei numeri. Per numeri reali, riuscite a capire la rappresentazione utilizzata? E per interi?

SCUOLA DI SPECIALIZZAZIONE IN FISICA MEDICA


Operazioni con vettori

•  Una delle forze di MATLAB è la potenza delle funzioni vettoriali •  Creazione di un vettore:

–  vettore riga: x = [ 1 2 5 4 ] oppure: x = [ 1, 2, 5, 4 ]

–  vettore colonna: y = [ 1; 2; 5; 4 ] oppure: y = x’

–  vettore equispaziato: x = inizio:passo:fine se passo è 1 si può omettere oppure: x = linspace(inizio,fine,numeropunti)

•  Molte operazioni vengono ripetute automaticamente su tutti gli elementi di un vettore: –  y = sqrt(2)*x –  z = y-x –  w = exp(x)


Costanti predefinite:

’: operatore di trasposizione N.B.: per vettori complessi fa il trasposto coniugato


Operazioni vettoriali

Operatore Descrizione

+ - somma e sottrazione elemento per elemento

* / prodotto e divisione per uno scalare prodotto e divisione tra matrici

^ esponenziazione (tra scalari e matrici quadrate)

.* ./ prodotto e divisione elemento per elemento

.^ esponenziazione elemento per elemento


•  Definito: –  x = 1:4

•  Confrontare i risultati di: –  x*x –  x*x’ –  x’*x –  x.*x

Esempio 1 Grafico di alcune funzioni Impariamo: •  utilizzo di indici per accedere ai dati •  funzione plot e proprietà di un grafico •  salvare istruzioni in una macro



Grafico di una funzione

•  Creare un vettore di ascisse: 101 punti tra -10 e 10 –  x = linspace(-10.,10.,101);

•  Calcolare un vettore di ordinate: –  y = sin(x)./x;

•  Verificare il valore del 51o elemento: –  y(51) corrisponde a x=0, ed il risultato è NaN = not a number –  y(51) = 1 assegnare il limite per x che tende a 0

•  Grafico di y(x) –  plot(x,y)

•  Titolo e assi –  title('Grafico di funzione') –  xlabel('x [rad]') –  ylim([-0.4,1.2]) –  ylabel('sin(x)/x')


Indice del vettore in (), inizia da 1

In caso di vettori molto grandi, ‘;’al termine del comando sopprime la stampa a schermo

Vettore di estremi dell’asse

Variabili di testo tra ‘...’


•  Sequenza di chiamata: –  plot(x,y,opzionibase,nomeopzione,valoreopzione,...)

•  Opzioni base (un singolo testo)

•  Opzioni avanzate:

Funzione plot

Carattere Colore Carattere Marker Carattere Linea

k black . punto - continua

r red * asterisco -- tratteggiata

g green o cerchio : punteggiata

b blue x croce :- dot-dash

c cyan

m magenta


‘LineWidth’ Spessore della linea, default: 1

‘MarkerFaceColor’ Colore di riempimento del marker, default: trasparente

‘MarkerSize’ Dimensione del marker, default: 10

Per documentazione: help plot


File macro

•  Si possono inserire i comandi in un file .m in modo da poterli richiamare più facilmente

•  Aprire l’editor di macro –  edit grafico

si apre un nuovo pannello ed il file grafico.m compare nella directory di lavoro

•  inserire nell’editor i comandi come sono dati •  salvare la macro (Save da GUI) •  eseguire la macro:

–  clear –  clc –  grafico


A. Andreazza – Introduzione a MATLAB SCUOLA DI SPECIALIZZAZIONE IN FISICA MEDICA 13

Commenti Il file appare nella directory

Esecuzione (GUI)

Esecuzione (linea di comando)


Altre funzioni

•  Per provare altre funzioni basta modificare le linee di definizione della y –  y = sqrt(2)*ones(1,length(x)); –  y = zeros(1,length(x));

•  Possiamo sfruttarle per vedere diversi modi di indicizzazione, usando vettori e condizioni logiche –  y = zeros(1,length(x)); y(51) = 1;

–  y = zeros(1,length(x)); y(26:76) = 1;

–  y = zeros(1,length(x)); y( x>0 ) = 1;

–  y = cos(x); y( y<0 ) = 0;


Funzione Delta

Funzione costante, length dà la lunghezza di un vettore ones dà una matrice di 1 zeros dà una matrice di 0

Rettangolo, vettore di indici

Gradino, vettore di elementi che soddisfano una condizione logica

Lobi positivi del coseno

Esempio 2 Generatori di numeri casuali Impariamo: •  generatori di distribuzioni uniformi e gaussiane •  funzioni statistiche •  istogrammi •  cicli for



Generatori di numeri casuali

•  MATLAB permette di riempire matrici con numeri pseudo-casuali: –  generatore uniforme in [0,1): ru = rand(1000,1); –  distribuzione normale: rn = randn(1000,1); –  interi equidistribuiti tra 1 e 6: ri = randi(6,1000,1);

•  Possiamo verificare visivamente che le distribuzioni siano quelle attese costruendo un istogramma: –  histogram(ru) –  è possibile anche dare un vettore che definisce gli estremi dei canali

dell’istogramma histogram(ru,linspace(0.,1.,21) per avere N canali servono N+1 punti

–  histogram(rn) –  histogram(ri,0.5:6.5)


Righe

Colonne


Funzioni statistiche

•  Possiamo verificare se i numeri generati possiedono le proprietà statistiche che ci attendiamo: –  valor medio: mean(ru) –  deviazione standard: std(ru)

•  I valori attesi sono: –  0.5 e 1/sqrt(12) per la distribuzione

uniforme –  0 e 1 pe la distribuzione normale –  e per gli interi equiprobabili tra 1 e 6? d = 1:6 mean(d) std(d)

•  Esistono molte altre funzioni che operano su un vettore (o sulle colonne di una matrice)


Funzione Descrizione

mean valor medio

std deviazione standard

var varianza

max valore massimo

min valore minimo

median mediana

sum somma degli elementi

prod prodotto degli elementi

diff vettore differenze tra el.

cumsum vettore somme cumlative


Lancio di più dadi

•  Come è fatta la distribuzione della somma di N dadi? –  il teorema del limite centrale dice che per N→∞ la forma tende ad essere

gaussiana –  2 dadi: ri = ri+randi(6,1000,1); histogram(ri) mean(ri) std(ri)

•  3 dadi: ri = ri+randi(6,1000,1); histogram(ri) mean(ri) std(ri)

•  Possiamo automatizzare?



Strutture di controllo

•  Cicli for for variabile=vettore

...codice in cui variabile assume i valori in vettore

end


•  Cicli while while condizione

...codice in cui può modificarsi la condizione

end

•  Esecuzione condizionata if condizione

codice se condizione è verificata

else if condizione2 codice se condizione2 è verificata

else codice se nessuna condizione è verificata

end


Lancio di più dadi

•  Costruiamo una macro sommadadi.m che guardi le distribuzioni fino a N dadi facendone gli istogrammi e mostrando l’andamento di valor medio e varianza.


Crea un nuovo pannello per i plot

I vettori vengono allungati a mano a mano che si inseriscono i dati

hold permette di sovrapporre diverse linee su uno stesso grafico

Esempio 3 Calcolo di fase ed ampiezza di una sinusoide Impariamo: •  matrici e loro indicizzazione •  operazioni tra matrici •  risoluzione di sistemi lineari



Dati con rumore

•  Supponiamo di avere una sinusoide misurata ad intervalli regolari, e con un certo errore di misura: –  x = (0:0.5:10)’; –  y = 1.5*sin(x+pi/4); –  yMis = y + 0.2*randn(length(y),1) –  sovrapporre i grafici di y e yMis

•  Vogliamo trovare i migliori parametri [a b] che descrivono i dati misurati:

•  Minimizzando il χ2:


y = Asin(x +φ) = asin x + bcos x

A = a2 + b2 , φ = tan−1 ba

χ 2 = (yi − asin xi − bcos xi )2

i∑


Fit di χ2

•  Il χ2:

•  può venire espresso in forma matriciale:

•  e la soluzione è data dal sistema


χ 2 = (yi − asin xi − bcos xi )2

i∑

χ 2 = (y −Ap)T (y −Ap)

y =

y1y2!yN

⎛

⎝

⎜⎜⎜⎜

⎞

⎠

⎟⎟⎟⎟

, p = ab( ), A =

sin x1sin x2!

sin xN

cos x1cos x2!

cos xN

⎛

⎝

⎜⎜⎜⎜

⎞

⎠

⎟⎟⎟⎟

ATy −ATAp = 0 p = ATA( )−1ATy


Fit di χ2

•  Soluzione in MATLAB: –  A = [ sin(x) cos(x)]; –  ATy = A’*yMis; –  ATA = A’*A; –  p = inv(ATA)*ATy –  ampiezza = sqrt(p’*p) –  fase = atan2(p(2),p(1))

•  Residui del fit: –  yFit = A*p; –  yRes = yMis-yFit; –  sigma2 = var(yRes)/(length(yRes)-1)

•  Matrice di covarianza dei parametri a e b:


Cp = ATA( )−1ATCyA ATA( )

−1Cy = σ y

21

Matrice da accostamento di colonne Usare ; per sovrapporre righe

Operazioni tra matrici

inv calcola la matrice inversa Per sistemi lineari esiste un metodo numericamente migliore: •  A*p=y ⇒ p = A\y •  p*A=y ⇒ p = y/A


Operazioni con matrici

•  Creazione: –  A = [ 1 2 3 ; 4 5 6 ; 7 8 9] –  A = eye(5) –  A = ones(Nrighe,Ncolonne) –  A = zeros(Nrighe,Ncolonne)

•  Indicizzazione: –  A(2,4) –  A(:,2) –  A(4,:)

•  Funzioni: –  det(A) –  inv(A) –  diag(A) –  size(A,dim)


; separa le righe

matrice identica 5x5

matrice di 1

matrice di 0

Elemento nella 2a riga 4a colonna

2a colonna

4a riga

determinante di A

Inversa di A

diagonale di A

lunghezza della dimensione di A


Fit di χ2

•  Completare l’esempio con: –  calcolo della matrice di covarianza di a e b –  calcolo della matrice di covarianza di A e φ

•  Vedere come si comportano le funzioni statistiche se applicate ad una matrice


ESERCIZI



Esercizio 1

•  Scrivere una macro che al suo interno: –  definisce un certo numero di punti su un intervallo –  calcola i valori di una funzione (definita all’interno della funzione)

su quei punti –  calcola numericamente l’integrale della funzione sull’intervallo



Esercizio 2

•  Prendere il file croce.mat: –  caricare i dati del file usando il comando load croce –  contiene due misure dei punti [x y] di una croce nelle matrici A e B –  fare la differenza delle misure e stimare la precisione in x e y –  la croce è ruotata: calcolare di quale angolo –  applicare una matrice di rotazione per “raddrizzarla”


Scuola di Specializzazione in Fisica Medica Introduzione a ... · A. Andreazza – Introduzione a...

Documents

Transcript of Scuola di Specializzazione in Fisica Medica Introduzione a ... · A. Andreazza – Introduzione a...