1 Interpolatie

1 Interpolatie en Approximatie

In dit hoofdstuk bespreken we methoden om een gegeven functie van een veranderlijke te be- naderen met een (gemakkelijk berekenbare) functie uit een voorgeschreven klasse, zoals polynomen (veeltermen) of trigonometrische polynomen (sommen van machten van sinus en cosinus). Bij ap- proximatie (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 { x1, · · · , xn} is het polynoom van graad n − 1, dat gelijk is aan f in deze punten:

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

π_f(x) =

n

X

i=1

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

n

Y

j=1, j6=i

x − x_j xi − xj

. (1.2)

Opgave 1.1: Bewijs dit.

Voor het Lagrange-interpolatiepolynoom geldt de volgende foutschatting:

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

f (x) = π_f(x) +

n

Y

i=1

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

n! , (1.3)

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

f (x) = π_f(x) + g(x)

n

Y

i=1

(x − x_i)

en beschouw de functie

ϕt(x) := f (x) − π_f(x) − g(t)

n

Y

i=1

(x − xi) .

Deze heeft (minstens) n+1 onderling verschillende nulpunten binnen het interval int(x, x₁, · · · , x_n), als t 6= xi, ∀ 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

dit interval ligt. Aan het functie-onafhankelijke deel van de restterm, Πⁿ_i=1 (x − x_i), 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 ni−1 afgeleiden, d.w.z.

π^(j)_f (xi) = f^(j)(xi) , voor j = 0 , · · · , ni− 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 !

n

Y

i=1

(x − x_i)ⁿⁱ (1.5)

voor zekere ξ_x ∈ int(x, x₁, · · · , x_n) .

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, · · · , x_i+m} door

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

, (1.6) voor m = 1, 2, · · · .

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

m+i

X

j=i

f (xj) Qm+i

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

en dat hieruit volgt, dat f (x_i, x_i+1, · · · , x_i+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 {x₁, · · · , xn} schrijven als

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

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

f (x) − πf(x) = f (x1, · · · , xn, x)

n

Y

j=1

(x − xj) . (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₀, · · · , x_k) = f^k(ξ)

k! , ξ ∈ int(x₀, · · · , x_k). (1.10)

Uit formule (1.6) zien we dat de gedeelde differenties bij een gegeven verzameling steunpunten {x₁, 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 − x_n−1) f_n⁽ⁿ⁻¹⁾ + f_n−1⁽ⁿ⁻²⁾) (x − x_n−2) + · · ·) (x − x₁) + f₁⁽⁰⁾ (1.12a) of in pseudoPASCAL:

som := fn⁽ⁿ⁻¹⁾ ;

forj := n − 1 downto 1 do som := som ∗ (x − x_j) + 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 steunpun- ten {xi} 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 som- meren, 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 − x₁| ≤ | x − x₂| ≤ | x − x₃| ≤ · · ·

