Karlsruhe Institute of Technology - math.kit.edu · Karlsruhe Institute of Technology 1...

Karlsruhe Institute of Technology1 Arithmetische Grundlagen

Am 4. Juni 1996 explodierte kurz nach dem Start dieerste Ariane 5 Rakete durch einen Softwarefehler.Die Horizontalgeschwindigkeit wurde durch eineGleitkommazahl

v ∈ [−10308,−10−308]∪{0}∪ [10−308,10308]

dargestellt.Innerhalb der Rechnung wurde die Zahl versehentlichin eine ganzzahlige Darstellung

i ∈ {0,1,2, ...,32767}

konvertiert.Als die Geschwindigkeit v > 32767 erreichte, verlor dieSoftware die Geschwindigkeitsinformation und damitschließlich die Orientierung.

http://ta.twi.tudelft.nl/users/vuik/wi211/disasters.html

http://www5.in.tum.de/ huckle/bugse.html

C. Wieners: Einführung in die Numerische Mathematik für Studierende der Fachrichtung Informatik und Ingenieurwissenschaften 1


Am 25. Februar 1991 während des ersten Golfkriegs inDharan, Saudi Arabien, verfehlte eine amerikanischePatriot–Rakete eine anfliegende irakische Scud Raketedurch eine falsche Zeitberechnung.Eine 1/10 Sekunde wurde ungenau dargestellt (durchRundungsfehler wurde die periodische Dualentwicklung

0.0001100110011001100110011001100....

in der Computerdarstellung zu

0.00011001100110011001100

abgeschnitten), so dass nach 100 Stunden Betriebszeiteine Zeitdifferenz von ca. 0.3 Sekunden entstand.Dieser Fehler wurde nicht in allen Teilen desBetriebsprogramms korrigiert.

http://ta.twi.tudelft.nl/users/vuik/wi211/disasters.html

http://www5.in.tum.de/ huckle/bugse.html



(1.3) a) Eine Gleitkommazahlen zur Basis B ∈ {2,3, ...} der Mantissenlänge M undExponentenlänge E ist die Menge

FL ={±Be

M

∑m=1

amB−m : e = e−+E−1

∑k=0

ck Bk , am,ck ∈ {0,1, ...,B−1}}

b) Eine Gleitkommaarithmetik wird durch eine Abbildung fl : R−→ FL mit fl(x) = x für x ∈ FLdefiniert. Sei bestimmt die Rundung: x ⊕y = fl(x + y), x �y = fl(x ·y), etc.Die zugehörige Maschinengenauigkeit ist

eps := sup{|x −fl(x)||x |

; 0 < |x |< M}.

(1.2) a) Ein Problem heißt sachgemäß gestellt, wenn es eindeutig lösbar ist und die Lösung stetigvon den Daten abhängt.b) Die Kondition eines Problems ist eine Maß dafür, wie stark die Abhängigkeit der Lösungvon den Daten ist.c) Die Stabilität eines numerischen Algorithmus ist eine Maß dafür, wie stark dieDaten-Abhängigkeit der numerischen Lösung im Vergeich zu der tatsächlichen Lösung ist.

(1.3) Sei f : RN −→ RK eine differenzierbare Funktion und x ∈ RN . Dann heißta) κkn

abs =∣∣∣ ∂

∂xnfk (x)

∣∣∣ absolute Konditionszahl.

b) κknrel =

∣∣∣ ∂

∂xnfk (x)

∣∣∣ |xn ||fk (x)| relative Konditionszahl.



Sei | · | eine Vektornorm, und sei ‖ · ‖ eine zugeordete Matrixnorm, d. h.,

|Ax | ≤ ‖A‖|x | , x ∈ RN , A ∈ RM×N .

Wir verwenden für x ∈ RN und A ∈ RM×N

|x |1 =N

∑n=1|xn| |x |2 =

√xT x |x |∞ = max

n=1,...,N|xn|

und die zugeordnete Operatornorm ‖A‖p = supx 6=0N

|Ax |p|x |p , d.h.

‖A‖1 = maxn=1,...,N

M

∑m=1|amn| , ‖A‖2 =

√ρ(AT A) , ‖A‖∞ = max

m=1,...,M

N

∑n=1|amn|

mit Spekralradius ρ(A) = max{|λ | : λ Eigenwert von A}.

(1.4) Sei A ∈ RN×N invertierbar. Dann heißt κp(A) = ‖A‖p‖A−1‖p die Kondition von A.

Sei b ∈ RN , b 6= 0N , 4b ∈ RN eine kleine Störung und b̃ = b +4b.Sei x ∈ RN Lösung von Ax = b, x̃ ∈ RN Lösung von Ax̃ = b̃, und 4x = x̃ −x der Fehler.

Dann gilt für den relativen Fehler|4x |p|x |p

≤ κp(A)|4b|p|b|p

.


Karlsruhe Institute of TechnologyAuslöschung bei Nullstellenberechnung

Wir betrachten die Gleichung

x2−2px + q = 0,

deren Nullstellen durch

x = p±√

p2−q

gegeben sind. Diese Berechnungsvorschrift ist aber für p� q nicht zu empfehlen, da dannAuslöschung bei der betragsmäßig kleineren Nullstelle auftritt. Wählt man beispielsweisep = 108 und q = 1, so berechnet Matlab

x1 = 2∗108, x2 = 0.

Die Auslöschung bei x2 kann man umgehen, indem man zuerst die größere Nullstelle durch

x1 = p + sign(p)√

p2−q

berechnet und dann (mit dem Satz von Vieta) die zweite Nullstelle durch

x2 =qx1

