Numpy - array multidimensionali¶
Nella lezioni precedenti abbiamo visto alcune nozioni base su numpy (ricordo che un tutorial base lo trovate qui).
- lavora per array
- gli array sono classi.
- somme, moltiplizazioni, funzioni lavorano sui singoli elementi dell'array, facendo dei cicli for impliciti.
- Permette da fare algebra lineare.
- È veloce.
In questa lezione approfondiamo un po' come funziona numpy, introducendo gli array multidimensionali
import numpy as np
import matplotlib.pyplot as plt
Array multidimensionali¶
Nelle lezioni precedenti aveamo visto come lavorare con array unidimensionali, in cui i singoli elementi erano numeri, in virgola mobile, o interi. Oggi introduciamo gli array multidimensionali, cioè ``vettori'' i cui singoli elementi possono essere altri vettori, matrici, o elementi di dimensione più ampia. Il modo di chiamare funzioni e operatori (cioè l'interfaccia) rimane lo stesso, ma bisogna afferrare alcuni concetti fondamentali su come sono organizzati i dati.
# gli array possono essere
# multidimensionali
a = np.array([0,1,2])
b = np.array([1,0,0])
c = np.array([0,-2,1])
print(f" array 1D: a={a},b={b},c={c}")
array 1D: a=[0 1 2],b=[1 0 0],c=[ 0 -2 1]
Possiamo costruire degli array 2D semplicemente chiamando np.array con una lista di array 1D.
Attenzione: devono essere tutti della stessa lunghezza!
ab = np.array([a,b])
ac = np.array([a,c])
print(f" array 2D:\n ab={ab},\n ac={ac}")
array 2D: ab=[[0 1 2] [1 0 0]], ac=[[ 0 1 2] [ 0 -2 1]]
no = np.array([[1,2],[3,4,5]]) # non funziona perché gli elementi [1,2] e [3,4,5] hanno lunghezza diversa
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[4], line 1 ----> 1 no = np.array([[1,2],[3,4,5]]) # non funziona perché gli elementi [1,2] e [3,4,5] hanno lunghezza diversa ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (2,) + inhomogeneous part.
Accedere agli elementi¶
Si accede agli elementi in 2 maniere diverse.
- come per le liste:
ab[0] = [0,1,2]eab[0][-1] = 2 - Usando lo stile di matlab, che permette di recuperare righe e colonne per intero:
ab[:,0] = [0,1]
Questo secondo stile è utile perché spesso i nostri array multidimensionali rappresentano serie di dati. Ad ogni riga associamo un indice temporale e le colonne sono i nostri osservabili. Da una misura vogliamo poi in genere verificare una qualche relazione tra le colonne, i.e. che una funzione abbia un certo valore o comportamento.
ab[:,0] + ac[:,-1] # sommare colonne di 2 array
array([2, 2])
ab[0,:] + ac[1,:] # sommare righe
array([ 0, -1, 3])
ab[:,:] # restituisce tutto l'array
array([[0, 1, 2],
[1, 0, 0]])
Vediamo un esempio: verificare che l'energia sia conservata in una simulazione¶
Abbiamo simulato un oscillatore armonico in assenza di attrito, salvando i dati sul file oscillatore_m1_K10.dat. I dati contengono posizione angolare $\theta(t)$ e velocità angolare $\omega(t)$. Se la simulazione è corretta, l'energia
$$E(t) = \frac{1}{2}m v^2 + \frac{1}{2} K x^2 $$
è conservata, cioè assume un valore costante.
Come facciamo a verificarlo?
- Carichiamo i dati con
np.loadtxt - Usiamo quanto visto sopra per calcolare un array che corrisponda all'energia $E(t)$.
- Plottiamo $E(t)$ per vedere che non abbia derive o andamenti regolari.
- Ci calcoliamo la media e la varianza di questa, e l'errore massimo.
data = np.loadtxt('./data/oscillatore_m1_K10.dat')
T = data[:,0] # tempi
X = data[:,1] # posizioni
V = data[:,2] # velocità
Prima di tutto diamo un'occhiata ai dati.¶
L'energia sembra conservata
f1,a1 = plt.subplots(1,2,figsize=(10,5))
a1[0].plot(T[:1000],X[:1000])
a1[0].set_xlabel('tempo')
a1[0].set_ylabel('posizione')
a1[1].plot(X[:1000],V[:1000])
a1[1].set_xlabel('posizione')
a1[1].set_ylabel('velocità')
Text(0, 0.5, 'velocità')
K = 10 # costante della molla
m = 1 # massa
In Python puro¶
Posso calcolare l'energia in Python, ma questo rallenta molto l'esecuzione ed è meno elegante. Lo faccio vedere per far capire cosa fa numpy dietro le quinte. Il comando "magico" %%timeit ci permette di stimare quanto tempo impiega l'esecuzione della cella. È un comando di Jupyter, non di python.
%%time
E = np.zeros(T.size) # creo un array di zeri
for i in range(T.size):
E[i] = 0.5*( m * V[i]**2 + K * X[i]**2)
CPU times: user 419 ms, sys: 3.47 ms, total: 423 ms Wall time: 423 ms
Ragiono per array¶
Ora, le mie colonne rappresentano, tempo, posizione e velocità. Vogliamo trovare l'energia ad ogni istante di tempo per cui abbiamo X e V, e sappiamo che gli operatori e le funzioni di numpy lavorano direttamente sugli array. Vediamo come fare in numpy
%%time
Enp = 0.5*(m*V**2 + K*X**2)
CPU times: user 4.77 ms, sys: 1.72 ms, total: 6.49 ms Wall time: 7.82 ms
Siamo passati da 388 ms a 3.94 ms, un fattore 100 di differenza nella velocità di esecuzione!
Vediamo se il risultato cambia.
(E - Enp).sum() # somma tutti i valori delle differenze
np.float64(-6.412648190234904e-13)
Il numero che troviamo dalla somma è sostanzialmente un errore dei floating point, il che ci ricorda come non sia il caso di usare l'operatore == lavorando con i numeri in virgola mobile.
f2,a2 = plt.subplots(1,1,figsize=(6,5))
a2.plot(T[:1000],E[:1000])
a2.set_xlabel('tempo')
a2.set_ylabel('energia')
Text(0, 0.5, 'energia')
?? L'energia non è esattamente costante, perché?¶
Ho due possibilità per questo errore. Quali sono? Se aprite il file di testo ve ne potete rendere conto.
Capire come sono organizzati gli array¶
Gli array di numpy possono essere di varie dimensioni. Raramente ne userete più di 3 o 4, ma è molto utile capire come sono organizzati e cosa sono gli "assi". L'immagine qui sotto dovrebbe essere esaustiva.
- Gli array sono salvati in celle contigue di memoria, quindi potete rappresentare solo forme "omogenee": gli elementi devono avere tutti le stesse dimensioni o forma (e.g. vettori $n$-dim, matrici $n\times m$, etc.)
- La forma si può recuperare tramite l'attributo
shape. - Gli array si possono scorrere lungo diversi assi. Lo zero è sempre quello che contiene gli elementi più "grandi", poi segue l'uno e così via.
È fondamentale aver chiaro come funzionano gli assi perché molte funzioni di numpy operano di default su tutti gli elementi, ma a noi spesso interessa applicarli solo ad un asse, come lo zero, e possiamo farlo passando un argomento opzionale a questi metodi. I metodi mean(), sum(), var(),... sono ottimi esempi in questo senso.
Esempio: posizione e velocità medie¶
Diciamo che vogliamo ottenere la posizione e la velocità medie nel nostro esempio di prima.
$$ \bar{X} = \frac{1}{N}\sum_{i=0}^N X(t_i) $$
Se uso data.mean() ottengo la media di tutto. Se invece chiamo data.mean(axis=0) ottengo la media colonna per colonna.
data.mean(axis=0) # l'ultima colonna è l'energia
array([ 5.00000000e+03, -1.40819546e-05, -1.12533021e-05, 5.00000000e+00])
data.std(axis=0)
array([2.88675423e+03, 7.07102383e-01, 2.23608189e+00, 0.00000000e+00])
Calcolare proprietà di un insieme di dati.¶
In fisica spesso ci troviamo ad analizzare serie dati riguardanti più osservabili, ad esempio le coordinate di un insieme di particelle interagenti tra loro, in funzione del tempo. Da queste vogliamo calcolare degli osservabili e studiare la loro evoluzione. Se usiamo numpy, questo tipo di sistema corrisponde ad un array tridimensionale: tempo, particelle, assi del sistema di riferimento. Se raccogliamo $T$ misure di $n$ particelle in 2 dimensioni, questo array avrà shape (T,n,2). Guardiamo ad esempio la figura riportata sopra. Ogni elemento dell'asse zero è un sotto-array di forma $(n=4,2)$, rappresentante le posizioni delle particelle in quell'istante di tempo. Se fissiamo anche l'asse 1, otteniamo un array con 2 elementi, le coordinate di una data particella. Fissare l'asse 2 vuol dire scegliere anche la coordinata $x$ o $y$ che ci interessa.
È importante notare che possiamo fissare uno o più assi a piacere, senza dover partire dall'asse zero. Ad esempio, se fissiamo il primo elemento dell'asse 1 e lasciamo liberi gli altri indici otteniamo un array di forma $(T=4,2)$ rappresentante l'evoluzione della prima particella.
Quindi
- fisso asse 0: fisso l'istante di tempo.
- fisso asse 1: fisso una particella.
- fisso asse 2: fisso una componente $x,y,z$.
Consideriamo per esempio un array the tiene conto della posizione $(x_i(t),y_i(t),z_i(t))$ di tre particelle di massa uguale a diversi istanti di tempo. Per semplicità diciamo che tutte e tre le particelle si muovono di moto uniforme sul piano $z=0$, con velocità parallela all'asse $y$ e partendo rispettivamente da $x_1(0) = -1$, $x_2(0) = 0$, $x_3(0) = 1$. Il baricentro del sistema quindi avrà coordinate $\vec{x}_{CM}(t) = (0,y_{CM}(t),0)$. Vediamo qui sotto come calcolare il baricentro del sistema ad ogni istante, la posizione di una data particella, e la proiezione del sistema su una data coordinata.
# costruisco il sistema desiderato: la coordinata x è la stessa ad ogni istante per tutte le particelle
# la coordinata y cambia, e la z è sempre zero.
T = 50
y_pos = np.linspace(0,10,T)
# np.column_stack crea una matrice le cui colonne sono date dagli array passati nella lista
# np.repeat(v,T) crea un array contenente T volte il valore v
# np.zeros(T) equivale a np.repeat(0,T)
p1 = np.column_stack([np.repeat(-1,T),y_pos,np.zeros(T)])
p2 = np.column_stack([np.repeat(0,T),y_pos,np.zeros(T)])
p3 = np.column_stack([np.repeat(1,T),y_pos,np.zeros(T)])
# metto assieme tutte le particelle in una serie dati
points = np.array([[p1[t],p2[t],p3[t]] for t in range(T)])
La serie dati points ora contiene T elementi ognuno dei quali contiene 3 array di 3 elementi, cioè le posizioni delle particelle. Verifichiamolo con shape.
print(points.shape)
(50, 3, 3)
Seguire una particella¶
Se fissiamo il secondo indice (cioè l'asse 1) troviamo le singole particelle! Verifichiamolo: se funziona la differenza tra gli array così ottenuti e p1, p2, p3 deve essere nulla.
q1 = points[:,0] # gli indici partono da zero!
q2 = points[:,1]
q3 = points[:,2]
print(f"Verifichiamo le differenze.\nPrima particella:\t{(q1 - p1).sum()}\nSeconda particella:\t{(q2-p2).sum()}\nTerza particella:\t{(q3-p3).sum()}")
Verifichiamo le differenze. Prima particella: 0.0 Seconda particella: 0.0 Terza particella: 0.0
Il plot dei tre array ci dà giustamente 3 linee verticali
f3,a3 = plt.subplots(figsize=(5,5))
a3.plot(q1[:,0],q1[:,1],label='p1')
a3.plot(points[:,1,0],points[:,1,1],label='p2') # equivalente a a3.plot(q2[:,0],q2[:,1],label='p2')
a3.plot(q3[:,0],q3[:,1],label='p3')
a3.set_xlabel('X')
a3.set_ylabel('Y')
a3.legend()
<matplotlib.legend.Legend at 0x11acebc80>
Posizioni del baricentro¶
La media sull'asse 1 mi calcola il baricentro del sistema ad ogni istante! Come atteso corrisponde con la linea centrale del grafico di sopra.
CM = points.mean(axis=1)
f4,a4 = plt.subplots(figsize=(5,5))
a4.plot(CM[:,0],CM[:,1],label='CM')
a4.set_xlabel('X')
a4.set_ylabel('Y')
a4.legend()
<matplotlib.legend.Legend at 0x11ad9fe60>
Proiezione del sistema lungo una coordinata¶
Se mi tengo solo un valore fisso dell'asse 2, ottengo una proiezione del mio sistema lungo una coordinata, ad esempio $y$.
points[:,:,1]
array([[ 0. , 0. , 0. ],
[ 0.20408163, 0.20408163, 0.20408163],
[ 0.40816327, 0.40816327, 0.40816327],
[ 0.6122449 , 0.6122449 , 0.6122449 ],
[ 0.81632653, 0.81632653, 0.81632653],
[ 1.02040816, 1.02040816, 1.02040816],
[ 1.2244898 , 1.2244898 , 1.2244898 ],
[ 1.42857143, 1.42857143, 1.42857143],
[ 1.63265306, 1.63265306, 1.63265306],
[ 1.83673469, 1.83673469, 1.83673469],
[ 2.04081633, 2.04081633, 2.04081633],
[ 2.24489796, 2.24489796, 2.24489796],
[ 2.44897959, 2.44897959, 2.44897959],
[ 2.65306122, 2.65306122, 2.65306122],
[ 2.85714286, 2.85714286, 2.85714286],
[ 3.06122449, 3.06122449, 3.06122449],
[ 3.26530612, 3.26530612, 3.26530612],
[ 3.46938776, 3.46938776, 3.46938776],
[ 3.67346939, 3.67346939, 3.67346939],
[ 3.87755102, 3.87755102, 3.87755102],
[ 4.08163265, 4.08163265, 4.08163265],
[ 4.28571429, 4.28571429, 4.28571429],
[ 4.48979592, 4.48979592, 4.48979592],
[ 4.69387755, 4.69387755, 4.69387755],
[ 4.89795918, 4.89795918, 4.89795918],
[ 5.10204082, 5.10204082, 5.10204082],
[ 5.30612245, 5.30612245, 5.30612245],
[ 5.51020408, 5.51020408, 5.51020408],
[ 5.71428571, 5.71428571, 5.71428571],
[ 5.91836735, 5.91836735, 5.91836735],
[ 6.12244898, 6.12244898, 6.12244898],
[ 6.32653061, 6.32653061, 6.32653061],
[ 6.53061224, 6.53061224, 6.53061224],
[ 6.73469388, 6.73469388, 6.73469388],
[ 6.93877551, 6.93877551, 6.93877551],
[ 7.14285714, 7.14285714, 7.14285714],
[ 7.34693878, 7.34693878, 7.34693878],
[ 7.55102041, 7.55102041, 7.55102041],
[ 7.75510204, 7.75510204, 7.75510204],
[ 7.95918367, 7.95918367, 7.95918367],
[ 8.16326531, 8.16326531, 8.16326531],
[ 8.36734694, 8.36734694, 8.36734694],
[ 8.57142857, 8.57142857, 8.57142857],
[ 8.7755102 , 8.7755102 , 8.7755102 ],
[ 8.97959184, 8.97959184, 8.97959184],
[ 9.18367347, 9.18367347, 9.18367347],
[ 9.3877551 , 9.3877551 , 9.3877551 ],
[ 9.59183673, 9.59183673, 9.59183673],
[ 9.79591837, 9.79591837, 9.79591837],
[10. , 10. , 10. ]])
Esempio: sistema a 3 corpi¶
Nella cartella data trovate i file
butterfly1.datperiodic_scarab.datyingyang.dat
Che corrispondono a traiettorie del sistema dei 3 corpi in due dimensioni. In questo famoso sistema abbiamo 3 corpi (ad esempio 3 pianeti) che interagiscono gravitazionalmente. Per semplicità si considerano corpi di massa $m=1$. Le orbite sono caotiche e si possono calcolare solo numericamente.
Gli header dei file sono i seguenti:
#format: time coords (X,Y), velocities (X,Y), accelerations (X,Y)
Dove le colonne si ripetono per il corpo 1, 2, e 3:
t p_1x p_1y p_2x p_2y p_3x p_3y v_1x v_1y v_2x v_2y ... etc
Vogliamo leggere il file, fare il plot delle orbite riportato qui sotto, e verificare poi che il centro di massa del sistema resta stabile.
Come facciamo?
#1. Carichiamo i dati da file
data3b_1 = np.loadtxt('./data/butterfly1.dat')
#2. salviamo le posizioni
pos = data3b_1[:,1:7]
b1 = pos[:,:2] # corpo 1
b2 = pos[:,2:4] # corpo 2
b3 = pos[:,4:6] # corpo 3
f5,a5 = plt.subplots(figsize=(7,7))
a5.plot(b1[:,0],b1[:,1],label='b1')
a5.plot(b2[:,0],b2[:,1],label='b1')
a5.plot(b3[:,0],b3[:,1],label='b1')
a5.set_xlabel('X')
a5.set_ylabel('Y')
a5.set_title('butterfly orbit')
a5.legend()
f5.savefig('./images/butterfly_orbit.png')
Calcoliamo l'evoluzione del baricentro del sistema¶
Ora il nostro array è un po' diverso: ad ogni istante non abbiamo 3 array di 2 colonne, bensì un array di 6 colonne. Possiamo quindi usare due modi per ottenere il baricentro che vogliamo.
- fare la media esplicitamente con i 3 array ottenuti in precedenza;
- cambiare la forma del nostro array con
np.reshape(più avanzato).
Vediamoli entrambi
1. media esplicita¶
Sommo le posizioni delle particelle e divido per il loro numero. Sono solo 3, per cui è un conto rapido.
CM_3b = (b1+b2+b3) / 3
2. reshape¶
numpy permette di cambiare la forma degli array per rendere più comodi i calcoli. Questo si fa con il metodo reshape. La docstring permette di verificare facilmente come si usa in caso di dubbi. Nel nostro caso vorremmo passare da un array di forma (T,6) ad uno di forma (T,3,2). Basta fare pos.reshape(10000,3,2). La cosa importante è che la moltiplicazione dei valori in shape prima e dopo la trasformazione deve rimanere la stessa. Nel nostro caso T*6 = T*3*2.
Questo metodo di calcolare il baricentro fa automaticamente il ciclo sulle 3 particelle e quindi funziona anche con un numero qualsiasi di particelle nell'asse 1.
pos_r = pos.reshape(10000,3,2)
CM_3b_r = pos_r.mean(axis=1)
print((CM_3b - CM_3b_r).sum())
0.0
CM_3b.std?
Stabilità del baricentro¶
Nella nostra simulazione dei 3 corpi il baricentro al tempo $t=0$ sta nell'origine, quindi a meno di errori numerici ci aspettiamo che resti lì, non essendoci forze esterne ad agire sul sistema. Possiamo verificarlo calcolando la media e deviazione standard delle serie dati che corrispondono a $x_{CM}, y_{CM}, z_{CM}$. Questo si fa sempre usando axis, solo che questa volta ci serve axis=0. I valori che otteniamo sono compatibili con errori numerici sui floating point.
print(f"Posizione media CM:\t {CM_3b.mean(axis=0)}\nStd dev coord CM:\t {CM_3b.std(axis=0)}")
Posizione media CM: [2.44048864e-15 2.49533499e-15] Std dev coord CM: [1.67926931e-13 1.64621411e-13]