Ciao,

con questo post completo il  tutorial scritto dai prof. L.A. Barba, G.F. Forsyth, C. Cooper sulla simulazione della diffusione del calore in una dimensione tramite python. Al seguente link potete trovare il tutorial in originale.

Per chi conosce poco o nulla del Python ho aggiunto una nota introduttiva per commentare alcuni tipi base di questo linguaggio ed introdurre i grafici e le animazioni in Python.

Nota introduttiva: qualche appunto su Python

Prima di iniziare la descrizione del codice con cui si simula l’evoluzione temporale del profilo termico monodimensionale, credo sia utile ed interessante discutere qualche elemento sui dati e la loro gestione in Python.

Python come molti sapranno è un linguaggio di programmazione ad oggetti. E’ ampiamente diffuso in molti ambiti applicativi. Nel campo scientifico è largamente utilizzato per l’elaborazione di dati e le simulazioni, rappresenta infatti un’ottima alternativa a programmi proprietari e potenti come Matlab. Python infatti è liberamente scaricabile ed ampiamente supportato nella documentazione e negli esempi: un sito da cui scaricare l’installer è questo link o in alternativa su Anaconda .

I temi che vorrei introdurre  sono:

  1. Le liste in Python;
  2. Numpy ed array;  
  3. Grafica in Python;
  4. Animazioni in Python;

 

  1. Tipi di Dato elementari in Python: le liste

Liste: cosa fa venire in mente la parola lista? Un elenco di “cose” diverse, in una lista posso mettere insieme la spesa da fare (pane, frutta…), numeri di telefono, nomi di amici …. un pò di tutto. In Python le liste sono esattamente uguali: un elenco ordinato di elementi diversi: stringhe alfanumeriche, numeri, liste o altri tipi di dato. Un esempio:

>>> elenco=[“mamma”,2,3,”pane”] 

La lista viene definita dalla variabile (elenco) e dalle parentesi quadre [ e ], che racchiudono gli elementi della lista ordinati. Ogni elemento della lista è separato da una virgola e la posizione all’interno della lista è data da un indice che comincia con 0 per il primo elemento a sinistra e cresce fino all’ultimo elemento. Tramite la posizione dell’elemento nella lista posso estrarre il singolo elemento ad esempio se scrivo:

>>> elenco[0]  l’out sarà  >>> ‘mamma’

>>> elenco[1]  l’out sarà  >>> 2

Un’importante caratteristica delle liste in python è la possibilità di selezionare gruppi di elementi molto velocemente tramite lo slicing, che consiste nell’indicare, fra parentesi quadre, l’indice iniziale e finale (escluso), separati da :, gli elementi che si vogliono considerare…..più facile da mostrare che spiegare:  

>>> a=[2,3,5,6,1,34,11,7]   %ho definito la varibile a tipo lista

>>> a[1:3]   l’out sarà >>> [3,5]   %ho prelevato gli elementi dal numero 1 fino al 3 escluso% 

>>> a[4:-1]  l’out sarà >>>[1,34,11]   % gli elementi di una lista da sx a dx si contano da 0 in poi mentre da dx a sx si contano con numeri negativi da -1 in poi; quindi la scrittura a[4:-1] indica di prelevare la lista dalla posizione 4 fino all’ultimo elemento escluso%

>>>a[:]        l’out sarà >>>[2,3,5,6,1,34,11,7] %l’ntera lista%

>>>a[3:]     l’out sarà >>>[6,1,34,11,7] % dalla quarta posizione in poi%

Per copiare una lista  >>>b=a. E cosi analogamente se voglio estrarre una porzione della lista ‘a’, ad esempio: >>> b=a[2:4] l’out sarà >>> [5,6] 

Sulle liste sono poi definite una serie di operazioni, per conoscerle basta definire una lista e scrivere dir(a). Fra le operazioni più comode ed usate ricordo: 

append: aggiunge un elemento in coda (esempio: >>>a=[‘si’,’no’,3] ; >>> a.append[‘forse’] out >>> [‘si’,’no’,3,’forse’])

pop: estrae ed elimina un elemento dalla lista

Infine la List comprehension, una modalità che Python offre per definire liste tramite delle operazioni che agiscono su altre liste:

>>> num=[1,2,3,4,5]

>>> doppi=[n*2 for in num]  >>> l’out sarà doppi=[1,4,6,8,10]

Ai fini delle simulazioni numeriche le liste sono un elemento importante ma serve una struttura più specifica, pensata solo per funzionare da vettore n dimensionale. Questa strutture sono gli array definiti nella libreria Numpy.

2. Numpy ed array

Numpy è la libreria matematica di python. Come tutte le librerie va caricara con il comando import Numpy as np (np diventa un alias per dire a python che ho bisogno di richiamare alcuni metodi ed oggetti di numpy). In Numpy è definita l’oggetto array che consiste in una tabella di numeri che possono essere indirizzati da una n-upla di numeri interi, l’esempio dipanerà ogni dubbio:

>>>> import Numpy as np

>>> a=np.array([[1,2,3],[4,5,6]])  %ho creato l’array a, in questo caso una matrice 2×3 di interi