erhält. Mit dieser Vorschrift berechnet Matlab die bessere Lösung

x1 = 2∗108, x2 = 0.5∗10−9.


Karlsruhe Institute of Technology2 Direkte Lösungsverfahren für lineare Gleichungen

(2.1) Sei L ∈ RN×N eine normierte untere Dreiecksmatrix und b ∈ RN . Dann ist L invertierbar unddas Lineare Gleichungssystem (LGS) Ly = b ist mit O(N2) Operationen lösbar.Entsprechend ist für eine invertierbare obere Dreiecksmatrix R ∈ RN×N das LGS Rx = y inO(N2) Operationen lösbar.

(2.2) Wenn eine Matrix A ∈ RN×N eine LR-Zerlegung A = LR mit einer normierten untereDreiecksmatrix L und einer invertierbaren obere Dreiecksmatrix R besitzt, dann ist Ainvertierbar und das LGS Ax = b ist mit O(N2) Operationen lösbar.

(2.3) Eine Matrix A ∈ RN×N besitzt genau dann eine LR-Zerlegung von A, wenn alleHauptuntermatrizen A[1 : n,1 : n] invertierbar sind. Die LR-Zerlegung ist eindeutig und lässtsich mit O(N3) Operationen berechnen.

(2.5) Eine Matrix A ∈ RN×N heißt strikt diagonal-dominant, falls |A[n,n]|>N∑

k=1k 6=n

|A[n,k ]|.

Sie heißt positiv definit, wenn xT Ax > 0 für x ∈ RN , x 6= 0.In beiden Fällen existiert eine LR-Zerlegung.

(2.6) Sei A ∈ RN×N symmetrisch und positiv definit.Dann existiert genau eine Cholesky-Zerlegung A = LLT mit einer unteren Dreiecksmatrix L.

Software: http://www.netlib.org/lapackC. Wieners: Einführung in die Numerische Mathematik für Studierende der Fachrichtung Informatik und Ingenieurwissenschaften 6

Karlsruhe Institute of Technology2 Direkte Lösungsverfahren für lineare Gleichungen

(2.7) a) Eine bijektive Abbildung π : {1, ...,N} −→ {1, ...,N} heißt Permutation.b) Sie ist eindeutig durch einen Permutationsvektor p ∈ RN mit p[n] = π(n) bestimmt.c) Die zugehörige Permutationsmatrix ist P = (eπ(1) · · · |eπ(N)) ∈ RN×N .

(2.8) Sei A ∈ RN×N invertierbar. Dann existiert eine Permutationsmatrix P, so dass PA eineLR-Zerlegung PA = LR besitzt, so dass |L[m,n]| ≤ 1 gilt.

(2.9) a) Zu v ∈ RN und k 6= n mit v [k ]2 + v [n]2 > 0 existiert eine Givens-Rotation G ∈ RN×N mit(G[k ,k ] G[k ,n]G[n,k ] G[n,n]

)=

(c s−s c

), c2 + s2 = 1 ,

und G[j][j] = 1 für j 6= k ,n und G[i][j] = 0 sonst, so dass für w = Gv gilt: w [n] = 0.Für |v [n]|> |v [k ]| setze τ = v [k ]

v [n] , s = 1√1+τ2

, c = sτ, sonst setze τ = v [n]v [k ] , c = 1√

1+τ2, s = cτ.

b) Zu v ∈ RN , v 6= 0, existiert eine Householder-Spiegelung H = IN −2

wT wwwT ∈ RN×N mit

w ∈ RN , w [1] = 1, sodass Hv = σe1 mit σ ∈ R.Falls v [1] > 0, setze σ =−|v |2, sonst setze σ = |v |2. Dann definierte w = 1

v [1]−σ(v −σe1).

Rotationen und Spiegelungen Q sind orthogonale Matrizen, d.h. QT Q = IN , Q−1 = QT ,‖Q‖2 = 1 und κ2(Q) = 1.

(2.10) Zu A ∈ RK×N existiert eine QR-Zerlegung A = QR mit einer orthogonalen Matrix Q ∈ RK×K

und eine oberen Dreiecksmatrixmatrix R ∈ RM×N .C. Wieners: Einführung in die Numerische Mathematik für Studierende der Fachrichtung Informatik und Ingenieurwissenschaften 7

Karlsruhe Institute of TechnologyLR-Zerlegung (ohne Vektorisierung)

N = size(A,1);

for n=1:N-1% Berechnung der n-ten Spalte von Lfor m=n+1:NA(m,n) = A(m,n)/A(n,n);

end;% keine Berechnung der n-ten Zeile von R erforderlich

% Berechnung der Restmatrixfor m=n+1:Nfor k=n+1:NA(m,k) = A(m,k) - A(m,n) * A(n,k);

end;end;

end;


Karlsruhe Institute of TechnologyLR-Zerlegung

N = size(A,1);

for n=1:N-1A(n+1:N,n) = A(n+1:N,n)/A(n,n);A(n+1:N,n+1:N) = A(n+1:N,n+1:N) - A(n+1:N,n) * A(n,n+1:N);

end;

x = b;for n=2:N

x(n) = x(n) - A(n,1:n-1) * x(1:n-1);end;

for n=N:-1:1x(n) = (x(n) - A(n,n+1:N)*x(n+1:N))/A(n,n);

end;


Karlsruhe Institute of TechnologyCholesky-Zerlegung

N = size(A,1);

for n=1:NA(n:N,n) = A(n:N,n) - A(n:N,1:n-1) * A(n,1:n-1)’;A(n:N,n) = A(n:N,n) / sqrt(A(n,n));

end;

