0 Titelblad en Inhoud en Foutenanalyse

Syllabus Numerieke Analyse I en II

P. de Groen

Abstract

Deze syllabus omvat hoofdstukken over een aantal onderwerpen die in de cursus behandeld wor- den: Afrondfouten, “Fast Fourier Transformatie”, interpolatie, approximatie door polynomen, spline-approximatie en numerieke integratie, Totale Kleinste Kwadraten, Geconjungeerde gradi- enten en Newton-methoden. Voor het hoofdstuk “Numerieke Lineaire Algebra” wordt verwezen naar de cursus van het tweede jaar Bachelor en het standaardwerk:

G.H. Golub & C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, Maryland, USA, 2^de druk, 1988.

Andere goede referentiewerken voor Numerieke Analyse zijn:

R. Bulirsch & J. Stoer, Introduction to Numerical Analysis, Springer Verlag, Berlin, 1977. (Ook verkrijgbaar in een goedkope duitstalige pocketeditie).

D. Kincaid & W. Cheney, Numerical Analysis, Brooks & Cole Publishing Company, Pacific Grove, California, USA, 1991; 2de druk, 1996.

0

CONTENTS 1

Contents

0 Foutenanalyse 3

0.a Elementaire definities . . . . 3

0.b Voorstelling van re¨ele getallen en “floating-point” aritmetiek . . . . 3

0.c Voorbeelden van een afrondfoutenanalyse . . . . 5

1 Interpolatie en Approximatie 8 1.a Lagrange interpolatie . . . . 8

1.b Alternatieven voor het representeren en uitrekenen van een interpolatiepolynoom . . . . 11

1.c Polynoomapproximatie . . . . 14

1.d Approximaties op deelintervallen . . . . 15

1.e Kubische Splines . . . . 16

1.f Praktisch rekenen met kubische splines . . . . 17

2 Fourier analyse en de “Fast Fourier Transform” 21 3 Numerieke Integratie 25 3.a Probleemstelling . . . . 25

3.b Gauss-integratie . . . . 27

3.c Samengestelde integratieformules . . . . 29

3.d Romberg integratie . . . . 32

3.e Voorbeeld van het Rombergschema . . . . 33

3.f Integratie met veranderlijke stapgrootte . . . . 35

3.g Numerieke stabiliteit . . . . 36

3.h De sommatieformule van Euler-McLaurin en de trapeziumregel . . . . 37

4 B-splines 40 4.a Definities en elementaire eigenschappen . . . . 40

4.b Het rekenen met B-splines . . . . 43

4.c B´ezier-polynomen en controlepunten . . . . 45

4.d B-spline-krommen en controlepunten . . . . 47

5 Totale kleinste kwadraten 49 5.a Beste benadering in IR . . . . 49

5.b Lineaire regressie in IR² . . . . 49

5.c Lineaire regressie van x op y . . . . 51

5.d De totale kleinste-kwadratenbenadering . . . . 52

5.e Regressie in meer dan twee dimensies . . . . 53

5.f De normaalvergelijkingen in IR^m . . . . 54

5.g Oplossing via de singuliere-waardenontbinding . . . . 55

5.h Totale kleinste kwadraten in IRⁿ . . . . 56

5.i Een alternatieve benadering voor totale kleinste kwadraten in IRⁿ . . . . 58

6 Geconjungeerde gradienten 59 7 Stelsels niet-lineaire vergelijkingen en minimalisatieproblemen 72 7.a Probleemstelling in ´e´en dimensie . . . . 72

7.b Intervalhalvering of binaire search . . . . 72

7.c Successieve substitutie . . . . 72

7.d Newton-Raphson . . . . 73

7.e Problemen in verscheidene dimensies . . . . 75

7.f Een aangepaste (gedempte) Newtonmethode . . . . 76

7.g De methode van de steilste helling (steepest descent) . . . . 79

8 Oefeningen Numerieke Analyse 82 8.a Fouriertransformatie . . . . 82

