Lezione del 16 marzo 2007

Metodo upwind

Programma che implementa un solutore a differenze finite centrate per un boundary layer (non funziona ovviamente):

File BLsolve.m

:open:
 1%
 2%  Risolve il problema BVP
 3%
 4%  y'' + (1/epsi) y' = (1/epsi)
 5%  y(0) = y(1) = 0
 6%  
 7%  con diferenze centrate al secondo ordine.
 8%  Il sistema risultante
 9% 
10%  +                                            + +        +   +        +
11%  |  alpha(1)  gamma(1)                        | | y(1)   |   | b(1)   |
12%  |  beta(2)   alpha(2)  gamma(2)              | |        |   |        |
13%  |              .......                       | |  ..    | = |        |
14%  |                                            | |        |   |        |
15%  |                       beta(N-1) alpha(N-1) | | y(N-1) |   | b(N-1) |
16%  +                                            + +        +   +        +
17%
18%  E' risolto con le primitive di MATLAB
19%
20function [x,y]=BLsolve(epsi,N)
21  h = 1.0/N ;
22  x = [h:h:1-h] ;
23  % riempimento del vettore alpha
24  al   = -2.*ones(N-1);
25  % riempimento del vettore beta
26  be   = (1 - h/(2*epsi)).*ones(N-1) ;
27  % riempimento del vettore gamma
28  ga   = (1 + h/(2*epsi)).*ones(N-1) ;
29  % riempimento del vettore r
30  bvec = h^2.*ones(N-1)./epsi ;
31  %
32  A = sparse(diag(al(1:N-1),0)) + ... 
33      sparse(diag(be(2:N-1),-1)) + ...
34      sparse(diag(ga(1:N-2),+1)) ;
35  y = A\bvec' ;
36end

Programma driver per verificare il codice:

File test.m

:open:
1epsi=1E-6;
2[x,y]=BLsolve(epsi,100);
3hold off ;
4plot(x,y,'.-b','LineWidth',1.5);
5hold on ;
6xe = [xi:h:xf] ;
7plot(xe,yesatta(xe,epsi),'-r','LineWidth',0.5);
8legend('Numerica');
9axis([0,1,-2,1]);

Programma che implementa un solutore upwind per un boundary layer:

File BLsolve1.m

:open:
 1%
 2%  Risolve il problema BVP
 3%
 4%  y'' + (1/epsi) y' = (1/epsi)
 5%  y(0) = y(1) = 0
 6%  
 7%  con diferenze in avanti per il termine di trasporto (y').
 8%  Il sistema risultante
 9% 
10%  +                                            + +        +   +        +
11%  |  alpha(1)  gamma(1)                        | | y(1)   |   | b(1)   |
12%  |  beta(2)   alpha(2)  gamma(2)              | |        |   |        |
13%  |              .......                       | |  ..    | = |        |
14%  |                                            | |        |   |        |
15%  |                       beta(N-1) alpha(N-1) | | y(N-1) |   | b(N-1) |
16%  +                                            + +        +   +        +
17%
18%  E' risolto con le primitive di MATLAB
19%
20function [x,y]=BLsolve1(epsi,N)
21  h = 1.0/N ;
22  x = [h:h:1-h] ;
23  % riempimento del vettore alpha
24  al = -(2+h/epsi).*ones(1,N-1);
25  % riempimento del vettore beta
26  be = ones(N-1) ;
27  % riempimento del vettore gamma
28  ga = (1 + h/epsi).*ones(1,N-1) ;
29  % riempimento del vettore r
30  bvec = h^2.*ones(1,N-1)./epsi ;
31  %
32  A = sparse(diag(al(1:N-1),0)) + ... 
33      sparse(diag(be(2:N-1),-1)) + ...
34      sparse(diag(ga(1:N-2),+1)) ;
35  y = A\bvec' ;
36end

Programma driver per verificare il codice:

File test1.m

:open:
1epsi=1E-6;
2[x,y]=BLsolve1(epsi,100);
3hold off ;
4plot(x,y,'.-b','LineWidth',1.5);
5hold on ;
6xe = [xi:h:xf] ;
7plot(xe,yesatta(xe,epsi),'-r','LineWidth',0.5);
8legend('Numerica');
9axis([0,1,-2,1]);

Soluzione esatta per fare confronti:

File yesatta.m

:open:
1%
2%
3function y=yesatta(x,epsi)
4  y = x - (1-exp(-x.*(1/epsi))) ./ ...
5          (1-exp(-1/epsi)) ;
6end