Esercitazioni di MATLAB: Lezione del 4 Marzo 2008¶
Esempio di programmazione
:open:
1%
2% Risolve il sistema lineare A*x = b
3%
4function x = linearsistemsolve( A, b )
5
6 [nr,nc] = size(A) ;
7
8 % controllo dimensioni
9 if nr ~= nc
10 error( sprintf( 'la matrice A (%d x %d) non e` quadrata ', nr,nc) ) ;
11 end
12
13 neq = length(b) ;
14
15 if nr ~= neq
16 error( sprintf( 'la matrice A (%d x %d) non e` compatibile\ncon la lunghezza del vettore b, length(b) = %d ', nr, nc, neq ) ) ;
17 end
18
19 M = [ A b ] ; % costruisco la matrice aumentata
20 for i=1:neq-1
21 % cerco l'elemento di massimo modulo colonna i-esima
22 [V,k] = max(abs(M(i:end,i))) ;
23
24 % rimetto l'indice k alla giusta posizione
25 k = k + i-1 ;
26
27 % scambio la riga i con la riga k
28 BF = M(i,:) ;
29 M(i,:) = M(k,:) ;
30 M(k,:) = BF ;
31
32 % controllo singolaritÃ
33 if V == 0
34 error( 'Matrice singolare' ) ;
35 end
36
37 % costruisco la matrice di Frobenius
38 L = eye(neq) ;
39 L(i+1:end,i) = -M(i+1:end,i)./M(i,i) ;
40
41 % azzero gli elementi della colonna i
42 M = L * M ;
43
44 end
45
46 % alloco il vettore soluzione
47 x = zeros(neq,1) ;
48
49 % a questo punto la matrice è in forma triangolare
50 % applichiamo le operazioni di ritorno
51 x(neq) = M(neq,neq) / M(neq,neq+1) ;
52 for i=neq-1:-1:1
53 x(i) = ( M(i,neq+1) - M(i,i+1:end-1) * x(i+1:end) ) / M(i,i) ;
54 end