Differenze finite per l`equazione di Laplace

Differenze finite per l’equazione di Laplace
Lucia Gastaldi
Dipartimento di Matematica,
http://dm.ing.unibs.it/gastaldi/
Laboratorio
Lucia Gastaldi
Differenze finite
Equazione di Laplace
Outline
1
Equazione di Laplace
Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Lucia Gastaldi
Differenze finite
Equazione di Laplace
Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Equazione di Laplace
Problema
Sia Ω =]0, a[×]0, b[. Dati f : Ω → R e g : Ω → R, si consideri il
seguente problema
−∆u = f in Ω
u=g
su ∂Ω
Dati Nx e Ny , poniamo hx =
a
b
e hy =
.
(Nx + 1)
(Ny + 1)
Notazioni
xi = ihx , per i = 0, 1, · · · , Nx + 1,
yj = jhy , per j = 0, 1, · · · , Ny + 1,
ui,j ≈ u(xi , yj ), per i = 0, 1, · · · , Nx + 1, j = 0, 1, · · · , Ny + 1.
Lucia Gastaldi
Differenze finite
Equazione di Laplace
Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Differenze finite
Discretizzazione delle derivate seconde
u(x − hx , y ) − 2u(x, y ) + u(x + hx , y )
hx2
u(x,
y
−
h
)
−
2u(x,
y ) + u(x, y + hy )
y
uyy (x, y ) ≈ δy2 u(x, y ) :=
2
hy
uxx (x, y ) ≈ δx2 u(x, y ) :=
Teorema
Se u è di classe C4 allora si dimostrano le seguenti stime dell’errore
di discretizzazione:
h2
uxx (x, y ) − δx2 u(x, y ) = − x uxxxx (x, y ) + o(hx2 )
12
uyy (x, y ) −
δy2 u(x, y )
hy2
= − uyyyy (x, y ) + o(hy2 )
12
Lucia Gastaldi
Differenze finite
Equazione di Laplace
Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Sistema lineare associato alla discretizzazione
Si ottengono Nx × Ny equazioni (per i = 1, · · · , Nx j = 1, · · · , Ny )
−
ui−1,j − 2ui,j + ui+1,j
ui,j−1 − 2ui,j + ui,j+1
−
= f (xi , yj ),
2
hx
hy2
alle quali si aggiungono le condizioni al bordo:
ui,0 = g (xi , 0)
ui,Ny +1 = g (xi , b)
u0,j = g (0, yj )
uNx +1,j = g (a, yj )
Lucia Gastaldi
per
per
per
per
i
i
j
j
= 0, · · ·
= 0, · · ·
= 0, · · ·
= 0, · · ·
Differenze finite
, Nx
, Nx
, Ny
, Ny
+ 1,
+ 1,
+ 1,
+ 1.
Equazione di Laplace
Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Sistema lineare
Vettore soluzione
u vettore di N = Nx × Ny componenti: u(k) = ui,j essendo
k = (j − 1) ∗ Nx + i;
ricostruzione della soluzione: ui,j = u(k) per
k −1
j = int(
) + 1, i = k − (j − 1) ∗ Nx .
Nx
Matrice delle differenze finite
La matrice del sistema A è una matrice
pentadiagonale con la struttura riportata in
figura.
Se hx = hy = h, ed il dominio è il quadrato
[−1, 1] × [−1, 1], si costruisce con il comando
A=delsq(numgrid(’S’,n+2))/hˆ2
Lucia Gastaldi
Differenze finite
Matrice delle differenze finite Nx=Ny=10
0
10
20
30
40
50
60
70
80
90
100
0
10
20
30
40
50
nz = 460
60
70
80
90
100
Equazione di Laplace
Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Termine noto I
b vettore di N = Nx × Ny componenti tale che b(k) = f (ihx , jhy )+
condizioni al bordo
Poniamo Nx = Nx , Ny = Ny , hx = hx , hy = hy




b(1:Nx) = 



1)
f (x1 , y1 ) + g (0,y
+ g (xh12,0)
hx2
y
···
f (xi , y1 ) + g (xhi2,0)
y
···
g (x ,0)
1)
f (xNx , y1 ) + g (a,y
+ hN2x
h2
x
Lucia Gastaldi
Differenze finite
y




,



Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Equazione di Laplace
Termine noto II
Per j = 2, Ny − 1




b(Nx*(j-1)+1:Nx) = 


Lucia Gastaldi
g (0,yj )
hx2
f (x1 , yj ) +
···
f (xi , yj )
···
f (xNx , yj ) +
Differenze finite
g (a,yj )
hx2







Equazione di Laplace
Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Termine noto III




b(Nx*(Ny-1)+(1:Nx)) = 



g (0,yN )
f (x1 , yNy ) + h2 y +
x
···
f (xi , yNy ) + g (xhi2,b)
y
···
f (xNx , yNy ) +
Lucia Gastaldi
Differenze finite
g (a,yNy )
hx2
g (x1 ,b)
hy2
+
g (xNx ,b)
hy2








Equazione di Laplace
Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Costruzione del termine noto in Matlab
Assegno la f : b(k) = f (ihx , jhy )
k=0;
for j=1:Ny
for i=1:Nx
k=k+1; f(k)=effe(x(i),y(j));
end
end
Impongo le condizioni al bordo
for i=1:Nx
f(i)=f(i)+g(i*hx,0)/hy^2;
f(Nx*(Ny-1)+i)=f(Nx*(Ny-1)+i)+g(i*hx,b)/hy^2;
end
for j=1:Ny
f((j-1)*Nx+1)=f((j-1)*Nx+1)+g(0,j*hy)/hx^2;
f(j*Nx)=f(j*Nx)+g(a,j*hy)/hx^2;
end
Lucia Gastaldi
Differenze finite
Equazione di Laplace
Problema e notazioni
Discretizzazione con differenze finite
Costruzione della matrice e del termine noto
Ricostruzione della soluzione
Ricostruzione della soluzione
Per ricostruire la soluzione si vuole ottenere una matrice Z di
dimensioni Ny × Nx che contiene i valori nei punti (xi , yj ).
Tenendo conto di come è stato costruito il vettore u si puo’
procedere a costruire la matrice Z riga per riga nel modo seguente:
for j=1:n
Z(j,:)=u((j-1)*n+1:j*n)’;
end
Come tenere conto delle condizioni al bordo per rappresentare
meglio la soluzione?
Lucia Gastaldi
Differenze finite