Zusammenfassung Numerische Mathematik Patrik Rohner ::

(1)

Zusammenfassung Numerische Mathematik

Patrik Rohner :: <rohnerpa@student.ethz.ch> :: FS’09

1 Grundlagen der Numerik

1.1 Gleitpunktarithmetik

Gleitpunktzahl: x ˆ = σM B

^E

σ: Vorzeichen B: Basis M: Mantisse E: Exponent

F¨ ur die Mantisse M gilt:

M = ∑

^l_j=1

m

_−j

B

^−j

0 ≤ m

_−j

≤ B − 1 : Ziffern der Mantisse

l : Mantissenl¨ ange (fest)

m

₋₁

≠ 0 f¨ ur x ˆ ≠ 0 : Normierung

F¨ ur den Exponenten E gilt:

E = τ ∑

^k_j=1

e

j

B

^k−j

τ : Vorzeichen

0 ≤ e

j

≤ B − 1 : Ziffern des Exponenten k : definierter Exponentenbereich Definition: M (B, l, k) = Menge der Maschinenzahlen

/ M (2,5,3) ⇒ M

max

=

³¹

32

≈ 1 M

min

=

¹

2

E

_max/min

= ±7 ˆ

x

max

= M

max

⋅ B

^E^max

=

³¹

32

⋅ 2

⁷

x ˆ

min

= M

min

⋅ B

^E^min

=

¹

2

⋅ 2

⁻⁷

/ M (B, l, k) ⇒ M

max

=

^B

l−1

B^l

M

min

=

¹

B

E

_max/min

= ±B

^k

− 1 ˆ

x

max

=

^B

l−1

B^l

⋅ B

^B^k⁻¹

ˆ x

min

=

¹

B

⋅ B

^−(B^k⁻¹⁾

Definition Maschinengenauigkeit: eps ist die kleinste Zahl, die zu 1 addiert noch eine Ver¨ anderung bewirkt.

∣ x − x ˆ

x ∣ ≤ B

^E−L

2B

^E−1

=

B

^1−L

2 =∶ eps

Pseudo-Arithmetik: (ˆ x + y) ⋅ ˆ z ˆ = z ˆ ⋅ (ˆ y + x) ˆ (ˆ x + y) ⋅ ˆ z ˆ ≠ ˆ x ⋅ ˆ z + ˆ y ⋅ z ˆ (ˆ x + y) + ˆ z ˆ ≠ x ˆ + (ˆ y + z) ⇒ ˆ von klein nach gross summieren

1.2 Fehlerfortpflanzung

Absoluter Fehler: ∆x = x ˆ − x

∆(x ± y) = (ˆ x ± y) − (x ˆ ± y) = ∆x ± ∆y

∆(x ∗ y) = ˆ xˆ y − xy ≅ x∆y + y∆x = xy(δx + δy)

∆(x/y) ≅ x/y(δx − δy) Relativer Fehler: δx =

^∆x

x

δ(x ± y) =

^∆(x±y)_x±y

=

_x±y^x

⋅

^∆x

x

±

_x±y^y

⋅

^∆y

y

=

_x±y^x

δx ±

_x±y^y

δy δ(x ∗ y) =

^∆(x∗y)_xy

≅ δx + δy

δ(x/y) =

^∆(x/y)_x/y

≅ δx − δy

Ausl¨ oschung: ∣x ± y∣ ≪ {∣x∣,∣y∣} ⇒ ∣δ(x ± y)∣ ≫ ∣δx∣ + ∣δy∣

Faustregel: wenn δx = 10

^−k

, dann sind ≈ k Ziffern richtig.

⇒ besser (x

³

− y

³

)/(x

²

+ xy + y

²

) als x − y rechnen

⇒ besser (x

⁶

− y

⁶

)/(x

⁴

+ x

²

y

²

+ y

⁴

) als x

²

− y

²

rechnen

1.3 Kondition eines Problems

κ

_H

= Konditionszahl = Verst¨ arkungsfaktor des relativen Fehlers von x Konditionszahl κ

H

: δH = ∣

^xH

′(x)

H(x)

∣ ∗ ∣δx∣ = κ

H

∗ ∣δx∣

1.3.1 Kondition einer Matrix κ ( A ) = ∣∣ A ∣∣ ⋅ ∣∣ A

⁻¹

∣∣ :

A ⋅ x ˆ = ˆ b ⇒ ∣∣δx∣∣ ≤ κ(A) ⋅ ∣∣δb∣∣

A ˆ ⋅ x ˆ = ˆ b ⇒ ∣∣δx∣∣ ≤

1−κ(A)∣∣δA∣∣^κ(A)

(∣∣δA∣∣ + ∣∣δb∣∣)

Faustregel: Bei d Dezimalstellen und κ(A) ≅ 10

^k

hat die L¨ osung im schlechtesten Fall noch d − k richtige Stellen.

in der 2-Norm: ∣∣A∣∣

2

= √

µ

max

(µ

max

= max. EW von A

^T

A)

A orthogonal: A

^T

= A

⁻¹

∣∣A∣∣

2

= 1

A symmetrisch: A

^T

= A ∣∣A∣∣

2

= ∣λ

max

∣ ∣∣A

⁻¹

∣∣

₂

=

_∣λ¹

min∣

A regul¨ ar: ∣∣A

⁻¹

∣∣

₂

=

√ ¹ µ_min

κ(A) =

√ µ

max

√ µ

min

und falls A

^T

= A gilt κ(A) = ∣λ

max

∣

∣λ

min

∣ κ ≈ 1 gut konditioniertes Problem κ ≫ 1 schlecht konditioniertes Problem

>> k_A=cond(A)

2 Lineare Gleichungssysteme

2.1 LRP-Zerlegung

LRP-Zerlegung von A aus Ax = b . . . . i LRP-Zerlegung von A: LR = P A >> [L,R,P]=lu(A) ii Vorw¨ artseinsetzen: Lc = P b nach c >> c=L\(P*b) iii R¨ uckw¨ artseinsetzen: Rx = c nach x >> x=R\c Der Algorithmus braucht

²₃

n

³

+ 2n

²

Operationen. O(n

³

) 2.1.1 Pivotstrategien (gegen Ausl¨ oschung)

Diagonalstrategie: Diagonalelement = Pivot (schlecht)

Spaltenmaximumstrategie: Das betragsm¨ assig gr¨ osste Element wird gew¨ ahlt.

relative Spaltenmaximumstrategie: Das betragsm¨ assig gr¨ osste Element, relativ zum zweitgr¨ ossten Element der Zeile, wird gew¨ ahlt.

2.2 Iterative Methoden

Alle hier ben¨ otigen pro Iterationsschritt grob eine Matrix-Vektor- Multiplikation (O(n

²

), wenn sparse: O(n)) und konvergieren linear (O(n)).

2.2.1 station¨ are Iterationsverfahren x

_k+1

= T ⋅ x

k

+ c

Der absolute Fehler ver¨ andert sich mit e

_k

= T

^k

⋅ e

0

, womit gilt:

Konvergenzbedingung: ∣∣T ∣∣ < 1 >> norm(T)

F¨ ur den Spektralradius ρ(T ) ∶= max ∣λ

i

∣ gilt ρ(T ) ≤ ∣∣T ∣∣ und so alternative Konvergenzbed.: ρ(T ) < 1 >> max(abs(eig(T))) Absch¨ atzung (worst case) der Anzahl Iterationen f¨ ur bestimmte Fehler:

∣∣e∣∣

2

= ∣∣T

^k

e

⁰

∣∣

₂

∣∣e

^k

∣∣

₂

≤ ∣∣T ∣∣

^k₂

∣∣e

⁰

∣∣

₂ ^∣∣e

k∣∣2

∣∣e⁰∣∣2

≤ ∣∣T ∣∣

^k₂

k ≥

ln(rel.Fehler) ln∣∣T∣∣2

Aufspaltung f¨ ur bestimmte Verfahren A = L + D + R

>> D=diag(diag(A)) >> L=tril(A)-D >> R=triu(A)-D Jacobi-Verfahren

T = −D

⁻¹

⋅ (L + R) c = D

⁻¹

b

Konsistenzbedingung: A regul¨ ar, d.h det(A) ≠ 0

hinreichende Konvergenzbedingung: A strikt diagonal dominant (max

∣.∣ in der Diagonalen)

notwendige Konvergenzbedingung (bei allen 3 Verfahren): ∣∣T ∣∣ < 1 Gauss-Seidel-Verfahren

T = −(D + L)

⁻¹

R c = (D + L)

⁻¹

b (D + L)x

_k+1

= −Rx

k

+ b >> x (k+1)=(D+L)\(-Rx k+b)

Falls A strikt diagonal dominant oder A symmetrisch positiv definit ist, konvergiert das Verfahren, und zwar schneller als das Jacobi-Verfahren.

SOR-Verfahren (Successive Over-Relaxation) x

_k+1

= (D + ωL)

⁻¹

[−ωR + (1 − ω)D]

´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶

=T(ω)

x

_k

+ (D + ωL)

⁻¹

ωb

´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶

=c(ω)

ω = Relaxationsparameter ω

optimal

=

²

1+√

1−[ρ(D⁻¹(L+R))]²

ρ: Spektralradius

(2)

2.2.2 Methode der konjugierten Gradienten (CG)

Zu l¨ osendes Gleichungssystem: Ax + b = 0

Voraussetzungen: A symmetrisch und positiv definit.

F(u) ∶= 1/2u

^T

Au + u

^T

b gradF =

⎛

⎜ ⎜

⎝

∂F

∂u₁

...

∂F

∂un

⎞

⎟ ⎟

⎠

= Au + b = r(u) (Residuen) Idee: Anstelle Ax + b = 0 zu l¨ osen, wird das Minimum von F(u) gesucht.

Es gilt F (x) = min

_u∈Rn

F (u).

CG-Verfahren Ax + b = 0 . . . .

Start: u

0

, r

0

= Au

0

+ b, T OL

1. Schritt: (Steilster Abstieg) i p

1

= −r

0

ii ρ

1

=

_(p^(r⁰^,r⁰⁾

1,Ap₁)

iii u

1

= u

0

+ ρ

1

p

1

iv r

1

= r

0

+ ρ

1

Ap

1

k. Schritt: (k ≥ 2 , konjugierter Gradient) i e

_k−1

=

^(r_(r^k−1^,r^k−1⁾

k−2,r_k−2)

ii p

k

= −r

_k−1

+ e

_k−1

p

_k−1

iii ρ

_k

=

^(r_(p^k−1^,r^k−1⁾

k,Ap_k)

iv u

k

= u

_k−1

+ ρ

k

p

k

v r

k

= r

_k−1

+ ρ

k

Ap

k

Abbruchkriterium: ∣∣r

_k

∣∣ ≤ T OL∣∣r

0

∣∣

Aufwand: O(n) Iterationen (typisch aber: √

n). Jeder Schritt O(n

²

) Operationen, O(n) bei sparse-Matrix.

Absch¨ atzung /obere Schranke: ∣∣u

k

−x∣∣ ≤ 2α

^k

∣∣u

0

−x∣∣ mit α =

√κ_A−1

√κ_A+1

/ Gesucht: obere Schranke der Anzahl Iterationen k, um den absoluten Anfangsfehler um den Faktor c zu reduzieren, mit κ

A

gegeben.

k =

^ln(0.5⋅c)_ln(α)

1

f u n c t i o n[ x , it , t i m e ]= m y _ c g ( A , b , u0 , maxit , tol )

2

tic; % s t a r t m e a s u r i n g t i m e

3

u = u0 ; % s t a r t w i t h i n i t i a l g u e s s

4

r0 = A * u0 + b ;

5

p = - r0 ;

6

r = r0 ;

7

it =0;

8

w h i l e ( it < m a x i t )

9

it = it +1;

10

rho =( r ’* r ) /( p ’( A p ) ) ;

11

u = u + rho * p ;

12

r _ o l d = r ;

13

r = r + rho ( A p ) ;

14

if ( n o r m ( r ) < tol *n o r m ( r0 ) )

15

b r e a k; % c o n v e r g e d

16

end

17

e =( r ’* r ) /( r_old ’* r _ o l d ) ;

18

p = - r + e * p ;

19

end

20

x = u ;

21

t i m e = toc; % s t o p m e a s u r i n g t i m e

22

if ( it == m a x i t )

23

e r r o r ([’ no c o n v e r g e n c e a f t e r ’ , i n t 2 s t r( it ) , ’ i t e r a t i o n s ’]) ;

24

end

>> [x,flag,relres,it]=pcg(A,b,tol,maxit) L¨ ost Ax = b!

Vorkonditionierung: Mit unvollst¨ andiger Cholesky-Zerlegung die Kondition von A verkleinern.

>> H=cholinc(A,1e-3)

>> [x,flag,relres,it]= pcg(A,b,tol,maxit,H’,H)

3 Nichtlineare Gleichungssysteme

3.1 Eine Gleichung in einer Unbekannten

Gegeben: Stetige Funktion y = f(x).

Gesucht: L¨ osung x

^∗

der Gleichung f(x) = 0 (Nullstelle).

F¨ ur ∣f

^′

(x∗)∣ ≪ 1 ist das Problem schlecht konditioniert:

x∗ = f

⁻¹

(0) = H(0) κ

H

= ∣

^H

′(0)

H(0)

∣ = ∣

_x∗f_′¹_(x∗)

∣

Die Konvergenz der Fixpunkt-Iteration und des Newton-Verfahrens werden mit dem Banachschen Fixpunktsatz untersucht. Sie konvergieren nur wenn der Fixpunktsatz erf¨ ullt ist.

3.1.1 Fixpunkt-Iteration

Gegeben: y = f (x), x ∈ [a, b]

Gesucht: L¨ osung der Fixpunktgleichung x = F(x) ⇒ x

_k+1

= F (x

k

)

Abbruchkriterium: ∣x

_k+1

− x

k

∣ ≤ ∣x

_k+1

∣ ⋅ RT OL + AT OL

Die Fixpunktgleichung erh¨ alt man durch Umformung (Addition, usw.) der Funktionsgleichung f(x) = 0. Die Fixpunktgleichungen sind nicht eindeutig, es gibt meistens mehrere verschiedene Arten die Gleichung darzustellen.

3.1.2 Newton-Verfahren (1D)

Gegeben: Stetige Funktion y = f(x)

Gesucht: L¨ osung x

^∗

der Gleichung f(x) = 0 (Nullstelle)

Idee: Ersetze Graph von f(x) durch Tangente in P (x

0

, f(x

0

))

Newton-Iteration in 1D . . . .

Start: x

0

, f

^′

(x), RT OL, AT OL

Vorgehen: x

_k+1

= F (x

k

) = x

k

−

^f(x^k⁾

f^′(xk)

Abbruchkriterium: ∣x

_k+1

− x

k

∣ ≤ ∣x

_k+1

∣ ⋅ RT OL + AT OL

⊞ schnell

⊟ lokal, teuer, f

^′

(x) muss bekannt sein

3.1.3 Banachscher Fixpunktsatz Folgende Bedingungen m¨ ussen erf¨ ullt sein:

F bidet das Intervall I = [a, b] in sich ab. Umkehrung m¨ oglich.

F ist eine Kontraktion: ∃ 0 < L < 1 L ≥ ∣

^F(x

1)−F(x²) x¹−x²

∣ Falls F

^′

stetig ist auf I gen¨ ugt ∣F

^′

∣ ≤ L < 1

▸ Falls F

^′′

in ganz I das gleiche Vorzeichen hat, nur ∣F

^′

(a)∣,∣F

^′

(b)∣

betrachten, sonst auch Maxima @ F

^′′

= 0 Daraus folgt dann:

F hat genau einen Fixpunkt x in [a, b]

Es gelten folgende Fehlerabsch¨ atzungen

∣x − x

_k

∣ ≤

^L

k

1−L

∣x

1

− x

0

∣ a priori

∣x − x

k

∣ ≤

_1−L^L

∣x

k

− x

_k−1

∣ a posteriori

Die a priori Absch¨ atzung wird vor allem dazu benutzt nach nur einer Iteration bereits eine obere Schranke f¨ ur die Zahl der Iterationen ang- zugeben. Die ben¨ otigt wird um x bis auf einen Fehler ε zu bestimmen.

∣x − x

k

∣ ≤

^L

k

1−L

∣x

1

− x

0

∣ ≤ ε ⇒ k ≥

ln_∣x^ε(1−L) 1−x0∣ lnL

Es gibt anziehende- und abstossende Fixpunkte

∣F

^′

(x)∣ < 1, so ist x anziehend

∣F

^′

(x)∣ > 1, so ist x abstossend

∣F

^′

(x)∣ = 1, so kann x anziehend, abstossend oder keines von beiden sein Je kleiner ∣F

^′

(x)∣ ist, desto schneller konvergiert das Verfahren. Die selben Bedingungen gelten f¨ ur die inverse Iteration f¨ ur (F

⁻¹

)

^′

.

3.1.4 Konvergenzordnung p

∣x − x

_n+1

∣ ≅ C∣x − x

n

∣

^p

CG-Verfahren: p = 1 Fixpunktiteration: p = 1 Newton-Verfahren: p = 2 Sekantenmethode: p = 1.6 Bisektionsverfahren: p = 1

3.1.5 Effizienzindex E

E = p

^1/m

wobei m die Zahl der Funktionsauswertungen in jedem Schritt bedeutet.

3.1.6 Die Sekantenmethode

Sekantenmethode . . . .

Start: x

0

, x

1

, RT OL, AT OL

Vorgehen: x

_k+1

=

^x^k−1_f(x^f^(x^k^)−x^k^f(x^k−1⁾

k)−f(xk−1)

Abbruchkriterium: ∣x

_k+1

− x

k

∣ ≤ ∣x

_k+1

∣ ⋅ RT OL + AT OL

⊞ braucht keine Ableitung, billiger als Newton

⊟ lokal, langsamer als Newton

(3)

3.1.7 Das Bisektionsverfahren

Gegeben: f(x), stetig im Intervall [a, b] und f (a) ⋅ f(b) < 0

Bisektionsverfahren (Divide and Conquer) . . . .

Start: a

0

, b

0

, x

0

=

^a⁰^+b⁰

2

, x

1

, RT OL, AT OL

Vorgehen:

i if f(x

k

) = 0 then break

ii if f(a

_k

) ⋅ f(x

_k

) < 0 then a

_k+1

= a

_k

, b

_k+1

= x

_k

else a

_k+1

= x

k

, b

_k+1

= b

k

x

_k+1

=

^a^k+1^+b^k+1

2

Abbruchkriterium: ∣x

_k+1

− x

k

∣ ≤ ∣x

_k+1

∣ ⋅ RT OL + AT OL

⊞ global

⊟ langsame Konvergenz

3.2 Mehrere Gleichungen mit mehreren Unbekannten

Gegeben: f(x) = 0 mit x =

⎛

⎜ ⎜

⎜

⎝ x

1

. . x

n

⎞

⎟ ⎟

⎟

⎠ f (x) =

⎛

⎜ ⎜

⎜

⎝

f

1

(x1, . . . , x

n

) . . f

n

(x

1

, . . . , x

n

)

⎞

⎟ ⎟

⎟

⎠

=

⎛

⎜ ⎜

⎜

⎝ 0 . . 0

⎞

⎟ ⎟

⎟

⎠ Jacobi-Matrix: J(x) =

^∂f_∂x

=

⎛

⎜ ⎜

⎝

∂f₁

∂x₁

∂f₁

∂x₂

. . .

_∂x^∂f¹

n

.. ..

∂f_n

∂x1

∂f_n

∂x2

. . .

_∂x^∂fⁿ

n

⎞

⎟ ⎟

⎠

⇒ Linearisierung von f in x

0

3.2.1 Newton Verfahren (nD)

Das Newton-Verfahren ist immer lokal konvergent, falls det(J(x

^∗

)) ≠ 0 Es konvergiert quadratisch. Es gilt der Banachsche Fixpunktsatz, ersetze ∣.∣

durch ∣∣.∣∣

Newtonverfahren in R

ⁿ

. . . .

Start: x

⁰

, RT OL, AT OL

Vorgehen:

i J(x

^k

)∆

^k

= −f (x

^k

) ⇒ ∆

^k

= J(x

^k

)/ − f (x

^k

) ii x

^k+1

∶= x

^k

+ ∆

^k

Abbruchkriterium: ∣∣∆

^k

∣∣ ≤ ∣∣x

^k+1

∣∣ ⋅ RT OL + AT OL

1

f u n c t i o n[ x , it ]= m y _ n e w t o n ( x0 , rtol , atol , m a x i t )

2

it =0; x = x0 ;

3

w h i l e ( it < m a x i t )

4

d e l t a = - J ( x ) \ f ( x ) ;

5

x _ o l d = x ;

6

x = x _ o l d + d e l t a ;

7

it = it +1;

8

if ( n o r m ( d e l t a ) <= n o r m ( x ) * r t o l + a t o l )

9

b r e a k; % c o n v e r g e d

10

end

11

end

12

if ( it == m a x i t )

13

e r r o r ([’ no c o n v e r g e n c e a f t e r ’ , i n t 2 s t r( it ) , ’ i t e r a t i o n s ’]) ;

14

end;

15

end

16

f u n c t i o n[ f_ ]= f ( x ) % f ( x ) =0

17

x1 = x (1) ; x2 = x (2) ;

18

f_ =[ x1 + x2 ^3 ; x1 ^2+ x2 ];

19

end

20

f u n c t i o n[ J_ ]= J ( x )

21

x1 = x (1) ; x2 = x (2) ;

22

J_ =[1 3* x2 ^2 ; 2* x1 1];

23

end

⊞ quadratische Konvergenz

⊟ Finden eines guten Startvektors x

⁰

.

⊟ Berechnen der Jacobi-Matrix in jedem Schritt.

3.2.2 Quasi-Newton Verfahren (nD)

In jedem Schritt wird die gleiche Jacobi-Matrix J(x

⁰

) verwendet.

1

f u n c t i o n[ x , it ]= m y _ q u a s i n e w t o n ( x0 , rtol , m a x i t )

2

it =0; x = x0 ;

3

J0 = J ( x0 ) ; % < - -

4

w h i l e ( it < m a x i t )

5

d e l t a = - J0 \ f ( x ) ; % < - -

6

[ . . . ] % s a m e as m y _ n e w t o n

⊞ J nur noch einmal berechnen.

⊟ nur noch lineare Konvergenz

3.2.3 ged¨ ampftes Newton Verfahren (nD) Wie Newton, aber mit D¨ ampfungsfaktor α

k

.

x

^k+1

∶= x

^k

+ α

k

∆

^k

α

k,start

≪ 1 α

_k,ende=1

⊞ gr¨ osserer lokaler Konvergenzbereich durch weniger ¨ Uberschiessen.

⊟ Durch D¨ ampfung aber langsamer.

4 Ausgleichsrechnung

Gesucht wird ein Vektor x, so dass die dazugeh¨ origen Residuen r minimal sind.

4.1 Lineare Ausgleichsrechnung

4.1.1 mit Normalengleichungen nach Gauss i Fehlergleichungen: Ax − c = r

ii Normalengleichungen: A

^T

A

²

A^∗

x = A

^T

c >> x=(A’A)\(A’c)

⊟ A

^∗

meist schlecht konditioniert ⇒ QR-Zerlegung

4.1.2 mit QR-Zerlegung

QR-Zerlegung von A aus Ax − c = r . . . . i QR-Zerlegung von A: A ↝ R c ↝ d

A

_k+1

= U

_pq^T

⋅ A

k

c

_k+1

= U

_pq^T

⋅ c

k

} bis A quadratisch, dann A = R und c = d U

pq

=

p → q →

⎛

⎜

⎝

cos φ 0 sinφ

0 1 0

− sinφ 0 cos φ

⎞

⎟

⎠

sinφ = −

√ ^A(q,p) A(p,p)²+A(q,p)²

cosφ =

√ ^A(p,p) A(p,p)²+A(q,p)²

ii R¨ uckw¨ artseinsetzen: Rx = d nach x

>> [Q,R]=qr(A) >> d=Q’c >> x=R0\d0 mit R0 quadratisch.

⊞ κ

Gauss

= κ(A

^T

A) =

^µ_µ^max

min

κ

QR

= κ(A) =

√ κ(A

^T

A)

4.1.3 mit Singul¨ arwertzerlegung A = U SV

^T

Gegeben: Fehlergleichungen Ax − c = r

1

f u n c t i o n[ x , y ]= m y _ s v d ( A , c )

2

[ u , s , v ]= svd( A ) ;

3

k = r a n k ( A ) ;

4

[ rows , c o l s ]= s i z e ( A ) ;

5

y =[0 0] ’; % y to be a c o l u m n v e c t o r

6

for i =1: k

7

y ( i ) =1/ s ( i , i ) * u (: , i ) ’* c ;

8

end

9

for i = k +1: c o l s

10

y ( i ) =0; % a r b i t r a r y , c h o s e n to be 0

11

end

12

x = v * y ;

4.2 Nichtlineare Ausgleichsrechnung

Wie bei der linearen Ausgleichsrechnung ist die Kondition von A

^T

A zu gross, weshalb man besser mit den Fehlergleichungen rechnet.

4.2.1 Gauss-Newton-Methode

Fehlergleichungen: f(x) − c = r

Linearisierte Fehlergln.: A

⁰

°

J@x0

ξ + f

⁰

− c

´¹¹¹¹¹¸¹¹¹¹¹¶

d⁰

= ρ

⁰

⇒ ξ

¹

= A

0

/(ρ

⁰

− d

⁰

)

x

¹

= x

⁰

+ ξ

¹

und wieder von vorne...

⊟ Je nach x

⁰

keine Konvergenz → Minimierung.

4.2.2 Gauss-Newton-Methode mit Minimierung

Gleiches Vorgehen bis ξ

^k+1

, dann mit ̃ x = x

^k

+ tξ

^k+1

und t = 1,

¹₂

,

¹₄

, ...

testen, ob S(̃ x) < S(x

^k

) erf¨ ullt ist und bei ∣S(̃ x) − S(x

^k

)∣ < T OL den Vorschlag ̃ x akzeptieren: x

^k+1

= ̃ x. (S(x) = r

^T

r)

⊞ Konvergenz garantiert, dank S(̃ x) < S(x

^k

)

⊟ Langsame Konvergenz am Anfang m¨ oglich, aber am Ende doch wieder

fast quadratisch.

(4)

5 Interpolation und Approximation

5.1 Das Interpolationspolynom

Gegeben: n + 1 St¨ utzstellen x

0

, x

1

, . . . , x

n

und die St¨ utzwerte f

0

= f (x

0

), f

1

= f (x

1

), . . . , f

n

= f (x

n

).

Gesucht: Polynom ≤ n-ten Grades P

n

(x) = a

0

x

ⁿ

+a

1

x

ⁿ⁻¹

+ ⋅ ⋅ ⋅ +a

_n−1

x+

a

n

mit P

n

(x

0

) = f

0

, P

n

(x

1

) = f

1

, . . . , P

n

(x

n

) = f

n

5.1.1 Darstellung nach Lagrange

P

n

(x) =

n

∑

i=0

f

i

l

i

(x) l

i

(x) =

n

∏

j=0,j≠i

x − x

j

x

i

− x

j

i = 0, 1, . . . , n / Gegeben: 4 St¨ utzstellen und ihre Werte ⇒ n = 3

P

n

(x) = ∑

ⁿ_i=0

f

i

l

i

(x) (Falls f

i

= 0 ⇒ l

i

nicht berechnen) l

0

(x) =

_(x^(x−x¹^)(x−x²^)(x−x³⁾

0−x1)(x0−x2)(x0−x3)

l

1

(x) =

_(x^(x−x⁰^)(x−x²^)(x−x³⁾

1−x0)(x1−x2)(x1−x3)

l

2

(x) =

_(x^(x−x⁰^)(x−x¹^)(x−x³⁾

2−x0)(x2−x1)(x2−x3)

l

3

(x) =

_(x^(x−x⁰^)(x−x¹^)(x−x²⁾

3−x0)(x3−x1)(x3−x2)

5.1.2 Baryzentrsiche Darstellung

P

n

(x) = ∑

ⁿ_i=0

µ

i

f

i

∑

ⁿ_i=0

µ

i

λ

i

∶=

1 ∏

ⁿ_j=0,j≠i

(x

i

− x

j

) µ

i

∶=

λ

i

x − x

i

mit St¨ utzkoeffizienten λ

i

(Nenner von l

i

(x)) und i = 0,1, . . . , n.

⊞ Bei verschiedenen Polynomen weniger aufwendig als Lagrange.

λ

i

→ O(n

²

) (1x); P

n

(x) → O(n); bei Lagrange immer O(n

²

).

5.1.3 Fehler des Interpolationspolynoms f(x) − P

n

(x) = f

⁽ⁿ⁺¹⁾

(ξ)

(n + 1)!

n

∏

i=0

(x − x

i

) = O(h

ⁿ⁺¹

) f(x) − P

n

(x) ≤

^f⁽ⁿ⁺¹⁾_(n+1)!^(ξ)

∏

ⁿ_i=0

(x − x

i

) mit ξ = x wo Fehler maximal

5.2 Splineinterpolation

Die Funktion f(x) wird st¨ uckweise mit Polynomen 3.Grades Q

i

(t) inter- poliert. Q

i

(t) ∶= P

i

(x

i

+ h

i

t) mit t =

^x−xⁱ

h_i

und h

i

= x

i+1

− x

i

.

Q

i

(t) = f

i

(1 − 3t

²

+ 2t

³

) + f

_i+1

(3t

²

− 2t

³

) + h

i

f

_i^′

(t − 2t

²

+ t

³

) +h

i

f

_i+1^′

(−t

²

+ t

³

)

Q ˙

i

(t) = f

i

(6t+6t

²

) +f

_i+1

(6t−6t

²

) +h

i

f

_i^′

(1−4t+2t

²

) +h

i

f

_i+1^′

(−2t+3t

²

)

Q ¨

i

(t) = f

i

(−6 + 12t) + f

i+1

(6 − 12t) + h

i

f

_i^′

(−4 + 6t) + h

i

f

_i+1^′

(−2 + 6t)

...

Q

_i

(t) = 12f

i

− 12f

i+1

+ 6h

i

f

_i^′

+ 6h

i

f

_i+1^′

Q

i

(0) = f

i

Q

i

(1) = f

_i+1

Q ˙

i

(0) = h

i

f

_i^′

Q ˙

i

(1) = h

i

f

_i+1^′

.

Sind die Ableitungen bekannt, handelt es sich um die Hermite- Interpolation und im Algorithmus beginnt man bei Schritt 2.

Sind die Ableitungen nicht bekannt, fehlen diese Gleichungen, man stellt daher zus¨ atzliche Forderungen:

1.Forderung: P

_i^′′

(x

_i+1

) = P

_i+1^′′

(x

_i+1

) ⇒

Q¨_i(1) h²_i

=

Q¨_i+1(0) h²_i+1

2.Forderung “not a knot“:

P

1

(x) = P

2

(x)

P

_n−2

(x) = P

_n−1

(x) } ⇒

⎧ ⎪

⎪ ⎪

⎪

⎨

⎪ ⎪

⎩

Q

...

₁(1) h³₁

=

Q

...

₂(0) h³₂ Q

...

_n−2(1)

h³_n−2

=

Q

...

_n−1(0) h³_n−1

2.Forderung “nat¨ urliche Randbedingung“:

P

₁^′′

(x

1

) = P

_n−1^′′

(x

n

) = 0 ⇒

Q¨1(0) h²₁

=

Q¨_n−1(1) h²_n−1

= 0

2.Forderung “periodische Funktion“:

P

₁^′′

(x

1

) = P

_n−1^′′

(x

n

) ⇒

Q¨₁(0) h²₁

=

Q¨_n−1(1) h²_n−1

periodische Funktion ⇒ { f

1

= f

n

f

₁^′

= f

_n^′

Spline-Interpolationsfunktion g(x). . . . i Finde L¨ osung aus Af

^′

= d

↘ not a knot:

⎛

⎜ ⎜

⎜

⎜ ⎜

⎝ b

1 a1

2

b

1

a

1

b

2

b

2

a

2

b

3

... ... ...

b

n−2

a

n−2

b

n−1 a_n−2

2

b

n−1

⎞

⎟ ⎟

⎟

⎟ ⎟

⎠

⎛

⎜ ⎜

⎜

⎜ ⎜

⎝ f

₁^′

f

₂^′

f

₃^′

...

f

_n−1^′

f

_n^′

⎞

⎟ ⎟

⎟

⎟ ⎟

⎠

=

⎛

⎜ ⎜

⎜

⎜ ⎜

⎝ d

1

d

2

d

3

...

d

n−1

d

n

⎞

⎟ ⎟

⎟

⎟ ⎟

⎠

↘ nat¨ urliche RB:

⎛

⎜ ⎜

⎜

⎜ ⎜

⎝

4 2

b

1

a

1

b

2

... ... ...

b

_n−2

a

_n−2

b

_n−1

2 4

⎞

⎟ ⎟

⎟

⎟ ⎟

⎠

⎛

⎜ ⎜

⎜

⎜ ⎜

⎝ f

₁^′

f

₂^′

...

f

_n−1^′

f

_n^′

⎞

⎟ ⎟

⎟

⎟ ⎟

⎠

=

⎛

⎜ ⎜

⎜

⎝

6(f1−f2) h₁

d

2

...

d

_n−1

6(fn−f_n−1) h_n−1

⎞

⎟ ⎟

⎟

⎠

↘ periodische RB, h = const.:

⎛

⎜ ⎜

⎜

⎝

8 2 2

b a b

... ... ...

b a b

b b a

⎞

⎟ ⎟

⎟

⎠

⎛

⎜ ⎜

⎜

⎝ f

₁^′

f

₂^′

...

f

_n−2^′

f

_n−1^′

⎞

⎟ ⎟

⎟

⎠

=

⎛

⎜ ⎜

⎜

⎝

6

h

(f

2

− f

n−1

) d

2

...

d

n−2

d

_n−1

⎞

⎟ ⎟

⎟

⎠ h

i

∶= x

_i+1

− x

i

, i = 1, . . . , n − 1

b

i

∶=

¹

h_i

, i = 1, . . . , n − 1 a

i

∶= 2(b

i

+ b

_i+1

), i = 1, . . . , n − 2 c

i

∶=

^fⁱ⁺¹^−fⁱ

h²_i

, i = 1, . . . , n − 1 d

1

∶= 2c

1

+

_h^h¹

1+h2

(c

1

+ c

2

) (nur bei “not a knot“) d

i

∶= 3(c

_i−1

+ c

i

), i = 2, . . . , n − 1

d

n

∶= 2c

_n−1

+

^hⁿ⁻¹

h_n−1+hn−2

(c

_n−1

+ c

_n−2

) (nur bei “not a knot“) ii Bestimme den Index i, f¨ ur den gilt x

i

≤ x ≤ x

i+1

iii Setze t =

^x−xⁱ

h_i

iv Berechne Qi(t) = α

0

+ [β

0

+ (γ

0

+ δ

0

t)(t − 1)]t , wobei β

−1

= h

i

f

_i^′

α

0

= f

i

γ

0

= β

0

− β

₋₁

β

0

= α

1

− α

0

δ

0

= γ

1

− γ

0

α

1

= f

_i+1

γ

1

= β

1

− β

0

β

1

= h

i

f

_i+1^′

v Setze g(x) = Q

i

(t).

Allgemeines Vorgehen beim L¨ osen eines Problems:

1 Entscheiden, welche 2.Forderung gew¨ ahlt wird ⇒ Gleichunssystem hin- schreiben.

2 Tabelle erstellen, um den Vektor der rechten Seite d

i

zu finden. Glei- chungssystem mit TR::rref() l¨ osen.

3 Aus f

_i^′

per Hermiteschen Algorithmus (ii bis v) die gesuchte Funktion g(x) finden, wobei g(x) = Qi(t) = α

0

+ [β

0

+ (γ

0

+ δ

0

t)(t − 1)]t mit TR::qi(α

0

, β

0

, γ

0

, δ

0

, t).

1

f u n c t i o n[ yy ]= m y _ s p l i n e _ n o t a k n o t ( x , fx , xx )

2

h = d i f f ( x ) ;

3

b = 1 . / h ;

4

a = 2 * ( b (1: end -1) + b (2: end ) ) ;

5

c = d i f f ( fx ) ./ h . ^ 2 ;

6

d = [ 2 * c (1) + h (1) /( h (1) + h (2) ) *( c (1) + c (2) ) ,...

7

3*( c (1: end -1) + c (2: end) ) ,...

8

2* c ( end) + h ( end ) /( h (end -1) + h ( end) ) * . . .

9

( c (end) + c ( end -1) ) ];

10

A = s p d i a g s ([[ b (1: end -1) ’; a ( end) / 2 ; 1 ] . . .

11

[ b (1) ; a ’; b (end ) ] . . .

12

[1; a (1) /2; b (1: end -1) ’]] ,...

13

-1:1 , l e n g t h( x ) , l e n g t h( x ) ) ;

14

f _ p r i m e = A \ d ’;

15

for k =1: l e n g t h( xx )

16

i = min( f i n d (( xx ( k ) - x ) <0) ) -1;

17

t =( xx ( k ) - x ( i ) ) / h ( i ) ;

18

a l p h a =[ fx ( i ) fx ( i +1) ];

19

b e t a =[ h ( i ) * f _ p r i m e ( i ) d i f f ( a l p h a ) ...

20

h ( i ) * f _ p r i m e ( i +1) ];

21

g a m m a = d i f f ( b e t a ) ;

22

d e l t a = d i f f ( g a m m a ) ;

23

yy ( k ) = a l p h a (1) +( b e t a (2) + . . .

24

( g a m m a (1) + d e l t a (1) * t ) ( t -1) ) t ;

25

end

>> [yy]= spline(x,fx,xx) ≡ my spline notaknot.m 5.2.1 Interpolationsfehler

Der Interpolationsfehler ist von der Ordnung O(h

⁴

). Verglichen mit dem Interpolationspolynom O(h

ⁿ⁺¹

) ist das viel besser, je gr¨ osser n ist. Beim Interpolationspolynom f¨ ur n gross hat es zwischen den St¨ utzstellen Extre- malwerte, die gegen den Rand und mit zunehmendem n immer gr¨ osser werden. Dieses Verhalten zeigt die Splineinterpolation nicht.

5.3 Trigonometrische Interpolation

Diskrete Fourier-Transformation (DFT) und Fast Fourier Transformation (FFT) spielen eine wichtige Rolle in der Signal- und Bildverarbeitung, im besonderen auch bei der Daten-Kompression.

Problemstellung: Es sind N St¨ utzstellen x

k

sowie deren St¨ utzwerte f

k

einer 2π-periodischen Funktion gegeben.

Gesucht: ein trigonometrisches Polynom

p(x) =

^a₂⁰

+ ∑

^N/2−1_j=1

[a

j

cos(jx) + b

j

sin(jx)] +

^a^N/2₂

cos (

^N

2

x)

(5)

5.3.1 DFT

Fourier-Reihe von f : f (x) = ∑

^∞_j=−∞

c

j

e

^ijx

c

j

=

¹

2π

∫

2π

0

f (x)e

^−ijx

dx Dieses Integral kann nat¨ urlich nicht gel¨ ost werden, da die Funtkion f(x) nur an N Stellen bekannt ist ⇒ Approximation des Integrals durch Tra- pezn¨ aherung c ˜

j

=

_N¹

∑

^N−1_l=0

f (x

_l

)e

^−ijx^l

, mit x

_l

=

^l2π_N

.

Analyse: ˜ c

^(N)

=

¹

N

W f

^(N)

Synthese: f

^(N)

= W ˜ c

^(N)

mit f

^(N)

=

⎛

⎜

⎝ f

0

⋮ f

_N−1

⎞

⎟

⎠

˜ c

^(N)

=

⎛

⎜

⎝

˜ c

0

⋮

˜ c

_N−1

⎞

⎟

⎠

und W =

⎛

⎜ ⎜

⎝

1 1 1 ⋯ 1

1 w w

²

⋯ w

^N−1

1 w

²

w

⁴

⋯ w

^2(N−1)

⋮ ⋮ ⋮ ⋱ ⋮

1 w

^N−1

w

^2(N−1)

⋯ w

^(N−1)²

⎞

⎟ ⎟

⎠

, wobei w = e

⁻ⁱ^2π^N

und W das konjugiert komplexe von W . trigonometrische Interpolationspolynome

komplex: p(x) ∶= ∑

^N/2+1_j=−N/2+1

˜ c

j

e

^ijx

+ ˜ c

_N/2

cos(

^N₂

x)

2π-periodisch:

⎧ ⎪

⎪ ⎪

⎨

⎪ ⎪

⎪ ⎩

p(x)

3

∶= ˜ c

0

+ c ˜

1

e

^ix

+ ˜ c

2

e

^−ix

p(x)

5

∶= ˜ c

0

+ c ˜

1

e

^ix

+ ˜ c

2

e

^−ix

+ ˜ c

3

e

^2ix

+ ˜ c

4

e

^−2ix

. . .

reell: p(x) ∶=

^a₂⁰

+ ∑

^N/2+1_j=1

[a

j

cos(jx) + b

j

sin(jx)] +

^a^N/2₂

cos(

^N₂

x) Allgemeines Vorgehen:

1 Diskrete Fourier-Analyse ˜ c

^(N)

=

_N¹

W f

^(N)

. 2 Auswertung des Polynoms p(x).

/ Bestimme das interpolierende trigonometrische Polynom von der 2π- periodischen Funktion mit x

i

0

^2π₃ ^4π₃

f

i

f

0

f

1

f

2

˜ c

⁽³⁾

=

¹

3

W f

⁽³⁾

=

¹

3

⎛

⎜

⎝

1 1 1

1 w w

²

1 w

²

w

⁴

⎞

⎟

⎠

⎛

⎜

⎝ f

0

f

1

f

2

⎞

⎟

⎠

=

⎛

⎜ ⎜

⎜

⎝

f0+f1+f2 3 f₀+wf1+w²f₂

3 f₀+w²f₁+w⁴f₂

3

⎞

⎟ ⎟

⎟

⎠ mit

⎧ ⎪

⎪ ⎪

⎨

⎪ ⎪

⎪ ⎩

w = e

⁻ⁱ^2π^N

= e

⁻ⁱ^2π³

= −

¹

2

−

√3 2

i w

²

= e

⁻²ⁱ^2π^N

= e

⁻ⁱ^4π³

= −

¹

2

+

√3 2

i w

⁴

= w = e

⁻ⁱ^2π³

= −

¹₂

−

√3 2

i

˜ c

⁽³⁾

=

⎛

⎜ ⎜

⎝

f₀+f1+f2 3

f₀+f1(−¹₂−^√₂³i)+f2(−¹₂+^√₂³i) 3

f₀+f1(−¹₂+^√₂³i)+f2(−¹₂−^√₂³i) 3

⎞

⎟ ⎟

⎠

p(x) = ˜ c

0

+ ˜ c

1

e

^ix

+ ˜ c

2

e

^−ix

= ˜ c

0

+ (˜ c

1

+ c ˜

2

) cos(x) + i(˜ c

1

− ˜ c

2

)sin(x) mit

⎧ ⎪

⎪ ⎪

⎨

⎪ ⎪

⎪ ⎩

˜

c

1

+ ˜ c

2

=

^2f⁰^−f¹^−f²

3

˜

c

1

− ˜ c

2

=

⁻ⁱ

√3f₁+i√ 3f₂ 3

i(˜ c

1

− ˜ c

2

) =

√3f₁−√ 3f₂ 3

p(x) =

^f⁰^+f₃¹^+f²

+

^2f⁰^−f¹^−f²

3

cos(x) +

√3f1−√ 3f2 3

sin(x)

5.3.2 FFT

Sei γ

^(N)

= W f

^(N)

= N˜ c

^(N)

und N = 2n , womit w

^N

= 1 und w

ⁿ

= −1.

gerade approximative Fourier-Koeffizienten: γ

2k

= ∑

ⁿ⁻¹_j=0

(w

²

)

^kj

(f

j

+ f

j+n

) ungerade approx. Fourier-Koeffizienten: γ

_2k+1

= ∑

ⁿ⁻¹_j=0

(w

²

)

^kj

(f

j

−f

_j+n

)w

^j

Die Anzahl komplexer Multiplikationen betr¨ agt µ =

^N

2

log

₂

(

^N

2

).

Wir betrachten den Fall N = 2

^m

, hier konkret: m = 3, N = 8, n = 4.

Es gilt: w = e

⁻ⁱ^2π^N

, w

⁸

= 1, w

⁴

= −1

W

^(N)

=

⎛

⎜ ⎜

⎜

⎜ ⎜

⎜

⎝ γ

0

γ

2

γ

4

γ

6

γ

1

γ

3

γ

5

γ

7

⎞

⎟ ⎟

⎟

⎟ ⎟

⎟

⎠

=

⎛

⎜ ⎜

⎜

⎜ ⎜

⎜

⎝

1 1 1 1 1 1 1 1

1 w

²

w

⁴

w

⁶

1 w

²

w

⁴

w

⁶

1 w

⁴

1 w

⁴

1 w

⁴

1 w

⁴

1 w

⁶

w

⁴

w

²

1 w

⁶

w

⁴

w

²

1 w w

²

w

³

w

⁴

w

⁵

w

⁶

w

⁷

1 w

³

w

⁶

w w

⁴

w

⁷

w

²

w

⁵

1 w

⁵

w

²

w

⁷

w

⁴

w w

⁶

w

³

1 w

⁷

w

⁶

w

⁵

w

⁴

w

³

w

²

w

⎞

⎟ ⎟

⎟

⎟ ⎟

⎟

⎠ Das ganze faktorisiert ergibt: W ̂

^(N)

= ( W

⁽ⁿ⁾

0 0 W

⁽ⁿ⁾

) D

^(N)

W ̂

^(N)

=

⎛

⎜ ⎜

⎜

⎜ ⎜

⎜

⎜ ⎜

⎝

1 1 1 1

1 w

²

w

⁴

w

⁶

1 w

⁴

1 w

⁴

1 w

⁶

w

⁴

w

²

1 1 1 1

1 w

²

w

⁴

w

⁶

1 w

⁴

1 w

⁴

1 w

⁶

w

⁴

w

²

⎞

⎟ ⎟

⎟

⎟ ⎟

⎟

⎟ ⎟

⎠

⎛

⎜ ⎜

⎝

1 1

1 w

⁴

w w

⁵

w

²

w

⁶

w

³

w

⁷

⎞

⎟ ⎟

⎠ mit

⎧ ⎪

⎪ ⎪

⎨

⎪ ⎪

⎪ ⎩ w

⁴

= −1 w

⁵

= −w w

⁶

= −w

²

w

⁷

= −w

³

Die Matrixmultiplikation z = D

^(N)

f

^(N)

l¨ asst sich einfach durchf¨ uhren. Die ubriggebliebene Rechnung ¨ ̂ γ

^(N)

= ( W

⁽ⁿ⁾

0 0 W

⁽ⁿ⁾

) z

^(N)

l¨ asst sich wie- derum auf zwei diskrete Fourier-Transformationen der Ordnung n = 4, im n¨ achsten Schritt weitere R¨ uckf¨ uhrung auf n = 2 und dann n = 1

1

f u n c t i o n []= m y _ f f t

2

N = 2 0 ; m =5; M = N *2^ m ;

3

% N b i g - - > m o r e i n i t i a l d a t a

4

% m , M b i g - - > m o r e i n t e r p o l a t i o n p o i n t s

5

x_k =2pi/ N [ 0 : N -1] ’; % i n i t i a l : n p o i n t s

6

f_k =sin (4* x_k ) . ^ 2 ; % g i v e n ( e x a m p l e )

7

c = fft( f_k / N ) ; % a n a l y s i s

8

c _ s t a r =[ c (1: N /2) ;z e r o s ( M - N ,1) ; c ( N / 2 + 1 : end ) ];

9

x _ s t a r =2pi/ M [ 0 : M -1] ’; % n e w : m > n p o i n t s

10

f _ s t a r = M *i f f t ( c _ s t a r ) ; % s y n t h e s i s

⊞ Aufwand fft: O(N log

₂

N) Aufwand dft: O(N

²

)

6 Numerische Integration und Differentiation

6.1 Numerische Integration - Quadratur

Bei der numerischen Integration soll das Integral I = ∫

_a^b

f (x)dx mit ̃ I approximiert werden. Es werden nun einige dieser Verfahren vorgestellt.

6.1.1 Trapezmethode

Idee: Unterteilung der Funktion in n Intervalle der L¨ ange h =

^b−a

n

und approximiere linear mit T (h) =

^h₂

[f

a

+ f

_b

]

I ̃ = T (h) = h

2 [f

0

+ f

1

] + h

2 [f

1

+ f

2

] + ⋯ + h

2 [f

_n−1

+ f

n

]

= h [ 1

2 f

0

+ f

1

+ ⋯ + f

_n−1

+ 1 2 f

n

]

Das Trapezverfahren approximiert mind. Polynome des 1. Grades exakt.

Absch¨ atzung: ∣I − T (h)∣ ≤

^M₁₂

(b − a)h

²

, M ≥ ∣f

^′′

(x)∣, x ∈ [a, b]

Das Trapezverfahren hat Ordnung 2.

Trapezverfahren mit Schritthalbierung . . . .

Start: a, b, RT OL, AT OL

1.Schritt:

i h

0

= b − a ii s

0

=

¹

2

(f(a) + f(b)) iii T

0

= s

0

h

0

iv N

0

= 1

weitere Schritte: (n = 0, 1,2, . . .) i h

n+1

=

^hⁿ

2

ii s

n+1

= s

n

+ ∑

^N_j=1ⁿ

f(a + (2j − 1)h

n+1

) iii T

_n+1

= h

_n+1

s

_n+1

iv N

n+1

= 2N

n

Abbruchbedingung: ∣T

n+1

− T

n

∣ ≤ ∣T

n+1

∣ ⋅ RT OL + AT OL

6.1.2 Simpson’sche Formel

Idee: Unterteilung der Funktion in 2n Intervalle der L¨ ange h =

^b−a

2n

und

approximiere quadratisch mit S(h) =

^h₃

[f

a

+ 4f

c

+ f

_b

]

(6)

I ˜ = S(h) = h

3 [f

0

+ 4f

1

+ f

2

] + h

3 [f

2

+ 4f

3

+ f

4

] + ⋅ ⋅ ⋅ +

h

3 [f

_2n−2

+ 4f

_2n−1

+ f

2n

]

= h

3 [f

0

+ 4f

1

+ 2f

2

+ 4f

3

+ 2f

4

+ ⋅ ⋅ ⋅ + 2f

_2n−2

+ 4f

_2n−1

+ f

2n

]

Die Simpson-Methode approximiert mindestens Polynome des 3. Grades exakt.

Absch¨ atzung: ∣I − S(h)∣ ≤

₁₈₀^M

(b − a)h

⁴

, M ≥ ∣f

⁽⁴⁾

(x)∣, x ∈ [a, b]

Die Simpson-Methode hat Ordnung 4.

Simpson-Methode mit Trapezwerten . . . .

Start: a, b, RT OL, AT OL

Berechnung der Trapezwerte: T

0

, T

1

, T

2

, . . .

Berechnung der Simpsonwerte:

S

1

=

^4T¹₃^−T⁰

, S

2

=

^4T²₃^−T¹

, S

3

=

^4T³₃^−T²

, . . .

Abbruchbedingung: ∣S

n+1

− S

n

∣ ≤ ∣S

n+1

∣ ⋅ RT OL + AT OL

6.1.3 Romberg-Verfahren

Idee: Durch Kombination von verschiedenen Trapezwerten kann man Fehlerterme tiefer Ordnung eliminieren und so die Ordnung des Verfah- rens erh¨ ohen → Simpsonwerte. Mit diesen verf¨ ahrt man genauso...:

R

0,0

= T (h) = I + c

1

h

²

+c

2

h

⁴

+ . . . T (

^h

2

) = I +

¹

4

c

1

h

²

+

¹

16

c

2

h

⁴

+. . . R

1,1

= S

1

= S (

^h

2

) =

^4T(

h 2)−T(h)

3

=

T(^h₂)−4⁻¹T(h) 1−4⁻¹

= I −

¹

4

c

2

h

⁴

+ O(h

⁶

) R

2,2

=

^16S(

h 4)−S(^h₂)

15

=

S(^h₄)−4⁻²S(^h₂)

1−4⁻²

= I + O(h

⁶

)

Romberg-Verfahren . . . .

Start:

h

0

= b − a, h

1

=

^h⁰

2

, h

2

=

^h¹

2

, . . . RT OL, AT OL

Romberg-Verfahren:

R

j,0

= T

j

R

j,k

=

^R^j,k−1⁻⁴

−kR_j−1,k−1 1−4^−k

Abbruchbedingung: ∣R

j,j

− R

j,j−1

∣ ≤ ∣R

j,j

∣ ⋅ RT OL + AT OL

6.1.4 Guass’sche Quadraturformeln

Die Gauss’sche Quadraturformel gilt f¨ ur das Integral I = ∫

₋₁¹

f(x)dx. Das allgemeine Integral I = ∫

_a^b

g(x)dx muss daher in diese Form transformiert werden.

Die n-Punkt Quadraturformel lautet:

I ˜ = Q

n

=

n

∑

j=1

w

j

f

j

, wobei die Funtionswerte f

j

und die Gewichte w

j

f(x) = b − a 2 g (

b − a 2 x +

a + b 2 ) w

j

= ∫

1

−1

⎡ ⎢

⎢ ⎢

⎢ ⎣

n

∏

k=1,k≠j

( x − x

k

x

j

− x

k

)

2

⎤

⎥

⎥ ⎥

⎥ ⎦

dx > 0 j = 1,2, . . . , n

⊞ Eine n-Punkt-Quadraturformel hat Genauigkeitsgrad von (2n-1), wenn die Nullstellen der Legende-Polynome verwendet werden.

⊞ Es werden viel weniger Funktionsaufrufe ben¨ otigt.

⊟ Aber daf¨ ur m¨ ussen viele Funktionswerte bekannt sein.

6.2 Numerische Differentiation

Notwendig ist numerische Differentiation, wenn nur einige Werte der Funk- tion bekannt sind. Mit der Talyorentwicklung um z:

f (z + h) = f(z) +

_1!^h

f

^′

(z) +

^h_2!²

f

^′′

(z) + ⋅ ⋅ ⋅ +

^h_k!^k

f

^(k)

(z) f (z − h) = f(z) −

_1!^h

f

^′

(z) +

^h_2!²

f

^′′

(z) − . . .

...lassen sich Approximationen f¨ ur f

^′

(z) herleiten:

Vorw¨ arts-Differenzenquotient:

^f^(z+h)−f_h ^(z)

= f

^′

(z) +

^h₂

f

^′′

+ O(h

²

) R¨ uckw¨ arts-Differenzenquotient:

^f^{(z)−f(z−h)}_h

= f

^′

(z) −

^h₂

f

^′′

+ O(h

²

) Symmetrische Differenzquotienten der Ordnung 1 , 2 und n:

∆

¹_h

(z) = f (z + h) − f(z − h)

2h = f

^′

(z) + h

²

6 f

^′′′

(z) + O(h

³

)

∆

²_h

(z) = f (z + h) − 2f (z) + f(z − h)

h

²

= f

^′′

+ O(h

²

)

∆

ⁿ_h

(z) = 1 (2h)

ⁿ

n

∑

j=0

(−1)

^j

a

j

f (z + b

j

h) a

j

= ( n

j ) b

j

= n − 2j mit dem Binomialkoeffizienten a

j

TR::ncr(n,j).

F¨ ur grosse h ist das Ergebnis schlecht, doch auch zu kleine h resultieren in schlechten Ergebnissen aufgrund Ausl¨ oschung.

h

optimal

≅

√

10

^−d

⋅ 2

_∣f^∣f(z)∣_′′_(z)∣

, d-stellige Gleitpunktarithmetik

7 Gew¨ ohnliche Differentialgleichungen

Gegeben: Die DGL x ˙ = f (t, x), sowie der Anfangswert x(t

0

) = x

0

. Gesucht: Die Approximation an der Endstelle ̃ x(t

f

).

7.1 Explizite Einschrittverfahren 7.1.1 Das explizite Eulerverfahren

Wir approximieren x(t + h) durch das Taylorpolynom vom Grad 1:

x(t) + x(t)h. Per Definition ist ˙ x(t) = ˙ f (t, x(t)) und somit F (t, x, h) ∶= x + hf (t, x).

Das Eulerverfahren ist wie folgt definiert:

̃ x

_j+1

= F (t

j

,̃ x

j

, h

j

) = ̃ x

j

+h

j

f (t

j

,̃ x

j

) x ˜

0

= x0 j = 0, 1, . . . , n − 1

Bei ¨ aquidistanten Schritten gilt: h ∶=

^t^f^−t⁰

n

t

j

∶= t

0

+ jh.

1

f u n c t i o n[ x _ t f ]= m y _ e e u l e r _ m a i n

2

f1 = @ ( x , t ) x ^2+ t ^2; % f u n c t i o n h a n d l e

3

h = 0 . 1 ; % ^ - - on f = x p r i m e

4

x0 =0;

5

t0 =0; tf =1;

6

n =( tf - t0 ) / h ;

7

x _ t f = m y _ e e u l e r ( f1 , t0 , x0 , h , n ) ;

1

f u n c t i o n[ x ]= m y _ e e u l e r ( f , t0 , x0 , h , n )

2

x = x0 ;

3

t = t0 ;

4

for it =1: n

5

x = x + h * f ( x , t )

6

t = t + h ;

7

end

Lokaler Fehler bei j + 1: l(t

j

,̃ x

j

, h) ∶= F (t

j

,̃ x

j

, h)

´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶

̃ x_j+1

−x(t

j

+ h)

wobei x(t

j

+ h) die L¨ osung an der Stelle t

j

+ h ist.

Globaler Fehler f¨ ur t = t

j

: ̃ x

j

− x(t

j

)

Die Fehlerordnung p (Ordnung des globalen Fehlers) eines Verfahrens mit Verfahrensfunktion F (t, x, h) ist gleich der Anzahl ¨ ubereinstimmender Ter- me in den Taylorentwicklungen nach h von x(t +h) und von F (t, x(t), h).

Der lokale Fehler des Verfahrens ist von der Ordnung O(h

^p+1

), beim expliziten Euler O(h

²

).

Der globale Fehler des Verfahrens ist von der Ordnung O(h

^p

), beim expliziten Euler O(h). Ein k-fach kleinerer Schritt h ergibt einen k-fach kleineren globalen Fehler.

7.1.2 Taylorverfahren h¨ oherer Ordnung Ordnung 2: Wir approximieren x(t + h) durch:

x(t) + x(t)h ˙ + x(t) ¨

^h₂²

(Taylorpolynom vom Grad 2)

¨

x(t) = f

t

(t, x(t)) + f

x

(t, x(t)) x(t) = ˙ f

t

(t, x(t)) + f

x

(t, x(t))f(t, x(t)) F (t, x, h) = x + hf (t, x) +

^h₂²

[f

t

(t, x) + f

x

(t, x)f(t, x)]

Das Taylor-Verfahren der Ordnung 2 ist also wie folgt definiert:

̃ x

j+1

=F (t

j

, ̃ x

j

, h

j

)

=̃ x

j

+ h

j

f (t

j

,̃ x

j

) + h

²_j

2 [f

t

(t

j

,̃ x

j

) + f

x

(t

j

, ̃ x

j

)f(t

j

,̃ x

j

)]

⊞ p = 2, lokaler Fehler somit O(h

³

) und globaler Fehler ist O(h

²

).

⊟ Die partiellen Ableitungen f

t

und f

x

m¨ ussen bekannt sein.

Grad p: Wir approximieren x(t + h) durch:

x(t) + x(t)h ˙ + ⋅ ⋅ ⋅ + x

^(p)

(t)

^h_p!²

(das Taylorpolynom vom Grad p)

⊞ p = n, lokaler Fehler somit O(h

ⁿ⁺¹

) und globaler Fehler ist O(h

ⁿ

).

⊟ Die partiellen Ableitungen von f bis zur Ordnung n−1 m¨ ussen bekannt

sein.

Zusammenfassung Numerische Mathematik Patrik Rohner ::