1

Iterative Losung großer Gleichungssysteme,Ubung 2, Aufgabe 1Lukas Kogler

Ein MATLAB-File mit der Implementation des Algorithmus und einem performance-Test im Vergleich zu fft bzw LU-Zerlegung befindet sich im Anhang.

Fur die Implementation der dst habe ich jedoch einen etwas anderen Weggewahlt, als er im Skript angegeben ist. Anstatt fft mit dem Vektor (0, v, 0, . . . , 0)der Lange 2n+ 1 zu machen, der einer Funktion auf [−1, 1], welche auf [−1, 0] 0ist entspricht, wende ich fft auf den Vektor (−vn,−vn−1, . . . ,−v1, 0, v1, v2, . . . , vn)an, der einer ungeraden Funktion entspricht. Dadurch erhalt man auf naturlicheWeise eine Darstellung durch reine sin-Terme und man kann sich theoretischsparen, am Ende den Imaginarteil nehmen zu mussen. Durch Rundungsfehlermusste ich (um die MATLAB-Warnung loszuwerden) dies allerdings trotzdemmachen.

C:\Users\wauzinger\Documents\a...\ue2_1.m Page 1

% --- first, verify that all produce the right result --- n = 10;%[xi,A,D,L] = make_pois_sparse(n);[xi,A,D,L] = make_pois_sparse(n);b = 8*ones(n,1); % solve -u''=8 -> u = 4*x*(1-x)xi = [0;xi;1]; uref = [4*xi.*(1-xi)];u_my_dst = [0;solve_my_dst(b,D);0];u_dst = [0;solve_dst(b,D);0];u_chol = [0;solve_chol(b,L);0];u_chol2 = [0;solve_chol2(b,A);0]; subplot(2,1,1);hold offplot(xi,uref, 'r');hold on;plot(xi, [0;A\b;0], 'c--');plot(xi, u_chol, 'm');plot(xi, u_dst, 'b');plot(xi, u_my_dst, 'g--');hold off;title('4*x*(1-x)');legend('ref','backslash','chol', 'dst','my-dst'); % --- performance tests ---%ns = []ns = [10:10:90,100:100:900,1000:1000:9000,10000:10000:90000,100000:100000:1000000];t = zeros(size(ns,2),4); nreps = 10;for k = [1:size(ns,2)] n = ns(k); [xi,A,D,L] = make_pois_sparse(n); b = ones(n,1); tic; for j=[1:nreps] y = solve_dst(b,D); end t(k,1) = toc; tic; for j=[1:nreps] y = solve_my_dst(b,D);

C:\Users\wauzinger\Documents\a...\ue2_1.m Page 2

end t(k,2) = toc; tic; for j=[1:nreps] y = solve_chol(b,L); end t(k,3) = toc; tic; for j=[1:nreps] y = solve_chol2(b,A); end t(k,4) = toc;endt = t./nreps;subplot(2,1,2)xi2 = 1:size(ns,1);cols = ['r','g','b','c'];hold offplot(ns,t(:,1), cols(1));hold onfor k=[2:4] plot(ns,t(:,k),cols(k));endxlabel('n')ylabel('t')title('timing');lgd = legend('dst','my-dst','chol', 'chol+factor');lgd.Location='northwest';hold off % --- functions --- function [xi,A,D,L] = make_pois_mat(n)h = 1/(1+n);xi = (h:h:1-h)';A = 1/(h^2)*(2*eye(n)-diag(ones(n-1,1),1)-diag(ones(n-1,1),-1));D = diag(4*sin(pi*xi/2).^2)/(h^2);L = chol(sparse(A));end function [xi,A,D,L] = make_pois_sparse(n)h = 1/(1+n);

C:\Users\wauzinger\Documents\a...\ue2_1.m Page 3

