Kapitel 2. L¨osung linearer Gleichungssysteme I Inhalt: 2 ... · Fehlerabsch¨atzungen 2.1...

Kapitel 2. Losung linearer Gleichungssysteme I

Inhalt:

2.1 Fehlerabschatzungen

2.2 Direkte Losungsverfahren

2.3 Spezielle Gleichungssysteme

Numerische Mathematik I 34

Prolog: Wiederholung aus der Linearen Algebra

Gegeben: A = (ajk ) j=1,...,mk=1,...,n

∈ Km×n und b = (bj )j=1,...,m ∈ Km.

Gesucht: x = (xk )k=1,...,n ∈ Kn mit Ax = b.

Satz: Das lineare Gleichungssystem Ax = b ist losbar genau dann, wenn

Rang(A) = Rang ([A, b]) = Rang

a11 . . . a1n b1...

. . ....

...am1 . . . amn bm

.

Im Fall m = n sind die folgenden Aussagen aquivalent:

(i) Ax = b ist fur jedes b ∈ Kn eindeutig losbar;

(ii) Rang(A) = n;

(iii) det(A) 6= 0;

(iv) Alle Eigenwerte von A sind ungleich Null;

(v) A ist invertierbar;

(vi) Ax = 0 hat als eindeutige Losung x = 0.


Fehlerabschatzungen

2.1 Fehlerabschatzungen

Wir behandeln die naturliche Stabilitat des mathematischen Problems

Zu A ∈ Kn×n, b ∈ Kn finde die Losung von Ax = b.

Dazu benotigen wir Normen auf dem Vektorraum Kn×n, die spezielleZusatzeigenschaften haben.

Beachte: Der Vektorraum ist endlich-dimensional, also sind alle Normenaquivalent.


Fehlerabschatzungen Definition (Matrixnorm):

2.1.1 Definition (Matrixnorm):

a) Eine Norm ‖ · ‖ auf Kn×n heißt Matrixnorm, wenn sie submultiplikativ ist,d.h.

(MN1) ‖A · B‖ ≤ ‖A‖ · ‖B‖ fur alle A,B ∈ Kn×n.

b) Eine Norm ‖ · ‖ auf Kn×n heißt vertraglich mit der (Vektor-)Norm ‖ · ‖′ aufKn, wenn

(MN2) ‖Ax‖′ ≤ ‖A‖ · ‖x‖′ fur alle A ∈ Kn×n, x ∈ K

n gilt.

Beispiele: Frobeniusnorm, Zeilensummennorm


Fehlerabschatzungen Definition (naturliche Matrixnorm):

2.1.2 Definition (naturliche Matrixnorm):

Zu einer gegebenen (Vektor-)Norm ‖ · ‖ auf Kn definieren wir

‖A‖op := supx∈Kn\0

‖Ax‖‖x‖ = sup

x∈Kn, ‖x‖=1

‖Ax‖.

‖ · ‖op heißt die naturliche Matrixnorm (oder Operatornorm) zur Norm ‖.‖ auf Kn.

2.1.3 Hilfssatz:

Die naturliche Matrixnorm ‖ · ‖op ist submultiplikativ und vertraglich mit dergegebenen Norm ‖.‖ auf Kn.

Speziell gilt ‖I‖op = 1 fur die Einheitsmatrix I ∈ Kn×n.


Fehlerabschatzungen Hilfssatz:

2.1.4 Hilfssatz:

a) Die naturliche Matrixorm zur Maximums-Norm ‖ · ‖∞ auf Kn ist

‖A‖∞ := max1≤j≤n

n∑

k=1

|ajk |. ( Zeilensummennorm)

b) Die naturliche Matrixorm zur Vektor-Norm ‖ · ‖1 ist

‖A‖1 := max1≤k≤n

n∑

j=1

|ajk |. ( Spaltensummennorm)

c) Die naturliche Matrixorm zur euklidischen Norm ‖ · ‖2 ist

‖A‖2 := max√λ : λ ist Eigenwert von A

TA. ( Spektralnorm)

2.1.5 Diskussion hermitescher und positiv-semidefiniter Matrizen


Fehlerabschatzungen Definition (Konditionszahl):

2.1.6 Definition (Konditionszahl):

Gegeben sei eine Matrixnorm ‖ · ‖ auf Kn×n. Fur eine regulare (=invertierbare)Matrix A ∈ Kn×n heißt

cond(A) = ‖A‖ · ‖A−1‖die Konditionszahl von A.

Die Fehleranalyse bei EXAKTER Matrix und gestorter rechter Seiteb = b +∆b ∈ Kn liefert den

2.1.7 Storungssatz (einfache Form):

Die Matrix A ∈ Kn×n sei regular und die Vektoren x , x ∈ Kn seien die Losungenvon

Ax = b bzw. Ax = b. (b 6= 0)

Weiter seien Vektor- und Matrixnorm vertraglich. Dann gilt fur den relativenFehler

‖x − x‖‖x‖ ≤ cond(A)

‖b − b‖‖b‖ .


Fehlerabschatzungen Hilfssatz:

2.1.8 Hilfssatz:

Es sei ‖ · ‖ eine Matrixnorm und B ∈ Kn×n mit ‖B‖ < 1. Dann ist die MatrixI+ B regular, und es gilt

‖(I+ B)−1‖ ≤ 1

1− ‖B‖ .

2.1.9 Korollar:

Es seien A,B ∈ Kn×n. Die Matrix A sei regular und die Matrix B erfulle‖B‖ < 1

‖A−1‖ . Dann ist auch A+ B regular.


Fehlerabschatzungen Storungssatz (vollstandige Form)

