Metodo trapezoidale (Crank-Nicolson) per equazione della diffusione - Python

Dopo aver visto l'algoritmo di Thomas per sistema tridiagonale, vediamo un altro caso di analisi numerica e possibile interesse pratico ingegneristico. Il metodo trapezoidale (anche noto come Crank-Nicolson) è la media aritmetica fra il metodo esplicito FTCS (Forward in Time, Central in Space) e il metodo implicito BTCS (Backward in Time, Central in Space). Il risultato di questo algoritmo è un'accuratezza di un ordine superiore rispetto ai precedenti!

Il caso in questione è stato risolto con Python (ovviamente a livello di efficienza, se fossero richieste più operazioni e quindi maggiore costo computazionale, ovviamente è più efficiente scriverlo in linguaggio C, vedi anche: Benchmark performance C, Python, MATLAB - GNU Octave, JavaScript, PHP).

Come condizioni iniziali per l'equazione della diffusione, abbiamo valore nullo ovunque tranne un impulso iniziale al centro; condizioni al contorno di flusso nullo alle estremità.

Nello specifico, i vettori A,x,b servono per il metodo BTCS (risoluzione del sistema Ax=b, tramite X=np.linalg.solve(A,b) e poi iterato per ogni step temporale), mentre FTCS è un vettore (del metodo omonimo, ovviamente) che si sovrascrive. Poi al termine (potrei farlo ad ogni iterazione ma è un'inutile spreco computazionale, non volendo memorizzare tutti i valori ad ogni tempo ma solo al tempo finale prefissato) sovrascrivo il vettore X, che diventerà quindi la soluzione con il metodo di Crank-Nicolson ovvero per ogni elemento del dominio: X[j]=0.5*(X[j]+FTCS[j]).

Ecco riportato tutto il codice Python:

# Diffusion Equation 1D
# Crank-Nicholson = 0.5*(FTCS+BTCS)
import numpy as np
from matplotlib import pyplot as plt
import time
L=20
V0=100
delta=0.05
dx=0.1
D=0.1
dt=delta/D*dx**2
T=20
N=(int)(L/dx)
A=np.zeros((N,N))
X=np.zeros(N)
b=np.zeros(N)
b[(int)(N/2)]=V0
X=b
FTCS=np.zeros(N)
FTCS[(int)(N/2)]=V0
# INIT
for j in range(N):
    for i in range(N):
        if i==j:
            A[i][j]=1+2*delta
        elif (i+1==j and i<N-1) or (i-1==j and i>0):
            A[i][j]=-delta
        else:
            A[i][j]=0
start=time.time()
# LOOP
for k in range((int)(T/dt)):
    # BTCS
    X=np.linalg.solve(A,b)
    b=X
    # FTCS
    for j in range(N):
            if j==0:
                FTCS[j]=FTCS[1]
            elif j==N-1:
                FTCS[j]=FTCS[N-2]
            else:
                FTCS[j]+=delta*(FTCS[j+1]-2*FTCS[j]+FTCS[j-1])
#C-N
for j in range(N):
    X[j]=0.5*(X[j]+FTCS[j])
end=time.time()
print("N iterazioni =",(int)(T/dt))
print("dt =",dt)
print("time =",end-start,"s")
plt.figure()
for i in range(N):
    plt.plot(i*dx,X[i],'o',color='blue')
plt.title("Diffusion Equation 1D\nCrank-Nicholson = 0.5*(FTCS+BTCS)")
plt.show()

Infine un'immagine che mostra l'andamento diffusivo, curva gaussiana.