xi = (h:h:1-h)';row = [2,-1,zeros(1,n-2)];e = 1/(h^2)*ones(n,1);A = spdiags([-e 2*e -e], -1:1, n, n);D = spdiags((4*sin(pi*xi/2).^2)/(h^2),0,n,n);L = chol(A);end function x = solve_dst(b,D)x = dst(b); % sine-transformx = D\x; % eigenvaluesx = idst(x); % inv. sin-transend function x = solve_my_dst(b,D)y2 = [0;b;0;-b(end:-1:1)];y2 = fft(y2);y3 = D\y2(2:size(b,1)+1);y2 = ifft([0;y3;0;-y3(end:-1:1)]);x = real(y2(2:size(b,1)+1)); %only necessary bc of rounding errors??%x = y2(2:size(b,1)+1);end function x = solve_chol(b,L)x = L'\b; %L^T y1 = bx = L\x; %L x = y1end function x = solve_chol2(b,A)L = chol(A); % L^T * L = Ax = L'\b; %L^T y1 = bx = L\x; %L x = y1end

2. Beispiel

a)

Geg.: b ∈ Rn, U ⊂ Rn m-dimensionaler UnterraumZ.z.: u = argminv∈U‖v − b‖2 ⇐⇒ u− b ⊥ UBeweis. Betrachte den m-dimensionalen Unterraum U als Bild einer Matrix A ∈ Rn×m. Es gilt mitu = Ax

Ax− b ⊥ Bild(A)

⇐⇒Ax− b ∈ ker(AT )

⇐⇒AT (AX − b) = 0

⇐⇒ATAx = AT b

Die Losung x der Normalengleichung minimiert aber genau ||Ax − b||2 (unter anderem SkriptumKapitel 7.2) also

argminv∈Bild(A)||v − b|| = Ax = u

Uu

b

b)

Geg.: Orthonormalbasis {q1...qm} von U , Q := (q1, ..., qm) ∈ Rn×m

� Z.z.: u = QQT bSei v ∈ U fest aber beliebig, dann gilt:

〈QQT b− b, v〉 = 〈(QQT − I)b, v〉 = 〈b, (QQT − I)v〉

Stelle v = Qw als Linearkombination in seiner Basis dar mit, so gilt weiters:

〈b, (QQT − I)v〉 = 〈b, (QQT − I)Qw〉= 〈b,QQTQw −Qw〉= 〈b,Qw −Qw〉 = 0

Gemeinsam mit a) folgt die Behauptung.

� Z.z.: P := QQT =⇒ P = PT = P 2

PT = (QQT )T = (QT )TQT = QQT = P

P 2 = (QQT )(QQT ) = QIQT = QQT = P

� Ges.: Koeffizienten von u = QQTx in der Form u =m∑i=1

ξiqi

QQT (x1...xm)T = Q(〈q1, x〉, ..., 〈qm, x〉)T =

m∑i=1

〈qi, x〉qi

1

c)

Geg.: {u1...um} beliebige Basis von UGes.: Orthogonalprojektor P auf U mittels einer Faktorisierung der Matrix U = (u1...um).P := QQT leistet gemaß b) das Gewunschte, wobei Q jene orthogonale Matrix aus der reduziertenQR-Zerlegung von U ist, also U = QR mit Q ∈ Rn×m und R ∈ Rm×m obere Dreiecksmatrix.

d)

Geg.: A ∈ Rn×m mit m < n und rgA = mGes.: Losung des linearen Ausgleichsproblems ‖Ax−b‖2 −→ minimal unter Verwendung der reduziertenQR-Zerlegung von A.Das lineare Ausgleichsproblem ist aquivalent zur Gauß’schen Normalengleichung ATAx = AT b. Einsetzender reduzierten QR-Zerlegung liefert:

(QR)TQRx = (QR)T b⇐⇒ RTQTQRx = RTQT b

⇐⇒ RTRx = RTQT b

⇐⇒ Rx = QT b

Es genugt also, das kleinere Gleichungssystem Rx = QT b zu losen.

e)

Im Folgenden fuhren wir einen Pseudocode fur die Aufgabenstellung an.

Algorithm 1 Reduced QR-Decomposition

