X Het geheel

(1)

P. de Groen

Abstract

Deze syllabus omvat hoofdstukken over een aantal onderwerpen die in de cursus behandeld worden: 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

(2)

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ê _X := X − Xê zodat X = X + Fê _X, relatieve fout in X : fê _X := X − Xê

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 groothedenX en^e Y geldt:^e FX+Y = FX + FY en f_X∗Y = fX + fY + fXfY .

Als we de (absolute of relatieve) fout in een (gemeten of berekende) grootheid X kennen, danê 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ëren we:ê

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

¯¯

¯

X − Xe X

¯¯

¯ ≤ δ^X. (0.2) Opgave 2: Bewijs de volgende rekenregels voor “de fouten” in de groothedenX 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

(5)

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 η := ¹₂β^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 + ε1) = x ⊙ y 1 + ε2

met |ε1| ≤ η en ε2| ≤ η . (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.

(6)

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 ε₁ en ε₂ ( met | εi| ≤ η ) zodat

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

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

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 ε₁, ε₂ en ε₃ met | εi| ≤ η , zodat

fl(√

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

= √

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

1 + ε₂(1 + ε₃) − 1 en ea := a (1 + ε₁) Bijgevolg zijn er ξ₂ en ξ₃, ( | ξi| ≤ η ) zodat:

fl(x) = −1 + √

1 + _ea c (1 + ξ₁)

c (1 + ξ₂) (1 + ξ₃)

= −1 + √

1 + _ea c

c (1 + ξ₂) (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.

(7)

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) = x₁y₁(1 + ξ₁) (1 + ε₂) · · · (1 + εn)

+ x₂y₂(1 + ξ₂) (1 + ε₂) · · · (1 + εn) + · · ·

+ x_n−2y_n−2(1 + ξ_n−2) (1 + ε_n−2) · · · (1 + εn) + 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 + ξi) (1 + ε_i) · · · (1 + εn) en | ζi| ≤ (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 k2k y k2

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

a = Xn i=1

xiyi, a, x1 · · · xn−1, y1 · · · yⁿ gegeven, en bepaal de afrondfout in de berekende waarde van x_n.

Algoritme:





S := a;

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

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

− x1y₁(1 + ξ₁) (1 + ε₁) · · · (1 + εn−1)

− x2y₂(1 + ξ₂) (1 + ε₂) · · · (1 + εn−1)

− · · ·

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

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

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

(8)

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

a = x₁y₁(1 + ξ₁) + x₂y₂ 1 + ξ₂

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

(1 + ε1) · · · (1 + εn−2) + x_e_ny_n 1 + ξ_n

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

=

n−1X

i=1

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

met

δi := 1 + ξ_i

(1 + ε₁) · · · (1 + εi−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_e_j, y_e_j := 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?

(9)

1 Interpolatie en Approximatie

In dit hoofdstuk bespreken we methoden om een gegeven functie van een veranderlijke te benaderen met een (gemakkelijk berekenbare) functie uit een voorgeschreven klasse, zoals polynomen (veeltermen) of trigonometrische polynomen (sommen van machten van sinus en cosinus). Bij approximatie (benadering) wordt alleen de eis gesteld, dat het verschil met de gegeven functie klein is in een of andere zin. Bij interpolatie moet de benaderende functie gelijk zijn aan de gegeven functie in een zeker aantal (vooraf voorgeschreven) punten.

1.a Lagrange interpolatie

Het Lagrange-interpolatiepolynoom π_f van een continue functie f op n steunpunten { x¹, · · · , xⁿ} is het polynoom van graad n − 1, dat gelijk is aan f in deze punten:

π_f(x_i) = f (x_i) , i = 1 · · · n. (1.1) Onder de aanname xi 6= x^j als i 6= j is πf uniek en er geldt:

π_f(x) = Xn i=1

f (x_i) L⁽ⁿ⁾_i (x), L⁽ⁿ⁾_i (x) :=

Yn j=1, j6=i

x − xj

xi − x^j . (1.2)

Opgave 1.1: Bewijs dit.

Voor het Lagrange-interpolatiepolynoom geldt de volgende foutschatting:

Stelling. Zij I ∈ IR een interval dat de steunpunten x1, · · · , xn bevat en zij f ∈ Cⁿ(I) (f is n keer continu differentieerbaar), dan geldt:

f (x) = π_f(x) + Yn i=1

(x − xi) f⁽ⁿ⁾(ξ_x)

n! , (1.3)

waar ξ_x ∈ int(x, x1, · · · , xn) , het interval opgespannen door de punten x , x₁, · · · , xn. Bewijs. Definieer de functie g door de vergelijking

f (x) = π_f(x) + g(x) Yn i=1

(x − xi)

en beschouw de functie

ϕ_t(x) := f (x) − πf(x) − g(t) Yn i=1

(x − xi) .

Deze heeft (minstens) n+1 onderling verschillende nulpunten binnen het interval int(x, x₁, · · · , xn), als t 6= xⁱ, ∀ i. Volgens de stelling van Rolle heeft de afgeleide dus nog minstens n verschillende nulpunten. Als we doorgaan met differenti¨eren en de stelling van Rolle toepassen, vinden we dat de n-de afgeleide nog minstens een nulpunt ξ_t binnen int(x, x₁, · · · , xn) heeft. Zo vinden we

0 = dⁿ

dxⁿ ϕ_t |x = ξt = f⁽ⁿ⁾(ξ_t) − n! g(t) , waaruit (1.3) volgt.

N.B. De rest in (1.3) heet de restterm van Lagrange. In deze restterm worden geen eisen gesteld aan de positie van het punt x t.o.v. de steunpunten. Als x tussen de steunpunten ligt, hebben we interpolatie in de eigenlijke betekenis van het woord, en we spreken van extrapolatie als x buiten

(10)

dit interval ligt. Aan het functie-onafhankelijke deel van de restterm, Πⁿ_i=1 (x − xi), kunnen we zien dat voor extrapolatie een veel grotere afbreekfout kunnen verwachten dan voor interpolatie.

Opgave 1.2: Een algemenere vorm van interpolatie is die, waarbij in de steunpunt xi niet alleen gelijkheid van de functiewaarde wordt ge¨eist, maar ook gelijkheid van eerste n_i−1 afgeleiden, d.w.z.

π^(j)_f (xi) = f^(j)(xi) , voor j = 0 , · · · , nⁱ− 1 , en i = 1 , · · · , n . (1.4) De graad van het interpolatiepolynoom is dan N − 1 met N := Σⁿi=1 n_i . Laat zien, dat analoog aan (1.3) geldt:

f (x) = π_f(x) + f^{(N )}(ξ_x) N !

Yn i=1

(x − xi)ⁿⁱ (1.5)

voor zekere ξ_x ∈ int(x, x1, · · · , xn) .

N.B. Als n_i = 2 ∀ i (in alle steunpunten stemmen dus de waarde van het interpolatiepolynoom en van zijn afgeleide overeen met die van de gegeven functie), dan spreken we van Hermite − interpolatie.

Formule (1.2) is niet erg geschikt om de waarde van het interpolatiepolynoom in een punt te berekenen, niet alleen vanwege de hoeveelheid rekenwerk, maar ook vanwege een mogelijke accumulatie van afrondfouten bij het optellen van de bijdragen van de verschillende steunpunten.

Een betere manier is de Newtonse algoritme met “gedeelde differenties”. We defini¨eren recursief de m-de gedeelde differentie van f op de m + 1 punten {xi, x_i+1, · · · , xi+m} door

f (x_i, x_i+1, · · · , xi+m) := f (x_i, x_i+1, · · · , xi+m−1) − f(xi+1, x_i+2, · · · , xi+m)

xi − x^i+m , (1.6)

voor m = 1, 2, · · · .

Opgave 1.3: Bewijs met volledige inductie f (x_i, x_i+1, · · · , xi+m) =

m+iX

j=i

f (x_j) Qm+i

l=i, l6=j (x_j − xl) (1.7)

en dat hieruit volgt, dat f (x_i, x_i+1, · · · , xi+m) een symmetrische functie is van haar argumenten, d.w.z. dat de functiewaarde invariant is onder verwisseling van argumenten.

Met behulp van deze gedeelde differenties kunnen we het Lagrange-interpolatiepolynoom op de punten {x1, · · · , xn} schrijven als

π_f(x) = f (x₁) + (x − x1)f (x₁, x₂) + (x − x1)(x − x2)f (x₁, x₂, x₃)+

+ · · · + (x − x1)(x − x2) · · · (x − xn−1)f (x₁, x₂, · · · , xn) (1.8) en voor de restterm geldt dan

f (x) − πf(x) = f (x1, · · · , xⁿ, x) Yn j=1

(x − x^j) . (1.9)

Opgave 1.4: Bewijs de formules (1.8) en (1.9) met volledige inductie. Laat ook zien, dat de k-de gedeelde differentie evenredig is met de k-de afgeleide:

f (x₀, · · · , xk) = f^k(ξ)

k! , ξ ∈ int(x0, · · · , xk). (1.10)

(11)

Uit formule (1.6) zien we dat de gedeelde differenties bij een gegeven verzameling steunpunten {x1, x₂, · · ·} een driehoekig tableau vormen, waarvan we de elementen alsvolgt kunnen berekenen:

fori := 1 to n do f_i⁽⁰⁾ := f (x_i) ; forj := 1 to n − 1 do

fori := n downto j + 1 do f_i^(j) := f_i^(j−1) − f_i−1^(j−1) xi− xi−j

;

(1.11)

Het tableau heeft dan de volgende vorm:

x 0 1 2 3 4 5

x1 f₁⁽⁰⁾ f₂⁽¹⁾ x₂ f₂⁽⁰⁾ f₃⁽²⁾

f₃⁽¹⁾ f₄⁽³⁾ x₃ f₃⁽⁰⁾ f₄⁽²⁾ f₅⁽⁴⁾

f₄⁽¹⁾ f₅⁽³⁾ · · · x₄ f₄⁽⁰⁾ f₅⁽²⁾ f₆⁽⁴⁾

f₅⁽¹⁾ f₆⁽³⁾ · · · x₅ f₅⁽⁰⁾ · · · ·

· · · ·

Vervolgens kunnen we de waarde van π_f in het gewenste punt x evalueren met een Horner-achtig schema:

π_f(x) := (· · · ( (x − xn−1) f_n⁽ⁿ⁻¹⁾ + f_n−1⁽ⁿ⁻²⁾) (x − xn−2) + · · ·) (x − x1) + f₁⁽⁰⁾ (1.12a) of in pseudoPASCAL:

som := fn⁽ⁿ⁻¹⁾ ;

forj := n − 1 downto 1 do som := som ∗ (x − xj) + f_j^(j−1) ; (1.12b) Op deze manier behoeven we de gedeelde differenties voor de benadering maar een keer te berekenen en hebben we vervolgens voor iedere evaluatie van π_f slechts n flops nodig.

Opgave1.5: Ga na, dat we bij de berekening van het differentietableau in de j-de slag f_i^(j)kunnen opbergen op dezelfde geheugenplaats als f_i^(j−1), zonder daarmee waarden te verliezen, die we later (bij de evaluatie van π_f) nog nodig hebben.

Bij exakt rekenen is de waarde van π_f(x) onafhankelijk van de volgorde, waarin we de steunpunten {xⁱ} ordenen, maar bij rekenen op een rekenmachine (met afrondende floating point operaties) is dit wel het geval. De totale afrondfout in het berekende resultaat van de som (1.8) zal i.h.a.

kleiner zijn, naarmate de termen van de reeks sneller dalen (mits we ze van klein naar groot sommeren, zoals in algoritme (1.12b). Aan de posities van de steunpunten en aan de co¨effici¨enten, de gedeelde differenties, kunnen we natuurlijk niets veranderen, maar we kunnen wel de volgorde van de steunpunten zo kiezen, dat x₁ het dichtst bij x ligt, en meer algemeen, dat geldt

| x − x1| ≤ | x − x2| ≤ | x − x3| ≤ · · ·

(12)

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

Stel dat we een functie f interpoleren op de equidistante punten 0, ±1, ±2, · · ·, en dat we de waarde van die functie exact kennen in alle punten, behalve in 0, waar we een fout van ´e´en eenheid hebben, dan kunnen we de doorwerking van deze fout zien in de volgende differentietabel:

x f out 1^stedif 2^dedif 3^dedif 4^dedif 5^dedif 6^dedif 7^dedif 8^stedif

−4 1

1

−3 1 − 8

1 − 7

−2 1 − 6 28

1 − 5 21

−1 1 −4 15 −56

1 −3 10 −35

0 1 −2 6 −20 70

−1 3 −10 35

1 1 −4 15 −56

−1 5 −21

2 1 − 6 28

− 1 7

3 1 − 8

− 1

4 1

Om de gedeelde differenties te verkrijgen, moeten we de k-de differenties nog delen door k!. We zien in deze tabel, dat de fouten in de hogere differenties ten gevolge van deze ene fout in f (0) snel aangroeien. Om de resulterende fout te vinden in de k-de term van formule (1.8), moeten we de fout in de k-de differentie delen door k! en vermenigvuldigen met (x − x¹) (x − x²) · · · (x − xk−1).

Als x₁ < x < x₂ < · · · < xn, dan zijn beide factoren ongeveer gelijk en geeft de k-de kolom de werkelijke fout in de k-de term weer; als · · · x4 < x₂ < x < x₁< x₃· · · (interpolatie in het midden), dan is de werkelijke fout ongeveer 2^−k kleiner. Maar ook dan blijft deze van orde 1 en kan de totale fout in π_f toenemen met de orde van de interpolatie. Het heeft dan ook weinig zin om via deze methode een interpolatiepolynoom te bepalen van een orde hoger dan ongeveer 10 bij de gebruikelijke machineprecisie van 10⁻¹⁵(reals van 64 bit met een mantisselengte van 53 bits).

Een alternatieve methode voor het bepalen van het interpolatiepolynoom zou het volgende kunnen zijn. We zoeken de coëfficiënten van π_f(x) = Σⁿ⁻¹_j=0 ajx^jbij gegeven functiewaarden π_f(xi) = f (x_i), i = 1, · · · , n. De coëfficiënten {aj | j = 0, · · · , n − 1} kunnen we dus eenvoudig bepalen door het oplossen van n lineaire vergelijkingen

n−1X

j=0

aj x^j_i = f (xi) , i = 1, · · · , n , (1.13a)

oftewel

A





 a₀ a1

..

a_n−1





 :=







1 x₁ · · · xⁿ⁻¹1

1 x2 · · · xⁿ⁻¹2

... ... · · · ... ... ... · · · ... 1 xn · · · xⁿ⁻¹n











 a₀ a1

..

a_n−1







=





 f (x₁) f (x2)

..

f (xn)







. (1.13b)

(13)

De matrix A in dit stelsel vergelijkingen (A is een V andermondematrix ) is in het algemeen helaas zeer slecht geconditioneerd.

Voorbeelden:

1. Met de steunpunten x_i = i/10, i = 0 , · · · , 10 vinden we een 11 × 11 matrix A met de conditiegetallen κ₂(A) = 1.16₁₀8 en κ_∞(A) = 3.4₁₀8.

2. Met de steunpunten x_i = cos(²ⁱ⁺¹₂₂ π), i = 0 , · · · , 10 (de nulpunten van het 10^de Chebyshev polynoom) T₁₀ vinden we een 11 × 11 matrix A met de conditiegetallen κ2(A) = 3.6₁₀3 en κ_∞(A) = 9.3103.

Het oplossen van het stelsel (1.13) is dus geen goede zaak. De reden is, dat de gekozen representatie

π_f(x) =

n−1X

j=0

a_j x^j (1.14)

nogal slecht is om mee te rekenen als n niet heel klein is.

Een betere representatie verkrijgen we door π_f te schrijven als een som van orthogonale polynomen in plaats van machten van x. Zij bij gegeven steunpunten {xⁱ | i = 1, · · · , n} de verzameling {Vj(x) | i = 0 , · · · , n − 1} een stelsel veeltermen met graad(Vj) = j en zo, dat

Xn i=1

V_j(x_i) V_l(x_i) =

½ 0 if j 6= l ,

1 if j = l , (j , l = 0 , · · · , n − 1), (1.15) d.w.z. de veeltermen zijn diskreet orthogonaal op de gegeven steunpunten. Als we nu het gezochte interpolatiepolynoom representeren als een lineaire combinatie van deze orthogonale veeltermen,

π_f(x) =

n−1X

j=0

b_jV_j(x) , (1.16a)

dan vinden we voor de co¨effici¨enten het stelsel vergelijkingen

B





 b₀ b1

..

b_n−1





 :=







V₀(x₁) V₁(x₁) · · · Vn−1(x₁) V0(x2) V1(x2) · · · Vn−1(x2)

· · · ·

· · · · V0(xn) V1(xn) · · · Vn−1(xn)











 b₀ b1

..

b_n−1







=





 f (x₁) f (x2)

..

f (xn)







. (1.16b)

De matrix B in dit stelsel heeft orthonormale kolommen en is dus orthogonaal (B^TB = BB^T = I) en is dus zeer goed geconditioneerd en zeer eenvoudig oplosbaar,

b_j = Xn i=1

V_j(x_i) f (x_i), i = 0, · · · , n − 1 .

Het grote probleem is echter om zo’n stelsel te vinden bij gegeven interpolatiepunten {x1, · · · , xn}.

Een mogelijkheid is om de veeltermen {1 , x , x², · · · } successievelijk te orthogonaliseren met behulp van de methode van Gram-Schmidt. Dit is echter equivalent met het oplossen van het stelsel (1.13b) met MGS (“modified Gram-Schmidt”, zie Golub & Van Loan) en dus met behoud van alle ellende wegens de slechte conditie van dit stelsel. Voor een willekeurig stel interpolatiepunten blijkt er niets beters te zijn! Als de interpolatiepunten gegeven zijn als de nulpunten van een “mooi” polynoom, zoals het n-de Chebyshev polynoom T_n, dan zijn er andere middelen om zo’n diskreet-orthogonale rij veeltermen te vinden.

(14)

Lemma:

n−1X

k=0

cos(2k + 1

2n mπ) =







0 , als m mod 2n 6= 0 ,

n , als m mod 2n = 0 en m/2n even

−n , als m mod 2n = 0 en m/2n oneven .

(1.17)

Bewijs.

n−1X

k=0

cos(2k + 1

2n mπ) = 1 2

n−1X

k=0

( e^imπ^2k+1²ⁿ + e^−imπ^2k+1²ⁿ )

= 1

2

2n−1X

k=0

e^imπ²ⁿ e^imπkⁿ ;

als m mod 2n 6= 0, dan kunnen we de meetkundige reeks in de laatste som sommeren en is het resultaat nul, en anders zijn alle termen +1 of −1 .

Stelling. Op de nulpunten {τ1, · · · , τn} van het n-de Chebyshev polynoom Tn geldt:

2 n

Xn k=1

T_i(τ_k) T_j(τ_k) =

( 0, als i 6= j,

1, als i = j, (0 < i, j < n) (1.18) Bewijs. Met de gelijkheid τ_k = cos ^2k+1_2n π en met de transformatie T_j(cos t) = cos jt kunnen we formule (1.17) toepassen.

Gevolg: Als π_f het Lagrange-interpolatiepolynoom van f is, op de nulpunten van T_n(π_f heeft dus graad n − 1), dan geldt:

π_f =

n−1X

j=0

c_jT_j met π_f(τ_k) =

n−1X

j=0

c_jT_j(τ_k) = f (τ_k)

en volgens(1.8) moet dan gelden

c_j = 2 n

Xn k=1

f (τ_k) T_j(τ_k).

Gebruik makend van de recurrente betrekking voor de Chebyshev polynomen, T_k+1(x) = 2xT_k(x) − Tk−1(x) voor k ≥ 1 ,

kunnen we de coëfficiënten {cj} uitrekenen in n²+ n flops en n evaluaties van de cosinus. Voor grote n kan dit echter veel efficiënter via een FFT (fast Fourier transform), immers:

cj = 2 n

Xn k=1

f_k Tj(τ_k) = 2 n

Xn k=1

f_k cos(2k − 1 2n jπ)

= 1

n Xn k=1

fk exp(2k − 1

2n jπi) + 1 n

Xn k=1

fk exp(−2k − 1 2n jπi)

= 1

n Xn k=1

f_k exp(2k − 1

2n jπi) + 1 n

Xn k=1

f_k exp(4n − 2k + 1

2n jπi)

= 1

n Xn k=1

f_k exp(2k − 1

2n jπi) + 1 n

X2n l=n+1

f_2n−l+1 exp(2l − 1 2n jπi)

= 1

n exp(j nπi)

2n−1X

k=0

f_k+1exp(jk nπi)

als f_2n−k := f_k voor 0 ≤ k < n. Behoudens een constante faktor is dit precies een (complexe) diskrete Fouriertransformatie.

(15)

Opgave 1.6: De Chebyshev polynomen van de tweede soort U_n worden gedefinieerd door U_n(cos t) := sin(n + 1)t

sin t , n = 0, 1, · · · .

a. Laat zien: U₀(x) = 1, U₁(x) = 2x en U_n+1(x) = 2 U_n(x) − Un−1(x) (dezelfde drieterms- recursierelatie als Tn maar met een ander begin).

b. Bewijs, dat U_k en U_l orthogonaal zijn t.o.v. het diskrete inprodukt Σⁿ_i=1 wi U_k(ξi) U_l(ξi) op de nulpunten ξ₁, · · · , ξn van U_nmet gewichten w_k := sin^{2 kπ}_n+1.

c. Bepaal de relatie tussen de ontwikkeling π_f = Σⁿ_i=1d_jU_j van het Lagrange-interpolatiepolynoom van f op de nulpunten van U_n en de complexe Fouriertransformatie.

1.c Polynoomapproximatie

Als we een functie willen benaderen (d.m.v. een polynoom), moeten we natuurlijk eerst afspreken op wat voor manier we de afwijking zullen meten. Veelgebruikte normen zijn de supnorm, die de maximale afwijking over een interval [a, b] meet,

k f k∞ := max

a ≤ x ≤ b | f(x) |, (1.21)

en de gewogen kwadraatnorm (of L²-norm), die het kwadraat van de afwijking maal een gewichts- functie w over het interval integreert,

k f k2,w := { Z _b

a |f(t)|² w(t) dt }¹². (1.22) Bij praktisch rekenen kennen we evenwel de te benaderen functie slechts in een eindig aantal punten (b.v. metingen) en zullen we afwijkingen moeten meten in een diskrete (semi-)norm, zoals het maximum of een gewogen som (zie b.v. 1.15) over de gegeven punten.

Volgens de stelling van Weierstrass kunnen we iedere continue functie op een kompakt interval willekeurig goed uniform benaderen met een polynoom van voldoend hoge graad. Deze stelling is constructief, d.w.z. we kunnen voor iedere functie expliciet een convergente rij polynomen aangeven.

De rij Bernstein-polynomen bijvoorbeeld, B_n(x) :=

Xn k=0

(n k) f (k

n) x^k (1 − x)^n−k , (1.20)

convergeert uniform naar f op het gesloten interval [0,1], als f continu is op dit interval. Voor numerieke doeleinden is een dergelijke polynoombenadering echter totaal onbruikbaar. Als alternatief voor de constructie van een goede polynomiale benadering kunnen we denken aan interpolatie. In het algemeen is dit een goede methode, maar we moeten voorzichtig zijn: niet ieder schema met interpolatie van steeds hogere orde convergeert¹! Een bekend voorbeeld is het volgende: Als we de functie 1/(1 − x²) op het interval [−5, 5] interpoleren in n equidistante punten x¹, · · · , xⁿ met x₁ = −5 en xn= 5, dan divergeert de rij interpolanten voor | x | ≤ 3.64 voor n → ∞, ondanks dat de gegeven functie oneindig vaak differentieerbaar is. Voor stellingen over het bestaan van “beste”

polynoombenaderingen verwijzen we naar de cursus “aanvullingen van de wiskunde”. Behalve met Fourier-Chebyshev approximaties (en dit voornamelijk wegens het verband met Fourier reeksen) wordt er in de praktijk weinig met polynoomapproximaties van hogere orde (hoger dan 5 `a 6) gewerkt. Een van de redenen is, dat interpolatie van hoge orde nogal slecht geconditioneerd is (zie

§1b) en niet noodzakelijk convergent. Een andere reden is, dat een convergerende rij interpolanten van steeds hogere graad n meestal vrij traag convergeert met een macht van 1/n als convergen- tieorde. Een verdeling van het interval in “kleine stukjes” en benadering met een polynoom van lage graad op ieder stukje apart levert meestal een veel betere banadering op met minder rekenwerk.

1Sterker: bij ieder interpolatieschema is er een C^∞-functie waarvoor het schema divergeert!

(16)

1.d Approximaties op deelintervallen

Laat f een voldoend gladde functie zijn op het interval [a, b]. We gaan deze nu benaderen door het interval op te splitsen in n deelintervallen,

a =: t₀ < t₁< · · · < tn−1 < t_n:= b , h := max

1<i<n t_i− ti−1, (1.23) met maximale maaswijdte h. We stellen on nu de vraag, hoe goed we f kunnen benaderen met stuksgewijze polynomen van de graad kleiner dan of gelijk aan k. D.w.z. op ieder deelinterval [t_i−1, t_i] benaderen we f met een polynoom van graad ≤ k . Het bepalen van een beste benadering kan zeer moeilijk zijn, maar een benadering, verkregen door f op ieder deelinterval te interpoleren blijkt goed genoeg voor het afleiden van een goede foutschatting voor dit soort benaderingen (de afbreekfout van de beste benadering kan alleen maar kleiner zijn). Om zo’n stuksgewijs polynomiale interpolant π te construeren kiezen we op het interval [0, 1] de referentiesteunpunten ξ0, ξ1, · · · , ξk. We defini¨eren de restriktie van π tot het deelinterval [t_i−1, ti] als het Lagrange- interpolatiepolynoom op de steunpunten

{ ti−1+ ξ_j(t_i− ti−1) | j = 0 , · · · , k} .

Volgens de restterm van Lagrange (1.3) geldt voor x ∈ [ti−1, t_i] de ongelijkheid

| f(x) − π(x) | = | f^(k+1)(ξ_x) (k + 1)!

Yk j=0

(k + 1)! . (1.24) De schatting van het produkt door de k+1-ste macht is vrij grof, maar geeft wel goed de orde van de (lokale) afbreekfout aan, deze is namelijk evenredig met de k+1-ste macht van de deelintervallengte.

Dus, als we deelintervallengte halveren (bij constante k), dan neemt de benaderingsfout af met een factor 2^−k−1. We zien dat een dergelijke manier van benaderen met stuksgewijze polynomen van vaste graad veel veiliger is dan benaderen met ´e´en polynoom van steeds hogere graad en niet kritisch afhangt van afgeleiden van hoge orde van f en van de verdeling van de steunpunten over het interval. Bovendien is deze wijze van benaderen flexibel; als de k-de afgeleide op een (klein) deel van het interval groot is en elders klein, dan kunnen we de verdeling (1.23) fijnmazig kiezen op plaatsen waar deze afgeleide groot is en veel grover elders, zodat de benaderingsfout in ieder punt ongeveer even groot is.

Laat ∆ het rooster zijn op [a, b], zoals gedefinieerd in (1.23). We defini¨eren de ruimte van stuksgewijze polynomen M^k,p(∆) voor k ≥ 0 en −1 ≤ p ≤ k−1 door

M^k,p(∆) := { f ∈ C^p([a, b]) | f|[ti−1, ti] is een polynoom van graad ≤ k ∀ i } . (1.25) Deze verzameling bestaat dus uit stuksgewijze polynomen van graad k op de deelintervallen, die tesamen met hun eerste p afgeleiden continu zijn op de deelintervalgrenzen. (Voor p = −1 is er dus geen continu¨ıteit op de punten van ∆. Voor p = k − 1 noemen we deze stuksgewijze polynomen splinesvan graad k). Een polynoom van graad k wordt bepaald door k + 1 co¨effici¨enten, zodat de dimensie van de ruimte van stuksgewijze polynomen op de verdeling (1.23) gelijk is aan n(k + 1).

De continu¨ıteitseisen op de interne deelpunten geven (n − 1)(p + 1) lineaire vergelijkingen voor deze co¨effici¨enten, zodat de dimensie van M^k,p(∆) gelijk is aan n(k + 1) − (n − 1)(p + 1). Voor benaderingen in ruimten van stuksgewijze polynomen met p < k/2 kunnen we volgens de methode, die hierboven geschetst is, het volgende resultaat bewijzen:

Stelling. Als f ∈ C^k+1([a, b]), dan geldt voor de beste benadering π ∈ M^k,p(∆) van f met p < k/2 de foutschatting:

k ( d

dx)^m (f − π) k∞ ≤ Cm h^k−m+1 k f^(k+1) k∞ , 0 ≤ m ≤ k, (1.26)

(17)

waar C_meen constante is die alleen van m afhangt, maar niet van de functie of de gekozen verdeling

∆.

Bewijs. Voor m = 0 en p = −1 behoeven we in (1.24) slechts het maximum te nemen over alle deelintervallen. Als p ≥ 0 moeten we de interpolatiepunten met enige zorg kiezen. Als p = 0, dan kiezen we ξ₀ = 0 en ξ_k = 0, zodat t_i zowel voor het polynoom op [t_i−1, t_i] als dat op [t_i, t_i+1] een interpolatiepunt is. Op deze manier is de stuksgewijs polynomiale benadering automatisch continu in t_i, ∀ i . Als p > 0, dan interpoleren we zo, dat op de uiteinden van ieder deelinterval de functiewaarde en de eerste p afgeleiden van het interpolerende polynoom met die van de functie overeenstemmen (zie 1.4-5). De graad van het interpolerende polynoom moet dan wel groter zijn dan 2p. Door het maximum van de restterm (1.5) te nemen over alle deelintervallen vinden we dan (1.26).

Voor de foutschattingen op de afgeleiden kunnen we een analoog bewijs geven. Als p = 1 en als alle referentiesteunpunten onderling verschillend zijn, dan heeft de restrictie van f − π tot een deelinterval k + 1 nulpunten volgens (1.3). Volgens de stelling van Rolle heeft de afgeleide f^′− π^′ nog k nulpunten; π^′ interpoleert f^′ dus op k onderling verschillende punten (onbekend, maar wel binnen het deelinterval). Hierop kunnen we dus de reststelling (1.3) weer toepassen. Met een schatting analoog aan (1.24) volgt het bewijs van (1.26) ook in dit geval. Voor hogere afgeleiden herhalen we deze redenering.

In het geval p > 1 doen we in principe hetzelfde, alleen moeten we dan een punt, waarin de functiewaarde en p afgeleiden overeenstemmen als een p+1-voudig nulpunt tellen; de afgeleide f^′−π^′ heeft in datzelde punt dan een p-voudig nulpunt, etc.

N.B. 1. De “beste benadering” π in formule (1.26) behoeft geen continue afgeleide te hebben in de punten van ∆, zodat deze formule eigenlijk niet correct is. We moeten dan de supnorm interpreteren als het maximum van de supnormen over de deelintervallen.

N.B. 2. De uitspraak (1.26) van deze stelling is meestal ook waar als p ≥ k/2, maar een bewijs daarvan is veel ingewikkelder, zoals we zullen zien in het geval van kubische splines (M^3,2(∆)).

1.e Kubische Splines

Laat ∆ een verdeling van het interval [a, b] zijn, zoals gedefini¨eerd in (1.23), en laat f een voldoend gladde functie zijn op dat interval. De kubische spline interpolant S_∆ van f is het stuksgewijze polynoom van graad 3 uit M^3,2(∆), dat voldoet aan de eisen:

S_∆(t_i) = f (t_i) , i = 0, · · · , n, (1.27) en aan ´e´en van de drie volgende randvoorwaarden:

(a) S_∆^′ (a) = f^′(a) , S_∆^′ (b) = f^′(b), (b) S_∆^′′(a) = 0 = S_∆^′′(b)

(c) S_∆(a) = S_∆(b), S_∆^′ (a) = S_∆^′ (b), S_∆^′′(a) = S_∆^′′(b) .

(1.28)

Voorwaarde (a) noemen we een vaste rand, voorwaarde (b) een vrije rand en (c) een periodieke rand. In het geval (c) moet de benaderde functie f natuurlijk wel periodiek zijn met periode b−a.

De spline-benadering van een functie f is een benadering met minimale tweede afgeleide, zoals volgt uit:

Stelling: Als f^′′ ∈ C⁴([a, b]), dan geldt

k f^′′ − S^′′∆ k²2 = k f^′′ k²2 − k S∆^′′ k²2 . (1.29)

(18)

Bewijs.

k f^′′ − S∆^′′ k²2 = k f^′′ k2² − 2 (f^′′, S_∆^′′) + k S∆^′′ k²2

= k f^′′ k²2 − 2 (f^′′ − S∆^′′, S_∆^′′) − k S^′′∆k²2 . Het inprodukt kunnen we parti¨eel integreren,

Z b

a (f^′′ − S∆^′′) S_∆^′′ dx = − Z b

a (f^′ − S∆^′ ) S_∆^′′′ dx + [(f^′ − S∆^′ ) S_∆^′′ ]^b_a .

De stokterm hierin is nul onder ieder van de drie voorwaarden (1.28). De integraal kunnen we splitsen in de integralen over de n deelintervallen en opnieuw parti¨eel integreren,

Rb

a (f^′ − S∆^′ ) S^′′′_∆ dx = ^Pⁿ_i=1 ^R_t^tⁱ

i−1 (f^′ − S∆^′ ) S_∆^′′′ dx

= ^Pⁿ_i=1 ^R_t^t_i−1ⁱ (f − S∆) S_∆⁽⁴⁾ dx − [(f − S∆) S_∆^′′′]^t_tⁱ_i−1.

De integraal is nul, omdat S_∆ een derde graads polynoom is en de vierde afgeleide dus nul is, en de stokterm is nul wegens (1.27).

Uit formule (1.29) zien we, dat onder alle C²-functies op [a, b], die f interpoleren op ∆, de spline- benadering van f de kleinste “gemiddelde tweede afgeleide” heeft (de funktionaal ϕ 7→ k ϕ^′′ k2

minimaliseert bij gegeven k ϕ^′′− f^′′ k²).

1.f Praktisch rekenen met kubische splines

Neem nu aan, dat van een functie f de functiewaarden in de n + 1 punten van ∆ en de afgeleiden in de randpunten a en b gegeven zijn. Welke berekeningen moeten we uitvoeren om bij gegeven x ∈ [a, b] de waarde van de spline-benadering (met vaste rand) in dat punt te vinden. Voor dit probleem zijn verschillende oplossingen mogelijk. Hier zullen we twee manieren beschrijven, waarbij de spline op ieder deelinterval op een speciale manier gerepresenteerd wordt. Hiertoe voeren we de volgende notaties in,

h_i := t_i − ti−1 , (deelintervallengte)

f_i := f (t_i) , (functiewaarden in de interpolatiepunten) d_i := ^fⁱ ^{− f}_h ⁱ⁻¹

i , (differentiequoti¨enten )

S_i := S_∆| [ti−1, ti] (restriktie tot het i-de deelinterval)

en voor x ∈ [ ti−1, t_i] schrijven we t := (x − ti−1)/h_i. De restriktie S_i van S_∆tot [ t_i−1, t_i] kunnen we dan alsvolgt schrijven:

S_i(x) = tf_i+ (1 − t)fi−1+ h_it(1 − t){(λi−1− di)(1 − t) − (λi− di)t} , (1.30) waarin λi en λ_i−1nog te bepalen constanten zijn. Deze representatie is zo gekozen, dat S_i−1 en Si

continu op elkaar aansluiten in t_i−1, hun waarden aldaar zijn immmers gelijk aan f_i:

S_i(t_i−1+ 0) = f_i−1= S_i−1(t_i−1− 0) ( limiet van rechts²= limiet van links). (1.31) Hetzelfde geldt voor de afgeleiden in dat punt,

S_i^′(t_i−1+ 0) = λ_i−1= S_i−1^′ (t_i−1− 0). (1.32) De tweede afgeleide in t_i−1 moet ook continu zijn. Met de identiteit

S_i^′′(x) = 1

hi(λ_i−1(6t − 4) + λⁱ(6t − 2) + dⁱ(6 − 12t)) (1.33)

2Met f (x + 0) duiden we de limiet van rechts aan, f (x+0) := lim

tց0f(x + t) en met f (x−0) dus de limiet van links.

X Het geheel

Contents

0 Foutenanalyse

1 Interpolatie en Approximatie