Dopo aver visto il caso Onda quadra, trasporto - metodi numerici a confronto - Python, sempre per l'equazione LAE (Linear Advection Equation), quindi un semplice trasporto, traslazione rigida (secondo x(t)=Ut), vediamo una serie di casi particolari di metodi numerici alle differenze finite, partendo da una condizione iniziale a gradino. Poi un programma scritto in Python che mostra il risultato grafico nei vari casi.
- soluzione analitica: appunto la semplice traslazione rigida,
x(t)=Ut
- metodo di Lax-Wendroff tradizionale, schema esplicito (con C=0.5): qn+1=0.5C(C+1)q[i-1]n+(1-C2)qn+0.5C(C-1)q[i+1]n
- metodo di Lax-Wendroff con instabilità dovuta al numero di Courant (C>1)
- metodo di Lax-Wendroff implicito, schema instabile: nota, alcuni metodi numerici possono essere usati alternativamente in forma esplicita e implicita (FTCS, BTCS e Crank-Nicolson la combinazione ovvero media dei due), altri metodi numerici invece possono dare origine ad instabilità, senza approfondire troppo volendo lo si può verificare con l'analisi di von Neumann; rispetto al caso precedente, schema di Lax-Wendroff esplicito, stabile (nel caso 0<C<1), la risoluzione di qn+1=0.5C(C+1)q[i-1]n+1+(1-C2)qn+1+0.5C(C-1)q[i+1]n+1 origina un metodo numerico instabile
Come dati del problema, U=1, V0=1, Tout=10, C=0.5, a=-10, b=15, N=100 punti.
Ecco tutto il codice in linguaggio Python.
from matplotlib import pyplot as plt
from numpy import zeros,linalg
N=100
V=zeros(N)
V0=1
U=1
a=-10 # Xinf
b=15 # Xsup
h=1.0*(b-a)/N # passo discretizzazione spaziale
Tout=10 # tempo finale
C=0.5 # numero di Courant
dx=1.0
dt=1.0*C*dx/U
fig,ax=plt.subplots(4,2)
def f(x): # INIT
if x<0.5*b:
return 1
else:
return 0
for i in range(N):
for j in range(4):
ax[j][0].plot(a+i*h,f(a+i*h),"o",color="blue") # T0
ax[j][0].set_title("Situazione iniziale (T=0)")
ax[0][1].plot(a+i*h,f(a+i*h),"o",color="green") # sol analitica Tout
ax[0][1].set_title("Soluzione analitica a Tout")
ax[1][1].set_title("Lax Wendroff standard, a Tout")
ax[2][1].set_title("Lax Wendroff instabile con C>1, a Tout")
ax[3][1].set_title("Lax Wendroff instabile implicito, a Tout")
# Lax-Wendroff esplicito standard
# INIT
for i in range(N):
V[i]=f(a+i*h)
# END
for time in range((int)(Tout/dt)):
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(a+i*h,V[i],'o',color='green')
# Lax-Wendroff instabile per numero di Courant
# INIT
for i in range(N):
V[i]=f(a+i*h)
# END
C=1.5 # Cmax=1, quindi instabile
for time in range((int)(Tout/dt)):
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[2][1].plot(a+i*h,V[i],'o',color='green')
# Lax-Wendroff implicito (instabile)
# INIT
for i in range(N):
V[i]=f(a+i*h)
# END
A=zeros((N,N))
for j in range(N):
for i in range(N):
if(i==j):
A[i][j]=C*C
elif i-1==j:
A[i][j]=-0.5*C*(C-1)
elif i+1==j:
A[i][j]=-0.5*C*(C+1)
else:
A[i][j]=0
for time in range((int) (Tout/(C*dx/U))):
V=linalg.solve(A,V) # Ax=b
for i in range(N):
ax[3][1].plot(a+i*h,V[i],'o',color='green')
plt.suptitle("Equazione LAE con condizione iniziale a gradino")
plt.tight_layout() # ottimizza gli spazi
plt.show()
Vediamo infine un'immagine che mostra il risultato grafico. Si nota chiaramente che gli schemi instabili portano ad una soluzione completamente sballata e l'andamento è diverso che si tratti di numero di Courant - passo di discretizzazione - troppo elevato (C>1) oppure uno schema intrinsecamente instabile (la versione implicita del metodo di Lax Wendroff che, a differenza di altri metodi numerici, in questo caso non è un metodo numerico stabile).
Cosa ne pensate? Avevate già un'idea dell'argomento? 🙂