8.b Approximatie . . . . 84

8.c Orthogonale polynomen . . . . 85

8.d Een spline-benadering . . . . 88

8.e Numerieke integratie . . . . 90

8.f Lineaire Algebra . . . . 95

CONTENTS 2

8.g Matrixalgoritmen . . . . 97

8.h De Cholesky-ontbinding . . . 100

8.i Givensrotaties etc . . . 102

8.j Niet-lineaire problemen . . . 105

9 Enige recente schriftelijke examens voor Numerieke Analyse I 106 9.a Examen Numerieke Analyse, juni 2001 . . . 106

9.b Examen Numerieke Analyse, april 2000 . . . 108

9.c Examen Numerieke Analyse, september 1999 . . . 110

0 FOUTENANALYSE 3

0 Foutenanalyse

0.a Elementaire definities

Gegeven is een grootheid X en haar benaderingX. De absolute en relatieve fouten in de benadering^e X worden gegeven door:e

absolute fout in X : F^e _X := X − X^e zodat X = X + F^e _X, relatieve fout in X : f^e _X := X − X^e

X zodat X = X(1 + f^e _X) (mits X 6= 0).

(0.1)

Het begrip “absolute fout” heeft in principe niets te maken met absolute waarden; absoluut staat slechts in tegenstelling tot relatief. De absolute fout heeft dezelfde dimensies (bv. lengte, gewicht, tijd) als de grootheid zelf, terwijl de relatieve fout dimensieloos is.

Opgave 1: laat zien, dat voor de absolute en relatieve fouten in de twee grootheden X en^e Y geldt:^e FX+Y = FX + FY en fX∗Y = fX + fY + fXfY .

Als we de (absolute of relatieve) fout in een (gemeten of berekende) grootheid X kennen, dan^e kennen we ook de grootheid zelf! Helaas zijn we bijna nooit in deze situatie en kennen we alleen een bovengrens voor de absolute waarde van de fout. In het gangbare spraakgebruik spreken we gewoonlijk over de “fout” in een grootheid terwijl we zo’n bovengrens bedoelen (of nog erger, terwijl we de spreiding in de stochastische fluctuaties rond de exacte waarde bedoelen). Dus, bij een gegeven benaderingX van een grootheid X defini¨eren we:^e

∆X is (een bovengrens voor) de absolute fout in X^e als |X − X | ≤ ∆^{e X}, δX is (een bovengrens voor) de relatieve fout inX^e als

¯¯

¯¯

¯

X − Xe X

¯¯

¯¯

¯ ≤ δ^X. (0.2) Opgave 2: Bewijs de volgende rekenregels voor “de fouten” in de grootheden X en^e Y :^e

∆X±Y ≤ ∆^X+ ∆Y ∆XY ≤ |Y |∆^X+ |X|∆^Y + ∆X∆Y , δX±Y ≤ |X|δX+ |Y |δY

|X ± Y | δXY ≤ δ^X+ δY + δXδY . (0.3) N.B. Lees deze regels alsvolgt: Als ∆X en ∆Y bovengrenzen zijn voor de fouten in X resp. Y , dan is er een bovengrens ∆_X±Y voor de fout in X ± Y waarvoor geldt ∆X±Y ≤ ∆X + ∆_Y . Hieruit volgt dus dat ∆X + ∆Y een bovengrens voor de fout in X ± Y is, etc.

Wat zijn de overeenkomstige rekenregels voor de absolute en relatieve fouten (bovengrenzen) in het quoti¨ent X/Y ?

0.b Voorstelling van re¨ele getallen en “floating-point” aritmetiek

Om een groot dynamisch bereik mogelijk te maken voor re¨ele getallen worden deze in een computer opgeslagen in de vorm mantisse maal exponent. Hiertoe wordt een grondtal β (meestal 2, soms 8 (vroeger op CDC) of 16 (IBM)) gekozen. Een x ∈ IR kan dan worden voorgesteld door een paar (m, e) met

