Wissenschaftliches Rechnen - Zusammenfassung Patrick Pletscher 5. Oktober 2004 1. Endliche Arithmetik

Double Precision (64 Bits)

1.1. Reelle Zahlen ↔ Maschinenzahlen

Double Precision ist heute der Normalfall. +-+-+-+-+-+-+-+-+-+-+-+-+ | S | E | M | +-+-+-+-+-+-+-+-+-+-+-+-+ 0 1 11 63

Jeder Computer ist ein endlicher Automat, d.h. er kann nur • endlich viele Operationen durchf¨ uhren • endlich viele Zahlen speichern

Die Zahl ist:

Die reellen Zahlen R werden durch die endliche Menge der V = (−1)S · 1.M · 2E−1023 falls 0 < E < 2047 Maschinenzahlen M approximiert. Dabei wird jeweils ein ganzes Intervall (z.B. alle reellen Zahlen deren ersten 10 Ziffern anders geschrieben: gleich sind) auf eine Maschinenzahl abgebildet. V = (−1)S · (1 + M/252 )· 2E−1023 Darstellung einer Maschinenzahl Ausnahmen: a ˜∈M a ˜ = ±m· B ±e E = 2047, M 6= 0 V = N aN (not a number) m = D.D · · · D Mantisse E = 2047, M = 0, S = 0 V = Inf (= ∞) e = D · · · D: Exponent E = 2047, M = 0, S = 1 V = Inf (= −∞) B: Basis E = 0, M 6= 0 V = (−1)s · 0.M · 2−1022 V = (−1)S · M/252 · 2−102 Wobei D eine Ziffer ∈ {0, 1, . . . , 9} (nicht normalisierte Zahl) E = 0, M = 0, S = 1 V = −0 Heute ist B meistens 2 (nach dem IEEE Standard for floating E = 0, M = 0, S = 0 V =0 Point operations), fr¨ uher oftmals 10. Es gibt eine gr¨osste und eine kleinste Maschinenzahl. Der Rechenbereich ist [−realmax, realmax], wobei realmax = 21024 . Es gibt auch eine kleinste normalisierte Zahl realmin = 2−1022 . Die kleinste nicht normalisierte Zahl ist 0.0 · · · 01· 2−1022 .

1.2. Gleitkommazahlen nach dem IEEE Standard Single Precision (32 Bits) +-+-+-+-+-+-+-+-+-+-+-+-+ | S | E | M | +-+-+-+-+-+-+-+-+-+-+-+-+ 0 1 9 31

Maschinengenauigkeit eps

S: Vorzeichen E: Exponent 8 Bits M: Mantisse 23 Bits

V = (−1)S · 1.M · 2E−127

Definition (Maschinengenauigkeit). Die Maschinengenauigkeit eps ist die kleinste pos. Maschinenzahl welche zu 1 addiert ein Resultat 6= 1 ergibt. Eine neuere Definition lautet: Der Abstand zwischen zwei aufeinanderfolgenden Floating Point Numbers zwischen 1 und 2.

falls 0 < E < 255

Ausnahmen: E E E E

= 255, M 6= 0 = 255, M = 0, S = 0 = 255, M = 0, S = 1 = 0, M 6= 0

E = 0, M = 0, S = 1 E = 0, M = 0, S = 0

In IEEE ist eps = 2−52 . V = N aN (not a number) V = Inf (= ∞) V = Inf (= −∞) V = (−1)s · 0.M · 2−126 (nicht normalisierte Zahl) V = −0 V =0

1.3. Rechnen mit Maschinenzahlen Rundungsfehler absoluter Rundungsfehler: ra = c˜ − c relativer Rundungsfehler: r = rca

1

Wie gross sind die Rundungsfehler? Bei heutigen Computern if (1/x)/x == 0 gilt: Noch besser ohne nichtnormalisierten Zahlen: ˜ = (a ⊕ b)(1 + r) a⊕b if (eps/x)/x == 0

Wobei r relativer Rundungsfehler mit |r| ≤ ε und ⊕ ∈ {+, −, ×, /}.

numerische Ausl¨ oschung Assoziativgesetz

Wenn 2 nahezu gleichgrosse Zahlen subtrahiert werden. Viele Gesetze der Mathematik gelten nicht, so z.Bsp. Assoziativge- Bits der Mantisse sind dann unbekannt und werden mit 0 ullt. setz. Man muss also u ur eine Sum- gef¨ ¨berlegen wie klammern. F¨ me ist es also von Vorteil mit dem kleinsten Element zu beginnen und das gr¨osste Element erst als Letztes zu addieren. 1.4. Maschinenunabh¨ angige Algorithmen These: Gute Algorithmen funktionieren dank der Rundungsfehler!

Kahan’s Summenalgorithmus

Um Summen auf dem Computer m¨ oglichst richtig zu berechnen, benutzt man den Algorithmus von Kahan. Beispiel: Exponentialfunktion s=

N X

function y= e2(x); % stable computation of the exponential % function using the series if x tol * abs(tnew), told = tnew; zh = 2 * zh; h = h / 2; s = s + sum(feval(f,a + [1:2:zh]*h)); tnew = h * s; end; T = tnew;

ξn yn

und durch Integration erhalten wir die Quadratur Regel: Z b n X Pn (x)dx = wi yi a

h · sneu 2

i=0

Wir benutzen die Lagrange Form von Pn : Pn (x) =

n X

Simpsonregel li (x)· yi

Ein Polynom zweiten Grades

i=0 n Y x − ξj li (x) = ξ − ξj j=0 i b

Z

Pn (x)dx = a

n X

Z yi

i=0

b

Z

f (x)dx ≈ a

b−a 6



 f (a) + 4f

a+b 2



 + f (b)

F¨ ur die zusammengesetze Simpsonregel m¨ ussen wir das Integrationsintervall [a, b] in 2n Unterintervalle aufteilen und es ergibt sich:

b

li (x)dx | a {z } wi

h (y0 + 4y1 + 2y2 + 4y3 + . . . + 2y2n−2 + 4y2n−1 + y2n ) 3

Newton-Cotes Regeln

S(h) =

Gegeben: St¨ utzstellen ξi Rb Berechne Gewichte wi = a li (x)dx

Der Integrationsfehler ist

Z

b

f (x)dx ≈ a

n X

Z b b−a f (x)dx − S(h) = h4 f (4) (ξ) a 180

wi f (ξi )

i=0

closedcotes := n -> factor(int(interp( [seq(i*h, i= 0..n)], [seq(f(i*h), i=0..n)], z),z=0..n*h));

Zur Programmierung (Einf¨ ugen von neuen ξi ”dazwischen”): F¨ ur h:

Trapezregel Z

b

f (x)dx ≈ a

ξ ∈ [a, b]

b−a (f (a) + f (b)) 2

Als zusammengesetzte Regel f¨ ur das Intervall [a, b] mit n gleich grossen Intervallen (xi , xi+1 ) der L¨ ange h = xi+1 −xi = (b − a)/n, wobei yi = f (xi ), ergibt sich: Dann f¨ ur

4

s1 s2 s4 h 2:

= y0 + y4 = y2 = y 1 + y3

sneu 1 sneu 2 sneu 4

∆f (x) = f (x + 1) − f (x) =

= s1 = s2 + s4 = neue Fkt. Werte

h 2

3

∆f (i) (x) =

a

n 1 2 4

wi 1 1 7

1 4 32

1 12

f (k+i) (x)

∆F (x) = F (x + 1) − F (x) = In Matrixschreibweise:   f (x)  f 0 (x)      ..   . f =   f (m) (x)    .. .  1 1!

n

32

7

Fehler 1 (2) h3 12 f (ξ) 1 (4) h5 90 f (ξ) 8 h7 945 f (6) (ξ)

Resultat von Steffenson Z b Z b P (x)dx − f (x)dx = hp+1 · k· f (p) (ξ),

∞ X

f (k−1)

k=1

b−aX wi fi Pn (x)dx = ns i=0 ns 2 6 90

1 k!

Und f¨ ur das Integral:

 A= Name Trapezregel Simpsonregel Milne-Regel



1 k!

∆F (x) ∆f (x) .. .

   ∆f =    ∆f (m−1) (x)  .. .  1 1 1 ... 2! 3! 4! 1 1 1 ...  1! 2! 3!  .. .. .. . . .

       

∆f = Af Nach f aufgel¨ost: f = A−1 ∆f Dabei ist A folgende Matrix:  b1 b 2  b1 A−1 = 

2.3. Fehler der Newton-Cotes Formeln

a

∞ X k=1

(sneu + 2· sneu + 4· sneu 1 2 4 )

¨ Uberblick b

1 k!

F¨ ur Ableitungen gilt:

function S = simpson(f,a,b,tol); h = (b-a)/2; s1 = feval(f,a) + feval(f,b); s2 = 0; s4 = feval(f,a + h); snew = h *(s1 + 4 * s4)/3; zh = 2; sold = 2*snew; while abs(sold-snew)>tol*abs(snew), sold = snew; zh = 2*zh; h=h/2; s2 = s2 + s4; s4 = sum(feval(f,a + [1:2:zh]*h)); snew = h*(s1 + 2*s2 + 4*s4)/3; end S = snew;

Z

f (k) (x)

k=1

Dann ergibt sich f¨ ur Simpson: S(h/2) =

∞ X

ξ ∈ (a, b)

b3 b2 .. .

 ... ...   .. .

a

Dabei h¨angen p und k nicht von f ab. Die Gleichung gilt Dabei sind die Koeffizienten bi aus den Bernoullizahlen Bk zu berechnen durch: nat¨ urlich nur f¨ ur gen¨ ugend oft differenzierbare f . Bk−1 bk = (k − 1)! 2.4. Euler-Mac Laurin Summationsformel Die ersten Bernoullizahlen:

Suchen einer geschlossenen Funktion f¨ ur β X

f (i)

Idee: Wir suchen eine Summenfunktion s(x) mit der Eigenschaft ∆s(x) = s(x + 1) − s(x) = f (x)

β X

Dann gilt:

Z f (i)

i=α

f (i) = s(β + 1) − s(α)

Wollen s(x) als formale (= k¨ ummern uns nicht um Konvergenzbereich) Potenzreihe berechnen. Wir betrachten die Taylorreihe f¨ ur f (x + h) f¨ ur h = 1: f (x + 1) =

k=0

f (k) (x)

β

=

f (x)dx + α

+

i=α

∞ X

B2 =

1 , 6

B3 = 0

Durch Umformen ergibt sich die Euler-Mac Laurin Summenformel

i=α

β X

1 B1 = − , 2

B0 = 1,

f (α) + f (β) 2

∞ X B2j (2j−1) (f (β) − f (2j−1) (α)) (2j)! j=1

Folgerung F¨ ur Summen von Polynomen Pm (x) werden alle Ableitungen gr¨osser m gleich 0 und somit ist es hinreichend die rechte Summe bei der Euler-Mac Laurin Formel nur bis zum Term m/2 zu bilden.

1 k!

5

Spezialisierung f¨ ur Trapezsumme

hi−1 h0 = i 2 2 Somit ergibt sich das Romberg Schema:

∞ X B2j (2j−1) (g (b) − g (2j−1) (a))h2j g(t)dt + (2j)! a j=1 | {z } | {z } exaktes Integral Fehler der Trapezregel b

Z

T (h) =

Spezialisierung des ANS f¨ ur Romberg Integration z = 0,

xi = h2i ,

Tij =

Wir haben die asymptotische Entwicklung der Trapezregel erhalten

hi =

4−j Ti−1,j−1 − Ti,j−1 4−j − 1

Satz 2.1. In der 2. Kolonne des Rombergschemas stehen die Simpson Werte. In der 3. Kolonne die Milne’s Regel. b

Z

g(t)dt + c1 h2 + c2 h4 + . . . + cm h2m + Rm

T (h) =

Integrale u ur eine Periode am ¨ber periodische Funktionen f¨ Besten mit Trapez. Romberg bringt hier keine Verbesserung.

a

Wobei f¨ ur das Restglied Rm gilt: Rm =

function [I, T] = romberg(f,a,b,tol); h = b - a; intv =1; s = (feval(f,a) + feval(f,b)) / 2; T(1,1) = s * h; for i = 2:15 intv = 2*intv; h = h/2; s = s + sum(feval(f,a+[1:2:intv]*h)); T(i,1) = s * h; vhj = 1; for j = 2:i vhj = vhj/4; T(i,j) = (vhj*T(i-1,j-1) - T(i,j-1)) / (vhj -1); end; if abs(T(i,i)-T(i-1,i-1)) < tol * abs(T(i,i)), I = T(i,i); return end end warning([’limit of extrapolation steps reached. ’, ... ’Required tolerance may not be met.’]); I = T(i,i);

B2m+2 2m+2 h (b − a)(2m+2) (ξ) (2m + 2)!

2.5. Romberg Integration b

Z T (h) = |a

g(t)dt +c1 h2 + c2 h4 + c3 h6 + . . . {z }

=c0

Somit T (h) ≈ Pm (h2 ) Wir wollen c0 ≈ Pm (0) berechnen. Vorgehen: • interpolieren durch Polynom vom Grade i

2.6. Gauss Quadratur • evaluieren f¨ ur x = 0 ergibt N¨ aherung f¨ ur Integral.

Wollen St¨ utzstellen optimieren.

F¨ ur die Interpolation k¨ onnte ein Lagrange-Polynom verwenR1 det werden, besser ist es aber das Aitken-Neville Schema zu Transformation von ”normalem” Integral zu −1 : ben¨ utzen.   Z Z b b−a a+b b−a 1 f t+ f (x)dx = dt 2 2 2 −1 a Aitken Neville Schema Nun versucht man die Gewichte wk und St¨ utzstellen ξk optimal zu w¨ a hlen: x y x0 T00 Z 1 n X x1 T10 T11 f (x)dx ≈ wk f (ξk ) (1) −1 x2 T20 T21 T22 k=1 .. . ... ... Satz 2.2. Die Ordnung der Gauss Formel 1 ist ≤ 2n − 1 xi

Ti0

Ti1

Ti2

...

...

...

...

...

Tii ...

..

Die wi und ξi sind symmetrisch. Ungerade Monome werden automatisch richtig integriert, Gleichungen nur f¨ ur gerade Monome aufstellen.

.

Wobei gilt

F¨ ur n ≤ 6 k¨onnen wir die L¨osung noch analytisch bestimmen. Ti0 Tij j

 = yi  (xi −z)Ti−1,j−1 +(z−xi−j )Ti,j−1 = i = 1, 2, 3, . . . xi −xi−j  = 1, 2, . . . , i

gauss := proc(n) local w, xi, i, j, firsteq, eqns, sols, m; global res;

z ist dabei eine feste Stelle, an dem die Polynome im ANS ausgewertet werden.

firsteq := 2 * sum(w[i], i=1..m); if irem(n,2)=1 then firsteq := firsteq+w[0] fi;

6

restart; lanczos := proc(p, alpha, beta, n) local k, q, x; p := array(0..n); alpha := array (1..n); beta := array (1..n-1); p[0] := 1/sqrt(2); alpha[1] := int(x*p[0]*p[0], x=-1..1); q := (x - alpha[1]) * p[0]; for k from 2 to n do beta[k-1] := sqrt(int(q*q, x=-1..1)); p[k-1] := expand(q / beta[k-1]); alpha[k] := int(x*p[k-1]*p[k-1], x= -1..1); q := (x - alpha[k])*p[k-1] - beta[k-1]*p[k-2]; od; RETURN (NULL); end;

eqns := {2 = firsteq}; for j from 2 by 2 to 2*(n-1) do eqns := eqns union{ int(x^j, x=-1..1) = 2*sum(w[i]*xi[i]^j, i=1..m)}; od; if irem(n,2)=1 then sols := {w[0]} else sols := {} fi; for j from 1 to m do sols := union {w[j], xi[j]}; od; res := solve(eqns, sols); evalf(res); end;

Charakterisieren der ξi , wi Z 1 n X P2n−1 (x)dx = P2n−1 (ξi )wi

Bez¨ uglich dem Skalarprodukt 3 ergeben sich f¨ ur die αk und βk : k (2) αk = 0, βk = √ −1 i=1 4k 2 − 1 Annahme: P2n−1 (x) = Hn−1 (x)· Qn (x) + Rn−1 (x) βk−1 l¨asst sich relativ bequem durch das Skalarprodukt (βk pk , βk pk ) berechnen, wobei βk pk = (x − αk )pk−1 − Somit ergibt sich f¨ ur die LHS von 2: βk−1 pk−2 , dies berechnet man am Besten immer am Anfang der neuen Iteration f¨ ur die vorhergehende, also eigentlich Z 1 Z 1 Z 1 β . k−1 P2n−1 (x)dx = Hn−1 (x)Qn (x)dx + Rn−1 (x)dx −1

−1

−1

Golub-Welsh Algorithmus 1. Wir w¨ahlen f¨ ur Qn das orthogonale Polynom zum Intervall [−1, 1] mit dem Skalarprodukt (f, g) = Die 3-gliedrige Rekursionsformel von R1 Matrixschreibweise f (x)g(x)dx −1 R1 Mit dieser Wahl ist −1 Hn−1 (x)Qn (x)dx = 0 2 0

2. W¨ahlenPξi = Nullstellen des orthogonalen Polynoms Qn ⇒ wi Hn−1 (ξi )Qn (ξi ) = 0

B B x· B @

3. Die ξi sind fest also kann man die Gewichte nach Idee von Newton-Cotes bestimmen R1 ⇒ wi = −1 li (x)dx (li = Lagrange Polynome) R1 Pn ⇒ −1 Rn−1 (x)dx = i=1 wi Rn−1 (ξi )

|

p0 (x) p1 (x) . . . pn−1 (x) {z p(x)

1 C C C A

=

α1 6 β1 6 6 6 6 6 6 4

β2 .. . ..

..

.

3

.

αn−1 βn−1

}

βn−1 αn

{z

|

Tn

0 B B +B @

Satz 2.3. Die Gauss-Quadraturregel 2 hat die Ordnung 2n− 1. Die ξi sind Nullstellen der orthogonalen Qn und die wi sind die Integrale der Lagrangepolynome.

β1 α2 .. .

2.4 in Vektor-

0 . . . 0 βn pn (x)

0 p0 (x) 7 7 7 B p1 (x) 7B 7B . . 7@ . 7 5 pn−1 (x)

1 C C C A

}

1 C C C A

x· p(x) = Tn · p(x) + en · βn · pn (x)

n := 20; X := sort([fsolve(orthopoly[P](n,x)=0,x)]); Q := int(interp(X, [seq(y[i], i=1..n)],z),z=-1..1); Orthogonale Polynome

Die ξi sind die Nullstellen von pn (x), d.h. ξi p(ξi ) = Tn · p(ξi ) Die ξi sind also die Eigenwerte von Tn , da Tn symmetrische Matrix ist, sind die EW sehr gut konditioniert und reell.

Da die Gaussquadraturregel sagt, dass Polynome bis Grad 2n − 1 korrekt integriert werden, im Speziellen auch die orthogonalen Polynome pi (x) f¨ ur i = 0, 1, . . . , n − 1 und alle Integrale u Polynomen Null ergeben mit Ausnahme ¨ber den √ ur die Gewichte wi : von p0 (x), welches 2 ergibt, gilt f¨ Satz 2.4. Die orthogonalen Polynome p0 (x), p1 (x), . . . √ gen¨ ugen der 3-gliedrigen Rekursionsformel Pw = 2· e1 Polynome Pk (vom Grade k) sind orthogonal auf [−1, 1] wenn  Z 1 0 k 6= i (Pi , Pk ) = Pi (x)Pk (x)dx = (3) 1 k=i −1

