Esercitazioni di MATLAB: Lezione del 17 Marzo 2008

Soluzione di problemi ellittici in dimensione 1

Solutore del BVP:

File bvp.m

:open:
 1%
 2% Genero la matrice dei coefficienti e il temine noto del 
 3% problema ellittico
 4%
 5%  u''(x) = f(x)     su (0,1)
 6%
 7%  u(0) = u0 u(1)=u1
 8%
 9%  N      = numero i suddivisioni dell'intervallo [0,1]
10%  f      = funzione termine sorgente
11%  u0, u1 = condizioni al contorno
12%
13%  A = matrice dei coefficienti
14%  b = termine noto
15%
16%  la soluzione è posta nei vettori colonna X e Y che contengono
17%  i valori nodali nella mesh compreso i termini al bordo
18%
19%  X = X(i) coodinata x del nodo (i)
20%  Z = Z(i) coodinata z della soluzione nel nodo (i)
21%
22%         +                 +
23%    1    |                 |
24%  -----  |  1    -2     1  |  u_i = f_i
25%   h^2   |                 |
26%         +                 +
27%
28function [A,b,X,Y] = bvp(N,f,u0,u1)
29  %
30  h   = 1/N ; % dimensione della mesh
31  neq = N-1 ;
32
33  % coefficienti dello stencil
34  h2 = h^2 ;
35
36  % costruzione della matrice dei coefficienti
37  A = sparse( neq, neq ) ;
38  b = zeros ( neq, 1   ) ;
39
40  for i=1:neq
41    x = i*h ;
42
43    % termine noto
44    b(i) = h2*feval( f, x ) ;
45
46    % coefficienti dello stencil
47    A(i,i) = -2 ;
48      
49    if i > 1
50      A(i,i-1) = 1 ;
51    else
52      b(i) = b(i) - u0 ;      
53    end
54
55    if i < neq
56      A(i,i+1) = 1 ;
57    else
58      b(i) = b(i) - u1 ;      
59    end
60  end
61  
62  % calcolo la soluzione
63  SOL = A\b ;
64
65  % riempio i vettori X e Y con la soluzione
66  Y = [ u0 ; SOL ; u1 ] ;
67  X = [ 0:h:1 ]' ;
68
69end

Script per testare il codice:

File testBVP.m

:open:
 1function testBVP
 2
 3N = 10 ;
 4[ A, b, X1, Y1 ] = bvp( N, @ff, ye(0), ye(1) ) ;
 5N = N*2 ;
 6[ A, b, X2, Y2 ] = bvp( N, @ff, ye(0), ye(1) ) ;
 7N = N*2 ;
 8[ A, b, X3, Y3 ] = bvp( N, @ff, ye(0), ye(1) ) ;
 9N = N*2 ;
10[ A, b, X4, Y4 ] = bvp( N, @ff, ye(0), ye(1) ) ;
11N = N*2 ;
12[ A, b, X5, Y5 ] = bvp( N, @ff, ye(0), ye(1) ) ;
13
14EX1 = ye(X1) ;
15EX2 = ye(X2) ;
16EX3 = ye(X3) ;
17EX4 = ye(X4) ;
18EX5 = ye(X5) ;
19
20E1 = norm( Y1 - EX1, inf ) ;
21E2 = norm( Y2 - EX2, inf ) ;
22E3 = norm( Y3 - EX3, inf ) ;
23E4 = norm( Y4 - EX4, inf ) ;
24E5 = norm( Y5 - EX5, inf ) ;
25
26disp( sprintf( 'E1 = %g', E1 ) ) ; 
27disp( sprintf( 'E2 = %g', E2 ) ) ; 
28disp( sprintf( 'E3 = %g', E3 ) ) ; 
29disp( sprintf( 'E4 = %g', E4 ) ) ; 
30disp( sprintf( 'E5 = %g', E5 ) ) ; 
31
32disp( sprintf( 'order = %g', log( E1/E2 ) / log( 2 ) ) ) ;
33disp( sprintf( 'order = %g', log( E2/E3 ) / log( 2 ) ) ) ;
34disp( sprintf( 'order = %g', log( E3/E4 ) / log( 2 ) ) ) ;
35disp( sprintf( 'order = %g', log( E4/E5 ) / log( 2 ) ) ) ;
36
37end
38
39function res = ff(x)
40  res = 100*sin(10*x) ;
41end
42
43function res = ye(x)
44  res = -sin(10*x) ;
45end