x = b;for n=1:Nx(n) = (x(n) - A(n,1:n-1) * x(1:n-1))/ A(n,n);

end;

for n=N:-1:1x(n) = (x(n) - A(n+1:N,n)’ * x(n+1:N))/ A(n,n);

end;


Karlsruhe Institute of TechnologyLR-Zerlegung mit Pivotsuche

N = size(A,1); p = (1:N)’;for n = 1:N-1

[r,m] = max(abs(A(n:N,n)));m = m+n-1;if abs(A(m,n))<eps

error(’*** ERROR *** LR-Zerlegung existiert nicht’);end;if (m ~= n)

A([n m],:) = A([m n],:); p([n m]) = p([m n]);end;A(n+1:N,n) = A(n+1:N,n)/A(n,n);A(n+1:N,n+1:N) = A(n+1:N,n+1:N) - A(n+1:N,n)*A(n,n+1:N);

end;x = b(p);for n=2:N

x(n) = x(n) - A(n,1:n-1)*x(1:n-1);end;for n=N:-1:1

x(n) = (x(n) - A(n,n+1:N)*x(n+1:N))/A(n,n);end;


Karlsruhe Institute of TechnologyBerechnung der Householder-Vektoren

function [v,beta] = householder(y)N = length(y);s = y(2:N)’ * y(2:N);if N == 1s = 0;

end;v = [1;y(2:N)];if s == 0beta = 0;

elsemu = sqrt(y(1)^2 + s);if y(1) <= 0v(1) = y(1) - mu;

elsev(1) = -s/(y(1) + mu);

end;beta = 2*v(1)^2/(s + v(1)^2);v = v / v(1);

end;return;


Karlsruhe Institute of TechnologyQR-Zerlegung

[M,N] = size(A);for m = 1:min(N,M-1)

[v,beta] = householder(A(m:M,m));if beta ~= 0w = beta * v’ * A(m:M,m:N);A(m:M,m:N) = A(m:M,m:N) - v * w;A(m+1:M,m) = v(2:M-m+1);

end;end;

for m = 1:min(N,M-1)v = [1;A(m+1:M,m)];beta = 2 / (v’ * v);if beta ~= 2

b(m:M) = b(m:M) - beta*(v’*b(m:M)) * v;end;

end;for n=min(N,M):-1:1x(n) = (b(n) - A(n,n+1:N) * x(n+1:N)) / A(n,n);

end;


Karlsruhe Institute of Technology3 Lineare Ausgleichsrechnung

(3.1) Sei A ∈ RK×N und b ∈ RK . Dann gilt:x ∈ RN minimiert |Ax −b|2 ⇐⇒ AT Ax = AT b.

(3.2) Zu A ∈ RK×N mit R = rang(A) existieren Singulärwerte σ1, ...,σR > 0 und eineSingulärwertzerlegung

A = V ΣUT

mit V ∈ RK×K , U ∈ RN×N orthogonal und Σ ∈ RK×N mit Σ[r , r ] = σr für r = 1, ...,R undΣ[k ,n] = 0 sonst.

(3.3) A+ = UΣ+V T ist die Pseudo-Inversemit Σ+ ∈ RN×K mit Σ+[r , r ] = 1/σr für r = 1, ...,R und Σ+[n,k ] = 0 sonst.

(3.4) x = A+b löst die Normalengleichung AT Ax = AT b.

(3.5) Sei A ∈ RK×N und b ∈ RK .Dann gilt für die Tikhonov-Regularisierung mit α > 0:x ∈ RN minimiert |Ax −b|22 + α |x |22 ⇐⇒ (AT A + αIN )x = AT b.

(3.6) Es gilt limα−→0

(AT A + αIN )−1AT b = A+b.


Karlsruhe Institute of TechnologyEin schlecht konditioniertes Gleichungssystem

Wir betrachten das lineare Gleichungssystem Ax = b mit der Hilbertmatrix

A =( 1

m + n + 1

)m,n=0,...,N

∈ RN+1×N+1

und der rechten Seite b =(

(−1)n(

log(2) +n

∑m=1

(−1)m))

n=0,...,N∈ RN+1.

Die exakte Lösung lautet:

N x0 x1 x2 x3 x4 x5 x61 0.93 −0.482 0.99 −0.80 0.333 1.00 −0.94 0.66 −0.224 1.00 −0.98 0.86 −0.53 0.155 1.00 −1.00 0.95 −0.77 0.42 −0.116 1.00 −1.00 0.98 −0.90 0.67 −0.32 0.07

(Beispiel aus Kress: Numerical Analysis)



Stören wir nun aber die rechte Seite geringfügig, indem wir log(2) nur bis auf5 Nachkommastellen auswerten, so erhalten wir folgende Lösung:

N x0 x1 x2 x3 x4 x5 x61 0.93 −0.482 0.99 −0.81 0.333 1.00 −0.96 0.70 −0.254 1.01 −1.16 1.63 −1.70 0.725 1.06 −2.74 12.68 −31.16 33.87 −13.266 1.39 −16.58 151.10 −584.81 1071.96 −926.77 304.50

Also: Eine geringfügige Störung der Daten führt zu einer großen Störung des Ergebnisses.Der Grund dafür liegt in der schlechten Kondition der Hilbertmatrix. Diese ist in derSpektralnorm:

N 2 3 4 5 6κ2(A) 19.28 524.06 1.55e + 04 4.77e + 05 1.495e + 07

(Beispiel aus Kress: Numerical Analysis)


Karlsruhe Institute of TechnologyInvertierung der Hilbert-Matrix (Matlab)

>> H = hilb(3); IH = inv(H); IH(1,1:3)

