Effetto Pivoting

Instabilità della fattorizzazione LU

P. Gervasio - Esercitazioni di Calcolo Numerico - A.A. 2003/2004 Partendo da (1) e definendo l'effetto benefico del pivoting (k): max |a_ij^(k)| = ρ_n max |a_ij| Problema: determinare ε e δ tali che: |L^(k) U^(k) - A|_ij ≤ n kAk O(kδAk u_∞ n_∞ Invece di risolvere il sistema Ax = b, si risolve il sistema: (A + δA)x̂ = b Attenzione: quando 1ρ · ≤ 1, si risolve il sistema un(A + δA)x̂ = b Se non si effettua pivoting o si effettua pivoting parziale, si ha: |L̂^(k) Û^(k) - A|_ij ≤ ρ_n²(3|A| + 5|δA|) O(|δA| ≤ nu_∞) Se si effettua pivoting totale si ha: dimensione del sistema n = 1/2 1/2 1/3 1/(n-1) 1/2(2 3 4 )≤ n · · . . . nρ n₁ = 1.1102e 16 in unità di roundoff (ε ≈ −u/2) MATLAB) N.B.In pratica si osserva che per la pivotazione= matrice con tutti gli elementi presi in|A| A 10.totale si ha <ρ nvalore assoluto. 1 2Risolvo con LU senza pivoting(k) n−1max = max = 2⇒ |a | |u |i,j iji,j,k ijNella maggior parte dei casi, quando si usa pi- max = 1|a |i,j ijn−12voting per righe si ha ρ <<nIl pivoting totale garantisce che resti limi-ρ n n−1= 2⇒ ρ n1.tato al crescere di e · <<n ρ un Guardiamo quando 1.ρ ≥unEsempio 1 1 1−tn−1 n−1= 2= 2 ε β· ·ρ u1 se = o =i j, j n n 2 2= seA −1 i > jij In MATLAB: = 2 e = 54 e risultaβ t0 altrove n−55= 2 .ρ un:Consideriamo come termine noto il vettore b 1 per 55.Segue che ≥ n ≥ρ un= 3 per = 1, La soluzione esatta− i, i ..., n.b i = è = [1; 1; 1; 1].del sistema ....;A x b x Per 54 il sistema dovrebbe essere risolton ≤senza problemi, per 55 la risoluzione deln

<p>≥sistema dovrebbe risultare inesatta.3 4Per svolgere la fattorizzazione LU senza pivoting, nonsi può usare il comando di Matlab, ma bisogna usarelu Proviamo a risolvere il sistema con = 55 con la fat-nuna function scritta appositamente: [L,U]=lufact(A). (usando la functiontorizzazione LU con pivot per righe[L,U,P]=lufactp(A)):Per la risoluzione del sistema triangolare inferiore:function [y]=risfor(L,b) [L,U,P]=lufactp(A);per la risoluzione del sistema triangolare superiore: y=risfor(L,P*b);function [x]=risback(U,y). x=risback(U,y);(Le function e si trovanolufact.m, risback.m risfor.m err=norm(x-xex,inf)/norm(xex,inf)in M:\matlab\calcolo). Ancora si ha un errore del 100% eclearn=input(’Dimensione matrice\n’) >> rho=max(max(U))/max(max(A))A=speye(n); A(:,n)=ones(n,1); rho =for i=2:n; A(i,1:i-1)=-1; end (1,1) 1.8014e+16xex=ones(n,1);b=(2:-1:3-n)’;b(n)=2-n; Per finire utilizziamo la pivotazione totale: (usando lafunction</p>

[L,U,P,Q]=lufactpt(A)):[L,U]=lufact(A);y=risfor(L,b); [L,U,P,Q]=lufactpt(A);x=risback(U,y); y=risfor(L,P*b);err=norm(x-xex,inf)/norm(xex,inf) x=risback(U,y);x=Q*x;err=norm(x-xex,inf)/norm(xex,inf)

Per = 54, si ha la soluzione esatta, mentre per = 55n n= [ ; ; ; ; ; ; ; ]si ottiene come soluzione il vettore: . . .x 1 1 1 1 1 0 1 In questo caso la risoluzione numerica è esatta: l’erroreè dello 0% econ un 100%.errore relativo sulla soluzione esatta del

Calcoliamo il fattore di crescita: >> rho=max(max(U))/max(max(A))rho =>> rho=max(max(U))/max(max(A)) (1,1) 2rho =(1,1) 1.8014e+16 5 6