x = m · β^e, (0.4)

waarin m de mantisse is en e de exponent. Omdat het paar (m · β , e − 1 ) hetzelfde getal voorstel kunnen we de mantisse normaliseren, b.v. door 1/β ≤ | m | < 1. Het spreekt vanzelf dat we in de

0 FOUTENANALYSE 4

praktijk een eindige representatie willen hebben en dus het aantal β-tallige cijfers in mantisse en exponent zullen beperken. De IEEE-standaard voor 64-bits REALs is een tweetallige representatie (β = 2) met 53 resp. 10 bits voor de absolute waarden van mantisse en exponent en twee tekenbits.

Omdat een genormaliseerde binaire mantisse altijd begint met een 1 (ga na!), hoeft dit eerste bit niet opgeslagen te worden. Met 10 bits is ook de grootte van de exponent aan een maximum gebonden.

Getallen die een exponent groter dan 2¹⁰ of kleiner dan 2⁻¹⁰ vragen (waarvan de absolute waarde dus kleiner dan (ongeveer) 10⁻³⁰⁰of groter dan 10³⁰⁰ is), kunnen dus niet gerepresenteerd worden;

we spreken dan van over- of underflow. De IEEE-standaard geeft de mogelijkheid om by underflow een getal op nul te zetten, en bij overflow een N aN (Not a Number) te genereren zodat er een soepele foutenopvang mogelijk is. Een re¨eel getal binnen het bereik zal in het algemeen niet exakt gerepresenteerd kunnen worden. Voor een gegeven x ∈ IR (binnen het bereik) noteren we met fl(x) het meest naburige wel repesenteerbare getal (machinegetal). Het verschil x − fl(x) is dan de afrondfout.

Stelling. Als voor een machinegetal een β-tallige representatie wordt gekozen met t bits in de mantisse, dan geldt voor de relatieve afrondfout bij afronding naar het dichtstbijzijnde machinegetal (behoudens over- en underflow):

¯¯

¯¯x − fl(x) x

¯¯

¯¯≤ η maar ook

¯¯

¯¯x − fl(x) fl(x)

¯¯

¯¯≤ η met η := ¹2β^1−t. (0.5) De grootheid η heet de machineprecisie.

Opgave 3: Bewijs deze stelling.

Ga ook na, dat er (behoudens over- en underflow) getallen ε1 en ε2 zijn bij iedere aritmetische operatie ⊙ ∈ {+, −, ∗, /} tussen twee machinegetallen x en y, zodat

fl(x ⊙ y) = (x ⊙ y)(1 + ε¹) = x ⊙ y 1 + ε2

met |ε¹| ≤ η en ε²| ≤ η . (0.6)

Opmerking. We kunnen η ook defini¨eren als het grootste re¨ele getal, zodat fl(1 + η) = 1, ga na!

Opgave 4: De reeksontwikkeling van de exponentiaal is: e^x = X∞

k=0

x^k k!

Hoeveel termen heb je nodig om e⁻⁵ te berekenen met een relatieve fout kleiner dan 10⁻³?

Kun je dit doen met een computer, waarin de variabelen van het type IR een mantisse van 4 decimalen hebben ?

Is er een betere manier om e⁻⁵ te berekenen met zo’n computer ?

Opgave 5: Laat f een voldoend gladde re¨ele funktie (b.v. f (x) = sin(x)) zijn met maxx |f^′′′(x)| ≤ M .

De afgeleide van f in x kunnen we dan benaderen met de centrale differentie D_hf (x) := f (x + h) − f(x − h))

2h .

Laat zien, dat voor de afbreekfout in D_hf geldt:

f (x + h) − f(x − h)

2h = f^′(x) + h²

