Vediamo questo esempio fisico-chimico, con programma e grafici creati con Python - Matplotlib. L'ho ritrovato come esercizio che avevo svolto in passato (ingegneria) e l'ho modificato leggermente.
- piscina, dimensioni 10x5x2 metri, quindi volume di 100 m3
- ad un'estremità viene versato un inquinante, concentrazione di carico organico 100.000 mg/l (valore altissimo, ad esempio i valori di BOD5/COD del percolato della discarica)
- diffusione costante pari a 103 m2/s (supponiamo quindi non ci siano altri fenomeni più complicati ed eterogenei per il mescolamento)
- equazione della diffusione:
dC/dt = D ∇^2(C)
, condizioni al contorno di flusso nullo sulle pareti (ovvero sono impermeabili) e soluzione numerica
- tempo 3600 secondi (un'ora)
- il metodo numerico che ho applicato è una variante di FTCS (prendendo spunto dall'estrapolazione di Richardson), passo di discretizzazione (adimensionale) pari a 0,1
Vediamo quindi cosa accade dopo un'ora, ho creato un grafico bello colorato tramite heatmap (rappresentando la situazione in superficie ovvero come la sostanza inquinante si diffonde). Prima il codice Python e poi l'immagine.
import numpy as np
from matplotlib import pyplot as plt
import time as timex
Tout=3600 #[s]
H=2 #[m]
L=10 #[m]
B=5 #[m]
dx=0.25 #[m]
delta=0.1 #[/]
D=0.001 #[m^2/s]
C0=100000 #[mg/l]
dt=dx**2*delta/D #[s]
M=np.zeros(((int)(L/dx),(int)(B/dx),(int)(H/dx))) #matrix
M[1][1][1]=1.0*C0 ##CC
#BCs: dC/dx=0, dC/dy, dC/dz=0
start=timex.time()
for time in range((int)(Tout/dt)):
for k in range((int)(H/dx)):
for j in range((int)(B/dx)):
for i in range((int)(L/dx)):
if(i>0 and i<(int)(L/dx)-1 and j>0 and j<(int)(B/dx)-1 and k>0 and k<(int)(H/dx)-1):
M[i][j][k]+=1.0*delta*(M[i+1][j][k]+M[i-1][j][k]+M[i][j+1][k]+M[i][j-1][k]+M[i][j][k+1]+M[i][j][k-1]-6.0*M[i][j][k])
temp=M[i][j][k]+2.0*delta*(M[i+1][j][k]+M[i-1][j][k]+M[i][j+1][k]+M[i][j-1][k]+M[i][j][k+1]+M[i][j][k-1]-6.0*M[i][j][k])
M[i][j][k]=2*M[i][j][k]-temp
elif i==0:
M[i][j][k]=M[i+1][j][k]
elif j==0:
M[i][j][k]=M[i][j+1][k]
elif k==0:
M[i][j][k]=M[i][j][k+1]
elif i==(int)(L/dx)-1:
M[i][j][k]=M[i-1][j][k]
elif j==(int)(B/dx)-1:
M[i][j][k]=M[i][j-1][k]
elif k==(int)(H/dx)-1:
M[i][j][k]=M[i][j][k-1]
end=timex.time()
print("tempo esecuzione codice:",round(end-start,2),"s")
print("step temporali:",(int)(Tout/dt))
print("celle in X:",(int)(L/dx))
print("celle in Y:",(int)(B/dx))
print("celle in Z:",(int)(H/dx))
color='plasma'
plt.imshow(M[:][:][1],color)
plt.show()
L'immagine dovrebbe rendere l'idea della concentrazione, in base ai diversi colori. L'estremo in alto a sinistra, come vedete, è quello "compromesso", dove è avvenuta la contaminazione.
Le considerazioni che possiamo fare, sono:
- se doveste trovare una sostanza inquinante, meglio stare lontani che più passa il tempo e più tende a diffondersi 😁
- Python è poco efficiente come linguaggio di programmazione (ma è comodo se vogliamo gestire anche la parte grafica; vedi anche: Python + C assieme con la libreria CFFI: efficienza top)
Cosa ne pensate, vi piace? 🙂 L'immagine è un po' semplice ma rende l'idea. La visualizzazione di questo tipo si chiama Heatmap (mappa di calore) e spesso è più intuitiva rispetto ad un grafico a linee.