k = 1, 2, . . . Wobei P die Eigenvektormatrix von T ist, dabei ist zu beachten, dass die Vektoren so normiert sind, dass die ersten KomMit αk = (xpk−1 , pk−1 ), βk = (xpk−1 , pk ) und √ mit der Initia- ponenten gleich √12 sind, da p0 (x) = √12 . Dies erreicht man lisierung β0 := 0, p−1 (x) := 0 und p0 (x) := 1/ 2. am besten in dem man von der EW-Zerlegung Tn = QDQT x· pk−1 (x) = βk−1 pk−2 (x)+αk pk−1 (x)+βk pk ,

7

1. Jede Kolonne von Q durch das erste Element dividieren. √ 2. Multiplizieren jede Kolonne mit 1/ 2.

Falls man das Integral nicht auf Maschinengenauigkeit eps integrieren m¨ochte, sondern nur bis zu einer Toleranz tol, kann man folgende Abbruchbedingung benutzen:

dies kann auch geschrieben werden als P = QV, Pw =

√

V = diag

1 1 √ ,... √ 2q1,1 2q1,n

2e1 kann nun gel¨ ost werden durch √ QVw = 2e1  w

 √ −1 T  2V Q e1 = 2  = 

2 q11 2 q12 .. . 2 q1n

tol tol · is + (i1 − i2) == · is eps eps

!

Um die Neu-Berechnung von Funktionswerten zu vermeiden, kann man fa , fm , fb als Parameter der rekursiv aufgerufenen Funktion u ¨bergeben. Als weitere Verbesserung kann man das Romberg-Schema auf die beiden mit der Simpsonmethode gerechneten i1, i2 anwenden 16· i2 − i1 i1neu = 15

    

function Q = simadpt(f,a,b,tol,trace,varargin) global warn1, warn2

Dadurch ergibt sich der Gauss-Legendre Algorithmus function [xi, w] = gausslegendre(n) beta = 0.5 ./ sqrt(1-(2*(1:n-1).^(-2))); [Q, D] = eig(diag(beta,1) + diag(beta,-1)); [x, i] = sort(diag(D)); w = 2 * Q(1,i).^2;

if (nargin < 4), tol = []; end; if (nargin < 5), trace = []; end; if (isempty(tol)), tol = eps; end; if (isempty(trace)), trace = 0; end; if tol a) a = a/2.0; end a = a*2.0; % kleinste pos. nicht-norm. b = 1.0; while ((b/2.0) > 0.0) b = b/2.0; end % gr"osste c = 2-e; while ((c*2.0) ~= Inf) c = c*2.0; end

I. Additionstheoreme sin2 ϕ + cos2 ϕ = 1 cos ϕ = sin(ϕ + π/2) sin(−ϕ) = − sin ϕ cos(−ϕ) = cos ϕ sin(α ± β) = sin α cos β ± cos α sin β cos(α ± β) = cos α cos β ∓ sin α sin β sin(2ϕ) = 2 sin ϕ cos ϕ cos(2ϕ) = cos2 ϕ − sin2 ϕ = 2 cos2 ϕ − 1 = 1 − 2 sin2 ϕ sin(3ϕ) = 3 sin ϕ − 4 sin3 ϕ cos(3ϕ) = 4 cos3 ϕ − 3 cos ϕ ϕ sin2 ( ϕ2 ) = 1−cos 2 ϕ cos2 ( ϕ2 ) = 1+cos 2 tan α±tan β tan(α ± β) = 1∓tan α tan β 2 tan ϕ tan(2ϕ) = 1−tan 2ϕ 3

ϕ−tan ϕ tan(3ϕ) = 3 tan 1−3 tan2 ϕ 1−cos ϕ 2 ϕ tan ( 2 ) = 1+cos ϕ ϕ sin ϕ = 1+cos tan ϕ2 = 1−cos sin ϕ ϕ

17