6 f^′′′(x + ϑh) met |ϑ| ≤ 1 . (0.7) Veronderstel, dat er voor de berekening van f een procedure beschikbaar is, die bij iedere waarde van x een resultaat aflevert met een relatieve fout kleiner dan of gelijk aan 2η . Geef dan een (goede) bovengrens voor de relatieve fout in de berekende waarde van D_hf als funktie van h en schets een grafiek van (een bovengrens voor) de totale fout (afbreek- plus afrondfout) in deze berekende waarde.

0 FOUTENANALYSE 5

0.c Voorbeelden van een afrondfoutenanalyse Gevraagd te berekenen x = ϕ(a).

Met een algoritme voor het berekenen van ϕ(a) vinden we ten gevolge van afrondfouten de berekende waarde: fl(x).

In een foutenanalyse proberen we fouten δx, δa of εa en εx te vinden, zodat fl(x) = x + δ_x voorwaartse foutenanalyse

= ϕ(a + δ_a) achterwaartse foutenanalyse

= ϕ(a + εa) + εx gemengde foutenanalyse

Definitie: We noemen de algoritme numeriek stabiel als we kunnen bewijzen:

δx of εx van de grootteorde van de onvermijdelijke fout, δa of εa van de grootteorde van de machineprecisie.

Voorbeeld 1: Er is een ε met | ε | ≤ η ( = machineprecisie ) zodat

fl(a + b) =





a + b + ε (a + b) voorwaarts

ea + ^eb met a := a(1 + ε)_e en ^eb := a(1 + ε) achterwaarts Voorbeeld 2: Er zijn ε1 en ε2 ( met | εⁱ| ≤ η ) zodat

fl(1 − x²) = (1 − x ∗ x ∗ (1 + ε¹)) ∗ (1 + ε²)

= (1 − xe²) (1 + ε2) met x := x_e √

1 + ε1 gemengd.

Voorbeeld 3: Geef een schatting van de afrondfout in de berekende waarde van de positieve wortel van de vierkantsvergelijking

a − 2x − c x² met a ≥ 0 en c ≥ 0 bij gebruik van de formule

x := −1 + √

1 + a c c

onder de aanname betreffende de afrondfout in de berekende waarde van de vierkantswortel fl(√

x) = √

x (1 + ε_x) met | εx| ≤ η voor iedere x . Antwoord: Er bestaan ε1, ε2 en ε3 met | εi| ≤ η , zodat

fl(√

1 + a c) = ^p(1 + a c (1 + ε1)) (1 + ε2) (1 + ε3)

= √

1 + _eac (1 + ξ1) met ξ1 := √

1 + ε2(1 + ε3) − 1 en ea := a (1 + ε1) Bijgevolg zijn er ξ2 en ξ3, ( | ξⁱ| ≤ η ) zodat:

fl(x) = −1 + √

1 + _ea c (1 + ξ1)

c (1 + ξ2) (1 + ξ3)

= −1 + √

1 + _ea c

c (1 + ξ2) (1 + ξ3) + ξ1

√1 + _ea c c De tweede term is groot t.o.v. x als | ac | ≪ 1.

Alternatieve (numeriek stabiele) rekenwijze voor deze wortel:

x := a

1 + √

1 + a c.

0 FOUTENANALYSE 6

Voorbeeld 4: Afrondfout in de berekende waarde van het inprodukt S :=

Xn

i=1

x_iy_i berekend met algoritme: S := 0;

for i := 1 to n do S := S + x_i ∗ yi

Voor de berekende waarde van S vinden we getallen ξi en ηi met | ξⁱ|, | εⁱ| ≤ η, i = 1 · · · n : fl(S) = x1y1(1 + ξ1) (1 + ε2) · · · (1 + εⁿ)

+ x2y2(1 + ξ2) (1 + ε2) · · · (1 + εn) + · · ·

+ xn−2yn−2(1 + ξn−2) (1 + εn−2) · · · (1 + εⁿ) + x_n−1y_n−1(1 + ξ_n−1) (1 + ε_n−1) (1 + ε_n) + xnyn(1 + ξn) (1 + εn)

zodat

S − fl(S) = Xn

i=1

x_iy_iζ_i met