ans = 9.00000000000003 -36.0000000000001 30.0000000000001

>> H = hilb(5); IH = inv(H); IH(1,1:3)

ans = 24.9999999999919 -299.999999999882 1049.99999999958

>> H = hilb(7); IH = inv(H); IH(1,1:3)

ans = 49.0000000578711 -1176.00000232576 8820.00002248178

>> H = hilb(9); IH = inv(H); IH(1,1:3)

ans = 80.9999332549633 -3239.99529943927 41579.919225348

>> H = hilb(11); IH = inv(H); IH(1,1:3)

ans = 120.91751059331 -7251.2942628623 141342.563141014



Wir können die Kondition verbessern, indem wir die Tikhonov-Regularisierung auf dieHilbertmatrix anwenden, d.h.

xα = (AT A + αIN )−1AT b .

Regularisieren wir mit einem Parameter α = 10−10, so erhalten wir

N x0 x1 x2 x3 x4 x5 x61 0.93 −0.482 0.99 −0.81 0.333 1.00 −0.95 0.69 −0.244 0.99 −0.89 0.47 0.06 −0.145 1.00 −0.91 0.52 0.02 −0.18 0.046 1.00 −0.94 0.58 0.08 −0.25 −0.17 0.20


Karlsruhe Institute of Technology4 Eigenwertberechung

(4.1) Eine Matrix H ∈ RN×N heißt Hessenberg-Matrix, wenn H[n + 2 : N,n] = 0N−n−1 fürn = 1, ...,N−2.

(4.2) Sei A ∈ RN×N . Dann existiert eine orthogonale Matrix Q ∈ RN×N , so dass H = QAQT eineHessenberg-Matrix ist. Die Berechnung benötigt O(N3) Operationen.Wenn A symmetrisch ist, dann ist H eine Tridiagonalmatrix.

(4.3) Sei A ∈ RN×N symmetrisch, tridiagonal, und irreduzibel, d.h.A[n + 2 : N,n] = A[n,n + 2 : N]T = 0N−n−1 und A[n−1,n] = A[n,n−1] 6= 0.Dann hat A paarweise verschiedene reele Eigenwerte λ1 < λ2 < · · ·< λN .

(4.4) Inverse Iteration mit variablem Shift

S0) Wähle z0 ∈ RN , z0 6= 0N , ε ≥ 0. Setze k = 0.S1) Setze wk = 1

|zk |2zk , µk = r(A,wk ).

S2) Falls |Awk −µk wk |2 ≤ ε STOP.S3) Berechne zk+1 = (A−µk IN )−1wk .S4) Setze k := k + 1, gehe zu S1).

Wenn der Startvektor z0 hinreichend nahe bei einem Eigenvektor vm mit isoliertemEigenwert λ = λm liegt, konvergiert die Iteration kubisch, d.h. |µk −λ | ≤ C |µk −λ |3 .


Karlsruhe Institute of Technology4 Eigenwertberechnung: QR-Iteration mit Shift

Sei A ∈ RN×N symmetrisch.S0) Berechne A0 = QAQT tridiagonal (Hessenberg-Transformation).Wähle ε ≥ 0. Setze k = 0.S1) Falls |Ak [n + 1,n]| ≤ ε für ein n:getrennte Eigenwertberechnung für Ak [1 : n,1 : n] und Ak [n + 1 : N,n + 1 : N].S2) Berechne dk = 1

2 (Ak [N−1,N−1]−Ak [N,N]) und

sk = Ak [N,N] + dk −sgn(dk )√

d2k + Ak [N−1,N]2 .

S3) Berechne QR-Zerlegung

Qk Rk = Ak −sk IN

und setze

Ak+1 = Rk Qk + sk IN .

S4) Setze k := k + 1, gehe zu S1).

Es gilt Ak+1 = QTk Ak Qk . Also haben Ak+1 und Ak disselben Eigenwerte.

Falls der shift sk = Ak [N,N] gewählt wird, entspricht die QR-Iteration der Inversen Iterationmit variablem Shift und Startvektor z0 = eN .


Karlsruhe Institute of TechnologyImpliziter QR-Algorithmus mit Wilkinson-Shift

Für symmetrische Matrizen (o.B.d.A. Hessenbergform) lässt sich der QR-Algorithmus somodifizieren, dass die QR-Zerlegung in jedem Iterationsschritt nur implizit durch N−1Transformationen mit Givensrotationen durchgeführt werden muss. Durch denWilkinson-Shift (siehe Golub, van Loan) erhalten wir außerdem kubische Konvergenz, imGegensatz zur linearen Konvergenz des QR-Algorithmus’ ohne Shift. Dabei können wir eineToleranz für die verschwindenden Nebendiagonalelemente vorgeben.

Wir betrachten

A = tridiag(−1,2,−1) ∈ RN×N .

Für die Toleranz tol = ε braucht der Algorithmus die folgende Anzahl von Iterationen:

N Iterationen100 281200 532500 1120

1000 2310

Das sind im Schnitt weniger als 3 Iterationen pro Eigenwert.Der maximale Fehler bei diesen Berechnungen beträgt zwischen 10−13 und 10−14.


Karlsruhe Institute of Technology5 Iterationsverfahren für lineare Gleichungen

(5.1) S0) Wähle x0 ∈ RN und ε > 0.Berechne r0 = b−Ax0. Setze k = 0.

S1) Falls |r0| ≤ ε STOP.

S2) Berechne

wk = Brk−1

xk+1 = xk + wk

r k+1 = rk −Awk

Setze k := k + 1 und gehe zu S1).

Sei A = L + D + R. Dann gilt B = D−1 für das Jacobi-Verfahren und B =(L + D

)−1 für dasGauß-Seidel-Verfahren.

(5.2) Sei A,B ∈ RN×N mit ρ(IN −BA) < 1. Dann ist A invertierbar, und es gilt für alle b ∈ RN undalle Startvektoren x0 ∈ RN konvergiert die Iteration

xk+1 = xk + B(b−Axk ) , k = 0,1,2, ...

gegen limk−→∞

xk = A−1b. Dann exisitiert eine Vektor-Norm | · | und dazu eine Matrix-Norm ‖ · ‖

mit ‖IN −BA‖< 1. Damit ergibt sich |x −xk | ≤ ‖IN −BA‖k |x −x0| (lineare Konvergenz).

Software: z.B. http://www.cise.ufl.edu/research/sparse/SuiteSparseC. Wieners: Einführung in die Numerische Mathematik für Studierende der Fachrichtung Informatik und Ingenieurwissenschaften 22

Karlsruhe Institute of Technology5 Iterative Lösungsverfahren: Krylov-Verfahren

Sei 〈·, ·〉V ein Skalarprodukt in V = RN .

S0) Wähle x0 ∈ RN .Berechne r0 = b−Ax0, z1 = Br0, h10 = |z1|V und v1 = 1

h10z1. Setze k = 1.

S1) Berechne

wk = BAvk

zk+1 = wk −k

∑j=1

hjk v j mit hjk = 〈v j ,wk 〉V

vk+1 =1

hk+1,kzk+1 mit hk+1,k = |zk+1|V

S2) Setzte k := k + 1 und gehe zu S1).

Dann ist v1, ...,vk eine Orthonormalbasis von dem Krylov-Raum

Vk = span{Br0,BABr0, ...,(BA)k−1Br0}= {Qk y : y ∈ Rk} , Qk =(v1|....|vk ) .

Es gilt BAvk =k+1∑

j=1hjk v j , also BAQk = Qk+1Hk mit Hk = (hjm) ∈ Rk+1,k .

GMRES-Verfahren: Wähle 〈v ,w〉V = vT w .cg-Verfahren (A,B symmetrisch positiv definit): Wähle 〈v ,w〉V = vT Aw .


Karlsruhe Institute of Technology5 Iterative Lösungsverfahren: GMRES-Verfahren

S0) Wähle x0 ∈ RN , ε > 0.Berechne r0 = b−Ax0, z1 = Br0, h10 = |z1|2 und v1 = 1

h10z1. Setze k = 1.

S1) Berechne

wk = BAvk

zk+1 = wk −k

∑j=1

hjk v j mit hjk = (v j )T wk

vk+1 =1

hk+1,kzk+1 mit hk+1,k = |zk+1|2

S2) Berechne yk ∈ Rk mit ρk = |Hk yk −h10e1|2 = min!Dabei ist Hk = (hjm)j=1,...,k+1, m=1,...,k ∈ Rk+1,k .

S3) Wenn ρk < ε, setze xk = x0 +k∑

j=1yk

j v j STOP.

S4) Setze k := k + 1 und gehe zu S1).

(4.4) Es gilt ρk = minz∈x0+Vk

|B(b−Az)|2.


Karlsruhe Institute of Technology5 Iterative Lösungsverfahren: cg-Verfahren

Seien A,B ∈ RN×N symmetrisch und positiv definit.S0) Wähle x0 ∈ RN , ε > 0.Berechne r0 = b−Ax0, w0 = Br0, ρ0 = (w0)T r0 und d1 = w0. Setze k = 0.S1) Falls ρk ≤ ε STOPS2) Setze k := k + 1 und berechne

uk = Adk

αk =ρk−1

(uk )T dk

xk = xk−1 + αk dk

rk = rk−1−αk uk

wk = Brk

ρk = (wk )T rk

dk+1 = wk +ρk

ρk−1dk

Gehe zu S1).

(4.5) Es gilt |xk −x |A ≤ 2(√

κ(BA)−1√κ(BA) + 1

)k|x0−x |A.


Karlsruhe Institute of Technology6 Iterationsverfahren für nichtlineare Gleichungen

Sei F : RN −→ RN differenzierbar. Sei x∗ ∈ RN eine Nullstelle von F (·), d.h. F (x∗) = 0.Dann gilt 0 = F (x) + DF (x)(x∗−x) + o(x∗−x).Falls DF (x) invertierbar ist, gilt x∗ = x −DF (x)−1F (x) + o(x∗−x).Falls DF (x∗) invertierbar ist, ist x∗ Fixpunkt von ΦF (x) = x−DF (x)−1F (x), d.h. ΦF (x∗) = x∗.

Newton-Verfahren: Wähle x0 ∈ RN und definiere xk+1 = ΦF (xk ) für k = 0,1,2, ...

(6.1) Sei DF (x∗) invertierbar. Dann ist das Newton-Verfahren für alle x0 hinreichend nahe bei x∗

wohldefiniert, und xk konvergiert gegen x∗. Wenn zusätzlich F (·) glatt genug ist, ist dieKonvergenz quadratisch, d.h. es existiert C > 0 mit |xk+1−x∗| ≤ C |xk −x∗|2.

Wenn DF (x)−1 durch B ∈ RN×N approximiert wird, erhalten wir eine einfacheFixpunktiteration mit ΦF ,B(x) = x −BF (x).

(6.2) Sei ρ(IN −B DF (x∗)) < 1. Dann gilt: Für alle x0 hinreichend nahe bei x∗ konvergiert dieeinfache Fixpunkt-Iteration xk+1 = ΦF ,B(xk ) linear gegen x∗.