>>> a[0] l’out sarà >>> [1,2,3] % la riga 0, tutte le colonne

>>> a[0,2] l’out sarà>>> [3] %elemento di riga 0 e colonna 2

>>> a[:,1] l’out sarà >>> [2,5] %slicing sulle righe, in questo caso tutte le righe, e solo la colonna 1

Come appena visto sugli array di numpy è valido lo slicing visto per le liste. Un altro metodo di selezione degli elementi di una matrice, più generale, si ottiene con gli operatori booleani, ovvero degli operatori che mi forniscono una condizione di vero o falso. Ad esempio se nella matrice a, prima definita, cerco in generale gli elementi maggiori di 4 (a>4), otterrò un matrice di vero o falso così fatta, ([False,False,False],[False,True,True]). La matrice booleana mi dice che la priam riga di a non ha nessun elemento che soddisfa la condizione, (tutto falso) mentre la seconda riga ha il primo elemento che non soddisfa la condizione, gli altri due si. Se ora passo ad a questa matrice boolena, python selezionerà gli unici elementi true e li restituirà in una nuova matrice, per cui:

>>> a[a>4]   l’out sarà >>> ([5,6])

Ancora un modo per selezionare tramite condizioni booleane è dato dall’operatore where:

>>> np.where(a>4) l’out sarà >>> ([1,1],[1,2]) % cioè una matrice in cui riporto gli indici che rispondono alla condizione a>4

Su gli array di Numpy sono definite le operazioni di somma,differenza , moltiplicazione e divisione. Queste ultime sono fatte elemento per elemento mentre il prodotto righe per colonne si ottiene con l’operazione .dot(array) che è un metodo della classe array: in poche parole, definite a e b due matrice opportune il prodotto axb  si fa scrivendo a.dot(b).

Ultimo punto di cui vorrei parlare è come costruire un array in numpy. Il metodo più semplice lo abbiamo già visto ed è la definizione diretta di tutti gli elementi, come fatto per a. Se vogliamo matrici particolari abbiamo i metodi zeros, ones, full ed eye:

>>> np. zeros((2,2)) l’out sarà >>> ([[0,0][0,0]]) % matrice 2×2 di zeri

>>> np. ones((2,2)) l’out sarà >>> ([[1,1][1,1]]) % matrice 2×2 di uno

>>> np. full((2,2),7) l’out sarà >>> ([[7,7][7,7]]) % matrice 2×2 di sette

>>> np. eye((2,2)) l’out sarà >>> ([[1,0][0,1]]) % matrice diagonale 2×2

 Per creare array di numeri i metodi che possiamo usare sono, random, arandom, linspace:

>>> np.random.rand(2,2) l’out sarà >>> ([0.23,0.56],[0.87,0.92]) % matrice 2×2 di numeri random fra 0 e 1

>>> np.arange(5) l’out sarà >>> ([0,1,2,3]) % array dei primi cinque numeri interi

>>> np.linspace(1,6,5) l’out sarà >>> ([1,2.25,3.5,4.75,6]) % prende l’intervallo da  start (1 in questo caso) a stop (6 in questo caso) e lo divide in  n-1 parti uguali dove n in questo caso è 5; il risultato è una sequenza di numeri equamente distribuiti da start a a stop (inclusi)

3. Grafica in Python

Perché scrivere codice che fa grafici quando esistono già molti programmi che lo fanno per noi? Potrebbe non essere la scelta più corretta, però : potreste includere i grafici direttamente nei programmi di simulazione o di elaborazione che avete sviluppato, invece di esportare i dati ed importarli successivamente in altri software; analogamente potreste analizzare interattivamente i vostri programmi dai grafici di out, o configurare i grafici come meglio vi serve. L’uso di librerie grafiche infatti non è cosi raro. Per fare questo in Python la libreria più utilizzata è Matplotlib ed in particolare pyplot (quindi il primo comando è import Matplotlib.pyplot as plt). Importata questa libreria per ottenere un grafico è sufficiente avere due liste X e Y di uguale dimensione con i dati e metterle insieme con il comando plot : 

>>> X=[1,2,3,4]

>>> Y=[8,9,10,11]

>>> plt.plot(X,Y)

>>> plt.show()

Python configurerà automaticamente la finestra all’interno della  quale riportare il grafico e dimensionerà automaticamente gli assi.

grafico_python_esempio

Esistono molti altri tipi di grafici disponibili in python. Per una specifica informazione vi rimando alla documentazione in linea.

Se volessimo controllare le dimensioni della finestra in cui verrà visualizzato il grafico possiamo definire l’oggetto figure: plt.figure(figsize=(12,8),dpi=100,facecolor=’1.0′). Con questa chiamata abbiamo definito una finestra di dimensioni 12×8 pollici con una risoluzione di 100 dpi, ed un colore di background bianco (le alternative vanno da 0 a 1 passando da nero al bianco con attraverso diverse tonalità di grigio)

4. Animazione in Python

