Dopo aver parlato dell'equazione della diffusione, risoluzione numerica tramite Crank-Nicolson, vediamo invece un'altra tipologia di diffusione.
Studiamo l'esempio di onda quadra, vale a dire, al tempo iniziale, ad esempio un segnale che presenta valore 1 in un certo range, 0 fuori da questo range. La soluzione analitica dell'equazione di trasporto è semplicemente del tipo x(t)=Ut
, una traslazione rigida senza fenomeni diffusivi ovvero fenomeni che smussano, attenuano l'onda (cosa che non dovrebbe avvenire), oltre ad altri "effetti collaterali" dei metodi numerici, come le oscillazioni.
Vediamo infatti il metodo di Lax-Wendroff alle differenze finite, uno schema numerico di secondo ordine che quindi pur avendo accuratezza migliore rispetto al primo ordine, presenta alcune controindicazioni. Il metodo numerico è il seguente: q[i]+=0.5*C*(C+1)*q[i-1]-C*C*q[i]+0.5*C*(C-1)*q[i+1]
, dove C è il numero di Courant, il passo di discretizzazione, approssimazione dello schema (più è piccolo e più il risultato è accurato, ma richiede maggiore costo computazionale).
In un secondo caso, dopo aver visto l'effetto indesiderato di Lax-Wendroff, oscillazioni in prossimità dei vertici dell'onda quadra, ho personalmente pensato di applicare questa correzione al metodo appena calcolato; in Python:
if V[i]>V[i-1] and V[i+1]<V[i] or V[i]<V[i-1] and V[i+1]>V[i]:
V[i]=0.5*(V[i-1]+V[i+1])
Quindi vale a dire che dev'esserci monotonia nella funzione e nel caso sia presente un'oscillazione, il valore centrale anziché "oscillare" viene impostato come pari alla media aritmetica dei due elementi adiacenti. Il concetto è molto più chiaro vedendo la successiva immegine, grafico realizzato tramite Matplotlib.
Ecco tutto il codice 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(3,2)
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(3):
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')
C=0.1
dx=1.0
dt=C*dx/U
Tout=10
# 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
#END
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[1][1].plot(i,V[i],'o',color='green')
# Lax-Wendroff aggiustato, correzione diffusione numerica
# 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
#END
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]
if V[i]>V[i-1] and V[i+1]<V[i] or V[i]<V[i-1] and V[i+1]>V[i]:
V[i]=0.5*(V[i-1]+V[i+1])
for i in range(N):
ax[2][1].plot(i,V[i],'o',color='green')
plt.suptitle("Onda quadra, LAE Equation: risoluzione analitica e numerica (Lax-Wendroff)")
for i in range(3):
ax[i][0].set_title("valori iniziali (T=0)")
ax[0][1].set_title("soluzione analitica a T=Tout")
ax[1][1].set_title("soluzione numerica Lax-Wendroff a T=Tout")
ax[2][1].set_title("soluzione numerica Lax-Wendroff corretto, a T=Tout")
plt.tight_layout() #ottimizza gli spazi
plt.show()
Infine l'immagine, che mostra prima la semplice soluzione analitica, traslazione rigida dell'onda quadra, poi il metodo di Lax-Wendroff tradizionale con l'effetto indesiderato delle oscillazioni, poi la mia correzione per risolvere questo problema.