Gleichungssystemeinstitute.unileoben.ac.at/amat/lehrbetrieb/num/vl-skript/... · 2010-03-17 ·...
Transcript of Gleichungssystemeinstitute.unileoben.ac.at/amat/lehrbetrieb/num/vl-skript/... · 2010-03-17 ·...
Funfte Vorlesung, 18. Marz 2010, Inhalt
GleichungssystemeDie Folien der vierten Vorlesung sind, weil großteils noch nicht behandelt,hier noch-mal enthalten:
• Matrizenrechnung
• Lineare Gleichungssysteme: Einfuhrung, Losbarkeit, MATLAB-Befehle, Fehler-Empfindlichkeit
Neue Folien (21–26):
• Iterative Verfahren: Jacobi, Gauß-Seidel, SOR
1
Lineares Gleichungssystem in n Gleichungen und Unbekannten.
a11 x1 + a12 x2+ . . . +a1n xn = b1a21 x1 + a22 x2+ . . . +a2n xn = b2... ... ...an1 x1 + an2 x2+ . . . +ann xn = bn
In Matrixschreibweise: Ax = b .
Gleichungssysteme lassen sich in Matrix-Schreibweise ubersichtlich
und pragnant formulieren.
Machen Sie sich mit den Bezeichnungen und Regeln der Matrizen-
rechnung vertraut!
2
Einheitsmatrix, Inverse und Transponierte
Die 3×3-Einheitsmatrix. Wichtige Ei-genschaft: A · I = I ·A = A
I =
1 0 00 1 00 0 1
Ein recht simples Beispiel:Das Magische Quadrat der Ordnung 3
M =
8 1 63 5 74 9 2
Die Inverse von M . Wichtige Eigen-schaft: M ·M−1 = M−1 ·M = I
M−1 =1
360
53 −52 23−22 8 38−7 68 −37
Die transponierte Matrix MT
MT =
8 3 41 5 96 7 2
MATLAB: eye(3), m=magic(3), inv(m), m’
3
Matrix-Algebra: Produkt von Matrizen
A ·B = C
Element cij ist (Skalar-)Produkt der i-ten Zeile von A mit der j-ten
Spalte von B.
A : (ℓ×m)-Matrix
B : (m× n)-Matrix
C : (ℓ×m) · (m× n) → (ℓ× n)-Matrix
Spaltenanzahl von A und Zeilenanzahl von B mussen ubereinstimmen!
Im Allgemeinen nicht kommutativ: A ·B 6= B ·A!
4
Matrizenrechnung in MATLAB
A’ Matrixtransposition AT
A+B, A-B Matrixaddition, -subtraktionA.*B, A./B elementweise Multiplikation, Division
A*B Matrixmultiplikation
A/B sehr schrage Schreibung:”Division“ A ·B−1
A\B noch schrager:”Division“ A−1 ·B
Wichtiger Befehl: x=A\b liefert Losung des Gleichungssystems
A · x = b
(aber nicht in der Form x = A−1 ·b, sondern durch Gauß-Elimination.
5
Gleichungssysteme: Aufgabenstellung
Eine der wichtigsten Aufgaben bei technisch-wissenschaftlichen Berechnungen istdas Losen linearer Gleichungssysteme.
Ax = b gegeben: A,b. gesucht: Vektor x
AX = B allgemeiner: B und X konnen auch Matrizen sein
XA = B diese Schreibweise kommt selten vor
Die Matrix A ist oft quadratisch und man sucht eine eindeutige Losung.
Ist A eine m×n-Matrix mit m > n (mehr Zeilen als Spalten, d.h. mehr Gleichungenals Unbekannte): uberbestimmtes System, man kann Losung mit kleinsten Fehler-
quadrat suchen.
Ist m < n: unterbestimmtes System. Gesucht: ein- oder mehrparametrige Losungs-
schar
6
Einige wichtige Themen bei großen linearen Gleichungssystemen
Rechenzeit und Speicherplatz: Speicherbedarf fur volle Matrix ∼ O(N2), Rechenzeit∼ O(N3)
Uber-, unterbestimmte oder schlecht konditionierte Systeme: Eindeutige Losungnicht moglich oder extrem fehlerempfindlich.
Dunn besetzte Matrizen: Oft sind nur wenige % der Matrixelemente 6= 0. SpezielleDatenstrukturen und Algorithmen benotigen signifikant weniger Speicherplatz undRechenzeit.
Direkte Verfahren sind Varianten des Gaußschen Eliminationsverfahrens (LR-Zerlegung,Cholesky-Zerlegung). Ublich bis n ≈ 1000.
Iterative Verfahren finden schrittweise verbesserte Naherungslosungen. Ublich nurfur dunn besetzte Matrizen, ab n ≈ 10.000.
7
Fallunterscheidungen bei Gleichungssystemen
Vorubung: eine Gleichung, eine Unbekannte
Die lineare Gleichung
ax = b a, b gegeben, x gesucht
• hat eine eindeutige Losung genau dann, wenn a 6= 0
• unendlich viele Losungen genau dann, wenn a = b = 0
• keine Losung genau dann, wenn a = 0, b 6= 0
Bei Systemen linearer Gleichungen treten diese Falle in ganz ahnlicher
Weise auf.
8
Zwei Gleichungen, zwei Unbekannte
[
1 23 4
]
·
[
xy
]
=
[
12
]
eindeutige Losungx = 0y = 1/2
[
1 23 6
]
·
[
xy
]
=
[
13
]
∞ viele Losungenx = 0, 1, 2, 3, . . .y = 1/2, 0 −1/2, −1, . . .
[
1 23 6
]
·
[
xy
]
=
[
12
]
keine Losung
Die Determinante determiniert: Im ersten Fall ist detA 6= 0,
in den anderen beiden Fallen ist detA = 0
9
Gleichungssysteme Ax = b mit quadratischer Matrix
Ein lineares Gleichungssystem mit quadratischer n× n-Matrix ist
• genau dann eindeutig losbar, wenn rangA = n (det(A) 6= 0).Andernfalls gibt es
• keine eindeutige (wenn rangA = rang[A,b]), oder
• uberhaupt keine Losung (wenn rangA 6= rang[A,b]).
Der Fall det(A) 6= 0 ist (algebraisch) gleichbedeutend mit
• Alle n Zeilen von A sind linear unabhangig (ebenso die Spalten).
• rang(A) = n Wissen Sie, was der Rang einer Matrix ist?
• A ist nichtsingular, A ist regular
MATLAB: x=A\b findet die eindeutige Losung, falls detA 6= 0
10
Sonderfall 1 bei singularer quadratischer Matrix , det(A) = 0
A =
1 2 34 5 67 8 9
,b =
123
Die Zeilen von A sind linear abhangig:A(3, :) = 2A(2, :)−A(1, :)
Dieselbe Abhangikeit gilt auch, wenn man die rechteSeite miteinbezieht
Wenn Matrix und erweiterte Matrix gleichen Rang haben,
rangA = rang[A,b] = k < n, dann hat das System eine
(n− k)-parametrige Losungsschar.
x =
−1/3+ z2/3− 2z
z
In diesem Beispiel: rangA = rang[A,b] = 2,daher (3− 2) = 1-parametrige Losungsschar;z ist freier Parameter
MATLAB: x=A\b gibt eine Warnung. Findet (manchmal) eine Losung, in der moglichstviele Komponenten = 0 sind. Fur dieses Beispiel x1 = −1
3, x2 = 2
3, x3 = 0
11
Sonderfall 2 bei singularer quadratischer Matrix , det(A) = 0
A =
1 2 34 5 67 8 9
,b =
124
Die Zeilen von A sind linear abhangig. Das gilt abernicht mehr fur die erweiterte Matrix [A,b].
Das Doppelte von Gl. 2, minus Gl. 1 liefert Gl. 3 mitanderer rechter Seite −→ Widerspruch!
Wenn der Rang der Matrix ungleich dem Rang der erweiterten Matrix
ist, rangA 6= rang[A,b], dann hat das System keine (exakte) Losung.
MATLAB: x=A\b gibt eine Warnung. Findet keine Losung, oder eine
”Losung“ mit vollig falschen Zahlenwerten.
12
MATLAB-Befehle fur quadratische Systeme Ax = b
x=A\b berechnet mit einem Eliminationsverfahren
• fur nichtsingulare Matrizen die eindeutige Losung.
• fur singulare Matrizen, wenn eine Losungsschar existiert, manchmal eine Losungmit moglichst vielen Null-Komponenten,
• fur singulare Matrizen, wenn keine Losung existiert, meistens Unsinn.
x=pinv(A)*b berechnet mit Hilfe der Pseudoinversen
• fur nichtsingulare Matrizen mit unnotig hohem Aufwand die eindeutige Losung.
• fur singulare Matrizen, wenn eine Losungsschar existiert, eine Losung mit mi-nimaler 2-Norm.
• fur singulare Matrizen, wenn keine Losung existiert, die”am wenigsten falsche
Losung“ : jene, fur die der Rest-Vektor r = Ax− b am kleinsten ist. Kleinste-
Quadrate-Losung ‖Ax− b‖2 → min!
13
MATLAB-Befehle fur quadratische Systeme Ax = b
rank(A) und rank([A,b]) geben MATLABs Meinung zum Rang von
Matrix und erweiterter Matrix.
Linear unabhangige Matrixzeilen konnen durch beliebig kleine Anderungen in den
Koeffizienten (z. B. aufgrund von Rundungsfehlern) abhangig werden. Rangbestim-
mung ist daher numerisch heikel und von einer Fehlertoleranz abhangig. MATLABs
rank-Funktion verwendet Singularwertzerlegung, ein aufwendiges, aber verlassliches
Verfahren.
rref([A,b]) berechnet mit einem speziellen Eliminationsverfahren (Gauß-
Jordan) eine reduzierte Treppenform des Systems (reduced row eche-
lon form), aus der alle Losungsfalle abgelesen werden konnen.
14
Beispiele zu rref
A =
1 1 11 2 31 3 6
, b =
2361114
, rref([A,b]) =
1 0 0 00 1 0 80 0 1 15
Eindeutige Losung; links die Einheitsmatrix, rechts der Losungsvektor.
A =
2 3 4 53 5 7 94 7 10 135 9 13 17
, b =
1111
, rref([A,b]) =
1 0 −1 −2 20 1 2 3 −10 0 0 0 00 0 0 0 0
Zwei Nullzeilen in Matrix und erweiterter Matrix: Rang = n-2, zweiparametrigeLosungsschar: wahle x3 und x4 beliebig, bestimme x1 und x2 aus den ersten beidenGleichungen.
x1 = 2+ x3 +2x4, x2 = −1− 2x3 +3x4
15
Mit derselben Matrix wie vorhin, aber anderer rechter Seite
A =
2 3 4 53 5 7 94 7 10 135 9 13 17
, b =
1110
, rref([A,b]) =
1 0 −1 −2 00 1 2 3 00 0 0 0 10 0 0 0 0
Zwei Nullzeilen in Matrix, aber nur eine in erweiterter Matrix: keine Losung.
Gleichungssysteme nicht durch Multiplikation mit Inverser losen!
Die Gleichung 7x = 13 losen Sie ja auch nicht, indem sie zuerst
7−1 = 0,1429 berechnen und dann x = 0,1429 ∗ 13 auswerten.
Sie formen die Gleichung naturlich um und rechnen x = 13/7.
MATLABs Befehl fur Gleichungslosen ist \ und nicht inv( ).
In MATLAB differieren inv(7)*13 und 7\13 um −2,2204 · 10−16.
Berechnung der Inversen und anschließende Multiplikation ist sehr viel
rechenaufwandiger und anfalliger gegen Rundungsfehler als direktes
Gleichungslosen!
16
Fehlerempfindlichkeit, Konditionszahl
Die Konditionszahl κ(A) einer Matrix A misst fur das Gleichungssysten
Ax = b, wie empfindlich der relative Fehler von x von kleinen relativen
Anderungen in A und b abhangt.
‖δx‖
‖x‖≤ κ(A)
(
‖δA‖
‖A‖+
‖δb‖
‖b‖
)
mit κ(A) = ‖A‖‖A−1‖
Der relative Fehler in x kann also κ(A) mal großer als der relative
Fehler in A und b sein.
17
Determinante, Rang und Konditionszahl
Das Kriterium detA 6= 0 fur die Losbarkeit eines Gleichungssystems ist nur beiMatrizen mit wenigen Zeilen und kleinen ganzzahligen Elementen angebracht. An-sonsten konnen Rundungsfehler das Ergebnis verfalschen. MATLAB-Beispiel:A=gallery(’rosser’) liefert 8× 8-Matrix, det(A)=-13017 (korrekt ware 0!).
Der Rang einer Matrix wird in MATLAB numerisch verlasslicher berechnet. DasKriterium rangA = n ist besser geeignet, um Losbarkeit festzustellen. Im obigenBeispiel rechnet MATLAB (korrekt): rank(A)=7
Die Konditionszahl lasst besser abschatzen, ob eine Matrix singular oder nahezusingular ist. MATLAB liefert cond(A)=4.7886e+016. Bei 16-stelliger Genauigkeit be-deutet das: A ist numerisch singular!
Der MATLAB-Operator \ zum Gleichungslosen macht automatisch eine Schatzungder reziproken Konditionszahl und liefert Warnung.
>> A\bWarning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 5.654224e-017.
18
Konditionszahl, Beispiel 4× 4 Hilbertmatrix H
H =
1 12
13
14
12
13
14
15
13
14
15
16
14
15
16
17
, H−1 =
16 −120 240 −140−120 1200 −2700 1680240 −2700 6480 −4200−140 1680 −4200 2800
Zeilensummennorm von Matrix und Inverser:
‖H‖∞ = 2,08; ‖H−1‖∞ = 13620; κ(H) = 28375
Storungen in den Daten werden ≈ 28000-fach verstarkt.
19
4× 4 Hilbertmatrix H, rechte Seite b = (1,1,1,1)T
Exakte Losung ist
x =
−460
−180140
Andere das Element h44: addiere 1/1000. Losung andert sich auf
x =
1,15789−1,89474−25,263236,8421
20
Losungsverfahren fur lineare Gleichungssysteme
Direkte Verfahren sind Varianten des Gaußschen Eliminationsver-
fahrens (LR-Zerlegung). Ublich bis n ≈ 10.000 Unbekannten.
Direkte Verfahren sind allgemeiner anwendbar und rechnen zu-
meist schneller, sofern die Matrix im schnell zuganglichen Speicher
des Rechners Platz hat.
Iterative Verfahren finden schrittweise verbesserte Naherungslosun-
gen. Ublich fur n ≫ 10.000.
Iterative Methoden sind nur fur spezielle Matrixtypen anwendbar,
die beispielsweise bei partiellen Differentialgleichungen auftreten.
21
Jacobi-Verfahren
salopp: Lose jede Gleichung nach ihrem Diagonal-Term auf, setze
Startwerte ein, iteriere.
Komponentenweise Schreibung (3× 3-Beispiel)
x1 = (b1 − a12x2 − a13x3)/a11
x2 = (b2 − a21x1 − a23x3)/a22
x3 = (b3 − a31x1 − a32x2)/a33
Matrix-Schreibweise: Jacobi-Verfahren entsteht aus Fixpunkt-Iteration
A = D + E, Ax = b umformen −→ x = D−1(b− Ex)
22
Weitere Schreibweisen des Jacobi-Verfahrens
Schreibung mit Indizes umstandlich und unubersichtlich; nicht merken!
fur i = 1, . . . , n
x(k+1)i =
bi −i−1∑
j=1
aijx(k)j −
n∑
j=i+1
aijx(k)j
/aii
Matlab ubersichtliche Matrix-Schreibweise
D=diag(diag(A)); E = A-D;
x = D\(b-E*x);
23
Konvergenz des Jacobi-Verfahrens
Das Jacobi-Verfahren konvergiert bei Gleichungssystemen mit stark
diagonaldominanter Matrix fur beliebige Startwerte zur exakten Losung.
Eine n× n-Matrix A = (aij) ist stark diagonaldominant, wenn
|aii| >n∑
j=1,j 6=i
|aij| fur i = 1,2, . . . , n
Weil Matrizen aus Anwendungen das oft nicht erfullen, gibt es allge-
meinere Fassungen der Konvergenzbedingung. Trotzdem:
Jacobi-Verfahren konvergiert in der Regel viel zu langsam fur prakti-
schen Einsatz!
24
Iterationsschritt des Gauß-Seidel-Verfahrens
salopp: Im Prinzip wie das Jacobi-Verfahren, nur: setze neue Werte,
soweit verfugbar, schon im aktuellen Schritt ein.
Matrix-Schreibweise (L ist Diagonale + Eintrage darunter in A)
setze A = L+ E, iteriere x(k+1) = L−1(b− Ex
(k))
Komponentenweise Schreibung umstandlich und unubersichtlich; nicht merken!
fur i = 1, . . . , n : x(k+1)i =
bi −
i−1∑
j=1
aijxj(k+1) −
n∑
j=i+1
aijx(k)j
/aii
Matlab ubersichtliche Matrix-Schreibweise
L = tril(A) ; E = A-L;x = L\(b-E*x);
25
Iterationsschritt des SOR-Verfahrens
salopp: Jeweils neuer Naherungswert zuerst als Zwischenresultat y(k+1)i aus Gauß-
Seidel-Schritt; endgultiger Naherungswert durch Extrapolation (Uberrelaxation) ausalter Naherung und Zwischenresultat.
Matrix-Schreibweise: Analog zu Jacobi- und Gauß-Seidel-Verfahren; ersetze dort Dbzw. L durch L+ (1
ω− 1)D.
Komponentenweise Schreibung umstandlich und unubersichtlich; nicht merken!
fur i = 1, . . . , n
y(k+1)i =
bi −
i−1∑
j=1
aijx(k+1)j −
n∑
j=i+1
aijx(k)j
/aii
x(k+1)i = ωy(k+1)
i + (1− ω)x(k)i
Theorie: 1 ≤ ω < 2, eher in der Nahe von 2. Fur ω = 1 reduziet sich SOR aufGauß-Seidel.
26