Differenze finite per le equazioni alle derivate parziali: calore, stabilità e CFL
Derivate come differenze, schema esplicito per l'equazione del calore, condizione di stabilità, Crank-Nicolson, Lax-Richtmyer, condizione di Courant e FDTD.
Le equazioni alle derivate parziali della lezione di Metodi matematici si risolvono a mano solo su domini semplici. Per un dissipatore di forma reale, per il campo in un dispositivo a semiconduttore, per la propagazione di un segnale su una linea con discontinuità, servono metodi numerici. Il più diretto è il metodo delle differenze finite: si sostituisce lo spazio e il tempo con una griglia di punti e le derivate con rapporti incrementali. In questa lezione lo applichiamo all’equazione del calore e vediamo i concetti di consistenza, stabilità e convergenza, e la condizione di Courant per le onde.
Derivate come differenze
Dallo sviluppo di Taylor, con passo 
:
![\[ u'(x) \approx \frac{u(x + \Delta x) - u(x)}{\Delta x}\quad\text{(errore } O(\Delta x)\text{)},\qquad u'(x) \approx \frac{u(x + \Delta x) - u(x - \Delta x)}{2\Delta x}\quad\text{(errore } O(\Delta x^2)\text{)} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-ed8cec81b5a5a16d99c33b415df07d30_l3.png "Rendered by QuickLaTeX.com")
La seconda, centrata, è più precisa a parità di passo, perché gli errori dei due lati si compensano. Per la derivata seconda, sommando gli sviluppi in 
e 
:
![\[ u''(x) \approx \frac{u(x + \Delta x) - 2u(x) + u(x - \Delta x)}{\Delta x^2} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-9247d18eadbf5cc91185039e9df633de_l3.png "Rendered by QuickLaTeX.com")
con errore 
. Queste formule bastano a discretizzare quasi tutte le equazioni della fisica.
L’equazione del calore: lo schema esplicito
su ![[0, L]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-1ffd428b4226e5d383c3fe072a154c75_l3.png "Rendered by QuickLaTeX.com")
. Si introduce la griglia 
, 
, e si indica con 
l’approssimazione di 
. Usando la differenza in avanti nel tempo e quella centrata nello spazio:
![\[ \frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha\,\frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-4d12720c02cb6a15265658ac80f05aa6_l3.png "Rendered by QuickLaTeX.com")
cioè, con 
,
![\[ u_i^{n+1} = u_i^n + r\big(u_{i+1}^n - 2u_i^n + u_{i-1}^n\big) \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-379b3cf40f202396e7cac23c526b7ac3_l3.png "Rendered by QuickLaTeX.com")
(rosso) si calcola dai tre valori vicini al passo precedente (blu): è lo “stencil” dello schema esplicito.Ogni nuovo valore si calcola dai tre vicini al passo precedente: si parte dalla condizione iniziale e si avanza riga per riga, applicando le condizioni al contorno agli estremi. Si può programmare in poche righe, anche in un foglio di calcolo.
La stabilità: 
Lo schema è stabile solo se
![\[ r = \frac{\alpha\,\Delta t}{\Delta x^2} \le \frac12 \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-a2d3b4f144db2a1d74455db98f7b08b7_l3.png "Rendered by QuickLaTeX.com")
Perché. Con 
il nuovo valore è una media pesata dei tre vicini con pesi positivi, 
: non può superare il massimo dei valori vecchi, quindi gli errori non crescono (principio del massimo discreto). Con 
il peso centrale diventa negativo, e una piccola oscillazione punto a punto viene amplificata a ogni passo: dopo poche decine di passi la soluzione esplode in uno zig-zag privo di senso.
La conseguenza pratica è pesante: dimezzando per avere più dettaglio, 
va diviso per quattro, e il lavoro totale cresce di otto volte. Come per le equazioni differenziali stiff, la soluzione sono gli schemi impliciti: lo schema di Crank-Nicolson, che fa la media tra il passo vecchio e quello nuovo, è stabile per qualsiasi e ha errore 
. A ogni passo richiede di risolvere un sistema lineare tridiagonale, che costa poco.
Consistenza, stabilità, convergenza
- Uno schema è consistente se, facendo tendere i passi a zero, l’equazione discreta tende a quella vera (l’errore di troncamento va a zero).
- È stabile se gli errori, di arrotondamento o di dati, non crescono senza limite durante i passi.
- È convergente se la soluzione numerica tende a quella vera quando i passi tendono a zero.
Il teorema di equivalenza di Lax-Richtmyer dice che per i problemi lineari ben posti un metodo consistente è convergente se e solo se è stabile. È il risultato centrale di tutta la teoria, ed è quello a cui si riferisce l’introduzione ai metodi alle differenze finite collegata a questa materia.
L’equazione delle onde e la condizione di Courant
Per 
, con differenze centrate nel tempo e nello spazio, lo schema esplicito è stabile se e solo se
![\[ C = \frac{c\,\Delta t}{\Delta x} \le 1 \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-970dc783d5b276fb1c25abbf4a07adbc_l3.png "Rendered by QuickLaTeX.com")
la condizione di Courant-Friedrichs-Lewy (CFL). Il significato è fisico: in un passo l’onda percorre 
, e lo schema usa solo i punti a distanza . Se l’onda va più lontano di un punto di griglia per passo, lo schema non “vede” da dove arriva l’informazione e la soluzione diventa instabile. Il dominio di dipendenza numerico deve contenere quello fisico.
È la condizione che regola i simulatori elettromagnetici FDTD (differenze finite nel dominio del tempo), usati per progettare antenne, guide d’onda e circuiti stampati ad alta velocità: con celle di 1 mm e la velocità della luce, il passo temporale deve stare sotto circa 2 picosecondi in tre dimensioni.
Due dimensioni e oltre
Per l’equazione di Laplace su una piastra, la formula delle differenze nelle due direzioni dà
![\[ u_{i,j} = \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1}}{4} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-9cade9b9e909d23e59c2f499d1fa5e02_l3.png "Rendered by QuickLaTeX.com")
ogni valore è la media dei quattro vicini, la versione discreta della proprietà della media delle funzioni armoniche. Ripetendo questa media su tutta la griglia finché i valori non cambiano più si ottiene il metodo di Jacobi della lezione sui sistemi lineari. Per geometrie complicate, con bordi curvi, si preferisce il metodo degli elementi finiti, che divide il dominio in triangoli o tetraedri e usa le formule di quadratura di Gauss: è il metodo dei programmi di simulazione strutturale, termica ed elettromagnetica usati nell’industria.
