Ho pensato a questo esempio, partendo dalle equazioni di Navier-Stokes di fluidodinamica, applico lo stesso concetto al flusso di traffico dei veicoli. Le equazioni di partenza sono:
Come argomento simile, vedi anche: Python - Matplotlib - esempio idrodinamica, correnti fluviali 1D
Consideriamo quindi velocità e densità in questo modo: la velocità media delle auto in un dato segmento (unità di misura, ad esempio un km) e la densità è il numero di auto presenti in quel segmento.
L'andamento ci dice che per ρmax abbiamo Umin e per ρmin abbiamo Umax (ovviamente, come intuibile, strada vuota la velocità è maggiore, strada intasata la velocità è minore).
Per semplicità ho creato quindi un programma che inizializza random la densità sulla strada mentre per la velocità applica un andamento lineare fra Umin e Umax, come detto prima, in funzione della densità.
Dopodiché, per ogni step temporale (ho fatto 5 iterazioni), ho stamppato due grafici con visualizzazione Heatmap, in particolare "Reds" (da rosso chiaro a rosso intenso) per quanto riguarda la densità e "Greys" (da bianco a nero) per quanto riguarda la velocità. Ad ogni step temporale questo sistema di equazioni differenziali 1D è stato risolto in modo esplicito, partendo dalle condizioni iniziali. Ovviamente poi dal punto di vista matematico si potrebbe approfondire ulteriormente (scelta di un metodo numerico più efficiente, stabilità, ecc) in questo caso ho semplicemente riportato delle condizioni di controllo all'interno del ciclo per evitare che vada oltre gli estremi (ρmax, ρmin, Umax, Umin).
Ecco tutto il codice in Python.
from matplotlib import pyplot as plt
import numpy as np
import random as rd
N = 20
T = 5
rho_min = 1
rho_max = 10
u_min = 1
u_max = 10
Rho = np.zeros((T, N)) # Matrice per salvare i valori di densità per ogni istante di tempo
U = np.zeros((T, N)) # Matrice per salvare i valori di velocità per ogni istante di tempo
# INIT densità
for i in range(N):
Rho[0][i] = rd.randint(rho_min, rho_max)
# INIT velocità, andamento lineare U(Rho)=f^-1(Rho)
for i in range(N):
U[0][i] = u_min + rho_min * (u_max - u_min) / (rho_max - rho_min) * (rho_max / Rho[0][i] - 1)
# LOOP
for k in range(1, T): # Parti da 1 poiché il primo istante (k=0) è stato inizializzato
for i in range(N):
if i > 0:
Rho[k][i] += -1.0 / T * Rho[k - 1][i] * (U[k - 1][i] - U[k - 1][i - 1]) - 1.0 / T * U[k - 1][i] * (Rho[k - 1][i] - Rho[k - 1][i - 1])
if Rho[k][i]!=0:
U[k][i] += U[k - 1][i] ** 2 * 0.5 / Rho[k][i] * (Rho[k][i] - Rho[k - 1][i])
if U[k][i] > u_max:
U[k][i] = u_max
elif U[k][i] < u_min:
U[k][i] = u_min
if Rho[k][i] > rho_max:
Rho[k][i] = rho_max
elif Rho[k][i] < rho_min:
Rho[k][i] = rho_min
fig, ax = plt.subplots(T, 2, figsize=(10, 10)) # T grafici
for k in range(T):
ax[k][0].imshow([Rho[k]], cmap='Reds', aspect='auto') # Visualizza la densità
ax[k][0].set_yticks([]) # rimuove etichetta inutile
ax[k][1].imshow([U[k]], cmap='Greys', aspect='auto') # Visualizza la velocità
ax[k][1].set_yticks([]) # rimuove etichetta inutile
plt.suptitle("Densità (sinistra) e velocità (destra) per i vari istanti temporali")
plt.tight_layout() # Ottimizza gli spazi
plt.show()
Ecco infine l'output grafico, visualizzazione di densità e velocità tramite Heatmap.
Come detto questa è una semplificazione, esempio indicativo, che però a livello concettuale può essere interessante l'analogia fra il moto dei fluidi e un modello per il traffico dei veicoli. Variando le condizioni e i dati di progetto, si studia poi un diverso andamento.