Die Fehleranalyse bei gestorter Matrix A = A+∆A und gestorter rechter Seiteb = b +∆b ∈ Kn liefert den

2.1.10 Storungssatz (vollstandige Form)

Die Matrix A ∈ Kn×n sei regular, die Storung ∆A = A− A erfulle ‖∆A‖ < 1‖A−1‖ .

Die Vektoren x , x ∈ Kn seien die Losungen von

Ax = b bzw. Ax = b. (b 6= 0)

Weiter seien Vektor- und Matrixnorm vertraglich. Dann gilt fur den relativenFehler

‖x − x‖‖x‖ ≤ cond(A)

1− cond(A) · ‖∆A‖‖A‖

(‖∆b‖‖b‖ +

‖∆A‖‖A‖

).



Einschub: Grundlagen zur Einfuhrung der Spektralnorm

Wir betrachten generell K = C.

1. Eigenwerte und Eigenvektoren, Jordan-Normalform:

Die Eigenwerte λj ∈ C, j = 1, . . . , n, einer Matrix A ∈ Cn×n sind die Nullstellen descharakteristischen Polynoms

p(λ) = det(A− λI);

zu jedem Eigenwert λj hat der Eigenraum

Eig(A, λj = v ∈ Cn : Av = λjv

die Dimension sj ≥ 1 (geometrische Vielfachheit von λj ).

Es existiert eine regulare Matrix S ∈ Cn×n (deren Spalten Eigen- und Hauptvektoren von A

sind) mit

S−1AS =

λ1

. . .

λn

+ N

mit einer (nilpotenten) Matrix N ∈∈ Cn×n, die uberall Nullen enthalt mit Ausnahme deroberen Nebendiagonalen, in der Nullen oder Einsen stehen konnen.



2. Spektralradius:

Fur A ∈ Cn×n ist(A) = max|λ| : λ ist Eigenwert von A

der Spektralradius von A: alle Eigenwerte liegen in einer Kreisscheibe in C mit MittelpunktNull und Radius (A).

Bei gegebener Norm ‖ · ‖ in Cn und vertraglicher Matrixnorm ‖.‖′ gilt

(A) ≤ ‖A‖′.

Denn: aus Av = λjv (mit ‖v‖ = 1) folgt

|λj | = |λj |‖v‖ = ‖λjv‖ = ‖Av‖ ≤ ‖A‖′ · ‖v‖ = ‖A‖′.

Speziell: mit der Zeilensummennorm ergibt sich

(A) ≤ ‖A‖∞ = max1≤j≤n

n∑

k=1

|ajk |.



3. Hermitesche Matrizen:A ∈ Cn×n ist hermitesch, wenn A = A

T, d.h. ajk = akj fur alle k, j = 1, . . . , n. (Reelle

Matrizen mit dieser Eigenschaft heißen symmetrisch.)

Hermitesche Matrizen haben nur reelle Eigenwerte; sie besitzen dazu eine Orthonormalbasis von

Eigenvektoren: mit einer Matrix S ∈ U(n) (der Gruppe der unitaren Matrizen, S−1 = ST) ist

STAS =

λ1

. . .

λn

∈ R

n×n.



4. Positiv (semi)-definite Matrizen:

Definition:

Die Matrix A ∈ Cn×n sei hermitesch.

a) A heißt positiv definit, wenn

xTAx =n∑

j,k=1

ajkxjxk > 0 fur alle x ∈ Cn \ 0.

b) A heißt positiv semi-definit, wenn

xTAx =n∑

j,k=1

ajkxjxk ≥ 0 fur alle x ∈ Cn.

A ist positiv (semi-)definit genau dann, wenn alle Eigenwerte von A positiv (bzw.nichtnegativ) sind.



5. Definition der Spektralnorm:Fur beliebiges A ∈ Cn×n ist die Matrix B := A

TA positiv semi-definit: mit der euklidischen

Norm ‖ · ‖2 in Cn ist

xTBx = xTATAx = (Ax)

T(Ax) = ‖Ax‖22 ≥ 0.

Also ist die Spektralnorm

‖A‖2 := max√λ : λ Eigenwert von A

TA =

√(A

TA)

wohldefiniert. Sie ist die naturliche Matrixnorm zur euklidischen Norm ‖ · ‖2 in Cn: furx ∈ Cn mit ‖x‖2 = 1 schreiben wir

x =n∑

j=1

αjvj ,

wobei (v1, . . . , vn) eine Orthonormalbasis des Cn ist, die aus Eigenvektoren von ATA

besteht. Dann ist

1 = ‖x‖22 = xT x =n∑

j,k=1

αjαk vTjvk =

n∑

j=1

|αj |2,

‖Ax‖22 = (Ax)T(Ax) = xTA

TAx =

n∑

j=1

αjvTj·

n∑

k=1

λkαkvk =n∑

j,k=1

αjαkλk vTjvk =

n∑

j=1

|αj |2λj ≤

Damit ist

max‖x‖2=1

‖Ax‖2 ≤√

(ATA),

und Gleichheit gilt fur einen Eigenvektor x zum großten Eigenwert λ = (ATA).



6. Speziell: Spektralnorm hermitescher Matrizen

Fur jede hermitesche Matrix A ∈ Cn×n ist

‖A‖2 = (A),

hier stimmen also Spektralnorm und Spektralradius uberein!

Denn die Eigenwerte von B = ATA = A2 sind die Quadrate der Eigenwerte von A.

Die Konditionszahl einer hermiteschen regularen Matrix A (bzgl. der Spektralnorm) istdaher

cond2(A) =max|λ| : λ Eigenwert von Amin|λ| : λ Eigenwert von A