Se il grafico non bastasse ma volessimo mostrare “l’evoluzione” di un grafico, Python offre la possibilità di creare delle animazioni. Per realizzare delle animazioni occorre importare la libreria dedicata di Matplotlib: import Matplotlib.Animation.

L’idea di base per costruire l’animazione e’ definire una funzione che aggiorna una parte del grafico. In python infatti i grafici sono degli  oggetti Line2D che possiamo parzialmente modificare con il comando set_xdata() o set_ydata(). Vediamo passo passo:

>>> x=[1,2,3] % lista di X

>>> y=[4,5,6] % primo gruppo di dati da mettere sul grafico

>>> z=[4,8,16] % secondo set di dati da mettere sul grafico

>>> lines=plt.plot(x,y)

>>> lines=plt.plot(x,z) % il comando plot crea il grafico generando un oggetto Line2D che abbiamo chiamato lines; in questo caso lines è fatto da due elementi in quanto abbiamo chiamato due volte il comando plot: nel primo elemento risiede l’associazione X,Y nel secondo risiede l’associazione X,Z.  I singoli elementi sono indirizzabili come se fossero delle liste (ma non lo sono):

>>> lines[0].set_ydata([5,6,7]) %con questo comando ho modificato solo la il dato y della primo elemento di lines (quello fatto dalla coppia X,Y)

Allora per creare un’animazione, definiamo un grafico vuoto con il comando lines=plt.plot([],[]) e poi definiamo una funzione che indichi cosa modificare del grafico e come modificarla:

def update(i):

     lines[0].set_xdata(X[:i])

     lines[0].set_ydata(Y[:i])

return lines

La funzione Funcanimation userà la funzione di aggiornamento per modificare il grafico e mostrare la sequenza di frame prodotte:

animation.FuncAnimation(fig,update,interval=50,blit=true,frames=20,repeat=True)

 Dove ho indicato con fig la figura da aggiornare, con update cose e come aggiornare e con gli altri parametri la frequenza di aggiornameto.

Con questi elementi proviamo a completare la simulazione già impostata nei post precedenti.

Ed ora la soluzione!

Siamo pronti per il number -crunch!

Le seguenti parti di codice inizializzano il programma importando le librerie necessarie, quindi definiamo i parametri necessari alla simulazione e definiamo le condizioni iniziali. Questa volta l’utente non può scegliere un  generico Δt; abbiamo pensato che questa scelta sarebbe stata non sicura: si rischierebbe di far saltare tutto. Invece, abbiamo costruito un codice che calcola un valore di Δt che si trovi nella regione stabile del problema, in accordo con la discretizzazione spaziali scelta! Si può provare a cambiare i parametri scelti per vedere come cambia la soluzione numerica, ma nulla salterà per aria!

insert1

insert2

insert3

Animazione

Poter vedere come cambia il profilo termico fra l’inizio e la fine della simulazione è utile ma è ancora più interessante poter vedere il suo cambiamento in modo dinamico!

A questo scopo, per prima cosa, importiamo il modulo di matplotlib che gestisce l’animazione ed un modulo speciale di IPython chiamato HTML (…).

Nota:

Dobbiamo inoltre installare un video encoder/decoder chiamato ffmpwg.

Se utilizzate Linux o OSx lo si può installare tramite il comando conda:

conda install -c conda-forge ffmpeg

Se invece utilizzate Windows potete trovare le istruzioni per l’installazione al link:

http://adaptivesamples.com/how-to-install-ffmpeg-on-windows/

insert4

Ora occorre resettare le condizioni iniziali!

insert5

 

Ora proviamo a creare un’animazione. Per farlo serve quale passo in più nelle programmazione, ma non è difficile da fare! Innanzitutto, dobbiamo creare una figura ed aggiungere gli assi. L’oggetto figura che abbiamo creato definisce il grafico che vogliamo animare. Gli assi ci permettono di creare gli elementi del grafico che noi cambieremo interattivamente per creare l’animazione. In fine creiamo una linea “blank” che noi aggiorneremo ad ogni iterazione. Il codice per tutto questo è:

insert6

Nota: il valore [0] come argomento di ax.plot() è necessario in quanto ax.plot può in modo opzionale restituire diversi valori. Poichè stiamo lavorando solo con una linea siamo interessati solo alla prima (indicata dall’indice zero).

Dopo aver configurato la figura, possiamo definire una funzione che descrive come deve evolvere il nostro grafico. Questo codice è semplicemente lo stesso che abbiamo già utilizzato, solo che ora vogliamo aggiornare i dati da graficare.

insert7

Inizializziamo l’animazione utilizzando animatio.FuncAnimation() e salviamo l’oggetto che così creiamo nella variabile anim. Alla funzione FuncAnimation() dobbiamo passare i seguenti parametri:

  • fig: il nome della figura che vogliamo animare (aggiornare)
  • diffusion: il nome della funzione che materliamente calcola il profilo termico per ogni iterazione
  • frames: il numero di frames da rappresentare (che impostiamo uguale ad nt)
  • interval: il numero di millisecondi per cui appare una frame.

insert8

Ok! Ora è tempo di guardare la nostra animazione.

IMG_0943

 

CC BY-NC-SA 4.0
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.