ζi := 1 − (1 + ξⁱ) (1 + εi) · · · (1 + εⁿ) en | ζⁱ| ≤ (n − i + 2) η als nη ≤ 0.1 . Bijgevolg geldt voor de voorwaartse fout:

|S − fl(S)

S | ≤ (n + 1) η

| S | Xn

i=1

| xiy_i| ≤ (n + 1) ηk x k²k y k²

| x^Ty| Voorbeeld 5: Bereken xn uit de vergelijking

a = Xn

i=1

xiyi, a, x1 · · · xⁿ⁻¹, y1 · · · yⁿ gegeven, en bepaal de afrondfout in de berekende waarde van xn.

Algoritme:





S := a;

for i := 1 to n − 1 DO S := S − xi ∗ yi; xn := S / yn

Voor de berekende waarden van S en xn vinden we voor zekere ξi en ηi met | ξⁱ|, | εⁱ| ≤ η : fl(S) = a (1 + ε1) · · · (1 + εn−1)

− x¹y1(1 + ξ1) (1 + ε1) · · · (1 + εn−1)

− x²y2(1 + ξ2) (1 + ε2) · · · (1 + εⁿ⁻¹)

− · · ·

− xn−2y_n−2(1 + ξ_n−2) (1 + ε_n−2) (1 + ε_n−1)

− xⁿ⁻¹yn−1(1 + ξn−1) (1 + εn−1) en

xe_n := fl(x_n) = fl(S) / ( y_n(1 + ξ_n) )

0 FOUTENANALYSE 7

Deling door (1 + ε1) · · · (1 + εn−1) geeft de achterwaartse foutschatting:

a = x1y1(1 + ξ1) + x2y2

1 + ξ2

1 + ε1 + · · · + x_n−1y_n−1 1 + ξ_n−1

(1 + ε1) · · · (1 + εⁿ⁻²) + x_enyn

1 + ξn

(1 + ε1) · · · (1 + εn−1)

=

n−1X

i=1

x_iy_i(1 + δ_i)x_eny_n(1 + δ_n)

met

δi := 1 + ξ_i

(1 + ε1) · · · (1 + εⁱ⁻¹) − 1 , zodat | δⁱ| ≤ (i + 1) η als n η < 0.1 . Conclusie: De berekende waardex_en is de oplossing van een naburige vergelijking

a = Xn

j=1

x_jy_ej, y_ej := y_j(1 + δ_j) .

Opgave 6: Voor de standaardafwijking S bestaan in de statistiek twee formules die wiskundig (in exacte re¨ele arithmetiek) gelijkwaardig zijn :

S² = 1 n − 1 (

Xn

i=1

x²_i − n g² ) en S² = 1 n − 1

Xn

i=1

(x_i− g)² met g het gemiddelde:

g := 1 n

Xn

i=1

xi .

Welke van de twee zou je gebruiken in een numeriek programma en waarom?

REFERENCES 112

References

[1] M. Hestenes & E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Research NBS, 49, pp. 409 – 436, 1952.

[2] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, J. Research NBS, 45, pp. 255 – 282, 1950.

[3] J.K. Reid, On the method of conjugate gradients for the solution of large sparse systems of linear equa- tions, Proc. Conf. on Large Sparse Sets of Linear Equations, Academic Press, New York, 1971.

[4] J.A. Meijerink and H.A. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Math.of Comp., 31, pp. 148 – 162, 1977.

[5] G.H. Golub & C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, Maryland, USA, 1^ste druk, 1983, 2^dedruk, 1988, 3^de druk, 1995.

[6] R. Bulirsch & J. Stoer, Introduction to Numerical Analysis, Springer Verlag, Berlin, 1977. (Ook verkri- jgbaar in een goedkope duitstalige pocketeditie).

[7] D. Kincaid & W. Cheney, Numerical Analysis, Brooks & Cole Publishing Company, Pacific Grove, California, USA, 1991; 2de druk, 1996.