Soluzione di problemi ellittici in dimensione 2

Solutore del problema ellittico:

File elliptic.m

:open:
  1%
  2% Genero la matrice dei coefficienti e il temine noto del 
  3% problema ellittico
  4%
  5%  u(x,y)    + u(x,y)      = f(x,y)     su (0,1) x (0,1)
  6%        xx          yy
  7%
  8%  u(x,y) = g(x,y)  su bordo (0,1) x (0,1)
  9%
 10%  N = numero i suddivisioni dell'intervallo [0,1]
 11%  f = funzione termine sorgente
 12%  g = funzione che implementa le condizioni al contorno
 13%
 14%  A = matrice dei coefficienti
 15%  b = termine noto
 16%
 17%  la soluzione � posta nelle metrici X, Y, e Z che contengono
 18%  i valori nodali nella mesh compreso i termini al bordo
 19%
 20%  X = X(i,j) coodinata x del nodo (i,j)
 21%  Y = Y(i,j) coodinata y del nodo (i,j)
 22%  Z = Z(i,j) coodinata z della soluzione nel nodo (i,j)
 23%
 24%
 25%         +                 +
 26%         |        1        |
 27%    1    |                 |
 28%  -----  |  1    -4     1  |  u_ij = f_ij
 29%   h^2   |                 |
 30%         |        1        |
 31%         +                 +
 32%
 33function [A,b,X,Y,Z] = elliptic(N,f,g)
 34  %
 35  h = 1/N ; % dimensione della mesh
 36  
 37  % numerazione delle equazioni
 38  EQ = zeros(N+1,N+1) ;
 39  
 40  % numero le equazioni (solo punti interni)
 41  EQ(2:end-1,2:end-1) = [1:N-1]'* ones(1,N-1) + ...
 42                        ones(N-1,1)*[0:N-2]*(N-1) ;
 43  neq = (N-1)^2 ;
 44  % questo corrisponde al ciclo
 45  % neq = 0 ;
 46  % for i=2:N
 47  %   for j=2:N
 48  %     neq     = neq + 1 ;
 49  %     EQ(i,j) = neq ;
 50  %   end
 51  % end
 52
 53  % ordinamento RED-BLACK
 54  % neq = 0 ;
 55  % for i=2:N
 56  %   for j=2:N
 57  %     if mod(i+j,2) == 0
 58  %       neq     = neq + 1 ;
 59  %       EQ(i,j) = neq ;
 60  %     end
 61  %   end
 62  % end
 63  % for i=2:N
 64  %   for j=2:N
 65  %     if mod(i+j,2) == 1
 66  %       neq     = neq + 1 ;
 67  %       EQ(i,j) = neq ;
 68  %     end
 69  %   end
 70  % end
 71
 72  % coefficienti dello stencil
 73  h2 = h^2 ;
 74
 75  % costruzione della matrice dei coefficienti
 76  A = sparse( neq, neq ) ;
 77  b = zeros ( neq, 1   ) ;
 78  for i=2:N
 79    x = (i-1)*h ;
 80    for j=2:N
 81      y = (j-1)*h ;
 82
 83      % numero equazioni vicine
 84      ic = EQ(i,j) ;
 85      il = EQ(i-1,j) ; % left
 86      ir = EQ(i+1,j) ; % right
 87      it = EQ(i,j+1) ; % top
 88      ib = EQ(i,j-1) ; % bottom
 89
 90      % termine noto
 91      b(ic) = h2*feval( f, x, y ) ;
 92
 93      % coefficienti dello stencil
 94      A(ic,ic) = -4 ;
 95      
 96      if il > 0
 97        % Punto interno aggiungo alla matrice il pezzo di stencil
 98        A(ic,il) = 1 ;
 99      else
