LU-faktorisering

Syfte: Matematik är ett av de verktyg som används inom numerisk analys. Detta innebär dock inte att man alltid bör använda matematiska  formler för praktisk problemlösning, ofta krävs en omformulering. Vi såg detta redan i den första laborationen. Här följer några fler exempel.

När man gör teori för linjära ekvationssystem är det bekvämt med matrisinverser; detta gäller även om man är numeriker. Detta faktum medför dock inte att matrisinverser är bra för praktiska beräkningar.

Kör kodavsnittet nedan. Klipp och klista. Kommentarer? Förklara tidsskillnaden! Spara körtiderna. 

n = 3000;
A = randn(n);
b = randn(n, 1);
tic; x = inv(A) * b; toc
tic; y = A \ b; toc

Det är vanligt att man arbetar med stora matriser som innehåller många nollor, så kallade glesa matriser. Glesa matriser uppkommer t.ex. när man löser olika typer av differentialekvationsproblem. Som Du strax kommer att se är det alldeles speciellt ineffektivt att arbeta med inverser av glesa matriser.

I följande övning kommer Du att använda kommandot spdiags för att skapa en så  kallad tridiagonal matris, dvs. en matris som har nollor överallt utom på (huvud)diagonalen och den första diagonalen under och över huvud-diagonalen. Man kan använda kommandot spy för att plotta matrisens gleshetsstruktur, dvs. man ritar en punkt för varje ickenolla, nollor ritas ej. Observera att Matlab använder ett annat lagringsformat för glesa matriser än för fulla matriser. När matrisen är gles så lagras endast ickenollorna och i Matlab (det finns andra datastrukturer) lagras väsentligen alla tripplar, (j, k, ajk), där ajk är skilt från noll.

Matrisen A är positivt definit. Studera gleshetsstrukturerna (help spy) hos A, chol(A) respektive inv(A). Hur många element krävs för att lagra inversen respektive Choleskyfaktorn? Titta i spy-graferna eller använd nnz. Spara körtider och uppgifter om minnesbehoven. Chol-kommandot beräknar den sk Choleskyfaktoriseringen av en matris, ett specialfall av LU-faktoriseringen, som i sin tur är en variant av Gausselimination. 

n = 6000;
e = ones(n, 1);
A = spdiags([e 2*e e], -1:1, n, n);
spy(A)
b = randn(n, 1);
tic; x = inv(A) * b; toc
tic; y = A \ b; toc

Du utför ett experiment som ger dig elementen i en matris, A. Matrisen används för att lösa ett Ax=b-problem. Högerledet b kommer från teori och utgörs inte av mätdata. Fär att få lite känsla för osäkerheten i lösningen x utför du fyra mätningar och får fyra matriser A1, A2, A3 och A4. Matriserna ger dig fyra lösningsvektorer x1, x2, x3 och x4. Matriser och högerled finns på följande fil, mat_data.mat (placera muspekaren på filen och håll nere höger musknapp, välj meny-alternativet Save Link As...). Ge Matlabkommandot load mat_data för att läsa in filen.

Beräkna och jämför de fyra lösningarna. Kommentarer?
Hur stämmer skillnaderna mellan lösningarna med teorin?
Visa hur du kan uppskatta osäkerheten i x om du bara hade utfört en mätning (A1 säg) och dessutom känner till att mätdata (matrisen) har ett maximalt relativt fel om 5e-5?