Lezione del 15 marzo 2007

Differenze finite per BVP

Funzione di test:

File p.m

:open:
1%
2%
3function y=p(x)
4  y = (x-1).*(x+1);
5  %y = x ;
6end

File q.m

:open:
1%
2%
3function y=q(x)
4  y = -1-x.^2;
5end

File r.m

:open:
1%
2%
3function y=r(x)
4  y = -exp(-x) .* (-3 .* sin(x) + 2 .* x.^2.*cos(x) + x.^2.*sin(x));
5  %y = -exp(-x) .* (-2 * sin(x) + x .* cos(x) + x .* sin(x) + cos(x) + x .^ 2 .* cos(x));
6end

Soluzione esatta del problema:

File yesatta.m

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

Solutore:

File BVPsolve.m

:open:
 1%
 2%  Risolve il problema BVP
 3%
 4%  y'' + p(x) y' + q(x) y = r(x)
 5%
 6%  con diferenze centrate al secondo ordine.
 7%  Il sistema risultante
 8% 
 9%  +                                            + +        +   +        +
10%  |  alpha(1)  gamma(1)                        | | y(1)   |   | b(1)   |
11%  |  beta(2)   alpha(2)  gamma(2)              | |        |   |        |
12%  |              .......                       | |  ..    | = |        |
13%  |                                            | |        |   |        |
14%  |                       beta(N-1) alpha(N-1) | | y(N-1) |   | b(N-1) |
15%  +                                            + +        +   +        +
16%
17%  E' risolto con le primitive di MATLAB
18%
19function [x,y]=BVPsolve(a,ya,b,yb,N)
20  h = (b-a)/N ;
21  %for i=1:N-1
22  %  x(i) = a + i*h ;
23  %end
24  x    = [a+h:h:b-h] ;
25  % riempimento del vettore alpha
26  al   = -2 + (h^2).*q(x) ;
27  % riempimento del vettore beta
28  be   = 1 - h.*p(x)./2 ;
29  % riempimento del vettore gamma
30  ga   = 1 + h.*p(x)./2 ;
31  % riempimento del vettore r
32  bvec = h^2.*r(x) ;
33  % condizioni al contorno
34  bvec(1)   = bvec(1)   - be(1)   * ya ; 
35  bvec(N-1) = bvec(N-1) - ga(N-1) * yb ; 
36  %
37  A = sparse(diag(al(1:N-1),0)) + ... 
38      sparse(diag(be(2:N-1),-1)) + ...
39      sparse(diag(ga(1:N-2),+1)) ;
40  %A = diag(al(1:N-1),0)  + ...
41  %	  diag(be(2:N-1),-1) + ...
42  %	  diag(ga(1:N-2),+1) ;
43  % soluzione del sistema
44  y = A\bvec' ;
45  spy(A)
46end

Programma driver per verificare il codice:

File test.m

:open:
1xi = 0 ;
2xf = 10 ;
3h  = (xf-xi)/1000;
4[x,y]=BVPsolve(xi,yesatta(xi),xf,yesatta(xf),10);
5plot(x,y,'.-b','LineWidth',1.5);
6hold;
7xe = [xi:h:xf] ;
8plot(xe,yesatta(xe),'-r','LineWidth',0.5);
9legend('Numerica','Esatta');

Esempio di output

../../_images/test.jpg

Programma per controllare l’ordine delle differenze finite:

File testordine.m

:open:
 1%
 2% controllo ordine della approssimazione alle
 3% differenze centrate
 4%
 5xi = 0 ;
 6xf = 10 ;
 7[x10,y10]=BVPsolve(xi,yesatta(xi),xf,yesatta(xf),10);
 8[x20,y20]=BVPsolve(xi,yesatta(xi),xf,yesatta(xf),20);
 9[x40,y40]=BVPsolve(xi,yesatta(xi),xf,yesatta(xf),40);
10[x80,y80]=BVPsolve(xi,yesatta(xi),xf,yesatta(xf),80);
11[x160,y160]=BVPsolve(xi,yesatta(xi),xf,yesatta(xf),160);
12[x320,y320]=BVPsolve(xi,yesatta(xi),xf,yesatta(xf),320);
13[x640,y640]=BVPsolve(xi,yesatta(xi),xf,yesatta(xf),640);
14
15err10  = max(abs(y10-yesatta(x10')));
16err20  = max(abs(y20-yesatta(x20')));
17err40  = max(abs(y40-yesatta(x40')));
18err80  = max(abs(y80-yesatta(x80')));
19err160 = max(abs(y160-yesatta(x160')));
20err320 = max(abs(y320-yesatta(x320')));
21err640 = max(abs(y640-yesatta(x640')));
22
23[s,err]=sprintf('log(E10/E20)/log(2)   = %f\n',log(err10/err20)/log(2)) ;
24disp(s);
25[s,err]=sprintf('log(E20/E40)/log(2)   = %f\n',log(err20/err40)/log(2)) ;
26disp(s);
27[s,err]=sprintf('log(E40/E80)/log(2)   = %f\n',log(err40/err80)/log(2)) ;
28disp(s);
29[s,err]=sprintf('log(E80/E160)/log(2)  = %f\n',log(err80/err160)/log(2)) ;
30disp(s);
31[s,err]=sprintf('log(E160/E320)/log(2) = %f\n',log(err160/err320)/log(2)) ;
32disp(s);
33[s,err]=sprintf('log(E320/E640)/log(2) = %f\n',log(err320/err640)/log(2)) ;
34disp(s);