1.b Alternatieven voor het representeren en uitrekenen van een interpolatie- polynoom

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 − x1) (x − x2) · · · (x − x_k−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 · · · x₄ < 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¨effici¨enten van πf(x) = Σⁿ⁻¹_j=0 ajx^jbij gegeven functiewaarden πf(xi) = f (xi), i = 1, · · · , n. De co¨effici¨enten {aj | j = 0, · · · , n − 1} kunnen we dus eenvoudig bepalen door het oplossen van n lineaire vergelijkingen

n−1

X

j=0

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

oftewel

A













 a₀ a1

..

..

an−1













 :=

















1 x₁ · · · xⁿ⁻¹₁ 1 x2 · · · xⁿ⁻¹₂ ... ... · · · ... ... ... · · · ... 1 xn · · · xⁿ⁻¹_n





























 a₀ a1

..

..

an−1















=













 f (x₁) f (x2)

..

..

f (xn)















. (1.13b)

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 xi = 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 xi = cos(²ⁱ⁺¹₂₂ π), i = 0 , · · · , 10 (de nulpunten van het 10^de Chebyshev polynoom) T₁₀ vinden we een 11 × 11 matrix A met de conditiegetallen κ₂(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−1

X

j=0

aj 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 poly- nomen in plaats van machten van x. Zij bij gegeven steunpunten {xi | i = 1, · · · , n} de verzameling {Vj(x) | i = 0 , · · · , n − 1} een stelsel veeltermen met graad(Vj) = j en zo, dat

n

X

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−1

X

j=0

bjVj(x) , (1.16a)

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

B













 b₀ b1

..

..

bn−1













 :=















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

· · · ·

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



























 b₀ b1

..

..

bn−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 =

n

X

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 {x₁, · · · , x_n}.

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.

Lemma:

n−1

X

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−1

X

k=0

cos(2k + 1

2n mπ) = 1 2

n−1

X

k=0

( e^imπ2k+12n + e^{−imπ2k+12n} )

= 1

2

2n−1

X

k=0

e^imπ2n e^imπkn ;

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 {τ₁, · · · , τn} van het n-de Chebyshev polynoom Tn geldt:

2 n

n

X

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 Tj(cos t) = cos jt kunnen we formule (1.17) toepassen.

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

π_f =

n−1

X

j=0

c_jT_j met π_f(τ_k) =

n−1

X

j=0

c_jT_j(τ_k) = f (τ_k)

en volgens(1.8) moet dan gelden

cj = 2 n

n

X

k=1

f (τk) Tj(τk).

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

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

cj = 2 n

n

X

k=1

fk Tj(τk) = 2 n

n

X

k=1

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

= 1

n

n

X

k=1

fk exp(2k − 1

2n jπi) + 1 n

n

X

k=1

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

= 1

n

n

X

k=1

f_k exp(2k − 1

2n jπi) + 1 n

n

X

k=1

f_k exp(4n − 2k + 1

2n jπi)

= 1

n

n

X

k=1

fk exp(2k − 1

2n jπi) + 1 n

2n

X

l=n+1

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

= 1

n exp(j nπi)

2n−1

X

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.

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) − U_n−1(x) (dezelfde drieterms- recursierelatie als Tn maar met een ander begin).

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

c. Bepaal de relatie tussen de ontwikkeling πf = Σⁿ_i=1djUj van het Lagrange-interpolatiepoly- noom 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 k_2,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, Bn(x) :=

n

X

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 nu- merieke 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 x1, · · · , xn 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!

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₁< · · · < t_n−1 < t_n:= b , h := max

1<i<n t_i− t_i−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, ti] 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 interpol- eren blijkt goed genoeg voor het afleiden van een goede foutschatting voor dit soort benaderin- gen (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 [ti−1, ti] als het Lagrange- interpolatiepolynoom op de steunpunten

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

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

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

k

Y

j=0

(x − t_i−1− ξ_j(t_i− t_i−1)) | ≤ | t_i− t_i−1|^k | f^(k+1)(ξx) |

(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)

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 ti, ∀ 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_∆(ti) = f (ti) , 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²₂ = k f^′′ k²₂ − k S_∆^′′ k²₂ . (1.29)

Bewijs.

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

= k f^′′ k²₂ − 2 (f^′′ − S_∆^′′, S_∆^′′) − k S^′′_∆k²₂ . 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 = ^Pn_i=1 ^R_t^ti

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

= ^Pn_i=1 ^R_t^ti

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 ϕ^′′ k₂ minimaliseert bij gegeven k ϕ^′′− f^′′ k2).

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,

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

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