100        % Punto di bordo aggiungo al termine noto la BC
101        b(ic) = b(ic) - feval( g, x-h, y ) ;      
102      end
103
104      if ir > 0
105        A(ic,ir) = 1 ;
106      else
107        b(ic) = b(ic) - feval( g, x+h, y ) ;      
108      end
109      
110      if it > 0
111        A(ic,it) = 1 ;
112      else
113        b(ic) = b(ic) - feval( g, x, y+h ) ;      
114      end
115
116      if ib > 0
117        A(ic,ib) = 1 ;
118      else
119        b(ic) = b(ic) - feval( g, x, y-h ) ;   
120      end
121    end
122  end
123  
124  % calcolo la soluzione
125  SOL = A\b ;
126
127  % riempimento delle matrici in uscita, inizializzazione
128  Z     = zeros(N+1,N+1) ;
129  [X,Y] = meshgrid(linspace(0,1,N+1),linspace(0,1,N+1)) ;
130  
131  % Cerco gli indici dei nodi interni, quelli numerati
132  IDX    = find( EQ > 0 ) ;
133
134  % copio la soluzione tenendo conto dell'ordine delle equazioni
135  Z(IDX) = SOL(EQ(IDX)) ;
136
137  % Cerco gli indici dei nodi di bordo
138  IDX    = find( EQ == 0 ) ;
139
140  % copio i valori al contorno sulla soluzione
141  Z(IDX) = feval( g, X(IDX), Y(IDX) ) ;
142
143end

Script per testare il codice:

File testElliptic.m

:open:
 1N = 5 ;
 2[ A, b, X1, Y1, Z1 ] = elliptic( N, @f, @g ) ;
 3N = N*2
 4[ A, b, X2, Y2, Z2 ] = elliptic( N, @f, @g ) ;
 5N = N*2
 6[ A, b, X3, Y3, Z3 ] = elliptic( N, @f, @g ) ;
 7N = N*2
 8[ A, b, X4, Y4, Z4 ] = elliptic( N, @f, @g ) ;
 9N = N*2
10[ A, b, X5, Y5, Z5 ] = elliptic( N, @f, @g ) ;
11
12EX1 = g(X1, Y1) ;
13EX2 = g(X2, Y2) ;
14EX3 = g(X3, Y3) ;
15EX4 = g(X4, Y4) ;
16EX5 = g(X5, Y5) ;
17
18E1 = max(max( abs(Z1 - EX1) )) ;
19E2 = max(max( abs(Z2 - EX2) )) ;
20E3 = max(max( abs(Z3 - EX3) )) ;
21E4 = max(max( abs(Z4 - EX4) )) ;
22E5 = max(max( abs(Z5 - EX5) )) ;
23
24disp( sprintf( 'E1 = %g', E1 ) ) ; 
25disp( sprintf( 'E2 = %g', E2 ) ) ; 
26disp( sprintf( 'E3 = %g', E3 ) ) ; 
27disp( sprintf( 'E4 = %g', E4 ) ) ; 
28disp( sprintf( 'E5 = %g', E5 ) ) ; 
29
30disp( sprintf( 'order = %g', log( E1/E2 ) / log( 2 ) ) ) ;
31disp( sprintf( 'order = %g', log( E2/E3 ) / log( 2 ) ) ) ;
32disp( sprintf( 'order = %g', log( E3/E4 ) / log( 2 ) ) ) ;
33disp( sprintf( 'order = %g', log( E4/E5 ) / log( 2 ) ) ) ;

Funzione per disegnare la soluzione:

File plotSol.m

:open:
1function plotSol(N,Z,EQ)
2  [X,Y] = meshgrid(linspace(0,1,N+1),linspace(0,1,N+1)) ;
3  surf(X,Y,Z) ;

Funzione termine sorgente:

File f.m

:open:
1function res = f(x,y)
2  res = 2*exp(x+y) ;
3  %res = 12*(x.^2+y.^2) ;
4  %res = exp(x.*y).*(x.^2+y.^2) ;
5end

Funzione valore al contorno:

File g.m

:open:
1function res = g(x,y)
2  res = exp(x+y) ;
3  %res = x.^4+y.^4 ;
4  %res = x.*y+exp(x.*y) ;
5end

Tutti i file MATLAB in un unico archivio

lez2.zip