Equazione di Poisson.

Alvise Sommariva

Università degli Studi di Padova

Dipartimento di Matematica

Equazione di Poisson, caso univariato.

20 maggio 2014 Consideriamo il problema ai valori al contorno

\{ \begin{matrix}−y′′(x) + q(x)y(x) = g(x), x ∈ (a, b), \\ y(a) = α, y(b) = β \end{matrix}  (1)

Supponiamo ▶ [a, b] sia un intervallo limitato, ▶ q ∈ C^{(0)}([a, b]), q(x) ≥ 0 per ogni x ∈ [a, b], ▶ g ∈ C^{(0)}([a, b]).

Per una dimostrazione di esistenza, unicità e regolarità della soluzione y, si consulti [1, p.671].

Di seguito consideriamo quale esempio il caso particolare di (1)

\{ \begin{matrix}−y′′(x) + cy(x) = f (x), x ∈ (a, b), \\ y(a) = 0, y(b) = 0 \end{matrix}  (2)

con ▶ [a, b] sia un intervallo limitato, ▶ c > 0, ▶ f ∈ C^{(0)}([a, b]).

Questo problema si dice di Poisson univariato.

Supponiamo che la soluzione u del problema di Poisson univariato, appartenga a C^{4}([a, b]). Allora dallo sviluppo di Taylor centrato in x, abbiamo

y(x + h) = y(x) + hy^{(1)}(x) + \frac{h^{2}}{2} y^{(2)}(x) + \frac{h^{3}}{6} y^{(3)}(x) + \frac{h^{4}}{24}y^{(4)}(ξ_{1}),

y(x − h) = y(x) − hy^{(1)}(x) + \frac{h^{2}}{2} y^{(2)}(x) − \frac{h^{3}}{6} y^{(3)}(x) + \frac{h^{4}}{24}y^{(4)}(ξ_{2}),

e quindi sommando membro a membro, dal teorema della media, per qualche ξ ∈ (a, b),

y(x + h) + y(x − h) = y(x) + h^{2}y^{(2)}(x) + \frac{h^{4}}{12}y^{(4)}(ξ),

da cui l’appross. della derivata seconda con le differenze centrate

y^{(2)}(x) ≈ \frac{y(x + h) − 2y(x) + y(x − h)}{h^{2}}.  (3)

Alvise Sommariva  Equazione di Poisson.  4/ 1

Consideriamo ora la discretizzazione equispaziata dell’intervallo [a, b]

x_{i} = a + ih, 0 ≤ i ≤ n + 1

con

h = \frac{b − a}{n + 1}.

Allora, il problema di Poisson univariato, può essere discretizzato utilizzando le differenze centrate, ricavando

−\frac{y(x + h) − 2y(x) + y(x − h)}{h^{2}} +cy(x) ≈ −y^{(2)}(x)+cy(x) = f (x)

con y(t) = 0 se t non appartiene ad (a, b).

Essendo, come visto,

−\frac{y(x + h) − 2y(x) + y(x − h)}{h^{2}} +cy(x) ≈ −y^{(2)}(x)+cy(x) = f (x)

con y(t) = 0 se t non appartiene ad (a, b), riscrivendo le equazioni relativamente ai punti x = x_{i}, intendiamo approssimare y(x_{k}) con u_{k} risolvente per k = 1, …, n,

− \frac{u_{k+1} − 2u_{k} + u_{k−1}}{h^{2}} + cu_{k} = f (x_{k}),  (4)

dove u_{0} = u_{n+1} = 0.

Equazione di Poisson, caso univariato.

Tali equazioni lineari possono essere riscritte come

\frac{1}{h^{2}}(−u_{k+1} + 2u_{k} − u_{k−1} + ch^{2}u_{k}) = f (x_{k}),  (5)

con u_{0} = u_{n+1} = 0, ed essere convenientemente definite matricialmente come Au = b, dove A è la matrice tridiagonale

A = \frac{1}{h^{2}} ( \begin{matrix}2 + ch^{2} & −1 & 0 & … & 0 \\ −1 & 2 + ch^{2} & −1 & 0 & … \\ 0 & −1 & 2 + ch^{2} & −1 & … \\ … & … & … & … & … \\ 0 & … & 0 & −1 & −2 + ch^{2} \end{matrix} )

e b_{k} = f (x_{k}) per k = 1, …, n.

Alvise Sommariva  Equazione di Poisson.  7/ 1

La matrice

A = \frac{1}{h^{2}} ( \begin{matrix}2 + ch^{2} & −1 & 0 & … & 0 \\ −1 & 2 + ch^{2} & −1 & 0 & … \\ 0 & −1 & 2 + ch^{2} & −1 & … \\ … & … & … & … & … \\ 0 & … & 0 & −1 & −2 + ch^{2} \end{matrix} )

è

▶ simmetrica, e definita positiva (per il primo ed il terzo teorema di Gershgorin), con lo spettro degli autovalori σ(A) = {λ_{i}} contenuto in [c, c + 4h^{−2}];

▶ a predominanza diagonale stretta, essendo c > 0.

Si prova che

▶ se y è la soluzione del problema di Poisson e u_{k} ≈ y(x_{k}) con x_{k} = a + kh, si ha

\frac{∥\{u_{i}\} − \{y(x_{i})\}∥_{2}}{\sqrt{n}} = O(h^{2})

▶ se la soluzione y ∈ C^{(4)}([a, b]), allora

max_{i}|y(x_{i}) − u_{i}| ≤ \frac{h^{2}}{24}∥y^{(4)}∥_{∞}∥(x − a)(b − x)∥_{∞}

(cf. [1, p.678]); ▶ essendo la matrice tridiagonale, il sistema lineare può essere

risolto con il metodo di eliminazione gaussiana con complessità O(n).

▶ se c = 0, si può ancora mostrare che A è def. positiva e che per ogni h che σ(A) = {λ_{i}} è contenuto in [δ, c + 4h^{−2}] con δ > 0 indipendente da h.

Alvise Sommariva  Equazione di Poisson.  9/ 1

Equazione di Poisson, caso bivariato.

Consideriamo l’equazione di Poisson

\{ \begin{matrix}Δu(x, y) = f (x, y), (x, y) ∈ Ω = (0, 1) × (0, 1) \\ u(x, y) = g(x, y), (x, y) ∈ ∂Ω \end{matrix}  (6)

dove

Δu(x, y) := \frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}}.  (7)
Alvise Sommariva  Equazione di Poisson.  10/ 1

Equazione di Poisson, caso bivariato.

Definita la griglia di punti G = {(x_{i}, y_{j})}_{i,j=0,...,n+1}

x_{i} = ih, y_{j} = jh, h = 1/(n + 1), i, j = 0, …, n + 1

risulta evidente che per i = 0 o j = 0, i = n + 1 o j = n + 1 abbiamo un punto del bordo e quindi in virtù delle condizioni di Dirichlet in (6), il valore della soluzione u^{∗} è determinato.

Vediamo cosa succede quando il punto della griglia G è interno al quadrato Ω = (0, 1) × (0, 1), caratterizzati dall’avere i, j = 1, …, n.

Si mostra facilmente, utilizzando la formula di Taylor bivariata, in cui vengono tralasciati termini di ordine O(h^{2}), che si può approssimare Δu(x, y) con

\frac{u(x + h, y) + u(x − h, y) + u(x, y + h) + u(x, y − h) − 4u(x, y)}{h^{2}}

e quindi essendo, x_{i} = ih, y_{j} = jh, Δu(x_{i}, y_{j}) con

\frac{u(x_{i+1}, y_{j}) + u(x_{i−1}, y_{j}) + u(x_{i}, y_{j+1}) + u(x_{i}, y_{j−1}) − 4u(x_{i}, y_{j})}{h^{2}}

(8)

da cui la discretizzazione dell’equazione di Poisson

u(x_{i+1}, y_{j})+u(x_{i−1}, y_{j})+u(x_{i}, y_{j+1})+u(x_{i}, y_{j−1})−4u(x_{i}, y_{j}) = h^{2} f (x_{i}, y_{j}),

(9)

per i, j = 1, …, n.

Alvise Sommariva  Equazione di Poisson.  12/ 1

Equazione di Poisson, caso bivariato.

con le condizioni al contorno

u(x_{i}, y_{j}) = g(x_{i}, y_{j}), i = 0, j = 1, …, n  (10)
u(x_{i}, y_{j}) = g(x_{i}, y_{j}), i = n + 1, j = 1, …, n  (11)
u(x_{i}, y_{j}) = g(x_{i}, y_{j}), i = 1, …, n, j = 0, j = n + 1. (12)
Alvise Sommariva  Equazione di Poisson.  13/ 1

Figura: La molecola della discr. del Laplaciano avente centro (0.5, 0.5) e h = 0.25. Si ricordi di dividere ogni valore nella molecola per h^{2}.

.

Figura: Una griglia avente passo h = 1/3 relativamente al quadrato [0, 1] × [0, 1].

.

Alvise Sommariva  Equazione di Poisson.  14/ 1

Equazione di Poisson, caso bivariato.

Purtroppo, la descrizione del sistema lineare non è troppo chiara. Vediamola scritta matricialmente. Sia B la matrice n × n

B = ( \begin{matrix}−4 & 1 & 0 & … & 0 \\ 1 & −4 & 1 & 0 & … \\ 0 & 1 & −4 & 1 & … \\ … & … & … & … & … \\ 0 & … & 0 & 1 & −4 \end{matrix} )

ed I la matrice identica di ordine n del tipo

I = ( \begin{matrix}1 & 0 & 0 & … & 0 \\ 0 & 1 & 0 & … & 0 \\ 0 & 0 & 0 & … & 1 \end{matrix} ) .

Allora se b è il vettore ottenuto dai contributi dei termini dovuti a f e g in (6) e (9), definita la matrice a blocchi

A = ( \begin{matrix}B & I & 0 & … \\ I & B & I & … \\ … & … & … & … \\ … & 0 & I & B \end{matrix} )

si ricava che il sistema da risolvere è Au = b, usando ad esempio il metodo di Jacobi, Gauss-Seidel, SOR o il gradiente coniugato.

Equazione di Poisson, caso bivariato.

Per una implementazione della matrice di Poisson A, utilizziamo la

funzione makefish

f u n c t i o n mat = makefish ( siz ) ;

% make a P o i s s o n m a t r i x

leng = siz∗siz ;

dia = z e r o s ( siz , siz ) ;

off = −eye ( siz , siz ) ;

f o r i=1:siz , dia ( i , i ) =4; end ;
f o r i=1:siz −1, dia ( i , i+1)=−1; dia ( i+1,i )=−1; end ;

mat = z e r o s ( leng , leng ) ;

f o r ib =1:siz ,
mat (1+(ib−1)∗siz : ib∗siz ,1+( ib−1)∗siz : ib∗siz ) = dia ; end ;
f o r ib =1:siz −1,

  mat (1+(ib−1)∗siz : ib∗siz ,1+ ib∗siz : ( ib+1)∗siz ) = off ;

mat(1+ib∗siz : ( ib+1)∗siz ,1+( ib−1)∗siz : ib∗siz ) = off ; end ;
r e t u r n ;

Vediamone un esempio dalla shell di Matlab/Octave:

>> makefish ( 2 )

ans =
4  −1  −1  0
−1  4  0  −1
−1  0  4  −1
0  −1  −1  4

>>

Si vede subito dal 4 sulla diagonale che makefish non calcola A ma −A e dovremo tener conto di questo dettaglio nell’implementazione. Osserviamo che non è proprio facile determinare, fissato i, j, quali siano i punti adiacenti a (x_{i}, y_{j}) che essendo sul bordo hanno valore della soluzione noto a priori e quindi tali da contribuire attivamente al termine noto.

Equazione di Poisson: un esempio bivariato.

Facciamo un esempio sulla risoluzione dell’equazione di Poisson via metodo alle differenze con 5 punti. Sia Ω = [0, 1] × [0, 1], h = 1/3 e siano

P_{i,j} = (ih, jh), i, j = 0, 1, 2, 3.

E’ chiaro che per ▶ per i = 0 i punti P_{0,j} sono sull’asse x = 0 (cioè l’asse y), ▶ per i = 3 i punti P_{3,j} sono sull’asse x = 1, ▶ per j = 0 i punti P_{i,0} sono sull’asse y = 0 (cioè l’asse x) ▶ per j = 3 i punti P_{i,3} sono sull’asse y = 1.

Date le condizioni al contorno, la soluzione in questi punti è nota ed è uguale a u_{i,j} = g(x_{i}, y_{j}).

I rimanenti punti P_{i,j}, con i, j = 1, 2 sono interni a Ω ed è u(x_{i+1}, y_{j})+u(x_{i−1}, y_{j})+u(x_{i}, y_{j+1})+u(x_{i}, y_{j−1})−4u(x_{i}, y_{j}) = h^{2} f (x_{i}, y_{j}),

(13)

Analizziamo caso per caso queste equazioni:

▶ Nel caso i = 1, j = 1 si ha

u(x_{2}, y_{1})+u(x_{0}, y_{1})+u(x_{1}, y_{2})+u(x_{1}, y_{0})−4u(x_{1}, y_{1}) = h^{2} f (x_{1}, y_{1}),

u(x_{0}, y_{1}) = g(x_{0}, y_{1}), u(x_{1}, y_{0}) = g(x_{1}, y_{0}).

Portando questi due termini a secondo membro otteniamo u(x_{2}, y_{1})+u(x_{1}, y_{2})−4u(x_{1}, y_{1}) = h^{2} f (x_{1}, y_{1})−g(x_{0}, y_{1})−g(x_{1}, y_{0}).

Equazione di Poisson: un esempio bivariato.

▶ Nel caso i = 2, j = 1 si ha

u(x_{3}, y_{1})+u(x_{1}, y_{1})+u(x_{2}, y_{2})+u(x_{2}, y_{0})−4u(x_{2}, y_{1}) = h^{2} f (x_{2}, y_{1}),

u(x_{3}, y_{1}) = g(x_{3}, y_{1}), u(x_{2}, y_{0}) = g(x_{2}, y_{0})

portando questi due termini a secondo membro otteniamo u(x_{1}, y_{1})+u(x_{2}, y_{2})−4u(x_{2}, y_{1}) = h^{2} f (x_{1}, y_{1})−g(x_{3}, y_{1})−g(x_{2}, y_{0}).

▶ Nel caso i = 1, j = 2 si ha u(x_{2}, y_{2})+u(x_{0}, y_{2})+u(x_{1}, y_{3})+u(x_{1}, y_{1})−4u(x_{1}, y_{2}) = h^{2} f (x_{1}, y_{2}), ed essendo

u(x_{0}, y_{2}) = g(x_{0}, y_{2}), u(x_{1}, y_{3}) = g(x_{1}, y_{3})

portando questi due termini a secondo membro otteniamo u(x_{2}, y_{2})+u(x_{1}, y_{1})−4u(x_{1}, y_{2}) = h^{2} f (x_{1}, y_{2})−g(x_{0}, y_{2})−g(x_{1}, y_{3}).

Alvise Sommariva  Equazione di Poisson.  21/ 1

▶ Nel caso i = 2, j = 2 si ha u(x_{3}, y_{2})+u(x_{1}, y_{2})+u(x_{2}, y_{3})+u(x_{2}, y_{1})−4u(x_{2}, y_{2}) = h^{2} f (x_{2}, y_{2}), ed essendo

u(x_{3}, y_{2}) = g(x_{3}, y_{2}), u(x_{2}, y_{3}) = g(x_{2}, y_{3})

portando questi due termini a secondo membro otteniamo u(x_{1}, y_{2})+u(x_{2}, y_{1})−4u(x_{2}, y_{2}) = h^{2} f (x_{2}, y_{2})−g(x_{3}, y_{2})−g(x_{2}, y_{3}).

Equazione di Poisson: un esempio bivariato.

Poniamo ora

b_{1} := h^{2} f (x_{1}, y_{1}) − g(x_{0}, y_{1}) − g(x_{1}, y_{0}),
b_{2} := h^{2} f (x_{2}, y_{1}) − g(x_{3}, y_{1}) − g(x_{2}, y_{0}),
b_{3} := h^{2} f (x_{1}, y_{2}) − g(x_{0}, y_{2}) − g(x_{1}, y_{3}),
b_{4}; = h^{2} f (x_{2}, y_{2}) − g(x_{3}, y_{2}) − g(x_{2}, y_{3}),

ordiniamo i punti da sinistra a destra, e dal basso verso l’alto (ordine lessicografico)

P_{1} = (x_{1}, y_{1}), P_{2} = (x_{2}, y_{1}), P_{3} = (x_{1}, y_{2}), P_{4} = (x_{2}, y_{2}), e infine poniamo

u_{1} = u(x_{1}, y_{1}), u_{2} = u(x_{2}, y_{1}), u_{3} = u(x_{1}, y_{2}), u_{4} = u(x_{2}, y_{2}),

Alvise Sommariva  Equazione di Poisson.  23/ 1

Equazione di Poisson: un esempio bivariato.

ottenendo così

u_{2} + u_{3} − 4u_{1} = b_{1},
u_{1} + u_{4} − 4u_{2} = b_{2},
u_{4} + u_{1} − 4u_{3} = b_{3},
u_{3} + u_{2} − 4u_{4} = b_{4},

da cui posto

A = ( \begin{matrix}−4 & 1 & 1 & 0 \\ 1 & −4 & 0 & 1 \\ 1 & 0 & −4 & 1 \\ 0 & 1 & 1 & −4 \end{matrix} )

basta risolvere il sistema Au = b per ottenere u(x_{1}, y_{1}), u(x_{2}, y_{1}), u(x_{1}, y_{2}), u(x_{2}, y_{2}).

Alvise Sommariva  Equazione di Poisson.  24/ 1

Notiamo che

>> makefish ( 2 )

ans =
4  −1  −1  0
−1  4  0  −1
−1  0  4  −1
0  −1  −1  4

>>

calcola proprio −A, mentre il termine noto b può essere facilmente calcolato dopo aver notato che

1. i termini −g(x_{i}, y_{0}) sono presenti nelle componenti b_{i};

2. posto n = 1/h, per i = 1, …, n − 1, i termini −g(x_{i}, y_{n}) sono presenti nelle componenti b_{(n−1)^{2}−(n−1)+i};

3. per j = 1, …, n − 1, i termini −g(x_{0}, y_{j}) sono presenti nelle componenti b_{s} con s ≡ 1 mod n − 1;

4. per j = 1, …, n − 1, i termini −g(x_{n}, y_{j}) sono presenti nelle componenti b_{s} con s ≡ 0 mod n − 1.

Equazione di Poisson: una implementazione.

Vediamo ora un’implementazione del metodo sopra descritto, detto per ovvi motivi a 5 punti (cf. (9)). Risulta importante ricordare la seguente stima dell’errore Teorema. Se u è soluzione dell’equazione di Poisson (6) ed è

almeno 4 volte differenziabile con continuità nel quadrato Ω := [0, 1] × [0, 1] ed u_{h} l’approssimazione ottenuta col metodo alle differenze con 5 punti, utilizzando una griglia G = {(x_{i}, y_{j})} con x_{i} = i h, y_{j} = j h, h = 1/(n + 1) allora

|u(x_{i}, y_{j}) − u_{h}(x_{i}, y_{j})| ≤ ch^{2}

con

c = (1/24) ( max_{(x,y)∈Ω} |\frac{∂^{4}u(x, y)}{∂x^{4}} | + max_{(x,y)∈Ω} |\frac{∂^{4}u(x, y)}{∂y^{4}} |)

Ci si aspetta quindi dai test numerici che effettueremo un errore dell’ordine di h^{2}.

Alvise Sommariva  Equazione di Poisson.  26/ 1

Equazione di Poisson: una implementazione.

Salviamo in poisson5pts.m la funzione

f u n c t i o n u=poisson5pts ( n , f , g_left , g_right , g_down , g_up )
A=−makefish ( n ) ;

h=1/(n+1) ;
x=(h : h:1−h ) ’ ; y=x ;
% SOLUZIONI IN BASSO .

x_loc=x ;

y_loc=z e r o s ( s i z e ( x_loc ) ) ;

b_down=f e v a l ( g_down , x_loc , y_loc ) ;

% SOLUZIONI IN ALTO.

y_loc=ones ( s i z e ( x_loc ) ) ;
b_up=f e v a l ( g_up , x_loc , y_loc ) ;

% SOLUZIONI A SINISTRA .

y_loc=x_loc ; x_loc=z e r o s ( s i z e ( x_loc ) ) ;
b_left=f e v a l ( g_left , x_loc , y_loc ) ;

% SOLUZIONI A DESTRA

x_loc=ones ( s i z e ( x_loc ) ) ;
b_right=f e v a l ( g_right , x_loc , y_loc ) ;

% COMPOSIZIONE TERMINE NOTO.

b1=b_down ;

% PRIMA RIGA IN BASSO .

b1 ( 1 )=b1 ( 1 )+b_left ( 1 ) ; b1 ( n )=b1 ( n )+b_right ( 1 ) ;

% PRIMA RIGA IN ALTO.

bn=b_up ; bn ( 1 )=bn ( 1 )+b_left ( n ) ; bn ( n )=bn ( n )+b_right ( n ) ;

% RIGHE INTERMEDIE .
bj = [ ] ;
f o r j =2:( n−1)
bjloc=z e r o s ( n , 1 ) ;

     bjloc ( 1 )=bjloc ( 1 )+b_left ( j ) ;
     bjloc ( n )=bjloc ( n )+b_right ( j ) ;
     bj=[bj ; bjloc ] ;

end

b=[b1 ; bj ; bn ] ;

% GRIGLIA LESSICOGRAFICA (METODO STANDARD) .
lunghezza_x=l e n g t h ( x ) ; X = [ ] ; Y = [ ] ;

 f o r index =1: lunghezza_x

X=[X ; x ] ;
Y=[Y ; y ( index ) ∗ones ( s i z e ( x ) ) ] ;

end

fXY=f e v a l ( f , X , Y ) ;

b_f=(h ˆ 2) ∗fXY ; b=b_f−b ; u=A\b ;

Equazione di Poisson: una implementazione.

Salviamo in demopoisson5pts.m la demo

% MODIFIED VERSION : MARCH 13 , 2008.

demo_example =2;
switch  demo_example
case 1

f=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_down=inline ( ’ ones ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_up=inline ( ’ ones ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_left=inline ( ’ ones ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_right=inline ( ’ ones ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
solution=inline ( ’ ones ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;

case 2

f=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_down=inline ( ’ exp ( p i ∗x ) ’ , ’ x ’ , ’ y ’ ) ;
g_up=inline ( ’−exp ( p i ∗x ) ’ , ’ x ’ , ’ y ’ ) ;
g_left=inline ( ’ cos ( p i ∗y ) ’ , ’ x ’ , ’ y ’ ) ;
g_right=inline ( ’ ( ( exp ( 1 ) ) ˆ p i ) ∗ cos ( p i ∗y ) ’ , ’ x ’ , ’ y ’ ) ;
solution=inline ( ’ ( exp ( p i ∗x ) ) . ∗ cos ( p i ∗y ) ’ , ’ x ’ , ’ y ’ ) ;

Equazione di Poisson: una implementazione.

Salviamo in demopoisson5pts.m la demo

case 3

f=inline ( ’ (−2∗( p i ˆ 2) ) ∗ s i n ( p i ∗x ) . ∗ s i n ( p i ∗y ) ’ , ’ x ’ , ’ y ’ ) ;
g_down=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_up=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_left=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_right=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
solution=inline ( ’ ( s i n ( p i ∗x ) ) . ∗ s i n ( p i ∗y ) ’ , ’ x ’ , ’ y ’ ) ;

otherwise

f=inline ( ’ ones ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_down=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_up=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_left=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
g_right=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ y ’ ) ;
solution=inline ( ’ ( s i n ( p i ∗x ) ) . ∗ s i n ( p i ∗y ) ’ , ’ x ’ , ’ y ’ ) ;

end

f o r index =2:5
n=2ˆ index ; h=1/(n+1) ; x=(h : h:1−h ) ’ ; y=x ;
[ X , Y]= meshgrid ( x , y ) ; X=X ’ ; Y=Y ’ ;
% VETT. SOL , NEI PUNTI DELLA GRIGLIA ORD. IN LESSICOGR . .

u=poisson5pts ( n , f , g_left , g_right , g_down , g_up ) ;

% USO RESHAPE COSI ’ LA SOLUZIONE HA LE STESSE DIMENSIONI DELLE MATRICI X, Y .
Z=( r e s h a p e ( u , n , n ) ) ;

i f demo_example <=3
     V=f e v a l ( solution , X , Y ) ;

err ( index )=norm ( V ( : )−Z ( : ) , inf ) ;

     i f index = 1

f p r i n t f ( ’\n \t [ n ] : %4.0 f [ ERR ] : %2.2 e ’ , n , err ( index ) ) ;

e l s e

f p r i n t f ( ’\n \t [ n ] : %4.0 f [ ERR ] : %2.2 e [ RATIO ] : %2.2 f ’ , …

             n , err ( index ) , err ( index −1)/ err ( index ) ) ;

end end

end

s u r f ( X , Y , Z ) ;

Equazione di Poisson: descrizione dell’implementazione.

Alcune osservazioni sui codici Matlab/Octave appena esposti.

▶ Posto h = \frac{1}{n+1}, allora la matrice B ∈ ℝ^{n×n} mentre il termine noto b e il vettore soluzione u apparterranno a ℝ^{n^{2}}. dobbiamo risolvere un sistema lineare Au = b con

A = ( \begin{matrix}B & I & 0 & … & 0 \\ I & B & I & … & 0 \\ … & … & … & … & … \\ … & 0 & I & B & I \end{matrix} )

e b un vettore i cui contributi dipendono dai valori che hanno sul bordo le funzioni f e g che definiscono l’equazione di Poisson

\{ \begin{matrix}− [\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}}] = f (x, y), (x, y) ∈ Ω \\ u(x, y) = g(x, y), (x, y) ∈ ∂Ω \end{matrix}  (14)

Equazione di Poisson: descrizione dell’implementazione.

▶ Dobbiamo tener conto che makefish non calcola A ma Ā = −A ed è per questo che scriviamo A=-makefish(n).

▶ (Per i più esperti) La funzione meshgrid crea a partire da due vettori di numeri reali x, y, le ascisse X e le ordinate Y dei punti facenti parte della griglia generata da x ed y. In particolare la matrice di punti

(x_{1}, y_{1}) (x_{2}, y_{1}) … (x_{n}, y_{1})
…  . . . . . . . . .
(x_{1}, y_{n}) (x_{2}, y_{n}) … (x_{n}, y_{n})

viene descritta tramite la griglia di ascisse e ordinate e cioè come

X = ( \begin{matrix}x_{1} & x_{2} & … & x_{n} \\ … & … & … \\ x_{1} & x_{2} & … & x_{n} \end{matrix} )

e

Y = ( \begin{matrix}y_{1} & y_{1} & … & y_{1} \\ y_{2} & y_{2} & … & y_{2} \\ y_{3} & y_{3} & … & y_{3} \end{matrix} )

Alvise Sommariva  Equazione di Poisson.  33/ 1

Equazione di Poisson: descrizione dell’implementazione.

Così ad esempio

>> h =1/3;
>> x=h : h:1−h ;
>> x

 x =

0.3333  0.6667
>> y=x ;
>> [ X , Y]= meshgrid ( x , y )

 X =

0.3333  0.6667
0.3333  0.6667

 Y =

0.3333  0.3333
0.6667  0.6667

>>

descrivendone le coordinate x, y.

Alvise Sommariva  Equazione di Poisson.  34/ 1

Equazione di Poisson: esempio 1.

Consideriamo l’equazione di Poisson nel quadrato unitario Ω = [0, 1] × [0, 1]

\{ \begin{matrix}\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}} = 0, (x, y) ∈ Ω \\ u(x, y) = 1, (x, y) ∈ ∂Ω \end{matrix}  (15)

la cui soluzione è u(x, y) = 1. Non è difficile osservare che

max_{(x,y)∈Ω} |\frac{∂^{4}u(x, y)}{∂x^{4}} | = 0, max_{(x,y)∈Ω} |\frac{∂^{4}u(x, y)}{∂y^{4}} | = 0

e quindi ci si aspetta che per qualsiasi h si abbia un errore dell’ordine della precisione di macchina.

Equazione di Poisson: esempio 1.

Lanciamo da shell il primo esempio della demo, ottenendo

>> demopoisson5pts

[ n ] :  4 [ ERR ] : 4 . 4 4 e−016 [ RATIO ] : 163181252362740.50
[ n ] :  8 [ ERR ] : 1 . 1 1 e−015 [ RATIO ] : 0 . 4 0
[ n ] :  16 [ ERR ] : 1 . 5 5 e−015 [ RATIO ] : 0 . 7 1
[ n ] :  32 [ ERR ] : 6 . 6 6 e−015 [ RATIO ] : 0 . 2 3

>>

Per ratio si intende il rapporto dell’errore tra due iterate successive.

Alvise Sommariva  Equazione di Poisson.  36/ 1

\{ \begin{matrix}\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}} = 0, (x, y) ∈ Ω \\ u(x, 0) = exp(π x), x ∈ [0, 1] \\ u(x, 1) = −exp(π x), x ∈ [0, 1] \\ u(0, y) = cos(π y), y ∈ [0, 1] \\ u(1, y) = exp(π) · cos(π y), y ∈ [0, 1] \end{matrix}  (16)

la cui soluzione è u(x, y) = exp(π x) · cos(π y).

Equazione di Poisson: esempio 2.

Figura: Soluzione del problema 2.

Alvise Sommariva  Equazione di Poisson.  38/ 1

Si nota subito che per x, y ∈ [0, 1] si ha

|\frac{∂^{4}u}{∂x^{4}}(x, y)| = π^{4} |exp(π x) · cos(π y)| ≤ π^{4} exp(π) ≈ 2254.1

|\frac{∂^{4}u}{∂y^{4}}(x, y)| = π^{4} |exp(π x) · cos(π y)| ≤ π^{4} exp(π) ≈ 2254.1

da cui

c ≤ (1/24) · 2254.1 · 2 ≈ 187.8428.

Quindi quale maggiorazione dell’errore assoluto in norma infinito, per n = 3, 7, 15, 31, avremo i valori immagazzinati qui sotto nel vettore err

Equazione di Poisson: esempio 2.

>> format short e
>> c =187.8428; err = [ ] ;
>> f o r n =2:5 , N=2ˆ n−1; h=1/(N+1) ; h2=h ˆ 2 ; err=[err ; c∗h2 ] ; end

>> err
 err =

1.1740 e+001
2.9350 e+000
7.3376 e−001
1.8344 e−001

>>

Lanciando la demo demopoisson5pts, per demoexample=2,

abbiamo

>> demopoisson5pts

[ n ] :  4 [ ERR ] : 9 . 7 5 e−002 [ RATIO ] : 0 . 7 4
[ n ] :  8 [ ERR ] : 3 . 2 0 e−002 [ RATIO ] : 3 . 0 4
[ n ] :  16 [ ERR ] : 9 . 0 5 e−003 [ RATIO ] : 3 . 5 4
[ n ] :  32 [ ERR ] : 2 . 4 5 e−003 [ RATIO ] : 3 . 6 9

>>

Alvise Sommariva  Equazione di Poisson.  40/ 1

1. Come ci si aspettava la maggiorazione è realizzata, ma purtroppo come stima è abbastanza conservativa.

2. Nella colonna [RATIO] abbiamo indicato il rapporto e_{2h}/e_{h} dove e_{h} è l’errore assoluto compiuto dal metodo a 5 punti con passo h (ovvero la quantità esposte nella colonna [ERR] nella stessa riga di h). Il fatto che la ratio sia 4 non è sorprendente. Infatti se l’errore decresce come h^{2} si può supporre che sia e_{h} ≈ ĉh^{2} per qualche ĉ indipendente da h e quindi

\frac{e_{2h}}{e_{h}} ≈ \frac{ĉ(2h)^{2}}{ĉh^{2}} ≈ 4.

Equazione di Poisson: esempio 3.

Consideriamo un la risoluzione dell’equazione di Poisson nel quadrato unitario Ω = [0, 1] × [0, 1] con un metodo alle differenze.

\{ \begin{matrix}\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}} = (−2 π^{2}) sin(π x) sin(π y), (x, y) ∈ Ω \\ u(x, 0) = 0, x ∈ [0, 1] \\ u(x, 1) = 0, x ∈ [0, 1] \\ u(0, y) = 0, y ∈ [0, 1] \\ u(1, y) = 0, y ∈ [0, 1] \end{matrix}  (17)

Equazione di Poisson: esempio 3.

Figura: Soluzione del problema 3.

Alvise Sommariva  Equazione di Poisson.  43/ 1

Equazione di Poisson: esempio 3.

la cui soluzione è u(x, y) = sin(π x) sin(π y). Ripetendo la stima basata sulle derivate quarte della soluzione u abbiamo facilmente

c ≤ (1/24) 2 π^{4} ≈ 8.1174

da cui

>> c =8.1174;
>> err = [ ] ; f o r n =2:5 , N=2ˆ n ; N=N−1; h=1/(N+1) ; h2=h ˆ 2 ; err=[err ; c∗h2 ] ; end
>> format short e

>> err
 err =

5.0734 e−001
1.2683 e−001
3.1709 e−002
7.9271 e−003

>>

Lanciando la demo demopoisson5pts, per demoexample=3,

abbiamo

>> demopoisson5pts

[ n ] :  4 [ ERR ] : 3 . 0 4 e−002 [ RATIO ] : 2 . 3 9
[ n ] :  8 [ ERR ] : 9 . 9 1 e−003 [ RATIO ] : 3 . 0 6
[ n ] :  16 [ ERR ] : 2 . 8 3 e−003 [ RATIO ] : 3 . 5 1
[ n ] :  32 [ ERR ] : 7 . 5 4 e−004 [ RATIO ] : 3 . 7 5

>>

Rispetto al caso precedente la stima è più precisa, e la ratio di circa 4 ci dice che la convergenza è ancora dell’ordine di h^{2}.

Risolvere l’equazione di Poisson nel quadrato unitario Ω = [0, 1] × [0, 1] con un metodo alle differenze. Quali passi scegliere h = 1/4, h = 1/8, h = 1/16, h = 1/32. Qual’è la soluzione esatta? Con quale fattore decrescono gli errori?

\{ \begin{matrix}\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}} = 4, (x, y) ∈ Ω \\ u(x, 0) = (x − 1)x, x ∈ [0, 1] \\ u(x, 1) = (x − 1)x, x ∈ [0, 1] \\ u(0, y) = (y − 1)y, y ∈ [0, 1] \\ u(1, y) = (y − 1)y, y ∈ [0, 1] \end{matrix}  (18)

Bibliografia

V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990.

Alvise Sommariva  Equazione di Poisson.  46/ 1