.


Fehlerabschatzungen Faustregel

2.1.11 Faustregel

Seien cond(A) · ‖∆A‖‖A‖ ≪ 1 und

cond(A) ≈ 10s ,‖∆A‖‖A‖ ≈ 10−k ,

‖∆b‖‖b‖ ≈ 10−k

mit k > s.Dann gilt ‖∆x‖‖x‖ ≈ 10s−k , d.h. man verliert circa s Stellen Genauigkeit in

jeder Komponente.


Fehlerabschatzungen Bemerkung:

2.1.12 Bemerkung:

(i) Falls A ∈ Kn×n regular und schlecht konditioniert ist, ist die Große desResiduums r = b − Ax ein schlechtes Maß fur die Gute von x : Es gilt

‖r‖‖b‖ ≤

‖A‖ · ‖x − x‖‖Ax‖ ≤ cond(A) · ‖x − x‖

‖x‖

und‖x − x‖‖x‖ ≤ ‖A‖ · ‖A

−1(Ax − Ax)‖‖b‖ ≤ cond(A) · ‖r‖‖b‖ .

Daraus folgt

cond−1(A) · ‖r‖‖b‖ ≤‖x − x‖‖x‖ ≤ cond(A) · ‖r‖‖b‖ .

(ii) Die Abschatzung in Satz 2.1.7 ist im Wesentlichen scharf; betrachte hierzudie Spektralnorm, positiv definites A, ∆A = 0, b = v1 und ∆b = ǫvn mitEigenvektoren v1, vn zum großten Eigenwert λ1 = (A) bzw. kleinstenEigenwert λn.


Direkte Losungsverfahren Satz: Elementare Umformungen des LGS

2.2. Direkte Losungsverfahren

Zur Losung des linearen Gleichungssystems Ax = b verwendet man einfacheAquivalenz-Umformungen des Gleichungssystems.

2.2.1 Satz: Elementare Umformungen des LGS

Die Menge der Losungen eines linearen Gleichungssystems bleibt unverandert,wenn man

E1 die Reihenfolge der Gleichungen vertauscht,

E2 beide Seiten einer Gleichung mit einer Zahl α 6= 0 multipliziert,

E3 eine Gleichung ersetzt durch die Summe dieser Gleichung und dem Vielfachen eineranderen Gleichung.

E4 Vertauscht man die Reihenfolge der Unbekannten x1, . . . , xn, setzt also

(y1, . . . , yn) := (xσ1 , . . . , xσn )

mit einer Permutation (σ1, . . . , σn) der Zahlen (1, . . . , n), so erhalt man die Losungsmengedes neuen Systems (bzgl. ~y) aus der Losungsmenge des alten (bzgl. ~x) durchentsprechende Vertauschung der Komponenten.


Direkte Losungsverfahren Die elementaren Umformungen in Matlab/Octave

2.2.2 Die elementaren Umformungen in Matlab/Octave:

MATLAB verwendet die NotationA(j,:) fur die j-te Zeile von A

A(:,k) fur die k-te Spalte von A

Damit lassen sich die elementaren Umformungen einer Matrix einfachschreiben:A([j1,j2],:)=A([j2,j1],:) Zeilentausch von Zeile j1 und j2A(:,[k1,k2])=A(:,[k2,k1]) Spaltentausch von Spalte k1 und k2A(j,:)=t*A(j,:) Multiplikation der Zeile j mit der Zahl tA(j2,:)=A(j2,:)+t*A(j1,:) Addition des t-fachen der Zeile j1 zur Zeile


Direkte Losungsverfahren Die elementaren Umformungen als Matrix-Multiplikation

2.2.3 Die elementaren Umformungen als Matrix-Multiplikation:

E1 Zeilentausch j1 ↔ j2: Pj1,j2 ∗ A mit einfacher Permutationsmatrix

Pj1,j2 =

1. . .

10 1

1. . .

11 0

1. . .

1

← j1

← j2


Direkte Losungsverfahren Die elementaren Umformungen als Matrix-Multiplikation

E3+ Fur j = j1 + 1, . . . , n: addiere das qj -fache der Zeile j1 von der Zeile j :(I + Lj1(q)) ∗ A mit strikter unterer Dreiecksmatrix

Lj1(q) =

0. . .

0qj1+1 0...

. . .

qn 0

← j1, q =

qj1+1

...qn

∈ K

n−j1 .


Direkte Losungsverfahren Hilfssatz:

2.2.4 Hilfssatz:

a) Die einfachen Permutationsmatrizen sind symmetrisch und orthogonal, alsoP−1j1,j2

= Pj1,j2 .

b) Fur 1 ≤ j1 ≤ j2 ≤ n ist

(I + Lj1(p)) · (I + Lj2(q)) = I + Lj1(p) + Lj2(q).

Insbesondere ist (I + Lj1(p))−1 = I − Lj1(p).

c) Fur 1 ≤ j1 < j2 < j3 ≤ n ist

Pj2,j3(I + Lj1(q)) = (I + Lj1(q))Pj2,j3 ,

mit q = Pj2,j3(0, . . . , 0, qj1+1, . . . , qn)T .


Direkte Losungsverfahren Definition: Dreiecksmatrizen

2.2.5 Definition: Dreiecksmatrizen

Eine Matrix der Form

R =

r11 r12 · · · r1n0 r22 · · · r2n...

. . .. . .

...0 · · · 0 rnn

heißt obere Dreiecksmatrix.Eine Matrix der Form

L =

ℓ11 0 · · · 0

ℓ21 ℓ22. . .

......

. . .. . . 0

ℓn,1 · · · ℓn,n−1 ℓnn