i , (differentiequoti¨enten )

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

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

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

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

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^′(ti−1+ 0) = λi−1= S_i−1^′ (ti−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) + λi(6t − 2) + di(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.

en de continu¨ıteitsvoorwaarde S_i^′′(t_i− 0) = S_i+1^′′ (t_i+ 0) vinden we voor de onbekenden λ_i het stelsel vergelijkingen

hi+1λi−1+ 2(hi+ hi+1)λi+ hiλi+1 = 3(hidi+1+ hi+1di) , (i = 1, · · · , n − 1) (1.34) Dit zijn n − 1 vergelijkingen voor n + 1 onbekenden. De onbekenden λ0 en λn vinden we uit de voorwaarden voor de vaste rand,

S₁^′(a + 0) = λ₀ = f^′(a) en S_n^′(b − 0) = λn= f^′(b). (1.35) Zo houden we een diagonaaldominant tridiagonaal stelsel van n−1 vergelijkingen over, dat gemakke- lijk en stabiel oplosbaar is met Gausseliminatie zonder rijverwisselingen. Om bij een gegeven rij functiewaarden de spline in een willekeurig punt uit te rekenen, moeten we dus eerst eenmalig een tridiagonaal stelsel oplossen (de rij λi berekenen). Daarna kunnen we met behulp van formule (1.30) de spline in ieder gewenst punt uitrekenen.

Opgave 1.7: Laat zien, dat een vrije randvoorwaarde de twee volgende vergelijkingen oplevert, 2λ₀+ λ₁ = 3d₁, en λ_n−1+ 2λn= 3dn, (1.36) en dat periodieke randvoorwaarden een extra vergelijking opleveren (behalve λ0 = λn),

hnλ1+ 2(h1+ hn)λn+ h1λn−1 = 3(d1hn+ dnh1). (1.37) In de representatie (1.30) hebben we de afgeleiden van de spline in de steunpunten als on- bekenden genomen en er vergelijkingen voor afgeleid uit de continu¨ıteitsvoorwaarden voor de tweede afgeleiden. Analoog kunnen we een representatie kiezen, waarin de tweede afgeleiden (de momenten Mi van S_∆ genaamd) de onbekenden zijn, waarvoor vergelijkingen worden afgeleid uit de continu¨ıteit van de eerste afgeleiden. Als we kiezen voor de representatie

S_i^′′(x) = M_i−1(1 − t) + M_it , (t := (x − t_i−1)/h_i) dan geeft twee maal integreren

Si(x) = 1

6h²_iM_i−1(1 − t)³+1

6h²_iMit³+ Ai(1 − t) + Bit , waarbij de integratieconstanten Ai en Bi bepaald worden uit de relaties

S∆(ti) = fi en S∆(ti−1) = fi−1. Zo vinden we op het deelinterval [ti−1, ti] de representatie

S_i(x) = (f_i−1

6h²_iM_i)t + (f_i−1−1

6h²_iM_i−1)(1 − t) +1

6h²_iM_i−1(1 − t)³+ 1

6h²_iM_it³. (1.38) waarin alleen de momenten als onbekenden voorkomen. Uit de continu¨ıteit van de eerste afgeleiden in t_i, (i = 1, · · · , n − 1) vinden we na enig rekenwerk de vergelijkingen

α_iM_i−1+ 2 M_i+ β_iM_i+1= 3m_i, (1.39) met

αi:= hi

h_i+ h_i+1 , βi := h_i+1

h_i+ h_i+1 , mi:= 2 d_i+1− di

h_i+ h_i+1.

voor i = 1, · · · , n − 1. Vanwege de vaste rand vinden we nog twee extra vergelijkingen, 2M₀ + M₁ = 3m₀, m₀ := _h²

1 (d₁ − f^′(a)),

M_n−1 + 2Mn = 3mn, mn := _h²_n(f^′(b) − dn) . (1.40)

De matrix in het zo verkregen stelsel vergelijkingen is diagonaaldominant, zodat ook dit stelsel een unieke oplossing heeft.

Opgave 1.8: Ga na dat in het periodieke geval (1.39) ook voor i = 0 geldt, als we de indices modulo n nemen. Ga ook na dat M₀ = Mn= 0 in het geval van een vrije rand.

Met behulp van de representatie (1.38) van de spline met de momenten als onbekenden kunnen we redelijk gemakkelijk een foutschatting afleiden. Hiervoor bewijzen we eerst het volgende lemma:

Lemma. Zij A de tridiagonale matrix behorende bij het stelsel vergelijkingen (1.39-40) voor de momenten,

A :=































2 1 0 · · · 0

α₁ 2 β₁ 0

0 α₂ 2 β₂ 0 ...

0 α₃ 2 β₃ 0

. .. ... . ..

... . .. . .. . .. 0

0 α_n−1 2 β_n−1

0 · · · 0 1 2































, (1.41)

dan geldt k x k∞≤ k A x k∞ voor elke vektor x ∈ IRⁿ⁺¹ (d.w.z. k A⁻¹k∞≤ 1 ).

Bewijs. Zij | x_k| := max_i | x_i| de grootste komponent van x , dan geldt voor de k-de komponent van de vector A x:

| (A x)k | = | 2xk + αkx_k−1 + βkx_k+1 | ≥ | xk | , omdat α_k+ β_k= 1. Hieruit volgt de bewering.

StellingAls f ∈ C⁴([a, b]) en als S_∆ de kubische spline interpolant van f met vaste rand is, dan voldoet deze tesamen met zijn eerste drie afgeleiden aan de foutschatting

k f⁽ⁱ⁾ − S_∆⁽ⁱ⁾ k∞ ≤ Ci h⁴⁻ⁱ k f⁽⁴⁾ k∞ , i = 0, 1, 2, 3, (1.42) waar h de maximale maaswijdte is van ∆ en Ci een constante is, die niet afhangt van f en van de verdeling ∆.

Bewijs. We zullen (1.42) eerst bewijzen voor i = 2 door te laten zien dat S_∆^′′ (een stuksgewijs lineair polynoom) dicht bij de lineaire interpolant van f^′′ligt. Daarna leiden we hieruit een schatting af voor de andere afgeleiden. Om de formules eenvoudig te houden, zullen we aannemen dat de verdeling ∆ equidistant is, ti− t_i−1= h ∀ i.

Voor het rechterlid van (1.39) geldt nu:

mi = fi−1 − 2fi + fi+1

h² = f^′′(ti) + h²

12 f⁽⁴⁾(ξi) . (1.43) Omdat αi= βi= ¹₂, vinden we uit (1.39) voor de verschillen vi := Mi− mi de relaties

v_i−1 + 4 v_i + v_i+1 = 2 m_i − m_i−1 − m_i+1 = h²

12 f⁽⁴⁾(η_i) . Met behulp van het lemma volgt hieruit

| vi | = | Mi − mi | ≤ h²

12 k f⁽⁴⁾ k∞ . (1.44)

Hieruit volgt dat het verschil tussen S_∆^′′ en de lineaire interpolant van f^′′ op de punten t_i en t_i−1 op het gehele deelinterval [ti−1, ti] begrensd is door (1.44). Het verschil tussen S_∆^′′ en f^′′ op dit deelinterval is dus hoogstens gelijk aan dit bedrag vermeerderd met de lineaire-interpolatiefout,

k S^′′_∆ − f^′′ k∞ ≤ ( 1 12 + 1

8) h² k f⁽⁴⁾ k∞ . (1.45)

Voor een schatting van f − S_∆gebruiken we de Lagrange-restterm (1.3) op de volgende manier.

Als g ∈ C² en g(x₁) = g(x₂) = 0, dan geldt

| g(x) | ≤ ¹₂ | (x − x₁)(x − x₂) g^′′(ξ_x) | ≤ 1

8(x₁− x₂)² k g^′′k_∞ , ∀ x ∈ [x₁, x₂]. (1.46) Het verschil f − Si is nul in de randpunten ti en t_i−1, zodat toepassing van (1.46) op dit verschil het volgende resultaat geeft:

k f − S_∆ k∞ ≤ h²

8 k S_∆^′′ − f^′′ k∞ ≤ 5

192 h⁴ k f⁽⁴⁾ k∞ .

Omdat het verschil f − Si nul is in de randpunten ti en t_i−1, heeft de afgeleide f^′− S_i^′ een nulpunt ζ binnen dat interval en kunnen we schrijven

f^′(x) − S_i^′(x) = Z x

ζ

f^′′(s) − S_∆^′′(s) ds.

Met behulp van (1.45) vinden we dus de schatting k f^′ − S_i^′ k∞ ≤ 5

24 h³ k f⁽⁴⁾ k∞ .

Tenslotte vinden we de schatting voor de derde afgeleide k f⁽³⁾− S_∆⁽³⁾ k_∞ door op het deelinterval [t_i−1, t_i] de derde afgeleide S_∆⁽³⁾ te vergelijken met de afgeleide van de lineaire interpolant van f^′′

op de steunpunten ti en t_i−1.

Deze stelling geldt ook voor een periodieke spline benadering van een periodieke functie (met een analoog bewijs), maar niet voor een spline benadering met een vrije rand. Bij een vrije rand is de tweede afgeleide in de randpunten a en b nul en kan dit i.h.a. onmogelijk een benadering geven van f^′′ van orde h². Het is echter wel te bewijzen, dat deze onnauwkeurigheid aan de rand naar het midden toe snel uitdempt en verwaarloosbaar klein wordt.

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.