1: %Erstelle 8x5 Zufallsmatrix:2: A = rand(8,5);3:

4: %Zufallsvektor:5: b = rand(8,1);6:

7: % Uberprufe, ob Matrix vollen Rang besitzt8: while rank(A) 6= 5 do9: A = rand(8,5);

10: end while11:

12: % Reduzierte QR-Zerlegung13: [Q,R] = qr(A,0);14:

15: %Lose Rx = QT b16: c = transpose(Q)*b;17: x = R\c;18:

19: %Vergleiche Ergebnisse20: display(’QR-Verfahren:’);21: display(x);22: display(’Matlab-Implementierung:’);23: display(A\b);

2

3. Beispiel

Geg.: A ∈ GLn(R), b ∈ Rn, U ein m-dimensionaler Unterraum des Rn

Ges.: Approximation fur die Losung x∗ des Ausgleichsproblems Ax − b in dem Unterraum U gemaßGalerkin-Approximation.

a)

Geg.: {u1...um} Basis von U mit Ansatz u = Uξ =m∑i=1

ξiui

Ges.: Lineares Gleichungssystem, dessen Losung die Koeffizienten ξi liefert.Es soll also fur alle j = 1...m gelten:

〈b−m∑i=1

ξiAui, uj〉 = 0⇐⇒ 〈m∑i=1

ξiAui, uj〉 = 〈b, uj〉

⇐⇒m∑i=1

ξi〈Aui, uj〉 = 〈b, uj〉

In Matrixschreibweise liefert das: 〈Au1, u1〉 · · · 〈Aum, u1〉...

...〈Au1, um〉 · · · 〈Aum, um〉

ξ1

...ξm

=

〈b, u1〉...〈b, um〉

⇐⇒ UTAUξ = UT b

wobei G := UTAU die Gramm-Matrix ist.Ist das Gleichungssystem eindeutig losbar? Nein, da zum Beispiel das Gleichungssystem(

1 0)( 0 1−1 0

)(10

)!=(1 0

)(10

)⇐⇒ 0 = 1

nicht einmal eine Losung besitzt.

b)

Sei im folgenden A symmetrisch und positiv definit.Z.z.: A SPD =⇒ G := UTAU SPDAufgrund der positiven Definitheit von A gilt:

〈Gx, x〉 = 〈UTAUx, x〉 = 〈A Ux︸︷︷︸y

, Ux︸︷︷︸y

〉 > 0

womit die Definitheit folgt. Die Symmetrie von G ist klar.Ges.: Die zu u gehorige Minimierungsaufgabe.Wir rufen folgende Resultate uber lineare Ausgleichsprobleme aus der Numerik in Erinnerung: Ein Vektorx∗ ∈ Rn lost genau dann das lineare Ausgleichsproblem

‖Ax− b‖22 = min!

wenn x∗ die Gauß’sche Normalengleichung

A∗Ax = A∗b

erfullt. Es existiert immer mindestens eine Losung des linearen Ausgleichsproblems, und die Menge allerLosungen ist gegeben durch

x∗ + ker(A).

1

Im Folgenden werden wir obiges Gleichungssystem in den Kontext der Normalengleichung setzen. Da Aund damit auch insbesondere A

12 symmetrisch sind, folgt

UTAUx = UT b⇐⇒ UT (A12 )TA

12︸ ︷︷ ︸

A

Ux = UT A12A− 1

2︸ ︷︷ ︸I

b

⇐⇒ ‖A 12Ux−A− 1

2 b‖22 = min!

⇐⇒ ‖A 12 (Ux−A−1b)︸ ︷︷ ︸

=x∗

‖22 = min!

⇐⇒ ‖Ux− x∗‖2A = min!

Uber diese Beziehung kann x als orthogonale Projektion von x∗ aufgefasst werden.

c)

Im Folgenden fuhren wir einen Pseudocode fur die Aufgabenstellung an.

Algorithm 1 Galerkin Approximation

