Neumann BC per l'equazione di Poisson (1D)
Ciao a tutti,
come descritto in questo post, nel caso dell'equazione di Poisson $u''(x)=-f(x)$ (prendo il caso semplice 1D), il problema con condizioni di Neumann poiché la soluzione è definita a meno di una costante.
Come esempio concreto prendo $f(x)=\cos(2 \pi x)$ e $u'(0)=-1$ e $u'(1)=0$.
Supponendo di utilizzare uno schema alle differenze finite centrate classico di ordine 2, i.e. $u''(x_i)=\frac{u(x_{i+1})-2u(x_i) + u(x_{i-1})}{h^2} + \mathcal{O}(h^2)$ la matrice risultante è singolare, come lecito aspettarsi.
Tuttavia, ho letto diverse risposte su scicomp che trattano l'argomento ma nessuna mi ha convinto sinceramente. Lasciando stare l'utilizzo di metodi iterativi o del GC, come posso rendere il problema associato ben posto?
come descritto in questo post, nel caso dell'equazione di Poisson $u''(x)=-f(x)$ (prendo il caso semplice 1D), il problema con condizioni di Neumann poiché la soluzione è definita a meno di una costante.
Come esempio concreto prendo $f(x)=\cos(2 \pi x)$ e $u'(0)=-1$ e $u'(1)=0$.
Supponendo di utilizzare uno schema alle differenze finite centrate classico di ordine 2, i.e. $u''(x_i)=\frac{u(x_{i+1})-2u(x_i) + u(x_{i-1})}{h^2} + \mathcal{O}(h^2)$ la matrice risultante è singolare, come lecito aspettarsi.
Tuttavia, ho letto diverse risposte su scicomp che trattano l'argomento ma nessuna mi ha convinto sinceramente. Lasciando stare l'utilizzo di metodi iterativi o del GC, come posso rendere il problema associato ben posto?
Risposte
Un modo per rendere il problema ben posto è quello di richiedere che valga
Detto $Ax=b$ il sistema con $A$ singolare a cui si giunge imponendo le condizioni di Neumann, si può aggiungere il vincolo qui sopra considerando il sistema lineare di matrice (non più singolare)
\( \begin{bmatrix} A & \vec{1} \\ \vec{1}^T & 0 \end{bmatrix} \begin{bmatrix} \vec{x} \\ \mu \end{bmatrix} = \begin{bmatrix} \vec{b} \\ 0 \end{bmatrix} \)
Il seguente codice Python riproduce quanto ho scritto sopra:
e produce il grafico seguente, dove visivamente le condizioni sulle derivate appaiono rispettate. Per validare il codice dovrei mostrare che questa è un'approssimazione del secondo ordine usando una soluzione di riferimento, tuttavia per ora mi basta sapere se fin qui tutto torna o se ci sono altre strade valide.
$ \int_{0}^{1} u(x)dx=0$
Detto $Ax=b$ il sistema con $A$ singolare a cui si giunge imponendo le condizioni di Neumann, si può aggiungere il vincolo qui sopra considerando il sistema lineare di matrice (non più singolare)
\( \begin{bmatrix} A & \vec{1} \\ \vec{1}^T & 0 \end{bmatrix} \begin{bmatrix} \vec{x} \\ \mu \end{bmatrix} = \begin{bmatrix} \vec{b} \\ 0 \end{bmatrix} \)
Il seguente codice Python riproduce quanto ho scritto sopra:
import numpy as np
import matplotlib.pyplot as plt
from math import pi
def rhs(x):
return np.cos(2*pi*x)
def setBoundary(A,f):
M = len(f) - 1
A[0, 0] = -2/h**2
A[0, 1] = 2/h**2
f[0] = -np.cos(2*pi*x[0]) - 2/h
A[M, M - 1] = 2/h**2
A[M, M] = -2/h**2
f[M] = -np.cos(2*pi*x[M])
return A,f
a = 0
b = 1
M = 40
h=1/M
x = np.linspace(a,b,M+1)
A = np.diag(np.ones(M),-1) -2*np.diag(np.ones(M+1)) + np.diag(np.ones(M),+1)
A = A/(h**2)
f = -rhs(x) #! u''(x) = -cos(2pix)
A,f = setBoundary(A,f)
G=np.zeros([M+2,M+2])
G[0:M+1,0:M+1]=A
G[-1,:]=np.ones(M+2)
G[:,-1]=np.ones(M+2)
G[-1,-1]=0
newf=np.append(f,0)
uRif=np.linalg.solve(G,newf)
uRif = uRif[0:-1] #exclude last because of Lagrange multipliers
plt.figure()
plt.plot(x,uRif,'-o')
plt.xlabel('x')
plt.show()
e produce il grafico seguente, dove visivamente le condizioni sulle derivate appaiono rispettate. Per validare il codice dovrei mostrare che questa è un'approssimazione del secondo ordine usando una soluzione di riferimento, tuttavia per ora mi basta sapere se fin qui tutto torna o se ci sono altre strade valide.
Un altro metodo possibile è quello di imporre la temperatura in un punto a tua scelta, anche interno al dominio. In generale, qualunque metodo che identifichi una, e una sola, delle soluzioni che differiscono di una costante renderà il problema ben posto.
Grazie per la risposta. Tale valore della soluzione va imposto solamente nell'interno del dominio da quello che ho letto in giro, non al bordo. In quanto se imponessi una condizione di Dirichlet a $x=1$ starei effettivamente risolvendo un altro problema.
Conosci per caso qualche testo dove viene trattato questo problema? Domande a parte, non ho trovato molto
Conosci per caso qualche testo dove viene trattato questo problema? Domande a parte, non ho trovato molto
Non vedo il problema nell'imporre il valore desiderato sul bordo :S
Non credo che ci sia molto da trattare su questo problema. Se hai un problema mal posto, buona parte della teoria cade; se questo argomento fosse in un libro, sarebbe in un paragrafo "trucchi da tenere a mente", perché alla fine si tratta di un trucco, nulla più.
Non credo che ci sia molto da trattare su questo problema. Se hai un problema mal posto, buona parte della teoria cade; se questo argomento fosse in un libro, sarebbe in un paragrafo "trucchi da tenere a mente", perché alla fine si tratta di un trucco, nulla più.
[strike]Imponendo $u(1)=1$ ad esempio, non avrei più il comportamente al bordo desiderato. Cioè sto proprio risolvendo un altro problema. Dove sto sbagliando?[/strike]
[strike]Per implementare questo mi basta cambiare l'ultima riga della matrice $A$? Per intenderci, quella che era singolare con prima riga $[-2,2,0,...]$ e ultima riga $[0,0,....,2,-2]$
L'ultima riga del sistema la ricavo così dunque. Dalla condizione di Neumann omogenea (usando un nodo fantasma) segue $u_{N+1}=u_{N-1}$ e dunque l'ultima riga è $\frac{2u_{n-1} - 2 u_n}{h^2}$. A questo punto, potrei fissare $u_n=a$. Però ho notato che facendo così non ottengo ciò che speravo[/strike]
L'ultima riga del sistema la ricavo così dunque. Dalla condizione di Neumann omogenea (usando un nodo fantasma) segue $u_{N+1}=u_{N-1}$ e dunque l'ultima riga è $\frac{2u_{n-1} - 2 u_n}{h^2}$. A questo punto, potrei fissare $u_n=a$. Però ho notato che facendo così non ottengo ciò che speravo[/strike]
Anzi no, sono riuscito
Ho scritto una risposta di getto ma mi sono accorto che non è così semplice come sembra a prima vista, quindi l'ho cancellata 
Questo il ragionamento che ho fatto:
il problema
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
\end{cases}
\]
ha infinite soluzioni, che differiscono di una costante.
Il problema, ad esempio,
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
u(0.5) = c \\
\end{cases}
\]
ha una soluzione unica, sia \(u_c\).
La stessa funzione \(u_c\) è soluzione del problema
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
u(\hat x) = u_c(\hat x) \\
\end{cases}
\]
per qualunque \(\hat x\), ma se prendo \(\hat x = 1\) ottengo il problema
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
u(1) = u_c(1) \\
\end{cases}
\]
e here be dragons, perché direi che questo problema ha la stessa soluzione \(u_c\), però sappiamo anche che il problema misto
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u(1) = u_c(1) \\
\end{cases}
\]
ha una soluzione unica tutta sua, e a priori non vedo garanzia che soddisfi la condizione che ho appena tolto sulla derivata al bordo destro. C'è qualcosa che sfugge al ragionamento.
Un errore che vedo è che non abbiamo controllato la condizione di compatibilità per il problema di Neumann. Se guardi bene, il grafico che hai messo nel secondo post ha derivata in \(x = 0\) pari a \(-1\), non \(1\) come scritto all'inizio.
Quindi non è vero che il problema
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
\end{cases}
\]
ha infinite soluzioni per ogni \(a\) e \(b\), ma solo per valori che soddisfano la condizione
\[
\int_{\partial\Omega} \partial_n u = \int_\Omega f.
\]
Bisogna fare più attenzione a questo dettaglio.
Questo il ragionamento che ho fatto:
il problema
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
\end{cases}
\]
ha infinite soluzioni, che differiscono di una costante.
Il problema, ad esempio,
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
u(0.5) = c \\
\end{cases}
\]
ha una soluzione unica, sia \(u_c\).
La stessa funzione \(u_c\) è soluzione del problema
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
u(\hat x) = u_c(\hat x) \\
\end{cases}
\]
per qualunque \(\hat x\), ma se prendo \(\hat x = 1\) ottengo il problema
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
u(1) = u_c(1) \\
\end{cases}
\]
e here be dragons, perché direi che questo problema ha la stessa soluzione \(u_c\), però sappiamo anche che il problema misto
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u(1) = u_c(1) \\
\end{cases}
\]
ha una soluzione unica tutta sua, e a priori non vedo garanzia che soddisfi la condizione che ho appena tolto sulla derivata al bordo destro. C'è qualcosa che sfugge al ragionamento.
Un errore che vedo è che non abbiamo controllato la condizione di compatibilità per il problema di Neumann. Se guardi bene, il grafico che hai messo nel secondo post ha derivata in \(x = 0\) pari a \(-1\), non \(1\) come scritto all'inizio.
Quindi non è vero che il problema
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b \\
\end{cases}
\]
ha infinite soluzioni per ogni \(a\) e \(b\), ma solo per valori che soddisfano la condizione
\[
\int_{\partial\Omega} \partial_n u = \int_\Omega f.
\]
Bisogna fare più attenzione a questo dettaglio.
"Raptorista":
Un errore che vedo è che non abbiamo controllato la condizione di compatibilità per il problema di Neumann. Se guardi bene, il grafico che hai messo nel secondo post ha derivata in x=0 pari a −1, non 1 come scritto all'inizio.
Difatti c'è un errore, il testo riportava derivata in $0$ pari a $-1$. Ho corretto, scusami.
Ad ogni modo, qui la condizione di compatibilità è soddisfatta poichè $f=\cos(2 \pi x)$ ha integrale nullo in $[0,1]$
Per implementare quanto hai detto sopra ho fatto così
ho aggiunto un grado di libertà alla mia matrice $A$ (quella del sistema singolare). Ossia ho aggiunto una colonna a sinistra e una riga in fondo. In colonna ho tutti $1$ eccetto per l'ultima componente. Mentre in riga ho tutti $0$ eccetto per la penultima.
Per $6$ nodi risulta così
Poi ho aggiungo una riga al termine noto, con il valore che volevo dare al bordo $x=1$, ad esempio $2$.
Mi è bastato cambiare il codice iniziale con
G[-1,:]=np.zeros(M+2) G[:,-1]=np.ones(M+2) G[-1,-2]=1 G[-1,-1]=0 newf=np.append(f,2)
e poi risolvere il sistema lineare e ottengo il risultato desiderato, almeno in termini di condizione ai bordi
Se l'integrale di \(f\) è nullo, i flussi non dovrebbero sommare a zero?
Ma nel caso 1D la condizione di compatibilità è proprio che l'integrale di $f$ sia nullo sul dominio, visto che il bordo è fatto da due punti. O forse non ho capito cosa intendi, scusami.
Non può essere: se fosse un'equazione del calore, significherebbe dire che stai cercando un equilibrio in una situazione in cui pompi calore [dentro o fuori] da un lato e non produci niente all'interno. Come fa a esistere questo equilibrio?
Sinceramente non mi ero posto da un punto di vista fisico, su cui comunque concordo con te.
Ho trovato questo, in particolare la seconda risposta scrive ancora che la condizione di compatibilità è l'integrale nullo su $[0,1]$.
Mi pare però di capire che invece porre $\int_0^1 u(x) dx=0$ sia accettabile "fisicamente", cioè sto richiedendo che $u$ sia a media nulla sul dominio. Sinceramente questo problema mi ha sorpreso
Ho trovato questo, in particolare la seconda risposta scrive ancora che la condizione di compatibilità è l'integrale nullo su $[0,1]$.
Mi pare però di capire che invece porre $\int_0^1 u(x) dx=0$ sia accettabile "fisicamente", cioè sto richiedendo che $u$ sia a media nulla sul dominio. Sinceramente questo problema mi ha sorpreso
"feddy":
Sinceramente non mi ero posto da un punto di vista fisico, su cui comunque concordo con te.
Bene che siamo d'accordo con questo. Non è il metodo più rigoroso di ragionare, ma di solito da delle buone linee guida.
"feddy":
Ho trovato questo, in particolare la seconda risposta scrive ancora che la condizione di compatibilità è l'integrale nullo su $[0,1]$.
Lo dice per un'equazione con condizioni al bordo omogenee.
Puoi vedere la cosa nel modo seguente: parti da
\[
\begin{cases}
-u'' = f \\
u'(0) = a \\
u'(1) = b
\end{cases}
\]
e fai un cambio di variabile \(u = \hat u + \bar u\) dove \(\hat u\) soddisfa \(\hat u'(0) = \hat u'(1) = 0\) e \(\bar u'(0) = a\) e \(\bar u'(1) = b\).
Allora il tuo problema diventa
\[
\begin{cases}
-\hat u'' = f + \bar u'' = \hat f\\
\hat u'(0) = 0 \\
\hat u'(1) = 0
\end{cases}
\]
nell'incognita \(\hat u\), che adesso ha condizioni omogenee.
Adesso imponi la condizione di compatibilità \(\int_\Omega \hat f = 0\), e vedi che questa incorpora le condizioni al bordo, come dev'essere.
"feddy":
Mi pare però di capire che invece porre $\int_0^1 u(x) dx=0$ sia accettabile "fisicamente", cioè sto richiedendo che $u$ sia a media nulla sul dominio.
Richiedere la media nulla è uguale cambiare unità di misura, per certi versi. In fluidodinamica incomprimibile è pratica comune risolvere l'equazione di Navier-Stokes con la pressione a media nulla, perché tanto la pressione è definita a meno di una costante, quindi la puoi traslare come vuoi. Non è tanto una questione di senso fisico.
Scusa per il ritardo nella risposta. Ciò che hai scritto mi torna, e ti ringrazio!
Giusto per non lasciare la discussione in sospeso. Ho chiesto a chi mi ha assegnato il problema delucidazioni (visto che non vi era alcuna condizione aggiuntiva se non quelle di Neumann) e di fatto non si era accorto di aver assegnato un problema mal posto.
Giusto per non lasciare la discussione in sospeso. Ho chiesto a chi mi ha assegnato il problema delucidazioni (visto che non vi era alcuna condizione aggiuntiva se non quelle di Neumann) e di fatto non si era accorto di aver assegnato un problema mal posto.
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo