1 Euklidische Approximation Sei V ein reeller euklidischer Vektorraum. Das Skalarprodukt in V wird mit h·, ·iV und die Norm mit k · kV bezeichnet. VN ⊂ V sei ein Teilraum der Dimension N < ∞ mit Basis {φn }n=1,...,N . (1.1)

Problemstellung: Sei v ∈ V . Bestimme vN ∈ V mit kv − vN kV = min kv − wN kV . wN ∈VN

(1.2)

Die Matrix   A = hφn , φk iV

n,k =1,...,N

∈ RN×N

ist symmetrisch und positiv definit. (1.3)

Problem (1.1) ist eindeutig lösbar. Es gilt N

vN =

∑ xn φn ,

n=1

wobei x ∗ ∈ RN die eindeutige Lösung des linearen Gleichungssystems Ax = b   mit b = hv , φk iV

k =1,...,N

∈ RN ist.

C. Wieners: Einführung in die Numerische Mathematik

1

2 Direkte Lösungsverfahren für lineare Gleichungen (2.1)

(2.2)

(2.4)

(2.5)

(2.6)

(2.7)

(2.8)

Sei L ∈ RN×N eine normierte untere Dreiecksmatrix und b ∈ RN , d.h. diag L = IN und L[1 : n, n + 1] = 0n für n = 1, ..., N − 1. Dann ist L regulär und das lineare Gleichungssystem Ly = b ist mit O(N 2 ) Operationen lösbar. Entsprechend ist für eine reguläre obere Dreiecksmatrix R ∈ RN×N (d.h. R[n, n] 6= 0 für alle n und R[n + 1, 1 : n] = 0Tn für n < N) das LGS Rx = y in O(N 2 ) Operationen lösbar. Die normierten unteren Dreiecksmatrizen bilden eine Gruppe. Die regulären oberen Dreiecksmatrizen bilden eine Gruppe. Wenn eine Matrix A ∈ RN×N eine LR-Zerlegung A = LR mit einer normierten untere Dreiecksmatrix L und einer regulären obere Dreiecksmatrix R besitzt, dann ist A regulär und das LGS Ax = b ist mit O(N 2 ) Operationen lösbar. Eine Matrix A ∈ RN×N besitzt genau dann eine LR-Zerlegung von A, wenn alle Hauptuntermatrizen A[1 : n, 1 : n] regulär sind. Die LR-Zerlegung ist eindeutig und lässt sich mit O(N 3 ) Operationen berechnen. n−1

N

k =1

k =n+1

Eine Matrix A ∈ RN×N heißt strikt diagonal-dominant, falls |ann | > ∑ |ank | +

∑

|ank |.

Wenn A strikt diagonal dominant ist, dann existiert eine LR-Zerlegung. Sei A ∈ RN×N symmetrisch und positiv definit. Dann existiert genau eine Cholesky-Zerlegung A = LLT mit einer regulären unteren Dreiecksmatrix L.

C. Wieners: Einführung in die Numerische Mathematik

2

2 Direkte Lösungsverfahren für lineare Gleichungen (2.9)


Sei π ∈ SN eine Permutation. Dann heißt Pπ = eπ −1 (1) |...|eπ −1 (N) ∈ RN×N Permutationsmatrix zu π. Es gilt (Pπ A)[n, k ] = A[π(n), k ] und (APπ )[n, k ] = A[n, π −1 (k )].

(2.10)

(2.11)

Die Permutationsmatrizen in RN×N bilden eine Gruppe. Es gilt Pσ Pπ = Pπ◦σ und Pπ −1 = PπT . Sei A ∈ RN×N regulär. Dann existiert eine Permutationsmatrix P, so dass PA eine LR-Zerlegung PA = LR besitzt und für die Einträge |L[m, n]| ≤ 1 gilt. Sie lässt sich mit O(N 3 ) Operationen berechnen. Die Lösung von Ax = b berechnet sich aus Ly = Pb und Rx = y . Sei | · | eine Vektornorm, und sei k · k eine zugeordete Matrixnorm, d. h., |Ax| ≤ kAk |x| ,

(2.12)

x ∈ RN , A ∈ RM×N .

Sei A ∈ RN×N regulär, und sei ∆A ∈ RN×N so klein, dass k∆Ak < kA−1 k−1 gilt. Dann ist die ˜ = A + ∆A regulär. Matrix A ˜ = b + ∆b. Sei b ∈ RN , b 6= 0N , ∆b ∈ RN klein und b ˜ x˜ = b ˜ Dann gilt für die Lösungen x ∈ RN von Ax = b und x˜ ∈ RN von A   |∆x| κ(A) |∆b| k∆Ak ≤ + . k∆Ak |x| |b| kAk 1 − κ(A) kAk Dabei ist ∆x = x˜ − x,

|∆x| |x|

der relative Fehler, und κ(A) = kAk kA−1 k die Kondition von A.

C. Wieners: Einführung in die Numerische Mathematik

3

2 Arithmetische Grundlagen (2.13)

a) Eine Gleitkommazahlen zur Basis B ∈ {2, 3, ...} der Mantissenlänge M und Exponentenlänge E ist die Menge n o M E−1 FL = ± B e ∑ am B −m : e = e− + ∑ ck B k , am , ck ∈ {0, 1, ..., B − 1} m=1

k =0

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

(2.14)

Sei f : RN −→ RK eine differenzierbare Funktion und x ∈ RN . Dann heißt kn = ∂ f (x) absolute Konditionszahl. a) κabs ∂ xn k kn = ∂ f (x) |xn | relative Konditionszahl. b) κrel ∂ xn k |f (x)| k

C. Wieners: Einführung in die Numerische Mathematik

4

2 Kondition und Stabilität (2.15)

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

(2.16)

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

|x|1 =

∑ |x[n]|

|x|2 =

p

xT x

|x|∞ = max |x[n]| n=1,...,N

n=1

und die zugeordnete Operatornorm kAkp = supx6=0N M

kAk1 = max

∑

n=1,...,N m=1

|A[m, n]| ,

kAk2 =

|Ax|p |x|p

q ρ(AT A) ,

, d.h. N

kAk∞ = max

∑ |A[m, n]|

m=1,...,M n=1

mit Spekralradius ρ(A) = max{|λ | : λ ∈ σ (A)} und Spektrum σ (A).

C. Wieners: Einführung in die Numerische Mathematik

5

3 Ausgleichsrechnung Sei Q ∈ RN×N orthogonal, d.h. Q T Q = IN . Dann ist κ2 (Q) = 1. (3.1)

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] c s = , c 2 + s2 = 1 , G[n, k ] G[n, n] −s c und G[j][j] = 1 für j 6= k , n und G[i][j] = 0 sonst, so dass enT Gv = 0 gilt: ] [n] √ 1 , c = sτ, sonst setze τ = vv [k √1 Für |v [n]| > |v [k ]| setze τ = vv [k [n] , s = ], c = 2

1+τ 2

1+τ

Es gilt G = c(ek ekT + en enT ) + s(ek enT − en ekT ) +

∑

, s = cτ.

ej ejT .

j6=k ,n (2.10)

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

2 ww T ∈ RN×N mit wT w

w ∈ RN , w[1] = 1, so dass Hv = σ e1 mit σ ∈ R und Hw = −w gilt: Falls v [1] > 0, setze σ = −|v |2 , sonst setze σ = |v |2 . Dann definierte w = (3.3)

(3.4)

1 v [1]−σ

(v − σ e1 ).

Zu A ∈ RK ×N existiert eine QR-Zerlegung A = QR mit einer orthogonalen Matrix Q ∈ RK ×K und eine oberen Dreiecksmatrix R ∈ RK ×N , d.h. QQ T = IK und R[k + 1 : K , k ] = 0K −k für k = 1, ..., K . Sei A ∈ RK ×N und b ∈ RK . Dann gilt: x ∈ RN minimiert |Ax − b|2

C. Wieners: Einführung in die Numerische Mathematik

⇐⇒

AT Ax = AT b . 6

Berechnung der Householder-Vektoren function [v,beta] = householder(y) N = length(y); s = y(2:N)’ * y(2:N); if N == 1 s = 0; end; v = [1;y(2:N)]; if s == 0 beta = 0; else mu = sqrt(y(1)^2 + s); if y(1) 0. R

R

r =1

r =1

Es gilt A = ∑ σr ur vrT und Ax = ∑ σr (vrT x)ur .

(2.22)

(2.22)

(2.23)

R

r =1

r =1

x = A+ b löst die Normalengleichung AT Ax = AT b. Sei A ∈ RK ×N und b ∈ RM . Dann gilt für die Tikhonov-Regularisierung mit α > 0: x ∈ RN minimiert |Ax − b|22 + α |x|22

(2.24)

R

A+ = UΣ−1 V T ist die Pseudo-Inverse. Es gilt A+ = ∑ σr−1 vr urT und Ax = ∑ σr−1 (urT x)vr .

⇐⇒

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

Es gilt lim (AT A + αIN )−1 AT b = A+ b. α−→0

Diskrepanzprinzip: Sei b ∈ Bild(A), x = A+ b und bδ eine Störung mit |b − bδ | < δ < |b|2 . Dann existiert ein α = α(δ ) > 0 mit |Ax α − bδ |2 = δ für x α = (AT A + αIN )−1 AT bδ . Es gilt α(δ ) −→ 0 für δ −→ 0. C. Wieners: Einführung in die Numerische Mathematik

9

4 Eigenwertberechnung (4.3)

(4.4)

H ∈ RN×N heißt Hessenberg-Matrix, wenn H[n + 2 : N, n] = 0N−n−1 für n = 1, ..., N − 2. Sei A ∈ RN×N . Dann existiert eine orthogonale Matrix Q ∈ RN×N , so dass H = QAQ T eine Hessenberg-Matrix ist. Die Berechnung benötigt O(N 3 ) Operationen. Wenn A symmetrisch ist, dann ist H eine Tridiagonalmatrix.

(4.8)

Sei A ∈ RN×N symmetrisch, tridiagonal, und irreduzibel, d.h. A[n − 1, n] = A[n, n − 1] 6= 0 und A[n + 2 : N, n] = A[n, n + 2 : N]T = 0N−n−1 .
Die charakteristischen Polynome Pn (t) = det A[1 : n, 1 : n] − tIn der Hauptuntermatrizen lassen sich durch eine Dreitermrekursion berechnen: Setze P0 ≡ 1. Dann gilt P1 (t) = A[1, 1] − t und Pn (t) = (A[n, n] − t)Pn−1 (t) − A[n − 1, n]2 Pn−2 (t) . Sie bilden eine Sturmsche Kette: Für die Nullstellen λ1n ≤ λ2n ≤ · · · ≤ λnn von Pn gilt n n−1 λkn−1 , −1 < λk < λk

(mit

λ0n

= −kAk∞ und

λkn

≤ λkn+1

0 (4.10) Sei |λ1 | = ρ(A) und |λn | < |λ1 | für n = 2, ..., N. Dann gilt für alle w ∈ R 1 1 k k A w = v1 . lim r (A, A w) = λ1 , lim k −→∞ k −→∞ |Ak w|2 (4.11) Sei |w|2 = 1 und µ = r (A, w). Dann gilt r (A, x) =

min |λn − µ| ≤ |Aw − µw|2 .

n=1,...,N

Eine konvergente Folge (d k ) in R mit Grenzwert d ∗ konvergiert a) linear, wenn c ∈ (0, 1) und k0 > 0 existieren mit |d k +1 − d ∗ | ≤ c |d k − d ∗ | für k ≥ k0 b) superlinear, wenn zu jedem ε > 0 ein k0 > 0 existiert mit |d k +1 − d ∗ | ≤ ε |d k − d ∗ | für k ≥ k0 c) von der Ordnung p > 1, wenn wenn C > 0 existiert mit |d k +1 − d ∗ | ≤ C |d k − d ∗ |p . C. Wieners: Einführung in die Numerische Mathematik

11

4 Eigenwertberechnung (4.12)

(4.13)

Inverse Iteration mit variablem shift S0) Wähle z0 ∈ RN , z0 6= 0N , ε ≥ 0. Setze k = 0. S1) Setze wk = |z1| zk , µk = r (A, wk ). k 2 S2) Falls |Awk − µk wk |2 ≤ ε STOP. S3) Berechne zk +1 = (A − µk IN )−1 wk . S4) Setze k := k + 1, gehe zu S1). Wenn der Startvektor z0 hinreichend nahe bei einem Eigenvektor vm mit isoliertem Eigenwert λm liegt, konvergiert die Iteration kubisch (d.h. von der Ordnung p = 3). QR-Iteration mit shift (A ∈ RN×N symmetrisch) S0) Berechne A0 = QAQ T 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 = 12 (Ak [N − 1, N − 1] − Ak [N, N]) und q sk = Ak [N, N] + dk − sgn(dk ) dk2 + 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 = QkT Ak Qk . Falls der shift µk = Ak [N, N] gewählt wird, entspricht die QR-Iteration der Inversen Iteration mit variablem shift und Startvektor z0 = eN .

C. Wieners: Einführung in die Numerische Mathematik

12

4 Eigenwertberechnung (4.14)

Gershgorin Zu A ∈ RN×N sind die Gershgorin-Kreise durch  Kn = λ ∈ C : |λ − A[n, n]| ≤ ∑ |A[n, k ]| ,

n = 1, ..., N

k 6=n

defininiert. Dann gilt σ (A) ⊂

N [

Kn .

n=1 (4.15)

Sei A ∈ RN×N symmetrisch mit Eigenwerten λ1 ≥ λ2 ≥ · · · ≥ λN . λn

=

λN+1−n

=

max

dim S=n

min

dim S=n

min r (A, x) ,

0N 6=x∈S

max r (A, x) .

0N 6=x∈S

C. Wieners: Einführung in die Numerische Mathematik

13

5 Iterative Lösungsverfahren für lineare Gleichungen (5.1)

Sei A, B ∈ RN×N mit ρ(IN − BA) < 1. Dann ist A invertierbar, und es gilt für alle b ∈ RN und alle Startvektoren x 0 ∈ RN konvergiert die Iteration x k +1 = x k + B(b − Ax k ) , gegen lim

k −→∞

(5.2)

xk

k = 0, 1, 2, ...

= A−1 b.

Sei K ∈ RN×N und ε > 0 Dann existiert eine Vektor-Norm | · | und eine zugeordnete Matrix-Norm k · k, so dass kK k ≤ ρ(K ) + ε gilt. Anwendung: Es gilt |x − x k | ≤ kIN − BAkk |x − x 0 | (lineare Konvergenz).

(5.3)

Konvergenz des Gauß-Seidel-Verfahrens Sei A = L + D + R ∈ RN×N symmetrisch positiv definit und sei B = (L + D)−1 . √ A Dann gilt bezüglich der Energienorm |x|A = x T Ax und kK kA = sup |Kx| |x| x6=0N

A

kIN − BAkA < 1 . k Anwendung der Neumannschen Reihe ergibt dann A−1 = ∑∞ k =0 (IN − BA) B.

C. Wieners: Einführung in die Numerische Mathematik

14

5 Iterative Lösungsverfahren für lineare Gleichungen (5.4)

Eine Matrix A ∈ RN×N heißt stark diagonal-dominant, wenn N

∑ |A[n, k ]| ≤ |A[n, n]| ,

n = 1, ..., N ,

k =1 k 6=n

N

und wenn ein j ∈ {1, . . . , N} existiert

∑ |A[j, k ]| < |A[j, j]|.

k =1 k 6=n

(5.5)

(5.6)

Eine Matrix A ∈ RN×N sei irreduzibel. Dann existiert zu jedem Paar j 6= n eine Folge j = j0 , j1 , j2 , . . . , jR = n mit A[j1 , j0 ] 6= 0, ..., A[jR , jR−1 ] 6= 0. Eine Matrix A ∈ RN×N sei stark diagonal-dominant und irreduzibel. Dann gilt: a) A ist regulär und das Jacobi-Verfahren x k +1 = x k + diag(A)−1 (b − Ax k ) konvergiert. b) Sei A[n, n] > 0 für alle n. Dann ist A positiv definit. c) Sei A[n, n] > 0 und A[n, k ] ≤ 0 für n 6= k . Dann ist A−1 [n, k ] ≥ 0 für alle n, k .

C. Wieners: Einführung in die Numerische Mathematik

15

5 Iterative Lösungsverfahren: Krylov-Verfahren (5.7)

(5.8)

Zu C ∈ RN×N und d ∈ RN ist k -te Krylov-Raum   Kk (C, d) = span d, Cd, ..C k −1 d = P(C)d : P ∈ Pk −1 . Zu einer regulären Matrix A ∈ RN×N , einer rechten Seite b ∈ RN und einem Startwert x 0 ∈ RN sei x ∈ RN die Lösung von Ax = b und r 0 = b − Ax 0 . Sei B ∈ RN×N eine reguläre Matrix (Vorkonditionierer). Wenn dim Kk (AB, r 0 ) < k für ein k gilt, dann ist x ∈ x 0 + Kk −1 (BA, Br 0 ). Sei h·, ·i ein Skalarprodukt in RN . Gram-Schmidt-Verfahren zur Berechnung einer Orthonormalbasis v 1 , ..., v k von Kk (BA, Br 0 ) = span{Br 0 , BABr 0 , ..., (BA)k −1 Br 0 } = {Vk y : y ∈ Rk } , S0) Wähle x 0 ∈ RN , setze r 0 = b − Ax 0 , z 1 = Br 0 , h10 = |z 1 |V und v 1 = S1) Für k = 1, 2, 3, ... berechne z k +1

wk


Vk = v 1 |....|v k . 1 1 h10 z .

= BAv k ,

k

=

w k − ∑ hjk v j mit hjk = hv j , w k iV j=1

v k +1

=

1 hk +1,k

z k +1 mit hk +1,k = |z k +1 |V

+1 Dann gilt BAv k = ∑kj=1 hjk v j , also BAVk = Vk +1 Hk mit Hk = (hjm ) ∈ Rk +1,k . C. Wieners: Einführung in die Numerische Mathematik

16

5 Iterative Lösungsverfahren: GMRES-Verfahren S0) Wähle x 0 ∈ RN , ε > 0. Berechne r 0 = b − Ax 0 , z 1 = Br 0 , h10 = |z 1 |2 und v 1 = S1) Berechne z k +1

wk

= BAv k

1 1 h10 z .

Setze k = 1.

und k

=

w k − ∑ hjk v j mit hjk = (v j )T w k j=1

v k +1

=

1 hk +1,k

z k +1 mit hk +1,k = |z k +1 |2

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

S3) Wenn ρk < ε, setze x k = x 0 + ∑ yjk v j

STOP.

j=1

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

Es gilt ρk =

min

z∈x 0 +span{v 1 ,...,v k }

|B(b − Az)|2 .

Das GMRES-Verfahren ist wohldefiniert, und wenn es abbricht, gilt |x − x k |2 ≤ k(BA)−1 k2 ε. (5.5)

Für C ≥ α > 0 gelte z T BAz ≥ αz T z und kBAk2 ≤ C.  k /2 α2 Dann gilt für das GMRES-Verfahren |x k − x|2 ≤ κ2 (BA) 1 − 2 |x 0 − x|2 . C

C. Wieners: Einführung in die Numerische Mathematik

17

5 Iterative Lösungsverfahren: CG-Verfahren S0) Wähle x 0 ∈ RN , ε > 0. Berechne r 0 = b − Ax 0 , y 0 = Br 0 , ρ0 = (y 0 )T r 0 und d 1 = y 0 . Setze k = 0. S1) Falls ρk ≤ ε

STOP

S2) Setze k := k + 1 und berechne

d

uk

=

Ad k ρk −1 (u k )T d k

αk

=

xk

=

x k −1 + αk d k

r

k

=

r k −1 − αk u k

y

k

=

Br k

ρk

=

k +1

=

(y k )T r k ρ yk + k dk ρk −1

Gehe zu S1). (5.6)

|x k − x|A

= ≤

min

z∈x 0 +span{d 1 ,...,d k }

|z − x|A

 pκ(BA) − 1 k |x 0 − x|A . max |P(λ )||x 0 − x|A ≤ 2 p P∈Pk , P(0)=1 λ ∈σ (BA) κ(BA) + 1 min

C. Wieners: Einführung in die Numerische Mathematik

18

5 Transformierte Chebychev-Polynome (5.6)

q min

max |P(t)| ≤ 2  q

P∈Pk , P(0)=1 t∈[a,b]

b a

−1

b a

+1

k 

P4 , P5 , P6 , P12 mit Pk (0) = 1 zu [a, b] = [0.1, 2.1] C. Wieners: Einführung in die Numerische Mathematik

19

6 Lösungsverfahren für nichtlineare Gleichungen (6.1)

(6.2)

Sei D ⊂ RN konvex und F : D −→ RN ein stetig differenzierbares Vektorfeld. N×N . Sei x ∗ ∈ D eine Nullstelle
von F , und sei B ∈ R Wenn ρ IN − BJF (x ∗ ) < 1 gilt, dann existiert ein δ > 0, so dass für alle x 0 ∈ B(x ∗ , δ ) die Fixpunktiteration x k +1 = Φ(x k ) mit Φ(x) = x − BF (x) linear gegen x ∗ konvergiert. Banachscher Fixpunktsatz Sei D ⊂ RN konvex und φ : D −→ D kontrahierend mit q ∈ (0.1), d.h. |φ (y ) − φ (z)| ≤ q |y − z| für y , z ∈ D . Dann existiert genau ein Fixpunkt x ∗ ∈ D von φ , d.h. φ (x ∗ ) = x ∗ , und es gelten für die Fixpunktiteration x k +1 = Φ(x k ) die Abschätzungen |x k − x ∗ | ≤

(6.3)

qk |x 0 − x 1 | 1−q

und

|x k − x ∗ | ≤

q |x k − x k −1 | . 1−q

Seien α, β , γ > 0 mit 2αγ < β 2 und P(t) = α − β t + 2γ t 2 . Dann konvergiert für t0 = 0 das Newton-Verfahren tk +1 = tk − P 0 (tk )−1 P(tk ) quadratisch gegen die kleinste Nullstelle t ∗ von P, und {tk } ist monoton steigend mit t0 = 0 < t1 < · · · < tk < tk +1 = tk +

C. Wieners: Einführung in die Numerische Mathematik

γ (tk − tk −1 )2 2α p ≤ t∗ = . 2 β − γtk β + β 2 − 2αγ 20

6 Lösungsverfahren für nichtlineare Gleichungen (6.3)

Sei D ⊂ RN offen und F : D −→ RN ein stetig differenzierbares Vektorfeld. Sei x 0 ∈ D mit a) |F (x 0 )| ≤ α b) |JF (x 0 )y | ≥ β |y | für y ∈ RN c) B(x 0 , 2α/β ) = {z ∈ RN : |z − x 0 | ≤ 2α/β } ⊂ D d) kJF (y ) − JF (z)k ≤ γ |y − z| für y , z ∈ B(x 0 , 2α/β ) e) 2αγ < β 2 . Dann ist das Newton-Verfahren x k +1 = x k − JF (x k )−1 F (x k ) wohldefiniert und konvergiert quadratisch gegen x ∗ ∈ D mit |x ∗ − x 0 | ≤

β+

2α p . β 2 − 2αγ

Gedämpftes Newton-Verfahren S0) Wähle x 0 ∈ D, ε > 0, θ ∈ (0, 1). Setze k = 0. S1) Falls |F (x k )| ≤ ε

STOP

S2) Löse JF (x k )d k = −F (x k ). S3) Bestimme tk ∈ {1, θ , θ 2 , ..., θ r } mit |F (x k + tk d k )| minimal. S4) Setze x k +1 = x k + tk d k , k := k + 1 und gehe zu S1). C. Wieners: Einführung in die Numerische Mathematik

21

7 Orthogonalpolynome Zu −∞ ≤ a < b ≤ ∞ und einem Gewicht W ∈ C(a, b), W > 0 definiere V = C(a, b) als Euklidischen Vektorraum mit dem Skalarprodukt (u, v )V =

Z b

u(t)v (t)W (t) dt .

a

Aufgabe: Zu f ∈ C(a, b) bestimme P ∈ PN mit kP − f kV = min kQ − f kV . Q∈PN

N

Sei Q0 , . . . , QN ∈ PN eine Orthogonalbasis. Dann gilt P(t) =

∑ (f , Qk )V

Qk (t).

k =0

(7.1)

Eine orthogonale Basis Qn , n = 0, 1, 2, ... berechnet sich mit Q−1 ≡ 0, Q0 ≡ ρn+1 Qn+1 (t) = (t − βn ) Qn (t) − ρn Qn−1 (t)

1 ρ0

und

n = 0, 1, 2, ...

mit ρ0 = k1kV , ρn = (Qn , tQn−1 )V > 0, und βn = (Qn , tQn )V . Es gilt grad Qn = n. (7.2)

Das Orthogonalpolynom Qn besitzt n verschiedene Nullstellen in (a, b). Die Orthogonalpolynome Q0 , Q1 , Q2 , ... bilden eine Sturmsche Kette.

C. Wieners: Einführung in die Numerische Mathematik

22

7 Orthogonalpolynome – Chebychev-Polynome 1 Wir betrachten den Fall (a, b) = (−1, 1) und W (t) = √ . 1 − t2 Definiere T0 ≡ 1, T1 (t) = t und Tn+1 ∈ Pn+1 mit Tn+1 (t) = 2tTn (t) − Tn−1 (t). Für t ∈ [−1, 1] gilt Tn (t) = cos(n arccos t) , n 1  p n 1 p 2 t + t −1 + t − t2 − 1 . (Tn , Tk )V = 0 für n 6= k , und für |t| ≥ 1 gilt Tn (t) = 2 2 (7.3)

Sei ξ 6∈ [−1, 1] und P(t) = Tn (t)/Tn (ξ ). Dann gilt max |P(t)| ≤ max |Q(t)| für alle Q ∈ Pn mit Q(ξ ) = 1. t∈[−1,1]

(7.4)

t∈[−1,1]

Sei 0 < α < β und κ = β /α. Dann existiert ein P ∈ Pn mit P(0) = 1 √ n κ −1 und max |P(x)| ≤ 2 √ . κ +1 x∈[α,β ]

Transformierte Chebychev-Polynome P4 , P5 , P6 , P12 mit Pn (0) = 1 zu [α, β ] = [0.1, 2.1] C. Wieners: Einführung in die Numerische Mathematik

23

7 Orthogonalpolynome (7.5)

Approximationssatz von Weierstraß Sei [a, b] ein kompaktes Intervall. Dann gilt: Jede stetige Funktion f ∈ C([a, b]) lässt sich gleichmäßig durch Polynome approximieren. Konstruktion: O.E. sei f ∈ C(R) mit f (t) = 0 für |t| > 0.5 und [a, b] = [−0.25, 0.25]. Dann wähle P2n (t) =

1 rn

Für |t| ≤ 0.25 gilt P2n (t) =

1 rn

Nun sei (v , w)0 =

R1

Z 1

(1 − t 2 )n f (s − t) ds ,

−1 n

∑ ankj t k

j,k =0

rn =

Z 1

(1 − t 2 )n dt .

−1

2n R1 j s f (s) dt ∈ P2n mit (1 − (s − t)2 )n = ∑ ankj t k sj . j,k =0

−1

v (t)w(t) dt das Skalarprodukt in L2 (−1, 1), seien Q0 , Q1 , Q2 , ... die

−1 n

Orthogonalpolynonme, und sei pn (f ) = ∑ (f , Qk )0 Qk die Orthogonalprojektion auf Pn . k =0

(7.6)

Für jede stetig Funktion f ∈ C[−1, 1] gilt lim kpn (f ) − f k0 = 0.

(7.7)

Sei f ∈ Ck [−1, 1]. Dann existiert C > 0 mit kpn (f ) − f k0 ≤ Cn−k

n−→∞

C. Wieners: Einführung in die Numerische Mathematik


d k dt

f . 0

24

8 Interpolation (8.1)

(8.2)

(8.3)

(8.4)

(8.5)

(8.6)

Lagrange-Interpolation Zu Stützstellen t0 < t1 < · · · < tN und Werten f0 , f1 , . . . , fN ∈ R bestimme P ∈ PN mit P(tn ) = fn . Die Interpolationsaufgabe (8.1) ist eindeutig lösbar. Hermite-Interpolation Zu t0 ≤ t1 ≤ · · · ≤ tN definiere dn = max{n − k : tk = · · · = tn } und d = max dn . Zu f ∈ C d (R) bestimme P ∈ PN mit  d dn  d dn P(tn ) = f (tn ) für n = 0, . . . , N. dt dt Die Interpolationsaufgabe (8.3) ist eindeutig lösbar.
1 d N Der eindeutig bestimmte höchste Koeffizient bN ≡ N! P(t) der Lösung der dt Interpolationsaufgabe (8.1) wird mit bN = f [t0 , . . . , tN ] bezeichnet. Zu t0 ≤ t1 ≤ · · · ≤ tN definiere die Newton-Basis von PN durch ω0 ≡ 1 und ωk +1 (t) = ωk (t)(t − tk ) ,

k = 0, ..., N − 1 .

N

a) P(t) = ∑ f [t0 , . . . , tk ]ωk (t) löst die Interpolationsaufgabe (8.3). k =0

b) Für f ∈ C(R) gilt f (t) = P(t) + f [t0 , . . . , tN , t] ωN+1 (t). C. Wieners: Einführung in die Numerische Mathematik

25

8 Interpolation (8.7)

Für tk < tn gilt f [tk , . . . , tn ] =

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

Sei f ∈ Cn−k (R). Für tk = tk +1 = · · · = tn gilt f [tk , . . . , tn ] =

1 (n − k )!



d dt

n−k f (tk ).

n

Sei Pkn (t) = ∑ f [tk , . . . , tj ]ωkj (t) mit ωkk ≡ 1 und mit ωk ,j+1 (t) = ωkj (t)(t − tj ) die Interpolation j=k

an tk ≤ · · · ≤ tn zu fk , ..., fn . Für tk < tn gilt Pkn (t) = (8.9)

t−tk tn −tk

−t Pk +1,n (t) + ttnn−t Pk ,n−1 (t). k

Für den Interpolationsfehler gilt f (t) − PN (t) =

 d N+1 1 f (τ) ωN+1 (t) (N + 1)! dt

mit einer Zwischenstelle τ = τ(t) ∈ I := [min{t, t0 } max{t, tN }]. (8.11)

Es gilt f [t0 , . . . , tN ] =

Z ΣN

 d N  N  f ∑ sn tn ds dt n=0

   s    0  N . mit ΣN =  ..  ∈ RN+1 : ∑ sn = 1, sn ≥ 0     n=0 sn  N 1 d (8.12) f [t0 , . . . , tN ] = f (τ) mit τ ∈ [min tn , max tn ]. N! dt C. Wieners: Einführung in die Numerische Mathematik

(Standard-Simplex)

26

9 Splines (9.1)

Zu einer Zerlegung ∆ : a = t0 < t1 < · · · < tN = b von [a, b] definiere den Spline-Raum von Grad k  Sk (∆) = S ∈ C k −1 [a, b] : Sn := S|[tn−1 ,tn ] ∈ Pk , n = 1, . . . , N . S ∈ Sk (∆) heißt interpolierender Spline zu f ∈ C 0 [a, b] wenn S(tn ) = f (tn ).

(9.2)

Es gilt dim Sk (∆) = k + N. Wir betrachten interpolierende Splines von Grad k = 2r + 1 mit einer der Randbedingungen (I) Natürliche Randbedingung  d m  d m S(a) = 0 und S(b) = 0 dt dt

für

r + 1 ≤ m ≤ 2r

∈ C r [a, b]

(II) Hermite-Randbedingung zu f  d m  d m S(a) = f (a) und dt dt (III) periodische Randbedingungen  d m  d m S(a) = S(b) für dt dt C. Wieners: Einführung in die Numerische Mathematik

 d m dt

S(b) =

 d m dt

f (b) für 1 ≤ m ≤ r

1 ≤ m ≤ 2r 27

9 Splines Zu g, h ∈ C m [a, b] definiere Z b  m  d m d (g, h)m = g(t) h(t) dt , dt dt a (9.3)

(9.4)

p (g, g)m .

 Zu a ≤ t0 < t1 < · · · < tN ≤ b und f0 , ..., fN definiere M = g ∈ C 2 [a, b] : g(tn ) = fn . Dann existiert genau ein S ∈ M mit |S|2 ≤ |g|2 für alle g ∈ M. Es gilt S 00 (t) = 0 für t ∈ [a, t0 ] ∪ [tN , b] und S|[tn−1 ,tn ] ∈ P3 . Sei S ∈ S2r +1 (∆) ein interpolierender Spline zu f und sei zusätzlich eine der Randbedingungen (I), (II) oder (III) erfüllt. Sei g ∈ C r +1 [a, b] mit g(tn ) = f (tn ) für n = 0, ..., N, und zusätzlich gelte im Fall (II)  d m  d m  d m  d m g(a) = f (a) und g(b) = f (b) für 1 ≤ m ≤ r , dt dt dt dt im Fall (III)  d m  d m g(a) = g(b) für 1 ≤ m ≤ r . dt dt Dann gilt (S − g, S)r +1 = 0

(9.5)

|g|m =

und

|S|r +1 ≤ |g|r +1 .

Die Spline-Interpolation S ∈ S2r +1 (∆) zu f mit einer der Randbedingungen (I), (II), (III) ist eindeutig lösbar.

C. Wieners: Einführung in die Numerische Mathematik

28

9 Splines (9.6)

Die Momente Mn = S 00 (tn ) von S ∈ S3 (∆) zu f sind eindeutig durch µn Mn−1 + Mn + λn Mn+1 = 3f [tn−1 , tn , tn+1 ] und µn =

hn , 2(hn + hn+1 )

Sn (t)

(9.7)

=

λn =

hn+1 mit hn = tn − tn−1 bestimmt und es gilt 2(hn + hn+1 )

Mn (t − tn−1 )3 + Mn−1 (tn − t)3 f (tn ) + f (tn−1 ) hn2 + − (Mn + Mn−1 ) 6hn 2 12  f (t ) − f (t ) h   tn + tn−1 n n n−1 + − (Mn − Mn−1 ) t − , t ∈ [tn−1 , tn ] . hn 6 2

Es gilt für m = 0, . . . , r und h = max hi , hn = tn − tn−1 i=n,...,N

|S − f |m ≤ 2−(r +1−m)/2 (9.8)

(r + 1) r +1−m h |f |r +1 m!

Es gilt für f ∈ C 2 [a, b] kS 00 k∞ ≤ 3 kf 00 k∞ .

(9.9)

Es gilt für f ∈ C 4 [a, b]  d 4 kS − f k∞ ≤ h4 k f k∞ . dt

C. Wieners: Einführung in die Numerische Mathematik

29

10 Numerische Integration – Newton-Cotes-Formeln (10.1)

Eine Quadraturformel IΞ : C[a, b] −→ R mit IΞ (f ) = ∑ ωξ f (ξ ) zu Stützstellen Ξ ⊂ [a, b] und ξ ∈Ξ

Gewichten ωξ ist eine Linearform zur Approximation des Integrals I(f ) =

Z b

f (t)dt.

a

Sie heißt exakt von der Ordnung p, wenn der Fehler I(P) = IΞ (P) für Polynome P ∈ Pp . (10.2)

Sei |Ξ| = N + 1. Wenn IΞ exakt von der Ordnung N ist, dann gilt Z b

ωξ =

a

Lξ (t)dt mit Lξ (t) = ∏η∈Ξ\{ξ }

t−η ξ −η

∈ PN .

Die Quadraturformeln von der Ordnung N zu äquidistanten Stützstellen ξn = a + nh, h = heißen Newton-Cotes-Quadraturformeln. Es gilt: Z N N b − a (−1)N−n ωn = ωN−n = ∏ (t − j)dt. N n!(N − n)! 0 j=0, j6=n (10.4) Die Newton-Cotes-Formeln sind für gerade N exakt von der Ordnung N + 1. (10.3)

(10.5)

Seien g, h ∈ C[a, b] mit g(t) ≥ 0 für t ∈ [a, b]. Dann gilt: Z b a

(10.6)

b−a N

g(t)h(t)dt = h(τ)

Z b

g(t)dt für eine Zwischenstelle τ ∈ (a, b).

a

Sei f ∈ C 4 [a, b]. Dann gilt für die Simpsonregel:  b − a 5 1 I2 (f ) − I(f ) = f 0000 (τ) für ein τ ∈ (a, b). 2 90

C. Wieners: Einführung in die Numerische Mathematik

30

10 Numerische Integration – Gauß-Quadatur (10.7)

Die maximale Ordnung einer Quadratur mit N Stützstellen ist 2N − 1.

(10.8)

Seien ξ1 , . . . , ξN die Nullstellen des Orthogonalpolynoms QN in (a, b) bzgl. des Gewichts W . N

Dann ist die zugehörige Gauß-Quadratur GN (f ) = ∑ ωn f (ξn ) exakt für P ∈ P2N−1 . n=1

Die Gewichte der Gauß-Quadratur sind positiv. Z b (10.9)

a

f (t)W (t) dt − GN (f ) =

N 1  d 2N f (τ) (RN , RN )W mit RN (t) = ∏ (t − ξn ). (2N)! dt n=1

Eine Folge {IΞN }N∈N von Quadraturen mit Stützstellen ΞN ⊂ [a, b] und Ge- wichten ωξN heißt konvergent, wenn lim |I(f ) − IΞN (f )| = 0 für alle f ∈ C[a, b]. N→∞

(10.10)

Sei {IΞN } für alle Polynome P konvergent (d.h. limN→∞ IΞN (P) = I(P)), und seien die Gewichte uniform beschränkt durch ∑ |ωξN | ≤ C für alle N. ξ ∈ΞN

Dann ist {IΞN } konvergent für alle f ∈ C[a, b]. (10.11)

Wenn alle Gewichte positiv sind und IΞN von der Ordnung N exakt ist, ist {IΞN } konvergent.

C. Wieners: Einführung in die Numerische Mathematik

31

10 Numerische Integration (10.12)

Sei Ih − I = Chp + O(hq ) mit q > p. Dann gilt für die Extrapolation ˆIh = Ih +


1 I −I 2p − 1 h 2h

a) |I − ˆIh | = O(hq )
b) |I − Ih | = 2p1−1 Ih − I2h + O(hq ) Sei g ∈ C 2m+2 [0, 1]. Es existieren Konstanten bj und eine Zwischenstelle ξ ∈ (0, 1) mit Z 1  d 2j−1  d 2j−1  m d 2m+2 1 1 g(0) − g(1) + bm+1 g(ξ ). g(t) dt = g(0) + g(1) + ∑ bj 2 2 dt dt dt 0 j=1 b−a 2m+2 (10.14) Für f ∈ C [a, b] und h = gilt n  d 2j−1   d 2j−1 m f (a) − f (b) . Th (f ) − I(f ) = ∑ cj h2j + O(h2m+2 ) mit cj = bj dt dt j=1

(10.13)

(10.15)

Zu Tk 0 = Thk (f ) mit hk = 2−k h für k = 0, 1, . . . , m definiere rekursiv 1 (Tk ,i−1 − Tk −1,i−1 ), i = 1, . . . , m und k = i, . . . , m . 4i − 1

 d 2m+2 Dann gilt: |Tmm − I| ≤ Ch2m+2 f ∞ mit C > 0 unabhängig von f . dt Tki = Tk ,i−1 +

(10.16)

a) |I − Tkk | = O(h2k +2 ) b) |I − Tk ,k −1 | = |Tkk − Tk ,k −1 | + O(h2k +2 ) .

C. Wieners: Einführung in die Numerische Mathematik

32

10 Numerische Integration – Romberg-Quadratur using namespace std; #include #include #include using __gnu_cxx::power; const int K = 6; const int N = power(2,K); double Romberg (double F(double), double a, double b, double Eps) { double f[N+1]; double h = (b-a) / N; for (int i=0; i