Cosa causa questo problema numerico?

[GianlucaN]

Salve a tutti, vorrei un consiglio su come trattare un problema numerico originato probabilmente da un comportamento dello schema di derìvazione ai bordi del dominio. Ho scritto un codice MATLAB per la risoluzione numerica dell'equazione di advezione monodimensionale:

$ \frac{\partial u}{\partial t}+ c \frac{\partial u}{\partial x}=0 $

dove $c$ è una costante rapresentativa della velocità con cui avanza la condizione iniziale, mantenendo intatto il suo profilo. Gli schemi utilizzati per la derivate sono i seguenti (tratti dall'articolo Lele 1992: "Compact finite differences schemes with spectral-like resolution")


per i punti interni (VI ordine), e



per i punti al bordo del dominio (IV ordine). L'avanzamento temporale è condotto con uno schema Runge Kutta al IV ordine temporale. La condizione iniziale è una funzione gaussiana posizionata come in figura:




Quando evolvo la soluzione numerica, trovo che la gaussiana avanza esattamente con velocità $c$, ma dal bordo sinistro nasce e cresce indefinitamente una coda che dopo un po' distrugge la soluzione:



sarei interessato ad indagare la causa di questo fenomeno ed avere consigli su come arginarlo. Grazie.

[feddy]

Ciao, scrivo di fretta due cose nella speranza che ti possano essere utili in qualche modo.

1. Per un'equazione così semplice non basta il metodo delle caratteristiche? Suppongo però che non sia la strada che tu debba seguire.

2. Guardando di fretta quello che hai fatto tu mi sembra il metodo delle linee: discretizzazione spaziale con differenze finite e poi integri nel tempo, in questo caso con un RK4.
Se così fosse, c'è da prestare attenzione: una volta discretizzata la derivata rispetto a $x$ (ad esempio $\frac{\partial u}{\partial x_i}=\frac{u_{i+1}(t)-u_{i-1}(t)}{2h}$, anche se tu hai fatto discretizzato diversamente) il problema è del tipo $y'(t)= A y(t)$, con $A$ matrice tridiagonale che realizza la discretizzazione della derivata prima rispetto a $x$, ed potrebbe essere stiff, e un metodo esplicito come RK4 di sicuro non è consigliato. Infatti $A$ è tale per cui i suoi autovalori sono reali negativi, e amplificati dal passo $h$.

Sempre ammesso che le condizioni ai bordi siano corrette

[feddy]

Adesso che ci penso avevo in mente che $A$ fosse la matrice che discretizza la derivata seconda... quella per la derivata prima è sempre tridiagonale, ma magari non da problemi con un metodo esplicito...

guardando la tua ultima sequenza di figure mi sa che non hai messo la condizione a $x=0$...

[Raptorista]

Che condizione al bordo metti a sinistra?

[GianlucaN]

Buongiorno feddy e raptorista, vi ringrazio per le risposte. Dunque, per rispondervi:
@feddy: certo, potrei procedere con metodo delle caratteristiche, ma sto cercando di ritrovare quella soluzione numericamente ed analizzare i problemi che sorgono. Mi è stato utile lo spunto per riflettere sull'instabilità dello schema RK esplicito.
@raptorista: in effetti mi sono reso conto che non ho imposto alcuna condizione al bordo sinistro. Lo schema per la derivazione al bordo così come l'ho utilizzato non impone alcun vincolo sul valore della derivata o della funzione al bordo sinistro. E riflettendo sul metodo delle caratteristiche, la regione dello spaziotempo nella quale si crea quell'effetto indesiderato è quella tra l'asse $t$ e la retta $t=x/c$, per $t>0$, nella quale non vi è informazione propagabile dal dato iniziale. Quindi immagino che una condizione al bordo sul valore di $u$ o della sua derivata potrebbe risolvere fornendomi dei valori "corretti" in tale regione, anche se non ne sono del tutto convinto, che valori dovrebbe produrre? C'è qualcosa che mi sfugge?

[Raptorista]

Le condizioni al bordo sono un dato del problema continuo: esistono prima di, e indipendentemente dal, problema discretizzato.
Generalmente mi aspetterei che quel problema abbia una condizione di dirichlet al bordo sinistro, che è un bordo cosiddetto di inflow.