heißt untere Dreiecksmatrix.

Bemerkung:

Das Produkt R1R2 oberer Dreiecksmatrizen ist eine obere Dreiecksmatrix.R ist genau dann regular, wenn alle Diagonalelemente ungleich Null sind.


Direkte Losungsverfahren Gauß-Algorithmus mit Pivotierung

2.2.6 Gauß-Algorithmus mit Pivotierung: Die Umformungen werden nach einemfestgelegten System durchgefuhrt. Dabei wird die LR-Zerlegung (LR = Produkteiner unteren und einer oberen Dreiecksmatrix) der zeilen-permutierten Matrix

PA = LR , P Permutationsmatrix,

erzielt. Die anschließende Losung des Gleichungssystems erfolgt durch Vorwarts-bzw. Ruckwartseinsetzen: die Losung von PAx = Pb (=: b) ergibt sich durch

Losen von Ly = b durch Vorwartseinsetzen

y1 =b1

ℓ11, y2 =

1

ℓ22(b2 − ℓ21y1), . . .

Losen von Rx = y durch Ruckwartseinsetzen

xn =yn

rnn, xn−1 =

1

rn−1,n−1(yn−1 − rn−1,nxn), . . .


Direkte Losungsverfahren Gauß-Algorithmus mit Pivotierung

1. Schritt:

Suche in der 1. Spalte das “Pivotelement” aj1, 1 ≤ j ≤ n, mit maximalemAbsolutbetrag und fuhre Zeilenvertauschung 1↔ j durch.

Fur 2 ≤ j ≤ n: Addiere −aj1/a11 · (1. Zeile) zur j-ten Zeile zum Zweck derElimination der 1. Spalte unterhalb von a11.

k-ter Schritt, 2 ≤ k ≤ n − 1:

Suche im unteren Teil der k-ten Spalte das “Pivotelement” ajk , k ≤ j ≤ m,mit maximalen Absolutbetrag und fuhre Zeilenvertauschung k ↔ j durch.

Fur k + 1 ≤ j ≤ n: Addiere −ajk/akk · (k-te Zeile) zur j-ten Zeile zum Zweckder Elimination der k-ten Spalte unterhalb von akk .

Falls im Laufe der Rechnung ein Pivotelement akk mit sehr kleinem Betragentsteht, wird eine Warnung oder Fehlermeldung “A fast singular” ausgegeben.


Direkte Losungsverfahren Beispiel: PA = LR in kompakter Schreibweise

2.2.7 Beispiel: PA = LR in kompakter Schreibweise A =

3 1 62 1 31 1 1

.

1. Schritt:

Pivot-Element a11 = 3 in 1. Spalte, also kein Zeilentausch erforderlich.