(6.3) Iterationsverfahren für nichtlineare Gleichungen mit DämpfungsstrategieS0) Wähle x0 ∈ RN , θ ∈ (0,1), smax ∈ N und ε > 0. Setze k = 0.

S1) Falls |F (xk )| ≤ ε STOP (Konvergenz)

S2) Wähle Bk ≈ DF (xk )−1 und berechne ck =−Bk F (xk ).

S3) Wähle tk ∈ {1,θ ,θ 2, ...,θ smax} mit |F (xk + tk ck )|< |F (xk )|.Falls |F (xk + tk ck )| ≥ |F (xk )| für tk = θ smax STOP (keine Konvergenz)

S4) Setze xk+1 = xk + tk ck , k := k + 1 und gehe zu S1).C. Wieners: Einführung in die Numerische Mathematik für Studierende der Fachrichtung Informatik und Ingenieurwissenschaften 26

Karlsruhe Institute of Technology7 Polynom-Interpolation

(7.1) Lagrange-Interpolation Zu Stützstellen t0 < t1 < · · ·< tN und Werten f0, f1, . . . , fN ∈ Rexistiert genau ein Polynom P ∈ PN mit P(tn) = fn .

Konstruktion: Definiere die Lagrange-Basis Ln(t) =N∏

k=0, k 6=n

t−tktn−tk

, setze P(t) =N∑

n=0fnLn(t).

(7.2) Hermite-Interpolation Zu t0 ≤ t1 ≤ ·· · ≤ tN definiere dn = max{n−k : tk = · · ·= tn}.

Zu Werten f0, f1, . . . , fN ∈ R existiert genau ein Polynom P ∈ PN mit 1dn!

(ddt

)dnP(tn) = fn.

Konstruktion: Zu t0 ≤ t1 ≤ ·· · ≤ tN definiere die Newton-Basis durch

ω0 ≡ 1, ω1(t) = t− t0, ωk (t) = (t− tk−1)ωk−1(t) ∈ Pk , k = 1, . . . ,N .

Nun definiere f [tn] = fn und f [tk , . . . , tn] = fn falls tk = tk+1 = · · ·= tn.Dann berechne rekursiv f [tk , . . . , tn] =

f [tk+1,...,tn ]−f [tk ,...,tn−1]tn−tk

.

Dann ist P(t) =N∑

k=0f [t0, . . . , tk ]ωk (t) das Interpolationspolynom.

Neville-Schema: Sei P0 ∈ Pn−k−1 Interpolation zu fk , ..., fn−1 an tk ≤ ·· · ≤ tn−1, und seiP1 ∈ Pn−k−1 Interpolation zu fk+1, ..., fn an tk+1 ≤ ·· · ≤ tn. Dann istP(t) = t−tk

tn−tkP1(t) + tn−t

tn−tkP0(t) Interpolation zu fk , ..., fn an tk ≤ ·· · ≤ tn.

(7.3) Sei f ∈ CN+1(a,b), t , tn ∈ (a,b), und fn = 1dn!

(ddt

)dnf (tn). Dann gilt für den

Interpolationsfehler f (t)−PN (t) = 1(N+1)!

(ddt

)N+1f (τ) ωN+1(t) mit τ ∈ (a,b).


Karlsruhe Institute of TechnologyRunge-Phänomen bei InterpolationWir betrachten das Interpola-tionsproblem für äquidistantverteilte Stützstellen. Wirwürden erwarten, dass füreine genügend glatte Funktiondie Folge der Interpolations-polynome, gehörig zu immerkleinerem Stützstellenabstand,in der Maximumsnorm gegendie Funktion konvergiert. Dassdem nicht so ist zeigt zumBeispiel das sogenannte Runge-Phänomen für die Funktionf : [−5,5]→ R,

f (x) :=1

1 + x2 .

Der maximale Fehler wird hier sogar immer größer, je höher der Grad desInterpolationspolynoms ist. Für den Grad 10 ist das zugehörige Interpolationspolynomeingezeichnet. Der Fehler wird offenbar durch die starke Oszillation verursacht, diesverschlimmert sich noch mit steigendem Polynomgrad.


Karlsruhe Institute of TechnologyDas Neville-Schema – Berechnung der Koeffizienten

function f = coeff(t,f)N = length(t);for k=1:N-1

for n=N:(-1):k+1if t(n) ~= t(n-k)f(n) = (f(n) - f(n-1)) / (t(n) - t(n-k));

endend

endreturn

function y = eval_newton(t,b,x)N = length(t);y = b(1);p = 1;for n=2:Np = p * (x - t(n-1));y = y + b(n) * p;

endreturn


Karlsruhe Institute of TechnologyNeville-Schema – Auswertung des Polynoms

function y = eval_neville(t,f,x)N = length(t);for k=1:N-1

for n=N:(-1):k+1if t(n) ~= t(n-k)f(n) = ((x-t(n-k))*f(n)+(t(n)-x)*f(n-1)) / (t(n)-t(n-k));

endend

endy = f(N);return


Karlsruhe Institute of Technology8 Spline-Interpolation

(8.1) a) Zu einer Zerlegung ∆: a = t0 < t1 < · · ·< tN = b von [a,b] definiere den kubischenSpline-Raum

S3(∆) ={

S ∈ C2[a,b] : Sn := S|[tn−1 ,tn ] ∈ P3, n = 1, . . . ,N}.

b) S ∈S3(∆) heißt interpolierender Spline zu f ∈ C0[a,b] wenn S(tn) = f (tn).