1: %Symmetrische 8x8 Integer-Zufallsmatrix2: A = randi(8,8);3: A = A*transpose(A);4:

5: %Losungsvektor:6: x = transpose([1 1 1 0 0 0 0 0]);7:

8: %Rechte Seite9: b = A*x;

10:

11: %Basis des Unterraums: R2-Ebene12: U = transpose([1 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0]);13:

14: %Reduziertes Gleichungssystem15: G = transpose(U)*A*U;16: c = transpose(U)*b;17:

18: %Koeffizienten:19: Xi = G\c;20:

21: %Linearkombination der Losung22: u = Xi(1)*U(:,1) + Xi(2)*U(:,2);23:

24: %Matlab-Implementiertung25: display(’Galerkin-Approximation:’);26: display(u);

2

Iteratives Losen linearer Gleichungssysteme - Ubung 2

Michael Neunteufel

11. April 2017

5. Beispiel

Geg.: Seien A und N symmetrisch und positiv definit mit der zugehorigen Iteration:

xk+1 = xk +N(b−Axk)

und der Iterationsmatrix M = I −W−1A mit W = N−1. Ferner gelte:

2W > A > 0

a)

Z.z.: Die adjungierte Matrix MA von M in Bezug auf das Skalarprodukt 〈., .〉A erfullt:

MA = A−1MTA

Beweis: Sei x beliebig aber fest gewahlt, dann gilt:

〈Mx, x〉A = 〈x,MTATx〉= 〈x,AA−1MTAx〉= 〈Ax,A−1MTAx〉= 〈x,A−1MTAx〉A = 〈x,MAx〉A

b)

Z.z.: ρ(M) = ‖M‖A < 1Beweis: Wir zeigen zuerst die Ungleichung. Da ρ(M) < 1 aquivalent zu σ(M) ⊂ (−1, 1) ist, gilt nachLemma 5.2 (iv)

ρ(M) < 1⇐⇒ I > I −W−1A und I −W−1A > −I

Nach Voraussetzung sind A und W , und damit auch W−1 positiv definit, womit die erste Ungleichheittrivialerweise erfullt ist und es genugt 2I −W−1A > 0 zu zeigen. Diese Ungleichung ist erfullt, da mit2W > A folgt:

2I −W−1A > 2I − 2W−1W = 0

Bleibt noch die Gleichheit ρ(M) = ||M ||A zu zeigen. Aufgrund der Symmetrie von A12MA−

12 gilt:

‖M‖AS.31= ‖A 1

2MA−12 ‖2 = ρ(A

12MA−

12 )

!= ρ(M)

Letztere Gleichheit ist erfullt, da unter Berucksichtigung des Determinantenmultiplikationssatzes folgen-de Identitaten gelten:

det(A12MA−

12 − λI) = det(A

12 (M − λI)A−

12 )

= det(A12 ) det(A−

12 )︸ ︷︷ ︸

=1

det(M − λI)

1

c)

Z.z.: λ ∈ (0,Λ] und 0 < λW ≤ A ≤ ΛW =⇒ σ(M) ⊂ [1− Λ, 1− λ]Beweis:

σ(M) ⊂ [1− Λ, 1− λ]⇐⇒ I − ΛI ≤M ≤ I − λI

Die erste Ungleichung liefert:

I − ΛI ≤M ⇐⇒ ΛI −W−1A ≥ 0

Diese Ungleichung ist erfullt, da:

ΛI −W−1A ≥ ΛI − ΛW−1W = 0

Analog gilt fur die zweite Ungleichung:

M ≤ I − λI ⇐⇒ −W−1A ≤ −λW−1W = −λI

womit die Behauptung folgt.

d)

Z.z.: ω ∈ (0, 1] =⇒ SSOR konvergiert fur alle SPD-Matrizen A. Es gilt:

MSSORω = I −W−1A

wobei W−1 = ω(2− ω)(D + ωU)−1D(D + ωL)−1

Beweis.: Wegen b) genugt es zu zeigen, dass 2W > A, oder dazu aquivalent 2W −A > 0, gilt.

