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 matrisen 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 alla tripplar, (j, k, a_jk), där a_jk ä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 = 3000;
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

Här följer en övning på illakonditionerade matriser, konditionstal och residualvektorer.

Den sk Hilbertmatrisen (help hilb) är mycket illakonditionerad för stora dimensioner. Vi vill lösa A x = b, där A är en Hilbertmatris av ordning 10. För att skaffa oss ett facit, fuskar vi lite och väljer x som består av idel ettor och genererar så b (känner vi x, så är ju b = A x). Vi skall nu räkna ut x och jämföra med den exakta lösningen (som vi ju känner).

Lös A x = b på vanligt sätt och jämför med den exakta lösningen. Kommentarer? Hur stämmer det med teorin?

Räkna nu ut residualvektorn, r = A x - b (x är alltså den beräknade lösningen och inte facit). Kommentarer? Hur stämmer det med teorin (i boken)?