Elimination in 1. Spalte: A(1) = (I − L1((23, 13)) ∗ A

3 1 62 1 31 1 1

3 1 62/3 1/3 −11/3 2/3 −1

“Buchfuhrung” uber die wesentlichen Eintrage von (I + L1((23, 13)) in der 1. Spalte


Direkte Losungsverfahren Beispiel: PA = LR in kompakter Schreibweise

2. Schritt:

Pivot-Element a32 = 2/3 in 2. Spalte, also Zeilentausch mit P2,3 =

1 0 00 0 10 1 0

.

Elimination in 2. Spalte: A(2) = (I − L2((12)) ∗ P2,3 ∗ A(1)

3 1 61/3 2/3 −12/3 1/3 −1

3 1 61/3 2/3 −12/3 1/2 −1/2

“Buchfuhrung” uber die wesentlichen Eintrage von I + L1((13, 23)) + L2((

12)) in der 2. Spalte

Ablesen des Ergebnisses:

P2,3A =

3 1 61 1 12 1 3

=

1 0 01/3 1 02/3 1/2 1

3 1 60 2/3 −10 0 −1/2

.


Direkte Losungsverfahren Satz zur LR-Zerlegung:

2.2.8 Satz zur LR-Zerlegung:

Zu jeder regularen Matrix A ∈ Kn×n existiert eine Permutationsmatrix P , eineuntere Dreiecksmatrix L mit Diagonalelementen 1 sowie eine obere DreiecksmatrixR mit

P · A = L · R .Die Matrizen P , L und R ergeben sich aus der Gauß–Elimination mitSpaltenpivotierung.


Direkte Losungsverfahren Komplexitat der LR-Zerlegung

2.2.9 Komplexitat der LR-Zerlegung: Zahlt man die Multiplikationen (mitanschließender Addition) als eine Rechenoperation, so erfordert Schritt k

die Berechnung von n − k Quotienten (=Eintrage von L in Spalte k),

die Berechnung der neuen Eintrage aj,ℓ im Bereich k + 1 ≤ j , ℓ ≤ n.

Der Gesamtaufwand der LR-Zerlegung ist also

n−1∑

k=1

(n−k+(n−k)2) =

n−1∑

k=1

(k+k2) =(n − 1)n

2+(n − 1)n(2n − 1)

6=

n3

3+O(n2).

Das anschließende Vorwarts- und Ruckwartseinsetzen erfordert zusammenn2 +O(n) Operationen.


Direkte Losungsverfahren Konditionierung der Gauß–Elimination mit Pivotierung:

Wenn man Ax = b mit der LR-Zerlegung lost, stellt sich die Frage, wie sichcond(A) zu cond(L)cond(R) verhalt.Wir geben nur eine grobe Auskunft uber L und R , die hierfur nicht ausreicht.

2.2.10 Konditionierung der Gauß–Elimination mit Pivotierung:

Sei A ∈ Kn×n regular und PA = LR die LR-Zerlegung nach demGauß-Algorithmus mit Pivotierung.

a) Die Eintrage der Matrix L erfullen |ℓjk | ≤ 1.

b) Die Eintrage der Matrix R erfullen

max1≤j,k≤n

|rjk | ≤ 2n−1 max1≤j,k≤n

|ajk |.

Vgl. Demmel, Applied Numerical Linear Algebra, SIAM Publ. 1997, S. 49. Der Faktor 2n−1 inder oberen Schranke wird angenommen fur Matrizen des Typs

A =

1 0 0 0 1−1 1 0 0 1−1 −1 1 0 1−1 −1 −1 1 1−1 −1 −1 −1 1

.


Direkte Losungsverfahren Satz: Numerische Stabilitat der LR−Zerlegung

Resultat von Wilkinson zur “Ruckwarts-Stabilitat” des numerischen Verfahrens:

2.2.11 Satz: Numerische Stabilitat der LR−Zerlegung:

Seien A ∈ Kn×n regular und b ∈ Kn. Das Gleichungssystem A · x = b werde mitGauß–Elimination mit Spaltenpivotierung gelost. Dann ist die unter dem Einflussvon Rundungsfehlern tatsachlich berechnete Losung

x = x +∆x

die exakte Losung des Gleichungssystems

(A+∆A) · (x +∆x) = b

mit‖∆A‖∞‖A‖∞

≤ 1.01 · 2n−1 · (n3 + 2n2) · eps.

Bemerkung: Diese Abschatzung ist fur praktische Falle oft zu pessimistisch. Insbesondere furdunnbesetzte Matrizen sind wesentlich bessere Resultate zu erwarten.


Direkte Losungsverfahren Erganzungen zur LR-Zerlegung:

2.2.12 Erganzungen zur LR-Zerlegung::

I. Totalpivotierung

II. Determinanten- und Rangbestimmung

III. Nachiteration

IV. Berechnung der Inversen: Austauschschritte

V. Direkte LR-Zerlegung (Verfahren von Crout)

Wir bezeichnen die Matrix nach dem k-ten Schritt mit A(k) =(a(k)ij

).



I. Totalpivotierung:

Im k-ten Schritt wird anstatt des betragsgroßten Elements a(k−1)j,k , j = k : n, der

k-ten Spalte in dem gesamten Bereich A(k−1)(k : n, k : n) nach dembetragsgroßten Element gesucht. Dieses wird mittels Zeilen- und Spalten-Tausch(E1 und E4) an die Stelle des Pivot-Elements ak,k gebracht.

Dadurch wird die maximal mogliche Vergroßerung der Elemente von A von2n−1 (siehe 2.2.10) auf den Faktor n reduziert.

Der gesamte Aufwand zur Suche der Pivot-Elemente steigt von n2/2 (beiSpaltenpivotisierung) auf n3/3.



II. Determinanten- und Rangbestimmung:

Determinante:

Wegen det(I + Lk(q)) = 1 (mit strikter unterer Dreiecksmatrix Lk(q)) folgtmit dem Determinanten-Multiplikationssatz det L = 1.

Fur einfache Permutationsmatrizen ist detPj,k = −1, 1 ≤ j < k ≤ n.

Also ist

± detA = det(PA) = det(LR) = detR =

n∏

k=1

rkk .

Das Vorzeichen richtet sich danach, ob eine gerade oder ungerade Anzahl vonPermutationen durchgefuhrt wurde.

Rangbestimmung:

Bei regularem A ist RangA = RangR = n, weil die elementarenUmformungen den Rang nicht verandern.

Im Fall RangA < n kann der Fall auftreten, dass alle Spalten-Eintrage

a(k−1)jk , j = k : n, im k-ten Schritt Null sind. Dann ist Spalten-Vertauschungerforderlich (siehe I.). Damit wird die Stufenform von A erzeugt, die eine odermehrere Nullzeilen enthalt. Der Rang von A ist die Anzahl der von Nullverschiedenen Zeilen in der Stufenform.



III. Nachiteration:

Problem: Lose Ax = b mit der LR-Zerlegung PA = LR .

Mit stark verfalschter LR-Zerlegung ∆A = LR − PA,

ǫ1 :=‖∆A‖‖A‖ ≫ eps,

(z.B. durch Abspeichern von L und R in “kurzem” Zahlenformat und

spateres Einlesen) erfullt die Losung von LRx = Pb nach Satz 2.1.10

‖x − x‖‖x‖ ≤ κ(A) · ǫ1, κ(A) =

cond(A)

1− cond(A) · ‖∆A‖‖A‖

.

Das Residuum

r0 = b − Ax (berechnet mit dem exakten A und Rundungsfehler eps)

ergibt den Korrekturvektor k1 als Losung von

LRk1 = Pr0.



Als verbesserte Losung setzen wir x1 = x + k1, denn

x1 − x = x − x + (LR)−1(Pb − PAx)

= (I− (LR)−1PA)(x − x)

= (LR)−1(LR − PA︸︷︷︸=∆A

)(x − x),

also

‖x1 − x‖‖x‖ ≤ ‖(LR)−1‖ ‖∆A‖

‖x − x‖‖x‖ ≤ ‖(LR)−1‖ ‖A‖ ‖∆A‖

‖A‖‖x − x‖‖x‖ .

Mit LR = A+∆A folgt wie im Beweis von Satz 2.1.10

‖x1 − x‖‖x‖ ≤ κ(A)

‖∆A‖‖A‖

‖x − x‖‖x‖ ≤ κ(A)2 ǫ21,

also eine Verdopplung der exakten Dezimalstellen! Weitere Schritte liefernVerbesserung bis zur Rechengenauigkeit eps. In der Praxis genugen 1− 3Schritte der Nachiteration.



IV. Berechnung der Inversen: Lose AX = I

1.Variante:

1. Aufstellen der LR-Zerlegung PA = LR .

2. Losen von LY = P mit Vorwartseinsetzen (simultan fur jede Spalte derrechten Seite).

3. Losen von RX = Y mit Ruckwartseinsetzen.

2.Variante:Vorwarts- und Ruckwarts-Elimination in A mit Skalierung der Zeilen

1. Starte mit (A | I).

2. Fuhre Elimination mit Pivotierung und Skalierung (E2) durch, um auf(I | A−1) zu kommen.

In Matlab/Octave: siehe Kommando rrefmovie([A,eye(n)])



3.Variante: Gauß–Jordan–Algorithmus (Austauschverfahren)Betrachte das LGS Ax = y mit variablen Vektoren x = (x1, . . . , xn)T , y = (y1, . . . , yn)T , undersetze nacheinander die Komponenten von xj durch die yk ’s:

Wenn apq 6= 0 gilt, kann Gleichung p nach xq aufgelost werden:

− ap1

apqx1 − · · · − ap,q−1

apqxq−1+

1

apqyp − ap,q+1

apqxq+1 − · · · − ap,n

apqxn = xq. (Pivotzeile)

Dies ist die neue p-te Zeile des aquivalenten Gleichungssystems

A(1)(x1, . . . , xq−1, yp, xq+1, . . . , xn)T = (y1, . . . , yp−1, xq, yp+1, . . . , yn)

T .

In den anderen Gleichungen ist xq ebenfalls zu ersetzen; dies liefert die Eintrage von A(1) zum 1.Austauschschritt:

Pivotelement : a(1)pq =

1

apq,

Pivotzeile : a(1)pk

= −apk

apq, k 6= q,

Pivotspalte : a(1)jq

=ajq

apq, j 6= p,

sonstige : a(1)jk

= ajk − ajqapkapq

, j 6= p, k 6= q

In den folgenden Austausch-Schritten wird nur außerhalb der p-ten Zeile und q-ten Spalte nachdem Pivotelement gesucht.

Bemerkung: Aufwand n3 +O(n2).



V. Direkte LR-Zerlegung (Verfahren von Crout)

Falls die LR-Zerlegung ohne Pivotierung durchgefuhrt werden kann, liefert der Ansatz

A = LR =

1ℓ21 1...

. . .

ℓn1 · · · ℓn,n−1 1

r11 r12 · · · r1nr22 · · · r2n

. . ....rnn

durch Betrachten der Eintrage von A nacheinander die Gleichungen

1. Zeile k = 1, . . . , n : a1k = r1k

1. Spalte j = 2, . . . , n : aj1 = ℓj1r11 ⇒ ℓj1 =aj1

r11...

p-te Zeile k = p, . . . , n : apk =

p−1∑

µ=1

ℓpµrµk + rpk

⇒ rpk = apk −p−1∑

µ=1

ℓpµrµ,k

p-te Spalte j = p + 1, . . . , n : ajp =

p∑

µ=1

ℓjµrµp

⇒ ℓjp =1

rpp

ajp −

p−1∑

µ=1

ℓjµrµp


Spezielle Gleichungssysteme Satz: LR-Zerlegung ohne Pivotierung

Spezielle Gleichungssysteme

Wir betrachten Gleichungssysteme mit spezieller Struktur, bei denen dieLR-Zerlegung ohne Pivotierung Vorteile hat.

2.3.1 Satz: LR-Zerlegung ohne Pivotierung

Es sei A ∈ Kn×n regular. Die LR-Zerlegung A = LR mit unterer Dreiecksmatrix L

und oberer Dreiecksmatrix R existiert genau dann, wenn jede Teilmatrix

Ak = (aij)i ,j=1,...,k , k = 1, . . . , n,

regular ist. Weiterhin ist dann die LR-Zerlegung von A eindeutig.

Beweis:1. Falls die LR-Zerlegung A = LR existiert, gilt fur 1 ≤ k ≤ n

Ak =

1ℓ21 1...

. . .

ℓk1 · · · ℓk,k−1 1

r11 r12 · · · r1kr22 · · · r2k

. . ....rkk

,

also detAk =∏k

j=1 rjj 6= 0. Deshalb ist Ak regular.

2. Fur die umgekehrte Schlussrichtung fuhren wir Induktion nach n.Numerische Mathematik I 73


Fur n = 1, also A ∈ K mit A = detA 6= 0, ist A = 1 · A die LR-Zerlegung.

Sei nun n ∈ N, n ≥ 2, A ∈ Kn×n regular und detAk 6= 0 fur k = 1, . . . , n.

Induktionsannahme: alle regularen Matrizen B ∈ KK (n−1)×(n−1) mit der Eigenschaft detBk 6= 0fur k = 1, . . . , n − 1 besitzen eine LR-Zerlegung.

Wegen a11 = detA1 6= 0 kann der 1. Eliminationsschritt ohne Zeilenvertauschung ausgefuhrtwerden. Dies liefert

A = (I + L1(q))

a11 a12 · · · a1n

0 B

.

mit q = 1a11

(a21, . . . , an1)T und B ∈ K(n−1)×(n−1). Die Teilmatrix Ak von A ergibt sich hieraus

als

Ak =

1q1 1q2 0 1...

.... . .

. . .

qk−1 0 · · · 0 1

a11 a12 · · · a1k

0 Bk−1

fur 2 ≤ k ≤ n, mit der entsprechenden Teilmatrix Bk−1 von B. Wegen0 6= detAk = a11 detBk−1 ist detBk 6= 0 fur k = 1, . . . , n − 1. Also besitzt die Matrix B nach

der Induktionsannahme die LR-Zerlegung B = LR, und damit ist

A = (I+L1(q))

a11 a12 · · · a1n

0 LR

= (I+ L1(q))

1 0 · · · 0

0 L

︸︷︷︸

a11 a12 · · · a1n

0 R



3. Die Eindeutigkeit der LR-Zerlegung folgt so: Aus A = L1R1 = L2R2 mit unterenDreiecksmatrizen L1,L2, deren Diagonalelemente 1 sind, und oberen Dreiecksmatrizen R1,R2

folgtL−12 AR−1

1 = L−12 L1 = R2R

−11 .

Die Matrix L−12 L1 ist untere Dreiecksmatrix mit Diagonalelementen 1, die Matrix R2R

−11 ist

obere Dreiecksmatrix. Aus der Gleichheit folgt

L−12 L1 = I = R2R

−11 .

Die Eindeutigkeit der Inversen liefert L2 = L1 und R2 = R1.


Spezielle Gleichungssysteme Satz: Positiv definite Matrizen

2.3.2 Satz: Positiv definite Matrizen

Eine symmetrische Matrix A ∈ Rn×n ist positiv definit genau dann, wenn fur dieTeilmatrizen Ak = (aij)i ,j=1,...,k gilt

detAk > 0, k = 1, . . . , n.

An Stelle der LR-Zerlegung A = LR wird eine ahnliche Zerlegung bevorzugt.

2.3.3 Folgerung: Cholesky-Zerlegung

Zu einer symmetrischen positiv-definiten Matrix A ∈ Rn×n existiert eine untereDreiecksmatrix L mit

A = LLT .

Achtung: L hat positive Diagonalelemente, die nicht alle 1 sein mussen (imGegensatz zur LR-Zerlegung).


Spezielle Gleichungssysteme Algorithmus zur Cholesky-Zerlegung:

Man berechnet den Cholesky-Faktor L ahnlich wie im Algorithmus von Crout:2.3.4 Algorithmus zur Cholesky-Zerlegung:Eingabe: A ∈ Rn×n symmetrisch

Berechnung: Fur k = 1, . . . , n

(1) Berechne ℓk,k =

(ak,k −

k−1∑

µ=1

ℓ2k,µ

)1/2

(2) Berechne fur j = k + 1, . . . , n ℓj,k =1

ℓk,k

(aj,k −

k−1∑

µ=1

ℓj,µℓk,µ

)

Bemerkung:

(i) Ob A positiv-definit ist, stellt man im Algorithmus fest: der Ausdruck unterder Wurzel muss positiv sein! Eine Fehlernachricht “A nicht positiv definit”

wird ausgegeben, falls diese Bedingung verletzt wird.(ii) Der Rechenaufwand fur die Berechnung der Cholesky–Zerlegung einer positiv

definiten Matrix (Wurzel zahlt als eine Operation wie Mult./Add.) ist

n∑

k=1

(k + (n − k)k) =

n∑

k=1

k(n − k + 1) =n3

6+O

(n2).


Spezielle Gleichungssysteme Algorithmus zur Cholesky-Zerlegung:

In Matlab/Octave (mit Tests der Eingangsdaten auf richtige Dimension, Symmetrie, positive Definitheit und Warnung bei sehr kleinem Diagonalelementvon L):

function L = mychol(A)

% teste Dimension von A

[m,n]=size(A);

if m = n

error(’Matrix ist nicht quadratisch’)

end

% teste auf Symmetrie

if norm(triu(A,1)-tril(A,-1)’)>n*eps

error(’Matrix ist nicht symmetrisch’)

end

L=zeros(n,n);

g=A(1,1);

for k=1:n-1

% teste auf positive Definitheit

if g<=0

error(’Matrix ist nicht positiv definit’)

end

L(k,k)=sqrt(g);

% teste auf numerische Invertierbarkeit

if L(k,k)<eps

fprintf(’Warnung: Matrix ist eventuell nicht invertierbar’)

end

L((k+1):n,k)=1/L(k,k)*(A((k+1):n,k)-L((k+1):n,1:(k-1))*L(k,1:(k-1))’);

g=A(k+1,k+1)-L(k+1,1:k)*L(k+1,1:k)’;

end

Der Test L=mychol(hilb(3)) ergibt

L =

1 0 012

1√

120

13

1√

121

√

180

Probe: norm(hilb(3)-L*L’)


Spezielle Gleichungssysteme Definition: Bandmatrizen

2.3.5 Definition: Bandmatrizen

Eine Matrix A ∈ Kn×n heißt Bandmatrix vom Typ (p, q) mit 0 ≤ p, q ≤ n − 1,wenn

aj,k = 0 fur alle Paare (j , k) mit k < j − p oder k > j + q.

Die Große p + q + 1 ist die Bandbreite von A.

Bemerkung:

Eine “kompakte” Speicherung einer Bandmatrix vom Typ (p, q) erfordert nurn(p+ q+1) Speicherplatze, jeweils n fur die Diagonale und Nebendiagonalen.Fur eine symmetrische Bandmatrix vom Typ (p, p) speichert man nur dieDiagonale und die unteren Nebendiagonalen, also nur n(p + 1) Eintrage.

Beispiel: symmetrische Tridiagonalmatrix (Typ (1,1)):

A =

a1 b1b1 a2 b2

. . .

bn−2 an−1 bn−1

bn−1 an

↔(a1 a2 · · · anb1 · · · bn−1 0

)


Spezielle Gleichungssysteme Satz:

2.3.6 Satz:

Es sei A ∈ Kn×n eine Bandmatrix vom Typ (p, q). Falls die LR−ZerlegungA = L · R existiert, dann ist L eine Bandmatrix vom Typ (p, 0) und R eineBandmatrix vom Typ (0, q).

Bemerkung:

Die LR-Zerlegung A = LR kann also in kompakter Speicherung durchgefuhrtwerden. Der Rechenaufwand in jedem Eliminationsschritt ist ungefahr gleich:

p Divisionen,

p∑

j=1

(q + j) Mult./Add.,

also insgesamt n(pq + p2

2

)+O (n · (p + q)) Operationen.

Fur positiv-definite Bandmatrizen ist der Cholesky-Faktor L eine Bandmatrixvom Typ (p, 0) ist.


Spezielle Gleichungssysteme Algorithmus: Cholesky-Zerlegung positiv-definiter Tridiagonalmatrizen

2.3.7 Algorithmus: Cholesky-Zerlegung positiv-definiter Tridiagonalmatrizen

Eingabe: Diagonale (a1, . . . , an) und untere Nebendiagonale(b1, . . . , bn−1)

Berechnung: Fur k = 1, . . . , n

(1) Berechne ℓk,k =(ak − ℓ2k,k−1

)1/2

(2) Berechne ℓk+1,k =1

ℓk,kbk

Bemerkung: Der Rechenaufwand fur die Berechnung der Cholesky–Zerlegung einerpositiv definiten Tridiagonalmatrix (Wurzel zahlt als eine Operation wieMult./Add.) ist 2n.


Spezielle Gleichungssysteme Beispiel: Positiv definite Bandmatrizen in den Anwendungen

2.3.8 Beispiel: Positiv definite Bandmatrizen in den AnwendungenDie Poisson-Gleichung (auch Potentialgleichung)

−∆u = f , (∆ =∂2

∂x2+

∂2

∂y2)

beschreibt das elektrostatische Potential u bei gegebener Ladungsdichte f in einem GebietΩ ⊂ R2. ∆ ist ein Differentialoperator 2. Ordnung, der Laplace-Operator. Mit Vorgabe von u aufdem Rand ∂Ω (Randbedingungen) stellt sich das Dirichlet-Problem:Finde u ∈ C2(Ω) mit

−∆u(x , y) = f (x , y), (x , y) ∈ Ω,u(x , y) = g(x , y), (x , y) ∈ ∂Ω.

Beispiel: Ω = [0, 1]2, f (x , y) = x2(x − 1)(2 − 6y) + y2(y − 1)(2 − 6x) in Ω, und homogeneRandbedingungen u(x , y) = 0 fur (x , y) ∈ ∂Ω ergibt

u(x , y) = x2(1− x)y2(1 − y).

Der Funktionsgraph in Matlab/Octave:

[xx,yy]=meshgrid(0:.05:1);

mesh(xx,yy,xx.∧2.*(1-xx).*yy.∧2.*(1-yy)) fur 3D-Plotcontour(xx,yy,xx.∧2.*(1-xx).*yy.∧2.*(1-yy)) fur Kontur-Plot



Differenzenverfahren: Eine numerische Naherungslosung wird durch Diskretisierung desLaplace-Operators als “5-Punkte-Stern” erzielt:

−∆u(x , y) ≈ 4u(x , y)− u(x + h, y)− u(x − h, y)− u(x , y + h) − u(x , y − h)

h2

an den Stellen (x , y) = (jh, kh) mit 1 ≤ j , k ≤ N und h = 1N+1

. Die Nummerierung dieser Stellenergibt den Losungsvektor

~u = (u1, . . . , uN , uN+1, . . . , u2N , . . . , uN2 )T

und den entsprechenden Vektor der rechten Seite ~f .

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1

5

9

13

2

6

10

14

3

7

11

15

4

8

12

16



Verwendung der homogenen Randbedingungen fuhrt dann zum linearen Gleichungssystem

B −I

−I B −I

. . .

−I B −I

−I B

~u = h2~f , B =

4 −1−1 4 −1

. . .

−1 4 −1−1 4

∈ RN×N ,

mit der N × N Einheitsmatrix I. Die gesamte Matrix A ∈ RN2×N2ist

symmetrische Bandmatrix vom Typ (N,N),

diagonaldominant, d.h.∑

k 6=j |ajk | ≤ |aj,j | fur alle j = 1, . . . ,N2,

unzerlegbar, d.h. keine Zeilen- und Spaltenpermutation fuhrt auf eine Form

(A1 A2

0 A3

)

mit quadratischen Matrizen A1,A3,

in mindestens einer Zeile sogar stark diagonaldominant.

Weil außerdem alle Diagonalelemente ajj positiv sind, ist A positiv-definit! Der Aufwand zur

Losung mit Cholesky-Zerlegung ist N4

2(im Vergleich zu N6

6bei voll besetzter Matrix, denn

n = N2).

Beachte: Der Cholesky-Faktor L ist Bandmatrix vom Typ (N, 0) mit von Null verschiedenenElementen in allen unteren Nebendiagonalen (sog. “Auffullen”).



Das Matlab/Octave-Skript poisson.m liefert die Losung u als N × N-Matrix.

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.01

0.02

0.03


Kapitel 2. L¨osung linearer Gleichungssysteme I Inhalt: 2 ... · Fehlerabsch¨atzungen 2.1...

Documents

Transcript of Kapitel 2. L¨osung linearer Gleichungssysteme I Inhalt: 2 ... · Fehlerabsch¨atzungen 2.1...