2W −A =2

ω(2− ω)(D + ωL)D−1(D + ωU)−A

=2

ω(2− ω)(I + ωLD−1)(D + ωU)−A

=2

ω(2− ω)(D + ωL+ ωU + ω2LD−1U)− (L+D + U)

=

(2

ω(2− ω)︸ ︷︷ ︸> 2

2−ω

−1

)D +

(2

2− ω− 1

)L+

(2

2− ω− 1

)U +

2ω

2− ωLD−1LT

>

(2

2− ω− 1

)︸ ︷︷ ︸

>1

A+2ω

2− ω︸ ︷︷ ︸>0

LD−1LT︸ ︷︷ ︸≥0

≥ A > 0

e)

Z.z.: W sei nur noch positiv definit und es gelte W +WT > A > 0. Dann folgt (ρ(M) ≤) ‖M‖A < 1.Beweis:

NT +N = W−T (W +WT )W−1 > W−TAW−1 = NTAN

Sei x mit ||x||A = 1 beliebig. Dann gilt mit a):

||Mx||2A = 〈Mx,Mx〉A= 〈x,MAMx〉A= 〈x,A−1(I −ANT )A(I −NA)x〉A= 〈x, (I −NTA−NA+NTANA)x〉A= 〈x, (I − (NT +N)A+NTANA)x〉A< 〈x, (I −NTANA+NTANA)x〉A= ||x||2A = 1,

womit sofort ||M ||A < 1 folgt.

2

6. Beispiel

Geg.: A SPD und MSORω = I − ω(D + ωL)−1A die Iterationsmatrix des vorwarts SOR-Verfahrens,

MSORω = I −ω(D+ωLT )−1A Iterationsmatrix des ruckwarts SOR-Verfahrens, MSSOR

ω = MSORω MSOR

ω

Iterationsmatrix des gedampften SSOR-Verfahrens mit:

MSSORω = I − ω(2− ω)(D + ωLT )−1D(D + ωL)−1A

a)

Z.z.: MSORω = (MSOR

ω )A

Beweis: Nach Beispiel 5 a) gilt (MSORω )A = A−1(MSOR

ω )TA und damit

A−1(MSORω )TA = A−1(I − ω(D + ωL)−1A)TA

= I − ωA−1A(D + ωLT )−1A

= MSORω

Damit kann unmittelbar auf die nicht-Negativitat des zugehorigen Spektrums geschlossen werden.Z.z.: σ(MSSOR

ω ) ⊂ R+0

〈MSSORω x, x〉A = 〈MSOR

ω MSORω x, x〉A

a)= 〈MSOR

ω x,MSORω x〉A = ‖MSOR

ω x‖2A ≥ 0

Damit ist MSSORω positiv semidefinit, womit die Behauptung folgt.

b)

Z.z.: ‖MSSORω ‖A = ‖MSOR

ω ‖2ABeweis. Nach der Definition der Matrixnorm, a) und der Symmetrie des euklidschen Skalarprodukts gilt:

‖MSORω ‖A = sup

‖x‖A=‖y‖A=1

〈MSORω x, y〉A = sup

‖x‖A=‖y‖A=1

〈x, MSORω y〉A = ‖MSOR

ω ‖A

Damit folgt mit der Submultiplikativitat der Matrixnorm:

‖MSSORω ‖A = ‖MSOR

ω MSORω ‖A ≤ ‖MSOR

ω ‖A‖MSORω ‖A = ‖MSOR

ω ‖2A

Andererseits gilt ebenso fur jene x mit ‖x‖A = 1 nach der Cauchy-Schwarz Ungleichung und der Sub-multiplikativitat:

‖MSORω x‖2A = 〈MSOR

ω x,MSORω x〉A

a)= 〈MSOR

ω MSORω x, x〉A ≤ ‖MSOR

ω MSORω ‖A‖x‖2A = ‖MSSOR

ω ‖A

womit die Gleichheit bewiesen ware.

3