Dopo aver visto Andamento modello di traffico - Python Matplotlib, vediamo questo esempio: trasporto di un'onda quadra. Ovvero partendo da una situazione iniziale in cui i valori sono 0 e 1 (0 prima e dopo una certa soglia, 1 nel mezzo, in modo da formare questa "onda quadra"), il modello fisico presuppone di avere un semplice trasporto (con velocità costante U), quindi semplicemente x(t)=Ut. Quindi in termini differenziali, abbiamo:
Dopodiché, ho creato un programmino Python con diversi schemi numerici alle differenze finite a confronto e la soluzione analitica.
I metodi presi in in considerazione sono:
- soluzione analitica: una semplice traslazione rigida, appunto data da
x(t)=Ut
- metodo di Godunov, per U>0 (tra l'altro ho scoperto che il matematico Sergei Godunov è morto nel 2023):
q[i]=q[i]-C*(q[i]-q[i-1])
- metodo di Lax-Wendroff:
q[i]=0.5*C*(C+1)*q[i-1]+(1-C*C)*q[i]+0.5*C*(C-1)*q[i+1]
- metodo di Lax-Friedrichs:
q[i]=0.5*(q[i+1]+q[i-1])-0.5*C*(q[i+1]-q[i-1])
Come dati del problema, U=1, V0=1, Tout=10, per i metodi numerici N=200 punti e C=0.1 (numero di Courant, passo di discretizzazione).
Ecco tutto il codice scritto in Python:
import numpy as np
from matplotlib import pyplot as plt
N=200
V=np.zeros(N)
V0=1
U=1
fig,ax=plt.subplots(4,3) #3x3 grafici
def f(x,U,T,V0):
if i<(int)(0.25*N+U*T) or i>(int)(0.75*N+U*T):
return 0
else:
return V0
# INIT
time=0
for k in range(4):
for i in range(N):
ax[k][0].plot(i,f(i,U,time,V0),'o',color='blue')
# END
time=10
for i in range(N):
if i<(int)(0.25*N+U*time) or i>(int)(0.75*N+U*time):
V[i]=0
else:
V[i]=V0
ax[0][1].plot(i,V[i],'o',color='green')
# Godunov upwind
C=0.1
dx=1.0
dt=C*dx/U
Tout=10
# INIT
for i in range(N):
if i<(int)(0.25*N+U*Tout) or i>(int)(0.75*N+U*Tout):
V[i]=0
else:
V[i]=V0
for time in range((int) (Tout/(C*dx/U))):
for i in range(1,N):
V[i]+=-C*(V[i]-V[i-1])
for i in range(N):
ax[0][2].plot(i,0,'o',color='red')
ax[1][1].plot(i,V[i],'o',color='green')
if f(i,U,Tout,V0)!=0:
ax[1][2].plot(i,(1-V[i]/f(i,U,Tout,V0)),'o',color='red') # err rel
else:
ax[1][2].plot(i,0,'o',color='red') # err rel
# Lax-Wendroff
# INIT
for i in range(N):
if i<(int)(0.25*N+U*Tout) or i>(int)(0.75*N+U*Tout):
V[i]=0
else:
V[i]=V0
for time in range((int) (Tout/(C*dx/U))):
for i in range(1,N-1):
V[i]+=0.5*C*(C+1)*V[i-1]-C*C*V[i]+0.5*C*(C-1)*V[i+1]
for i in range(N):
ax[0][2].plot(i,0,'o',color='red')
ax[2][1].plot(i,V[i],'o',color='green')
if f(i,U,Tout,V0)!=0:
ax[2][2].plot(i,(1-V[i]/f(i,U,Tout,V0)),'o',color='red') # err rel
else:
ax[2][2].plot(i,0,'o',color='red') # err rel
# Lax-Friedrichs
# INIT
for i in range(N):
if i<(int)(0.25*N+U*Tout) or i>(int)(0.75*N+U*Tout):
V[i]=0
else:
V[i]=V0
for time in range((int) (Tout/(C*dx/U))):
for i in range(1,N-1):
V[i]=0.5*(V[i+1]+V[i-1])-0.5*C*(V[i+1]-V[i-1])
for i in range(N):
ax[3][1].plot(i,V[i],'o',color='green')
if f(i,U,Tout,V0)!=0:
ax[3][2].plot(i,(1-V[i]/f(i,U,Tout,V0)),'o',color='red') # err rel
else:
ax[3][2].plot(i,0,'o',color='red') # err rel
print("Numero iterazioni =",(int)(Tout/(C*dx/U)))
plt.suptitle("Onda quadra, LAE Equation: risoluzione analitica e numerica")
for i in range(4):
ax[i][0].set_title("valori iniziali (T=0)")
ax[0][1].set_title("soluzione analitica a T=Tout")
ax[0][2].set_title("err. rel. soluzione analitica")
ax[1][1].set_title("soluzione numerica Godunov a T=Tout")
ax[1][2].set_title("err. rel. Godunov")
ax[2][1].set_title("soluzione numerica Lax-Wendroff a T=Tout")
ax[2][2].set_title("err. rel. Lax-Wendroff")
ax[3][1].set_title("soluzione numerica Lax-Friedrichs a T=Tout")
ax[3][2].set_title("err. rel. Lax-Friedrichs")
plt.tight_layout() #ottimizza gli spazi
plt.show()
La cosa interessante, che possiamo commentare, la vediamo poi nell'immagine che segue.
- soluzione analitica: come è ovvio, semplice traslazione rigida, errore pari a zero
- metodo di Godunov: metodo numerico di primo ordine, non distorce il risultato ma non è molto accurato (tende a renderlo "gaussiano")
- metodo di Lax-Wendroff: metodo numerico di secondo ordine, si vede un migliore tentativo nell'approssimare l'onda quadra anche se in prossimità dei vertici potete osservare una "cosa strana", questa è la diffusione numerica! Ovvero un effetto indesiderato di uno schema numerico di ordine superiore (equazione di primo grado, semplice trasporto rigido, non dovremmo avere diffusione), questo caso meriterebbe una discussione a parte
- metodo di Lax-Friedrichs: metodo numerico di primo ordine, rispetto a Godunov risulta più approssimativo in generale, non rappresenta bene l'andamento, quanto all'errore massimo l'effetto è comunque più mediato
Ecco infine l'immagine con tutti i risultati.
Cosa ne pensate? Vi chiedo se avete già sentito parlare di questi metodi numerici oppure semplicemente, a vederli, quale vi piace di più? 🙂
