Vediamo un caso pratico di applicazione di analisi numerica, il metodo di Crank-Nicolson (noto anche come metodo trapezoidale), per risolvere l'equazione differenziale della diffusione. Nello specifico, vediamo un'applicazione matematica di questo metodo numerico, che è semplicemente la media dei metodi FTCS (esplicito) e BTCS (implicito), che abbiamo già trattato. Ne risulta un metodo numerico semi-implicito, incondizionatamente stabile e con accuratezza di secondo ordine nel tempo e nello spazio.
Nota: Crank-Nicolson fa la media di FTCS e BTCS ad ogni step temporale, non le soluzioni finali (anche se il risultato può essere molto simile, formalmente è importante questa distinzione). Più precisamente ciò che avviene è la media dell'operatore spaziale (∂²U/∂x²) valutato ai tempi n (FTCS, esplicito) e n+1 (BTCS, implicito).
Per il caso di studio (Diffusion Equation 1D, schema di Crank-Nicolson), ci concentriamo sul metodo matematico, non il significato fisico. L'equazione della diffusione può applicarsi alla diffusione di temperatura, diffusione di una sostanza chimica (soluto non reattivo, per semplicità), la struttura matematica di base è la stessa, con l'equazione differenziale ∂U/∂t = D * ∂²U/∂x² (il significato delle variabili, i relativi valori e le condizioni al contorno, variano poi a seconda dell'applicazione specifica). Condizioni al contorno di Neumann, flusso nullo (∂U/∂x=0, ai bordi).
Per il codice Python, riportato in seguito, è stato fatto uso dell'IA (Qwen3.8-Max, Qwen AI), per una revisione e migliorie che rendono più efficiente il codice (ad esempio vettorizzazione anziché classici cicli for).
Importante: prima di riportare il codice, commentiamo un errore piuttosto comune nell'implementazione (di cui avevamo già discusso nello studio del metodo FTCS). Il seguente codice di esempio
while(time<=Tout):
for i in range(1,N-1):
U[i]=U[i]+Fo*(U[i+1]-2*U[i]+U[i-1])
sarebbe SBAGLIATO, anche se può apparire intuitivamente corretto come implementazione di FTCS. Questo perché U viene aggiornato in-place, quindi non è più il metodo FTCS ma uno schema misto. In pratica il metodo FTCS prevede di calcolare U[ i ] al tempo n+1 usando i valori U[ i ], U[ i+1 ] e U[ i-1 ] al tempo n; tuttavia U[ i-1 ] è già stato calcolato, questo codice aggiornerebbe U[ i ] con un valore U[ i-1 ] del tempo n+1 anziché n, quindi appunto sbagliato.
Vediamo dunque la versione completa e corretta del codice Python, (vettorizzazione e slicing di Numpy, per evitare l'errore precedentemente illustrato, in alternativa si può usare esplicitamente un comando come u_copy = u.copy() per una copia temporanea). A seguire lo screenshot del risultato grafico.
import numpy as np
from matplotlib import pyplot as plt
import time
# --- Parametri Fisici e di Griglia ---
L = 20
V0 = 100
r = 0.05 # Numero di diffusione (spesso chiamato 'r' o 'lambda' invece di 'delta')
dx = 0.1
D = 0.1
dt = r / D * dx**2
T = 20
N = int(L / dx)
steps = int(T / dt)
# --- Inizializzazione ---
x = np.linspace(0, L, N)
u = np.zeros(N)
u[N // 2] = V0 # Impulso iniziale al centro
# --- Costruzione Matrice (Vero Crank-Nicolson) ---
# Schema: (I - r/2 * L) u^{n+1} = (I + r/2 * L) u^n
# La matrice A_CN rappresenta il termine implicito (I - r/2 * L)
main_diag = np.full(N, 1 + r)
off_diag = np.full(N - 1, -r / 2)
A_CN = np.diag(main_diag) + np.diag(off_diag, k=1) + np.diag(off_diag, k=-1)
# Correzione Condizioni al Contorno: Flusso Nullo (Neumann)
# u_{-1} = u_1 -> il coefficiente raddoppia
A_CN[0, 1] = -r
A_CN[N-1, N-2] = -r
# Nota: Per massima efficienza in Python su matrici tridiagonali costanti,
# qui si potrebbe pre-fattorizzare A_CN (es. scipy.linalg.lu_factor)
# o usare l'Algoritmo di Thomas (O(N)) invece di np.linalg.solve (O(N^3)).
start_time = time.time()
# --- Loop Temporale ---
for k in range(steps):
# 1. Calcolo del Laplaciano (parte esplicita al tempo n)
laplacian = np.zeros(N)
laplacian[1:-1] = u[2:] - 2*u[1:-1] + u[:-2]
# Applicazione flusso nullo ai bordi per il termine esplicito
laplacian[0] = 2*u[1] - 2*u[0]
laplacian[N-1] = 2*u[N-2] - 2*u[N-1]
# 2. Costruzione del Termine Noto (RHS)
rhs = u + (r / 2) * laplacian
# 3. Risoluzione del sistema lineare
u = np.linalg.solve(A_CN, rhs)
end_time = time.time()
# --- Output e Grafico ---
print(f"N iterazioni = {steps}")
print(f"dt = {dt:.5f}")
print(f"Time = {end_time - start_time:.4f} s")
plt.figure(figsize=(8, 5))
plt.plot(x, u, 'o-', color='blue', markersize=3, label=f't = {T}')
plt.title("Equazione della Diffusione 1D\nMetodo di Crank-Nicolson")
plt.xlabel("Spazio (x)")
plt.ylabel("Concentrazione / Temperatura")
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()
Vediamo infine la rappresentazione grafica del risultato, tramite la libreria Matplotlib.
Un'ultima precisazione: come descritto nei commenti al codice, un ulteriore miglioramento di efficienza (tanto più rilevante quanto più è grande la dimensione della matrice, dunque il carico computazionale richiesto) si avrebbe dalla sostituzione del semplice metodo np.linalg.solve(A,b) che ha complessità O(N3), con l'Algoritmo di Thomas, trattandosi di un sistema tridiagonale (vedi BTCS), che ha complessità O(N).
