Ytterligare ett problematiskt exempel och tillväxtfaktorer i LU-faktorisering 1 Introduktion

På föreläsning IV LU-faktoriserade vi matrisen A = 1 2 53 +2 1 1 1 .

Vi fick att L = 2 53 1 + 2 0 1 , U = 1 2 53 +2 0 − (2 1 53 + 1) .

Detta fick förödande konsekvenser då på grund av avrundningsfel element U inte kunde representeras i datorn och istället lagrades som åtgärda detta problem introducerade vi LU-faktorisering med pivotering. Ett annat problem som också uppstod är att både element (2 , − 1) (2 i 53 L +2) . För att och element (2 , 2) i U är mycket stora relativt elementen i (2 , 2) A . Till exempel vid lösning av i systemet U y = b uppstår division med element genom att pivotera raderna i A (2 , 2) i U flera decimalers nogranhet i vår lösning. Även detta problem tillintetgjordes . Betrakta nu matrisen detta gör att vi förlorar A =  1      − 1 − − 1 1 − 1 0 1 − 1 − 1 − 1 En LU-faktorisering med pivotering ger 0 0 1 − 1 − 1 0 0 0 1 1  1 1 1      − 1 1 U =  1 0      0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1  2 4 8 16      (1) (2) Testa gärna i Matlab. Om vi låter struktur d.v.s., − 1 A under diagonalen, vara en 1 m × m matris med samma generella på diagonalen, 1 ’or i sista kollumnen och 0 ’or i övrigt, kommer element ( m, m ) katastrofalt vid lösning av systemet i U y U vara = b 2 m − 1 oavsett pivotering. Detta är då division med 2 m − 1 motsvarar en förlust av m − 1 bitars nogganhet. I praktiska beräkningar är det inte ovanligt med mycket stora koefficientmatriser A .

1

Vi befinner oss nu i en obekväm situation, vi vill använda en algoritm för att lösa lineära system för vilken det finns probleminstanser som ger oanvändbara resultat. Vi kan dock tillåta oss att vara optimistiska, exemplet ovan är artifi ciellt konstruerat och kommer inte från någon tillämpning. Om detta exempel har sagts (L.N.Trefethen and D.Bau III, Numerical Linear Algebra, pp 166): Despite examples like (1) Gaussian elimination with partial pivoting (LU faktorisering med pivotering) is utterly stable in practice. Large factors U like (2) never seem to appear in real applications. In fifty years of computing, no matrix problems that excite an explosive instability are known to have arisen under natural circumstances.

2 Vår uppgift

Vi vill undersöka om att det är säkert att använda LU-faktorisering med piv otering för att lösa lineära system. Först definierar vi ett mått på hur mycket större U är än A genom ρ = max i,j | U i,j | max i,j | A i,j | d.v.s., det till belopp största elementet i U dividerat med det till beloppet största elementet i A . Vi kallar ρ för tillväxtfaktor. I exemplet ovan är med pivotering och räkna ut tillväxtfaktorn ρ ρ = 16 . Nu är vi redo att börja på vår undersökning. Uppgiften är att konstruera, LU-faktorisera för en miljon slumpmässiga m × m matriser, ρ för varje matris skall sparas i en vektor och illustreras i ett histogram.

Detta skall göras för m = 8 , 16 och 32 . Instruktionerna för fallet m = 8 kan skrivas i Matlab liknande kod som rho_8 = zeros(10^6,1); for n = 1:10^6 A = randn(8,8); [L,U,P] = lu(A); rho_8(n) = max(max(abs(U))/max(max(abs(A)))); end figure(1) centers = min(rho_8):0.2:max(rho_8); hist(rho_8,centers) I kodsnutten ovan har kommentarer utelämnats använd Matlabs hjälp funktion om det behövs för att kommentera på vad koden utför. Gör experimentet, redovisa histogrammen och kommentera följande i en rapport: • Tänk på att teoretiskt så skulle vi kunna få ett 32769 , 2 32 − 1 = 2147483648 ρ , vad får vi i praktiken?

= 2 8 − 1 = 128 , 2 16 − 1 = • Vad säger histogrammen? Vilken tillväxtfaktor är vanligast för För m = 16 ? För m = 32 ?

m = 8 ?

• Statistiskt sett, är det säkert att använda LU-faktorisering med pivotering för att att lösa lineära system?

2