(8.2) Es gilt dimS3(∆) = 3 + N. Mit einer der Randbedingungen(I) Natürliche Randbedingung S′′(a) = 0 und S′′(b) = 0(II) Hermite-Randbedingung zu f ∈ C1[a,b] S′(a) = f ′(a) und S′(b) = f ′(b)(III) Periodische Randbedingungen S′(a) = S′(b) und S′′(a) = S′′(b)

ist die Spline-Interpolation S ∈ S3(∆) zu f eindeutig lösbar.

(8.3) Die Momente Mn = S′′(tn) von S ∈S3(∆) zu f sind eindeutig durch

µnMn−1 + Mn + λnMn+1 = 3f [tn−1, tn, tn+1]

und µn =hn

2(hn + hn+1), λn =

hn+1

2(hn + hn+1)bestimmt und es gilt

Sn(t) =Mn(t− tn−1)3 + Mn−1(tn− t)3

6hn+

f (tn) + f (tn−1)

2− h2

n12

(Mn + Mn−1)

+( f (tn)− f (tn−1)

hn− hn

6(Mn−Mn−1)

)(t− tn + tn−1

2

).

(8.4) Sei f ∈ C4[a,b]. Dann gilt für die Spline-Interpolation ‖S− f‖∞ ≤ h4 ‖f ′′′′‖∞.C. Wieners: Einführung in die Numerische Mathematik für Studierende der Fachrichtung Informatik und Ingenieurwissenschaften 31

Karlsruhe Institute of TechnologySpline-Interpolation mit Matlab

n = 7;x = rand(n,1);y = rand(n,1);plot(x,y,’.’)axis([-0.1 1.1 -0.1 1.1]);t = 1:n;ts = 1:1/10:n;xs = spline(t,x,ts);ys = spline(t,y,ts);hold on;plot(xs,ys,’r’);hold off;


Karlsruhe Institute of Technology9 Trigonometrische Interpolation und FFT

(9.1) Sei f : R−→ C eine 2π-periodische stetige Funktion und N = 2P . Sei tn = n2π

Nund fn = f (tn).

Dann existiert genau ein trigonometrisches Interpolationspolynom

TN (t) =N−1

∑n=0

cn exp(int) , cn ∈ C

mit TN (tn) = fn für n = 0, ...,N−1. Es gilt

ck =1N

N−1

∑n=0

fnω−kn , ω = exp(i2π/N) .

Die Koeffizienten lassen sich rekursiv mit O(N logN) Operationen berechnen:

N−1

∑n=0

fnω−kn =

N/2−1

∑m=0

f2mξ−km + ω

−1N/2−1

∑m=0

f2m+1ξ−km , ξ = ω

2 .

Für s-mal stetig differenzierbare Funktionen konvergieren die Koeffizienten cn gegen dieFourier-Koeffizienten

f̂k =1

2π

∫ 2π

0f (t)exp(ikt)dt ,

und es existiert eine Konstante Cs > 0 mit

‖f −TN‖L2(0,2π) ≤ CsN−s‖( d

dt

)sf‖L2(0,2π) , ‖g‖L2(0,2π) =

(∫ 2π

0|g(t)|2 dt

)1/2

.


Karlsruhe Institute of TechnologyFFT mit Matlab

Fs = 1000; % Sampling frequencyT = 1/Fs; % Sample timeL = 1000; % Length of signalt = (0:L-1)*T; % Time vector% Sum of a 50 Hz sinusoid and% a 120 Hz sinusoidx = 0.7*sin(2*pi*50*t)+sin(2*pi*120*t);y = x+2*randn(size(t)); % plus noiseplot(Fs*t(1:50),y(1:50))title(’Signal Corrupted Random Noise’)xlabel(’time (milliseconds)’)

NFFT = 2^nextpow2(L);% Next power of 2 from length of y

Y = fft(y,NFFT)/L;f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.plot(f,2*abs(Y(1:NFFT/2+1)))title(’Single-Sided Amplitude Spectrum of y(t)’)xlabel(’Frequency (Hz)’)ylabel(’|Y(f)|’)


Karlsruhe Institute of Technology10 Numerische Integration

(10.1) Sei [a,b]⊂ R ein Intervall, und sei Ξ⊂ [a,b] eine endliche Menge mit #Ξ = N.Dann existiert genau eine Quadratur

QΞ : C[a,b]−→ R , QΞ(f ) = ∑ξ∈Ξ

ωξ f (ξ )

mit Gewichten ωξ ∈ R zu den Stützstellen ξ ∈ Ξ, die für Polynome vom Grad N−1 exakt ist:

QΞ(P) =∫ b

aP(t)dt , P ∈ PN−1 .

Konstruktion: Sei Lξ (t) = ∏η∈Ξ\{ξ}t−η

ξ−η. Dann gilt ωξ =

∫ ba Lξ (t)dt .

(10.2) Es gibt keine Quadratur QΞ mit #Ξ = N, die für Polynome P ∈ P2N exakt ist.Es gibt genau eine Quadratur GN mit N Stützstellen, die für Polynome P ∈ P2N−1 exakt ist.Wähle dazu die Nullstellen des Orthogonalpolynoms LN ∈ PN bzgl. (f ,g) =

∫ ba f (t)g(t)dt .

(10.3) Sei QΞ eine Quadratur, die für Polynome P ∈ PK−1 exakt ist. Dann existiert C > 0 mit∣∣∣QΞ(f )−∫ b

af (t)dt

∣∣∣≤ C∥∥∥( d

dt

)Kf∥∥∥

∞

(10.4) Für die summierte Trapezregel

TM (f ) = h(1

2f (a) +

M−1

∑m=1

f (a + mh) +12

f (b))

mit h =b−a

M

gilt∣∣TM (f )−

∫ ba f (t)dt

∣∣≤ C h2 ‖f ′′‖∞.C. Wieners: Einführung in die Numerische Mathematik für Studierende der Fachrichtung Informatik und Ingenieurwissenschaften 35

Karlsruhe Institute of Technology10 Extropolierte Trapezsummen (Romberg)

Die Trapezsummen lassen sich extrapolieren: Setze Tk0 = T2k M (f ) und definiere rekursiv

Tki = Tk ,i−1 +1

4i −1(Tk ,i−1−Tk−1,i−1), i = 1, . . . ,m und k = i , . . . ,m .

Dann ist Tmk ist für Polynome P ∈ P2k+1 exakt.

Beispiel: Approximiere∫ 1

0

11 + t2 dt =

π

4

Neville-Schema zur Extrapolation von Tk0 = T2k+4T00 = 0.785235403010347T01 = 0.785357473293744 T11 = 0.785398163388209T02 = 0.785387990871414 T12 = 0.785398163397304 T22 = 0.785398163397910T03 = 0.785395620265938 T13 = 0.785398163397446 T23 = 0.785398163397456 T33 = 0.785398163397449

Konvergenzabschätzung durch Vergleich von TrapezsummenT01−T00 = 0.00012207028T02−T01 = 0.00003051757 T12−T11 = 0.0000000000090942T03−T02 = 0.00000762939 T22−T12 = 0.0000000000001426

FehlerT00− π

4 =−0.00016276038T01− π

4 =−0.00004069010 T11− π

4 =−9.2390539e−12T02− π

4 =−0.00000101725 T12− π

4 =−1.4477308e−13 T22− π

4 = 4.6151971e−13T03− π

4 =−0.00000025431 T13− π

4 =−2.1094237e−15 T23− π

4 = 7.4384942e−15 T33− π

4 = 2.2204460e−16


Karlsruhe Institute of Technology11 Numerische Integration von Differentialgleichungen

(11.1) Sei t0 ∈ R und T > 0. Zu einem Anfangswert u0 ∈ RM und f ∈ C([t0, t0 + T ]×RM ,RM ) suchenwir eine Lösung u ∈ C1([t0, t0 + T ],RM ) der Anfangswertaufgabe

u̇(t) = f(t ,u(t)

)für t ∈ [t0, t0 + T ] und dem Anfangswert u(t0) = u0 .

(11.2) Wenn L > 0 existiert mit |f (t ,w)− f (t ,z)| ≤ L |w −z|, dann gilt für zwei Lösungenu̇(t) = f

(t ,u(t)

)und v̇(t) = f

(t ,v(t)

)|u(t)−v(t)| ≤ exp(L(t− t0)) |u(t0−v(t0)| für t ≥ t0 .

Daher ist die Lösung der Anfangswertaufgabe eindeutig.

(11.3) Klassisches Runge-Kutta-Verfahren: Sei tn = t0 + nτ, u0 = u0 und setze

k1 = f (tn−1,un−1)

k2 = f (tn−1 +τ

2,un−1 +

τ

2k1)

k3 = f (tn−1 +τ

2,un−1 +

τ

2k2)

k4 = f (tn−1 + τ,un−1 + τk3)

un = un−10 +

τ

6(k1 + 2k2 + 2k3 + k4) .

Dann gilt |u(tn)−un| ≤ C exp(L(tn− t0)τ4.

(11.3) Implizites Euler-Verfahren: In jedem Schritt bestimme un als Lösung der nichtlinearenGleichung un = un−1 + τf (tn,un). Dann gilt |u(tn)−un| ≤ C exp(L(tn− t0)τ.


Karlsruhe Institute of Technology11 Simulation: Radioaktiver Zerfall u̇ =−λu

N exakt explizites Eulerverfahren implizites Eulerverfahren Mittelpunktsregelun = un−1 + τf (un−1) un = un−1 + τf (un) un = un−2 + 2τf (un−1)

4 0.50000 0.46711 0.52770 0.512668 0.50000 0.48431 0.51440 0.50323

16 0.50000 0.49233 0.50735 0.5008132 0.50000 0.49621 0.50371 0.5002064 0.50000 0.49811 0.50187 0.50005

O(1/N) O(1/N) O(1/N2)

Table: Vergleich der Konvergenzordnung |c(tn)−cn|= O(N−β ) = O(hβ

N ) im Zeitintervall [0,5730].

N = 5 N = 100Figure: Stabilität der numerischen Approximation. Vergleich im Zeitintervall [0,57300] für N = 5,100.


Karlsruhe Institute of TechnologyODE mit Matlab

function vdpode(MU)tspan = [0; max(20,3*MU)];y0 = [2; 0];options = odeset(’Jacobian’,@J);[t,y] = ode15s(@f,tspan,y0,options);figure;plot(t,y(:,1));title([’Solution of van der Pol

Equation, \mu = ’ num2str(MU)]);xlabel(’time t’);ylabel(’solution y_1’);axis([tspan(1) tspan(end) -2.5 2.5]);

function dydt = f(t,y)dydt = [ y(2)

MU*(1-y(1)^2)*y(2)-y(1) ];endfunction dfdy = J(t,y)dfdy = [ 0 1-2*MU*y(1)*y(2)-1 MU*(1-y(1)^2)];

endend


Karlsruhe Institute of Technology - math.kit.edu · Karlsruhe Institute of Technology 1...

Documents

Transcript of Karlsruhe Institute of Technology - math.kit.edu · Karlsruhe Institute of Technology 1...