Implementazione delle condizioni al contorno per un'equazione non stazionaria
Salve a tutti,
Sono uno studente di ingegneria aerospaziale alle prese con l'esame di metodi numerici.
In questo esame si tratta la risoluzione di PDE con il metodo delle differenze finite, sia stazionarie che non stazionarie.
Purtroppo la materia, che trovo interessante, è piuttosto complessa per me in quanto manca l'adeguato materiale didattico; esistono solo dei codici compilati da qualche collega, senza un po' di commento.
La teoria l'ho studiata qua e la da qualche libro, e pertanto trovo difficoltà, non tanto nel campo dello stazionario, quanto in quello della non stazionarietà ed in particolare nell'implementare le opportune condizioni al contorno.
Anche perché molti testi sono molto specifici, mentre il mio esame consta di soli 6 CFU.
Fatto questo preambolo, vengo al problema.
Ho, ad esempio, questa equazione di diffusione non stazionaria
$ (d\varphi)/dt=k(d^2\varphi)/dx^2 $
Il testo chiede di risolverla con un metodo FTCS, ovvero Eulero esplicito.
A questo punto so che la derivata temporale la si svolge tra due intervalli di tempo mentre quella spaziale viene effettuata per i soli punti interni del mesh.
In particolare, se suppongo un mesh eqispaziato di passo h avrò la seguente discretizzazione
$ \varphi_i^(n+1)=\beta\varphi_(i-1)^n+(1-2\beta)\varphi_i^n+\beta\varphi_(i+1)^n $
dove $ \beta=(k\Deltat)/(h^2) $.
Cioè in questo caso, essendo il mesh spaziale uniforme, la matrice che discretizza la derivata seconda nei punti interni del mesh è tridiagonale ed ha sulle diagonali 1, -2, 1.
A questo punto, però, se ho una condizione alla Neumann non so come procedere.
Il testo mi consiglia una discretizzazione della derivata prima sui due nodi estremanti asimmetrica e del primo ordine; inoltre mi suggerisce di applicarla solo ai nodi all'intervallo di tempo n+1.
Qualcuno sa darmi qualche suggerimento sul da farsi?
Grazie a tutti.
Sono uno studente di ingegneria aerospaziale alle prese con l'esame di metodi numerici.
In questo esame si tratta la risoluzione di PDE con il metodo delle differenze finite, sia stazionarie che non stazionarie.
Purtroppo la materia, che trovo interessante, è piuttosto complessa per me in quanto manca l'adeguato materiale didattico; esistono solo dei codici compilati da qualche collega, senza un po' di commento.
La teoria l'ho studiata qua e la da qualche libro, e pertanto trovo difficoltà, non tanto nel campo dello stazionario, quanto in quello della non stazionarietà ed in particolare nell'implementare le opportune condizioni al contorno.
Anche perché molti testi sono molto specifici, mentre il mio esame consta di soli 6 CFU.
Fatto questo preambolo, vengo al problema.
Ho, ad esempio, questa equazione di diffusione non stazionaria
$ (d\varphi)/dt=k(d^2\varphi)/dx^2 $
Il testo chiede di risolverla con un metodo FTCS, ovvero Eulero esplicito.
A questo punto so che la derivata temporale la si svolge tra due intervalli di tempo mentre quella spaziale viene effettuata per i soli punti interni del mesh.
In particolare, se suppongo un mesh eqispaziato di passo h avrò la seguente discretizzazione
$ \varphi_i^(n+1)=\beta\varphi_(i-1)^n+(1-2\beta)\varphi_i^n+\beta\varphi_(i+1)^n $
dove $ \beta=(k\Deltat)/(h^2) $.
Cioè in questo caso, essendo il mesh spaziale uniforme, la matrice che discretizza la derivata seconda nei punti interni del mesh è tridiagonale ed ha sulle diagonali 1, -2, 1.
A questo punto, però, se ho una condizione alla Neumann non so come procedere.
Il testo mi consiglia una discretizzazione della derivata prima sui due nodi estremanti asimmetrica e del primo ordine; inoltre mi suggerisce di applicarla solo ai nodi all'intervallo di tempo n+1.
Qualcuno sa darmi qualche suggerimento sul da farsi?
Grazie a tutti.
Risposte
Non ho molto tempo, comunque si chiama metodo delle linee: prima discretizzi gli operatori differenziali spaziali con differenze finite (in questo caso centrate del secondo ordine), così da ridurti ad una ODE, poi integri nel tempo, stando attento al metodo numerico utilizzato, in quanto il problema è in generale "stiff", a causa della matrice che hai menzionato sopra.
Per semplicità, chiamo $\varphi=u$, e supponi di avere una condizione tale PDE con condizione al bordo $x=b$ data da $\frac{\partial u}{\partial x}(t,b)=0$.
Introduci la variabile $y_{m+1}(t) \approx u(t,x_{m+1})$, dove $x_{m+1}=b+h$, e imponi $y_{m+1}(t)=y_{m-1}(t)$.
Quindi, quando vai ad approssimare la derivata seconda nel tuo problema, avrai al bordo:
$\frac{\partial^2 u}{\partial x^2}(t,b) \approx \frac{u(t,x_{m+1})-2*u(t,x_m)+u(t,x_{m+1} )} {h^2}=\frac{y_{m+1}(t)-2y_m(t)+y_{m-1}(t)}{h^2}=\frac{2*y_{m-1}(t)-2y_m(t)}{h^2}$
Quindi nella matrice tridiagonale dovrai cambiare le entrate $(m,m-1)$ e $(m,m)$ con rispettivamente $2$ e $-2$, ricordandoti di dividere per $h^2$!
Per semplicità, chiamo $\varphi=u$, e supponi di avere una condizione tale PDE con condizione al bordo $x=b$ data da $\frac{\partial u}{\partial x}(t,b)=0$.
Introduci la variabile $y_{m+1}(t) \approx u(t,x_{m+1})$, dove $x_{m+1}=b+h$, e imponi $y_{m+1}(t)=y_{m-1}(t)$.
Quindi, quando vai ad approssimare la derivata seconda nel tuo problema, avrai al bordo:
$\frac{\partial^2 u}{\partial x^2}(t,b) \approx \frac{u(t,x_{m+1})-2*u(t,x_m)+u(t,x_{m+1} )} {h^2}=\frac{y_{m+1}(t)-2y_m(t)+y_{m-1}(t)}{h^2}=\frac{2*y_{m-1}(t)-2y_m(t)}{h^2}$
Quindi nella matrice tridiagonale dovrai cambiare le entrate $(m,m-1)$ e $(m,m)$ con rispettivamente $2$ e $-2$, ricordandoti di dividere per $h^2$!
Ciao, ti ringrazio per la celere risposta.
Il problema che mi ritrovo ad affrontare è il modo con cui implementare questa condizione in MatLab.
Cioè, per i nodi interni non c'è problema... il problema mi si pone quando vado a discretizzare per i nodi di bordo.
Ad esempio, se ho N nodi che discretizzano il mesh, devo allocare una matrice NXN?
Il problema che mi ritrovo ad affrontare è il modo con cui implementare questa condizione in MatLab.
Cioè, per i nodi interni non c'è problema... il problema mi si pone quando vado a discretizzare per i nodi di bordo.
Ad esempio, se ho N nodi che discretizzano il mesh, devo allocare una matrice NXN?
Rileggi l'ultima frase nella mia risposta... è quello che c'è da fare nella pratica per nua condizione di Neumaann omogenea. Insomma, come per quelle di Dirichlet, vanno cambiate la prima e l'ultima riga della matrice che discretizza il problema.
Ovvio che devi allocare una matrice NxN, e tale matrice sarà chiaramente sparsa, con, a priori, $N+2*(N-1)=3N-2$ elementi diversi da 0
Ovvio che devi allocare una matrice NxN, e tale matrice sarà chiaramente sparsa, con, a priori, $N+2*(N-1)=3N-2$ elementi diversi da 0
Ti ringrazio per la risposta, ho risolto il quesito!
Nel caso avessi una derivata prima (quindi dovessi discretizzare Neumann) da discretizzare non su due punti ma, ad esempio, su quattro?
Ad esempio, se mi metto in x = b, posso effettuare una discretizzazione su quattro punti asimmetrica?
Il procedimento è lo stesso, solo che devo cambiare gli ultimi quattro posti della matrice, giusto?
Grazie!
Nel caso avessi una derivata prima (quindi dovessi discretizzare Neumann) da discretizzare non su due punti ma, ad esempio, su quattro?
Ad esempio, se mi metto in x = b, posso effettuare una discretizzazione su quattro punti asimmetrica?
Il procedimento è lo stesso, solo che devo cambiare gli ultimi quattro posti della matrice, giusto?
Grazie!
Felice di esserti stato utile! 
Sostanzialmente vuoi usare differenze finite non centrate. Sì, sostanzialmente userai i punti a sinistra. Di solito ci sono già le formule, ad esempio, per discretizzare la derivata prima al bordo di destra usando due punti a sinistra, o cose così.
Sostanzialmente vuoi usare differenze finite non centrate. Sì, sostanzialmente userai i punti a sinistra. Di solito ci sono già le formule, ad esempio, per discretizzare la derivata prima al bordo di destra usando due punti a sinistra, o cose così.
Ti ringrazio ancora per l'aiuto ed il prezioso supporto...
Il fatto è che all'esame il professore di diverte coi più disparati metodi di interpolazione; stencil asimmetrici, con punti iniziali particolari e condizioni al contorno periodiche la cui implementazione faccio fatica a comprendere...
Non riesco a trovare nemmeno un testo che tratti, in maniera chiara, le equazioni del calore (diffusione e convezione).
Il fatto è che all'esame il professore di diverte coi più disparati metodi di interpolazione; stencil asimmetrici, con punti iniziali particolari e condizioni al contorno periodiche la cui implementazione faccio fatica a comprendere...
Non riesco a trovare nemmeno un testo che tratti, in maniera chiara, le equazioni del calore (diffusione e convezione).
Ciao, perdona ancora il disturbo..
Ho provato a svolgere questa traccia, ma i risultati sono deludenti...
Questo è il mio codice...
Non riesco a capire dove sbaglio, ho seguito le indicazioni di qualche collega ma il mio risultato oscilla...fa schifo insomma
Ti ringrazio
Ho provato a svolgere questa traccia, ma i risultati sono deludenti...
Questo è il mio codice...
Non riesco a capire dove sbaglio, ho seguito le indicazioni di qualche collega ma il mio risultato oscilla...fa schifo insomma
Ti ringrazio
Ho fatto una "versione" senza quella strana forzante e con condizione iniziale di tipo sinusoidale...e pare venga.
Ho, inoltre, rivisto alcune cosette...
Se puoi buttarci un occhio te ne sarei infinitamente grato!
Ho, inoltre, rivisto alcune cosette...
Se puoi buttarci un occhio te ne sarei infinitamente grato!
Onestamente non ho guardato bene il codice, anche perché ora non ho molto tempo. L'imposizione della condizione di Neumann con FD decentrate va bene: hai cambiato la prima riga della matrice di discretizzazione.
Ti ringrazio, l'importante è sapere di aver ben appreso come si implementa la matrice!
Un commento: per metodi espliciti in tempo allocare una matrice è uno spreco di tempo e memoria, faresti meglio a fare un ciclo su tutti gli elementi del vettore e applicare direttamente la formula in questione, cioè \(u_i^{n+1} = u_i^n + dt\cdot(\dots\).
I nodi di bordo vanno trattati a parte, chiaramente, ma non è difficile.
I nodi di bordo vanno trattati a parte, chiaramente, ma non è difficile.
Da quel che ho letto ci sarebbero una barcata di commenti da fare.
Cominciamo con qualche piccolo accorgimento sul problema: è chiaro che se discretizzi con passo h in N intervalli il tuo dominio, avrai, una matrice di N+1xN+1 nodi, dove "la cornice esterna" di 2 righe e 2 colonne saranno i nodi di boundary, cioè dove sono definite le equazioni di condizioni al contorno; mentre nella matrice interna N-1xN-1 le equazioni alle differenze finite che hai ottenuto a partire dalla PDE.
Per quanto riguarda "dove saperne di più". Mi sembra di capire che il tuo sia il primo corso di calcolo scientifico. Mai capito il senso di proporre questi temi senza partire dalle basi. Mi viene da segnalarti due testi: "Iterative Methods for Sparse Linear Systems" di Y. Saad, è opne source e lo trovi facilmente digitando il titolo dle libro su Google. Il capitolo 2 ti spiega passo dopo passo come costruire diverse stencils, le differenze finite, come costruirti la matrice dei coefficienti del tuo dominio discretizzato e un sacco di altre cose. E' un po' avanzato in certi passaggi, ma è senz'altro una fonte di rilievo. "Calcolo Scientifico" di Quarteroni, Saleri e Gervasio è un libro di introduzione al Calcolo Scientifico; hai due capitoli dedicati alle Equazioni Differenziali (uno alle ODE e uno alle PDE), esempi ed esercizi, codici in Matlab riportati in modo completo.
Cominciamo con qualche piccolo accorgimento sul problema: è chiaro che se discretizzi con passo h in N intervalli il tuo dominio, avrai, una matrice di N+1xN+1 nodi, dove "la cornice esterna" di 2 righe e 2 colonne saranno i nodi di boundary, cioè dove sono definite le equazioni di condizioni al contorno; mentre nella matrice interna N-1xN-1 le equazioni alle differenze finite che hai ottenuto a partire dalla PDE.
Per quanto riguarda "dove saperne di più". Mi sembra di capire che il tuo sia il primo corso di calcolo scientifico. Mai capito il senso di proporre questi temi senza partire dalle basi. Mi viene da segnalarti due testi: "Iterative Methods for Sparse Linear Systems" di Y. Saad, è opne source e lo trovi facilmente digitando il titolo dle libro su Google. Il capitolo 2 ti spiega passo dopo passo come costruire diverse stencils, le differenze finite, come costruirti la matrice dei coefficienti del tuo dominio discretizzato e un sacco di altre cose. E' un po' avanzato in certi passaggi, ma è senz'altro una fonte di rilievo. "Calcolo Scientifico" di Quarteroni, Saleri e Gervasio è un libro di introduzione al Calcolo Scientifico; hai due capitoli dedicati alle Equazioni Differenziali (uno alle ODE e uno alle PDE), esempi ed esercizi, codici in Matlab riportati in modo completo.
Come motivazione a sostegno dell'ottimo libro citato di @anonymous_40e072 (altrimenti l'OP potrebbe non capire il senso dell'utilizzo di metodi iterativi per sistemi lineari sparsi), pensa a questo esempio:
Hai l'equazione di Poisson su un griglia bidimensionale (con condizioni omogenee) con 40 nodi su ogni lato, e la matrice di stiffness ( o rigidezza) del problema (di taglia 1600) ha circa 7840 elementi non nulli. Vista la sua simmetria, uno può pensare ad una fattorizzazione di Cholesky di $A$ in $A=R^{T}R$ per la risoluzione del sistema lineare corrispondente. Ma questo produce la matrice $R$ con 68135 elementi diversi da zero, e questo non si è approfittato della "sparsità" della matrice $A$. Da qui la necessità dell'utilizzo di metodi iterativi
Hai l'equazione di Poisson su un griglia bidimensionale (con condizioni omogenee) con 40 nodi su ogni lato, e la matrice di stiffness ( o rigidezza) del problema (di taglia 1600) ha circa 7840 elementi non nulli. Vista la sua simmetria, uno può pensare ad una fattorizzazione di Cholesky di $A$ in $A=R^{T}R$ per la risoluzione del sistema lineare corrispondente. Ma questo produce la matrice $R$ con 68135 elementi diversi da zero, e questo non si è approfittato della "sparsità" della matrice $A$. Da qui la necessità dell'utilizzo di metodi iterativi
"feddy":
Come motivazione a sostegno dell'ottimo libro citato di @anonymous_40e072 (altrimenti l'OP potrebbe non capire il senso dell'utilizzo di metodi iterativi per sistemi lineari sparsi), pensa a questo esempio:
Hai l'equazione di Poisson su un griglia bidimensionale (con condizioni omogenee) con 40 nodi su ogni lato, e la matrice di stiffness ( o rigidezza) del problema (di taglia 1600) ha circa 7840 elementi non nulli. Vista la sua simmetria, uno può pensare ad una fattorizzazione di Cholesky di $A$ in $A=R^{T}R$ per la risoluzione del sistema lineare corrispondente. Ma questo produce la matrice $R$ con 68135 elementi diversi da zero, e questo non si è approfittato della "sparsità" della matrice $A$. Da qui la necessità dell'utilizzo di metodi iterativi
Sì, confermo. Tendenzialmente i motivi a supporto dei metodi iterativi sono due: quello che hai spiegato tu, ossia quello che in gergo è chiamato fill-in, e la scalabilità dei metodi iterativi. In altri termini, la complessità computazionale dei metodi diretti (i.e. basati su fattorizzazioni varie, come di Cholesky o LU) aumenta più velocemente rispetto ai metodi iterativi. Quando i problemi si fanno grossi e tosti (e quasi sempre, anche sparsi) i metodi iterativi diventano l'unica scelta possibile.
Tutto senz'altro molto interessante, ma probabilmente un "po' troppo" per l'OP
Wow, vi ringrazio per i preziosi commenti e suggerimenti!
Purtroppo i metodi iterativi non sono presenti nel programma stilato dal prof, in effetti ho già letto di queste metodologie proprio sul Quaternoni che ho usato come "base", ma che si discosta molto dall'approccio del mio corso.
E' una tematica molto interessante quella del calcolo scientifico, sicuramente la approfondirò in un corso di laurea magistrale!
Grazie ancora a tutti
Purtroppo i metodi iterativi non sono presenti nel programma stilato dal prof, in effetti ho già letto di queste metodologie proprio sul Quaternoni che ho usato come "base", ma che si discosta molto dall'approccio del mio corso.
E' una tematica molto interessante quella del calcolo scientifico, sicuramente la approfondirò in un corso di laurea magistrale!
Grazie ancora a tutti
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo