4a lezione - laboratorio
Corso di Laurea Ingegneria
MECCANICA
a.a 2004-2005
1
Esercizio 1 (Esame 02/12/2002)
Si considerino i sistemi lineari Ai xi=bi,
i=1,2,3 con i vettori dei termini noti bi,
i=1,2,3 scelti in modo che la soluzione dei
sistemi sia i=[1,1,1,1]T, i=1,2,3.
Supponiamo che:
15 6 8 11
6 6 5 3
i


A1 
, A i   A1  per i  2, 3.
8 5 7 6


11 3 6 9 
2
Quesito 1
Si determini con MATLAB, il
condizionamento K2(Ai ), i=1, 2, 3 e si
verifichi che K2(Ai ) = (K2 ( A1 ) )i , i=2, 3.
Si spieghi il motivo di tale relazione e se
ne prevedano le conseguenze.
3
Quesito 2
Si costruisca un file MATLAB:
Cognome_Nome.m, che una volta
avviato:
a) faccia visualizzare una schermata con i
dati personali (Cognome, nome,
matricola, corso di laurea) ed una
breve presentazione del problema;
4
Quesito 2b
b) mediante un ciclo for, determini i dati Ai , bi
, i=1, 2, 3;
risolva quindi i sistemi Ai x i  b i
applicando il metodo di Gauss con pivoting
parziale;
calcoli l’errore relativo in norma 2;
5
Quesito 2c
c) faccia visualizzare una tabella riassuntiva che riporti:
intestazione: indice iter soluzione errore;
e su ogni riga il valore dell’indice i della
matrice, il numero di iterazioni iter effettuate nel
raffinamento, la soluzione corrispondente x i scritta
come vettore riga e l’errore relativo err;
Si utilizzino i seguenti formati di stampa:
1 cifra intera per il valore di i;
2 cifre intere per il valore di iter;
10 cifre decimali e formato virgola fissa per le
componenti di x i ;
2 cifre decimali e formato esponenziale per err.
6
Soluzione teorica del Quesito 1
Proprietà della matrice A e conseguenze.
1
A1  A1T  A11  A1T
A1
2


  A A1  
T
1
A11
K 2  A1   A1
2
2


 A1 
2

  2  A1     A1 


  A1T A11   A11
1
1 2
A

   A1   A
1
1


max  A1 
min  A1 
7
Soluzione teorica del Quesito 1
Proprietà delle matrici A i , i  2, 3
A 2   A1   A 2  A per la simmetria di A 1 .
2
A2
A
T
2
  A 2   A1  
2
2
1
2 2

 A
K 2  A2   A 2
1
2
2
2
A 
1
    A  quindi:
A
    A    A 
2
1
2 2
Analogamente per i = 3.
1
1
1
K2(A1)
1
1
8
2
Istruzioni relative al Quesito 1
clear all
disp('Numero di condizionamento delle matrici Ai')
A1=[15 6 8 11 ; 6 6 5 3 ; 8 5 7 6; 11 3 6 9];
cond_Ai=[];cond_A=[];
for i =1:3
Ai=A1^i;
cond_Ai=[cond_Ai,cond(Ai)]; % vettore dei cond(Ai)
cond_A=[cond_A,cond(A1)^i]; % vettore dei
% cond(A1^i)
end
disp('cond(Ai)')
disp(num2str(cond_Ai,'%13.3e'))
disp(’(cond(A1))^i')
disp(num2str(cond_A,'%13.3e'))
Nome file script: punto1
9
Output punto1
>> punto1
Numero di condizionamento delle
matrici Ai
cond(Ai)
6.499e+003
4.224e+007
2.746e+011
(cond(A1))^i
6.499e+003
4.224e+007
2.745e+011
>>
Conseguenze del numero
di condizionamento grande?
10
Calcolo della soluzione di
 446
199
A2  
 272

 330
199 272 330 
106 131 141 
,
131 174 199 

141 199 247 
A2 x  b 2
1247 
 577 

b2  
 776 


 917 
>> x1=A2\b2 % Operatore \
Numero di cifre significative perse
x1 =
>> nc=log10(K2)
1.00000000154209
nc =
7.6257
0.99999999870797
1.00000000095304
0.99999999790943
>> err1=norm(x1-alfa)/norm(alfa)
9 cifre signific.,
err1 = 1.53e-009=.153e-008

11
8 decimali corretti
Come incide il valore di K2(A2)?
>> A2m=A2;
>> A2m(2,2)=A2(2,2)+1e-3; % perturbazione data
% sulla matrice
>>pert=norm(A2-A2m)/norm(A2)
pert =
1.0546e-006 % entità della perturbazione
>> x1m=A2m\b2 % soluzione perturbata
x1m =
2.06050813486122
0.11143302854021
1.65544907526039
-0.43770893726797
>> err1m=norm(x1m-alfa)/norm(alfa)
12
err1m = 1.0501 >100% !!!
Istruzioni relative al Quesito 2a
clear all
disp('Cognome e nome studente: XXXX XXX')
disp('Numero di matricola: XXXX')
disp(‘Corso di Laurea: XXXX’)
disp(' ')
disp('Questo programma calcola e
visualizza la soluzione dei ')
disp('sistemi lineari A_i x_i=b_i,
i =1,2,3, con i vettori b_i tali che sia')
disp('alfa=[1,1,1,1]'',essendo: ')
A1=[15 6 8 11;6 6 5 3;8 5 7 6;11 3 6 9];
disp('A1=');disp(A1)
13
disp( 'e A_i= A1^i per i=2,3.')
Output file Cognome_Nome.m
>> Cognome_Nome.m
Cognome e nome studente: XXXX XXX
Numero di matricola: XXXX
Corso di Laurea:XXXX
Questo programma calcola e visualizza
la soluzione dei sistemi lineari
A_i x_i=b_i, i =1,2,3, con i vettori
b_i tali che sia alfa=[1 1 1 1]',
essendo:
A1= 15
6
8
11
6
6
5
3
8
5
7
6
11
3
6
9
14
e A_i= A1^i per i=2,3.
Istruzioni di gausspv_r.m
Istruzioni relative al Quesito 2b
tab=[];
alfa=ones(4,1); % soluzione
for i =1:3
A=A1^i;
b=A*alfa; % vettore termini noti
[L,U,P] = lu(A);
%[L,U]=lu(A);
y = forwsub(L,P*b); % y = L\(P*b);
x = backsub(U,y);
% x=U\y;
R=eye(size(A));
[x,iter]=raff_iter(A,b,L,U,P,R,x);
err=norm(alfa-x)/norm(alfa);
tab=[tab;[i,iter,x',err]];
end
15
Risultati File Punto 2b
Tabella Punto 2c
fprintf('i iter \t\t\t soluzione
\t\t \t\t
errore \n')
fprintf('%1d %2d %14.10f %14.10f %14.10f
%14.10f %12.2e \n',tab');
i iter
soluzione
errore
1
0
1.0000000000 1.0000000000 1.0000000000 1.0000000000
2.42e014
2
1
1.0000000016 0.9999999987 1.0000000010 0.9999999978
1.59e009
3
5
0.9999992309 1.0000006444 0.9999995247 1.0000010427
7.62e007
16
Esercizio 2
Sia dato il sistema lineare avente la matrice
dei coefficienti ed il vettore dei termini noti
così assegnati:
4 1 0 0 0 0
 3
1 5 2 0 0 0 
 2 


 
0 2 6 3 0 0
 1
A
 b 
 0 0 3 5 1 0 
 1
 0 0 0 1 3 1 
 3


 
0 0 0 0 1 2
 4
17
Problema I°: quesiti 1, 2
1. Si studino le caratteristiche della matrice A
e si stabilisca quale tipo di fattorizzazione è
possibile effettuare.
2. Si risolva il sistema usando ogni tipo di
fattorizzazione possibile e con l’operatore \
applicato alla matrice A ed al vettore b.
18
Caratteristiche di A: conseguenze
Caratteristiche di A:
a) simmetrica
b) tridiagonale;
c) diagonalmente dominante;
d) definita positiva perché è a), c) e tutti gli
elementi in diagonale principale sono >0.
Conseguenze:
a) È possibile la fattorizzazione A=L*U;
b) È possibile la fattorizzazione di Cholesky
A=R*R’
19
Verifica MATLAB e soluz. sistema
I0=[4 5 6 5 3 2];I1=[1 2 3 -1 1];[n,m]=size(A);
A=diag(I0)+diag(I1,-1)+diag(I1,1); % matrice A
b=[3 –2 1 –1 3 4]; % termine noto
A==A’; % controllo simmetria
for i=1:size(A,1)
minore(i)=det(A(1:i,1:i));
end
% minori principali di testa
disp(minore)
[L U P]=lu(A) % fattorizzazione LU
Rt=chol(A) % fattorizzazione di Cholesky
y=Rt’\b;alpha=Rt\y; %soluzione sistema utilizzando Chol
y=L\b;alpha1=U\y; %soluzione sistema utilizzando A=Lu
alpha2=A\b
% soluzione sistema (\ utilizza A=LU)
20
Problema II°: quesiti 1 e 2
1. Si studi la convergenza dei metodi di
Jacobi, Gauss-Seidel e Rilassamento in
serie (SOR) per il sistema assegnato.
2. Si dica quale di questi metodi è il più
veloce giustificando teoricamente la
risposta e calcolando, mediante Matlab, i
raggi spettrali delle rispettive matrici di
iterazione ed il valore ottimale del
parametro  per il metodo SOR.
21
Quesito 3a
3. Si costruisca un file MATLAB che:
a) calcoli la soluzione numerica del
problema assegnato applicando il metodo
con convergenza migliore e fissando una
precisione non inferiore a 1.e-4 ed un
vettore di innesco pari a x0=[1 1 1 1 1 1]T;
22
Quesito 3b
b) faccia visualizzare una tabella riassuntiva che riporti:
intestazione: iterazioni soluzione errore;
e su ogni riga il numero di iterazione iter , la
soluzione corrispondente x i e la norma del residuo
err;
Si utilizzino i seguenti formati di stampa:
2 cifre intere per il valore di iter;
5 cifre decimali e formato virgola fissa per la
soluzione approssimata x i ;
1 cifra decimale e formato esponenziale per err.
23
Soluzione quesito 1:Convergenza dei
metodi
Caratteristiche di A:
a) diagonalmente dominante  Jacobi conv.
b) tridiagonale  anche Gauss-Seidel
converge, inoltre:
  BGS     BJ   R  BGS   2R  BJ 
2
c) simmetrica
d) definita positiva  anche SOR converge
per
 0 <  < 2;
24
Qual è il metodo più veloce?
Il metodo più veloce è SOR se si assume:
   ott 
2
1  1    BGS 
.
25
Istruzioni relative ai quesiti 1, 2
% file Punto2
I0=[4 5 6 5 3 2];I1=[1 2 3 -1 1];[n,m]=size(A);
A=diag(I0)+diag(I1,-1)+diag(I1,1)
D=diag(diag(A));
B_J=eye(n)-inv(D)*A; % metodo di Jacobi
rho_J=max(abs(eig(B_J)))
R_J=-log(rho_J)
omega=1; % metodo di Gauss-Seidel
OE=-omega*tril(A,-1);
B_GS=eye(n)-omega*inv(D-OE)*A;
rho_GS=max(abs(eig(B_GS)))
R_GS=-log(rho_GS)
omega_ott=2/(1+sqrt(1-rho_GS)) % metodo SOR
OE=-omega_ott*tril(A,-1);
B_r=eye(n)-omega_ott*inv(D-OE)*A;
26
rho_r=max(abs(eig(B_r))) R_r=-log(rho_r)
Risultati file Punto2
rho_J
=
0.7196
R_J
=
0.3290
rho_GS =
0.5179
R_GS =
0.6580
omega_ott =
rho_r
=
1.1804
0.1804
R_r
=
1.7126
27
Istruzioni relative ai quesiti 3a, 3b
b=[3 -2 1 -1 3 4]';
K=cond(A,inf);
precisione=input('precisione = '); % 1.e-4
toll=precisione/K
% toll = 9.3310e-006
x0=ones(n,1);omega=omega_ott;
[x,iter,err,rho]=Gauss_Seidel_ril(A,b,x0,omega,nmax,toll);
it=[1:iter]';
tab=[it x err];
s='-------------------------------------------------';
disp(s)
fprintf('iter
soluzione
errore\n')
disp(s)
28
fprintf('%2d %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %9.1e\n‘…
,tab');
Risultati quesiti 3a, 3b
iter
soluzione
errore
1
0.40980 -1.22145 -0.09326 -0.11434
0.56155
1.84897
1.1e+000
2
1.17181 -0.48443
0.47164 -0.41692
0.18755
1.91655
4.2e-001
3
0.81686 -0.80030
0.67260 -0.59295
0.15916
1.92112
1.3e-001
4
0.97411 -0.87533
0.76977 -0.63672
0.14526
1.92850
3.4e-002
5
0.96788 -0.90620
0.79022 -0.64659
0.14099
1.92969
8.7e-003
6
0.97811 -0.91270
0.79491 -0.64914
0.14028
1.92989
1.1e-003
7
0.97819 -0.91376
0.79599 -0.64961
0.14015
1.92993
3.4e-004
8
0.97849 -0.91415
0.79622 -0.64972
0.14011
1.92995
5.1e-005
9
0.97855 -0.91421
0.79627 -0.64974
0.14010
1.92995
1.3e-005
10
0.97855 -0.91422
0.79628 -0.64975
0.14010
1.92995
2.5e-006
Si invitano gli studenti a determinare la soluzione con
precisione = 1.e-8 e riportare le componenti della
29
soluzione approssimata con 10 cifre e virgola fissa.