Neumann BC per l'equazione di Poisson (1D)

[feddy]

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?

[feddy]

Un modo per rendere il problema ben posto è quello di richiedere che valga

$ \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:

Codice:

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.

[Raptorista]

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.

[feddy]

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

[Raptorista]

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ù.

[feddy]

Imponendo $u(1)=1$ ad esempio, non avrei più il comportamente al bordo desiderato. Cioè sto proprio risolvendo un altro problema. Dove sto sbagliando?

[feddy]

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

[feddy]

Anzi no, sono riuscito

[Raptorista]

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.

[feddy]

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

Codice:

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