3 Numerieke Integratie

3 Numerieke Integratie

3.a Probleemstelling

Gegeven een (voldoend gladde) functie f op een begrensd interval [a, b], bepaal een benadering voor de integraal

I :=

Z b

a f (t) dt (3.1)

en geef een foutschatting voor de verkregen benadering.

Een voor de hand liggende aanpak van dit probleem is het benaderen van f met een eenvoudig integreerbare functie, zoals een polynoom (zie hoofdstuk 1), en deze benadering exact te integreren.

Voorbeelden:

(1) Benader f met een lineair polynoom, f (x) = f (a) x − b

a + f (b) x − a

b + 1

2(x − a) (x − b) f^′′(ξ_x) en integreer,

Z _b

a

f (x) dx = b − a

2 (f (a) + f (b)) + 1 2

Z _b

a

(x − a) (x − b) f^′′(ξ_x) dx . (3.2) Omdat (x − a) (x − b) definiet is op [a, b], mogen we op de restterm de middelwaardestelling van de integraalrekening toepassen: Er is een punt η in [a, b] zodat

Z b a

(x − a) (x − b) f^′′(ξ_x) dx = f^′′(η) Z b

a

(x − a) (x − b) dx . Zo vinden we de trapeziumregel (met restterm)

Z b

a f (x) dx = b − a

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

12 f^′′(η) . (3.3)

(2) Benader f met een stukje Taylorreeks op het midden van het interval, f (x) = f (b + a

2 ) + (x − b + a

2 ) f^′(b + a 2 ) + 1

2 (x − b + a

2 )² f^′′(ξ_x) ,

integreer en gebruik opnieuw de middelwaardestelling (omdat (x − ^b+a₂ )² definiet is). Zo vinden we de middelpuntsregel (met restterm):

Z _b

a

f (x) dx = (b − a) f (b + a

2 ) + (b − a)³

24 f^′′(η) (3.4)

voor zekere η ∈ [a, b].

In het algemeen vinden we zo een (lineaire) integratief ormule Σⁿ_i=1 w_i f (t_i) als benadering van een integraal,

Z _b

a f (x) dx =

n

X

i=1

wi f (ti) + rest , (3.5)

waar {t₁, · · · , t_n} de steunpunten en {w₁, · · · , w_n} de gewichten heten. Bij gegeven steun- punten kunnen we het gewicht w_i bepalen door het i-de Lagrange-polynoom L⁽ⁿ⁾_i , cf. (1.2), te integreren (hetgeen i.h.a. veel werk is). Een eenvoudige methode om bij gegeven steunpunten de

gewichten te bepalen is de methode der onbepaalde co¨effici¨enten. Deze methode berust op het feit, dat interpolatie en dus ook integratie op n steunpunten exact is voor alle polynomen van graad kleiner dan n en dat het voldoende is gelijkheid in (3.5) te testen voor een basis in deze ruimte van polynomen:

n

X

i=1

w_i (t_i − a)^k = Z _b

a

(x − a)^k dx = (b − a)^k+1

k + 1 , k = 0, · · · , n − 1 . (3.6) Als de steunpunten onderling verschillend zijn, t_i 6= t_j (i 6= j), dan is de matrix in het linkerlid van (3.6) een reguliere Vandermondematrix, cf (1.13b), en heeft het stelsel vergelijkingen dus een eenduidige oplossing. Natuurlijk geldt hier bij het uitrekenen van de gewichten eenzelfde bezwaar als bij interpolatie, n.l. dat deze matrix meestal zeer slecht geconditioneerd is (tenzij n heel klein is).

Op grond van de restterm bij Lagrange-interpolatie (1.3) vinden we voor de rest in formule (3.5) de volgende foutschatting:

| rest | ≤ 1

n! max

a ≤ x ≤ b | f⁽ⁿ⁾(x) | Z _b

a n

Y

j=1

| x − t_i | dx ≤ C (b − a)ⁿ⁺¹ k f⁽ⁿ⁾ k_∞ , (3.7)

waar C een constante is, die alleen van n en de keuze van de steunpunten afhangt maar niet van f . Het kan voorkomen dat de integratieformule, gevonden bij een zekere verzameling van n steun- punten, ook exact is voor sommige polynomen van hogere graad. Als we bijvoorbeeld op het interval [0 , 1] de drie volgende steunpunten kiezen, t₁ = 0, t₂ = ¹₂ en t₃ = 1, dan vinden we met behulp van de methode der onbepaalde co¨effici¨enten de vergelijkingen:

w₁ + w₂ + w₃ = 1 , ¹₂ w₂ + w₃ = ¹₂ , 1

4 w₂ + w₃ = 1 3 .

De gewichten zijn dus w1 = w3 = 1/6 en w2 = 2/3. De gevonden integratieformule is de formule van Simpson. Deze formule is ook exact voor polynomen van graad 3 en zij is van orde 4. Hiervoor is dus een betere foutschatting dan (3.7) beschikbaar. Om in zo’n geval een betere foutschatting te vinden, gebruiken we de Reststelling van Peano:

Stelling. Zij L een lineaire functionaal op C^m([a, b]) zo, dat

L f = 0 , voor alle polynomen van graad kleinerdan m , (3.8) dan geldt

L F = Z _b

a

f^(m)(t) K(t) dt , K(t) := 1

(m − 1)! L_x (x − t)^m−1 Y (x − t) . (3.9) Hierin staat L_x(x − t)^m−1Y (x − t) voor de functionaal L toegepast op de functie

x 7→ (x − t)^m−1Y (x − t) (als functie van x).

Y is de Heaviside functie, Y (z) = 1 als z ≥ 0 en Y (z) = 0 als z < 0.

Bewijs. Pas L toe op de Taylorontwikkeling van f in a f (x) =

m−1

X

j=0

f^(j)(a) (x − a)^j

j! +

Z x a

(x − t)^m−1

(m − 1)! f^(m)(t) dt .

Aangezien L verdwijnt op het polynomiale gedeelte van deze ontwikkeling, volgt de bewering on- middellijk hieruit.

Opmerking. Een integratieformule, die exact is voor alle polynomen van graad kleiner dan m heet een formule van orde m.

We kunnen deze stelling toepassen op de boven afgeleide formule van Simpson. Deze formule is ook exact voor polynomen van graad 3. We vinden:

rest = L f = ^R₀¹ f (x) dx − ¹₆(f (0) + 4f (¹₂) + f (1))

= ₂₄^{1 R}₀¹ f⁽⁴⁾(t)^hR₀¹(x − t)³Y (x − t) dx − ²₃(¹₂ − t)³Y (¹₂− t) − ¹₆(1 − t)³ⁱdt . Aangezien

Z ₁

0 (x − t)³Y (x − t) dx = Z ₁

t (x − t)³dx = 1

4(1 − t)⁴, vinden we

K(t) := [ · · · ] = ( ₁

4(1 − t)⁴ − ¹₆(1 − t)³ = ₁₂¹ (1 − 3t) (1 − t)³, t ≤ ¹₂,

1

4(1 − t)⁴ − ²₃(¹₂ − t)³ − ₆¹(1 − t)³ = ₁₂¹ (3t − 2) t³, t ≥ ¹₂. Daar de “kern” K(t) negatief is voor alle t in het interval [0 , 1], mogen we de middelwaardestelling toepassen. Er is dus een ξ ∈ [0 , 1] zodat

rest = 1

24f⁽⁴⁾(ξ) Z 1

0

K(t) dt = − 1

2880 f⁽⁴⁾(ξ) . (3.10)

Opmerking. In dit geval is het toevallig eenvoudiger de restterm voor de formule van Simpson af te leiden uit de restterm voor het kubische interpolatiepolynoom π van f , waarvoor geldt

π(0) = f (0) , π(¹₂) = f (¹₂) , π^′(¹₂) = f^′(¹₂) , π(1) = f (1) . Op grond van formule (1.5) voor de restterm van dit interpolatiepolynoom vinden we

Z 1

0 f (x) − π(x) dx = 1 24

Z 1

0 f⁽⁴⁾(ξx) x (x − ¹₂)²(x − 1) dx = − 1

2880 f⁽⁴⁾(ξ) , omdat x (x −¹₂)²(x − 1) negatief is op het gehele interval [0 , 1].

3.b Gauss-integratie

We vragen ons nu af wat de maximale orde is, die we met een integratieformule van de vorm Σⁿ_i=1w_if (t_i) kunnen bereiken. Voor de eenvoud van de formules nemen we aan dat het integratiein- terval [−1 , 1] is. De integratieformule bevat 2n parameters, de gewichten w_i en de steunpunten ti . Door substitutie van de polynomen 1, x, x², · · · , x²ⁿ⁻¹ (met de methode der onbepaalde co¨effici¨enten) vinden we 2n niet-lineaire vergelijkingen voor de 2n parameters van de gezochte in- tegratieformule. We zullen via een omweg laten zien dat deze vergelijkingen precies ´e´en oplossing hebben en dat de maximale orde van een n-punts formule dus 2n is. Zo’n formule van maximale orde heet een Gauss-formule.

Bij een gegeven verzameling steunpunten {t₁, · · · , t_n} kunnen we via het lineaire stelsel (3.6) de gewichten {w₁, · · · , w_n} bepalen zo, dat de formule minstens van orde n is,

Z 1

−1 f (x) dx =

n

X

j=1

wj f (tj) als f een polynoom van graad kleiner dan n is.

Hiermee is de formule exact voor polynomen van graad kleiner dan n. Als de formule ook exact is voor polynomen met graad tussen n en 2n, dan moet dit in het bijzonder gelden voor de polynomen x^k Pn(x), k = 0, · · · , n − 1 ,

Z ₁

−1

P_n(x) x^k dx =

n

X

i=1

w_i P_n(t_i) t^k_i . (3.11)

Maar we weten ook dat het n-de Legendre polynoom loodrecht staat op alle polynomen van lagere graad, d.w.z.

Z 1

−1

P_n(x) x^k dx = 0 , k = 0, · · · , n − 1 .

Dit betekent, dat het rechterlid van (3.11) nul is voor alle k, 0 ≤ k < n. Omdat de Vandermon- dematrix met elementen t^k_i niet singulier is, impliceert dit, dat P_n(t_i) = 0, i = 1, · · · , n. Zonder expliciet aan een zeer slecht geconditioneerde Vandermondematrix te moeten rekenen, kunnen we de gewenste steunpunten dus vinden als de nulpunten van het n-de Legendre polynoom. Deze nulpunten zijn enkelvoudig en liggen alle binnen het interval (−1, 1). Bovendien zien we hieruit, dat er op dit interval precies ´e´en integratieformule bestaat van orde 2n met n steunpunten.

Om de restterm van de gevonden n-punts Gaussformule te vinden voor een willekeurige functie f ∈ C²ⁿ([−1 , 1]) interpoleren we f met een polynoom π van graad 2n − 1 d.m.v. Hermite- interpolatie op de nulpunten van {t1, · · · , tn} van Pn. Hiervoor geldt volgens (1.5)

f (x) = π(x) +

n

Y

i=1

(x − t_i)² f⁽²ⁿ⁾(ξx) (2n)! .

Als we deze vergelijking integreren, de integraal over het polynoom π vervangen door de inte- gratieformule en op de restterm de middelwaardestelling toepassen, vinden we:

Z ₁

−1

f (x) dx =

n

X

j=1

w_j f (t_j) + f⁽²ⁿ⁾(ξ_x) (2n)!

Z ₁

−1 n

Y

i=1

(x − t_i)² dx .

Omdat het polynoom Πⁿ_i=1 (x − t_i) van gelijke graad is als P_n en dezelfde nulpunten heeft, zijn beide polynomen gelijk op een constante factor na. Deze vinden we uit de Rodrigues relatie voor P_n ,

P_n(x) = 1 2ⁿ n!

dⁿ

dxⁿ (x² − 1)ⁿ = (2n)!

2ⁿ (n!)² xⁿ + termen van lagere orde = (2n)!

2ⁿ (n!)²

n

Y

i=1

(x − t_i) De integraal in de restterm kunnen we hiermee alsvolgt uitrekenen:

Z ₁

−1 n

Y

i=1

(x − t_i)² dx = 2²ⁿ(n!)⁴ ((2n)!)²

Z ₁

−1

P_n(x)² dx = 2²ⁿ⁺¹(n!)⁴ (2n + 1) ((2n)!)² . Zo vinden we dus de n-punts Gauss formule

Z ₁

−1

f (x) dx =

n

X

j=1

w_j f (t_j) + f⁽²ⁿ⁾(ξ) 2²ⁿ⁺¹(n!)⁴

(2n + 1) ((2n)!)³ . (3.12) De gewichten van een Gauss formule zijn positief, immers

wi = Z ₁

−1 n

Y

j=1, j6=i

(x − tj

t_i− t_j)² dx > 0 , (3.12a) daar de formule exact is voor alle polynomen van graad ≤ n. Overigens zal men nooit de gewichten met behulp van deze expressie uitrekenen. Via de recurrente betrekking voor de Legendre poly- nomen kan men voor de gewichten de volgende formule vinden:

w_i = 2

(1 − t²_i) [P_n^′(t_i)]² . (3.13) In de literatuur zijn ook aan een aantal andere typen integratieformules een naam verbonden:

Een n-punts Newton-Cotes formule (n ≥ 2) op het interval [a, b] is de formule met maximale orde op n equidistante steunpunten {t₁ := a, t₂, · · · , t_n := b}. Voor n > 6 zijn dit soort formules numeriek minder aantrekkelijk, omdat sommige gewichten dan negatief zijn. Voorbeelden zijn de trapeziumregel (n = 2) en de regel van Simpson (n = 3).

Opgave 3.1: Laat zien dat de orde van zo’n formule n is voor even n en n + 1 voor oneven n.

Een n-punts Lobatto formule (n ≥ 2) op het interval [a, b] is de formule met maximale orde op de n steunpunten {t₁ := a, t₂, · · · , t_n := b}, waarbij de steunpunten t₂, · · · , t_n−1 zo gekozen worden, dat de orde maximaal is. De orde van zo’n n-punts formule is 2n − 2. De steunpunten van deze formules zijn verbonden met de nulpunten van P_n−1^′ (als [a, b] = [−1, 1]) en ze liggen alle binnen het interval. De gewichten zijn altijd positief. De restterm heeft de vorm C f⁽²ⁿ⁻²⁾(ξ).

Voorbeelden zijn opnieuw de trapeziumregel en de regel van Simpson.

Een n-punts Radau formule (n ≥ 2) op het interval [a, b] verkrijgen we door ofwel t₁ = a ofwel t_n = b te kiezen en vervolgens de overige steunpunten en gewichten zo bepalen, dat de orde maximaal is. De orde van zo’n formule is 2n − 1 (in feite hebben we er twee, die elkaars gespiegelde zijn). Ook hier liggen de steunpunten binnen het integratie-interval en zijn de gewichten positief.

De restterm heeft de vorm C f⁽²ⁿ⁻¹⁾(ξ).

3.c Samengestelde integratieformules

Als we de afbreekfout, die ontstaat door toepassing van een integratieformule op een functie, te groot vinden, kunnen we overgaan op een formule van hogere orde, maar ook kunnen we het integratie-interval in stukken verdelen en op ieder stuk afzonderlijk de integratieformule toepassen.

Als voorbeeld kiezen we de samengestelde middelpuntsregel M_n. We verdelen het integratie-interval in n gelijke stukken van lengte h := (b − a)/n,

Z _b

a

f (x) dx =

n

X

j=1

h f (a + jh − h 2) +

n

X

j=1

h³

24 f^′′(ξ_j) = M_n(f ) + R_n(f ) , (3.14) waar ξj een punt is in het deelinterval [a + jh − h, a + jh]. We kunnen de afbreekfout Rn dus willekeurig klein maken, door h voldoend klein te kiezen. Evenals bij interpolatie is het ook hier vaak gunstiger de orde van de integratieformule beperkt te houden en de nauwkeurigheid tot de gewenste waarde op te voeren door de lengte van de deelintervallen, waarop deze wordt toegepast, te verkleinen. Een dergelijke strategie convergeert voor een grotere klasse van functies, dan integratie d.m.v. een serie (enkelvoudige) formules van steeds hogere orde. Bovendien heeft de restterm van een samengestelde formule op ieder deelinterval dezelfde vorm en kunnen we deze gebruiken om de afbreekfout te schatten. In het voorbeeld (3.14) is de restterm een Riemannsom voor de integraal van de tweede afgeleide,

n

X

j=1

h f^′′(ξ_j) − Z b

a

f^′′(x) dx → 0 , voor n → ∞ ,

zodat we de rest ook kunnen schrijven als³ R_n =

n

X

j=1

h³

24 f^′′(ξ_j) = h²

24 (f^′(b) − f^′(a)) + o(h²) , (h → 0) . (3.15) We concluderen hieruit dat de restterm in eerste benadering gelijk is aan een constante maal h², met een constante, die niet van h (maar wel van f ) afhangt,

R_n = C h² + o(h²) , (h → +0).

Als de o(h²)-term voldoend klein is t.o.v. C h², dan kunnen we C schatten door M te berekenen voor twee verschillende waarden van h (of n) en hieruit de onbekende integraal te elimineren:

I = M_n + C h² + o(h²)

= M_2n + C(^h₂)² + o(h²) . Aftrekken geeft:

M_n − M_2n = 3

4 C h² + o(h²) We zien hieruit,

R_2n = 1

4 C h² + o(h²) = 1

3(M_n − M_2n) + o(h²). (3.16) De belangrijke conclusie, die we hieruit kunnen trekken, is dat we de afbreekfout in M_2n kunnen schatten door naast M_2n ook M_n te berekenen. Dit geeft ons de mogelijkheid om een automatisch afbreekcriterium te formuleren, dat een benadering levert met vooraf voorgeschreven precisie onder de voorwaarde, dat de o(h²)-term in (3.15) klein is:

< kies startwaarde voor n > ; < bereken M_n > ; repeat

n := 2 ∗ n;

< bereken M_n >;

until| M_n − M_n/2| < eps .

(3.17)

Als formule (3.16) een goede schatting van de fout oplevert, dan kunnen we deze ook gebruiken om de verkregen benadering te verbeteren:

I = M_2n + 1

3(M_2n − M_n) + o(h²). (3.18)

Dit proces heet extrapolatie.

Opgave3.2: Laat zien dat de algoritme (3.17) niet altijd het gewenste resultaat aflevert. Kies als voorbeeld f = (x −¹₄)(x − ¹₂)(x −³₄) en start met n = 1.

3Het grote O-symbool en het kleine o-symbool zijn gedefinieerd alsvolgt.

f(x) = O(g(x)) (x → a + 0) ,

als er een halfopen interval (a, a + δ] bestaat (met δ > 0) en een positieve constante C zodat | f (x) | ≤ C | g(x) | voor alle x ∈ (a, a + δ].

f(x) = o(g(x)) (x → a + 0) , als lim

xa, x→a

f(x)

g (x) = 0 ,

of anders gezegd, als er bij iedere ε > 0 een δ > 0 bestaat, zodat | f (x) | ≤ ε | g(x) | voor alle x ∈ (a, a + δ].

Opgave 3.3: Definieer P als het lineaire polynoom, dat voldoet aan P (h²) = M_n en P (h²/4) = M_2n. Laat zien, dat P (0) precies de extrapolatie van formule (3.18) levert,

P (0) = M2n + 1

3(M2n − Mn) .

Conclusie: De afbreekfout is een functie van h². Deze functie benaderen we met een lineair poly- noom (in h²). De waarde van dit polynoom in nul is de waarde van de extrapolatie. Extrapolatie is dus niets anders als een speciale vorm van polynoombenadering. Omdat de afbreekfout een functie van h² is, spreken we in dit geval over h-kwadraat extrapolatie.

Extrapolatie is een veelgebruikt hulpmiddel in de numerieke analyse om een berekende bena- dering te verbeteren (niet alleen voor integralen!). Extrapolatie kunnen we altijd toepassen als er een (asymptotisch) verband is tussen de fout en een diskretisatie-parameter (b.v. de lengte van het integratie-interval). Een ander voorbeeld is het volgende:

Uit de Taylorontwikkeling van f (x ± h) van een voldoend gladde functie f vinden we voor de centrale differentie als benadering van de eerste afgeleide de formule

D_hf := f (x + h) − f (x − h)

2h = f^′(x) + 1

3h²f⁽³⁾(x) + 1 12h

Z _h

−h

f⁽⁵⁾(t) (x − t)⁴dt , zodat

f^′(x) = D_hf + C h² + O(h⁴) . (3.19) Van dit verband tussen de afbreekfout en de stapgrootte h kunnen we gebruik maken on de schatting van de afgeleide te verbeteren in het volgende getallenvoorbeeld. Voor de functiewaarden van log(x) in de buurt van x = 2 hebben we de volgende waarden, korrekt afgerond tot 5 decimalen (en dus met een afrondfout kleiner dan .000005):

x f(x)

1.8 .58779 1.9 .64185 2.0 .69315 2.1 .74194 2.2 .78846

Voor de afgeleide f^′(2) vinden we met centrale differenties de benaderingen (met tussen haakjes de uiterste grenzen voor de afrondfouten in de berekende waarden)

D_0.1 f_|x=2 = .50045 (±.00005), D_0.2 f_|x=2 = .501675 (±.000025).

De geschatte afbreekfout in D_0.1 f_|x=2 is dus 1

3(D_0.1 f_|x=2 − D_0.2 f_|x=2 ) = −.00040633 (±.000025) en als we D_0.1 f_|x=2 hiermee verbeteren vinden we

f^′(2) ≈ .50004 (±.000075) .

Opgave 3.4: Ga na, dat de grenzen voor de afrondfouten in bovenstaande formules kloppen.

Opgave 3.5: Stel dat we in bovenstaand voorbeeld ook f (2.3) en f (1.7) kennen, dan zouden we f^′(2) ook op basis van D_0.1 en D_0.3 kunnen benaderen. Laat zien, dat de formule hiervoor is:

f^′(2) = D_h f + 1

8 (D_h f − D_3h f ) + O(h^4).

Aangezien we (3.19) voor voldoend gladde f verder kunnen ontwikkelen, f^′(x) = D_hf + C₁h² + C₂h⁴ + O(h⁶) , kunnen we uit beide schattingen voor f^′(2) de h⁴-term elimineren. Leid de formule hiervoor af.

Overigens merken we op, dat de nauwkeurigheid van de boven gegeven waarden van ln(x) geen hogere orde extrapolatie rechtvaardigen.

3.d Romberg integratie

Systematisch gebruik van extrapolatie bij numerieke integratie kan leiden tot een effici¨ent integrati- eschema van hoge orde voor een grote klasse van functies. We zullen dat laten zien aan de hand van het integratieschema van Romberg. Dit gaat uit van de asymptotische-reeksontwikkeling van de samengestelde trapeziumregel. Als we het integratie-interval [a, b] verdelen in n deelintervallen van lengte h := (b − a)/n en op ieder deelinterval de trapeziumregel (3.3) toepassen, dan vinden we de samengestelde trapeziumregel

T_h := h

2 (f (a) + f (b)) + h

n−1

X

j=1

f (a + jh) . (3.20)

Als f ∈ C^2k+1([a , b]), dan heeft T_h een (asymptotische) reeksontwikkeling in even machten van h, T_h = I + c₁h² + c₂h⁴ + · · · + c_kh^2k + O(h^2k+1) , I :=

Z b

a f (x) dx , (3.21) waarin c₁, · · · , c_kconstanten zijn, die alleen afhangen van f (maar niet van h). Het bewijs hiervan zullen we geven in §3h over de sommatieformule van Euler-McLaurin. Op grond van het bestaan van deze reeks (waarvan we de co¨effici¨enten i.h.a. niet kennen) kunnen we herhaald extrapoleren en zo successievelijk de termen van orde h², h⁴, etc. in de afbreekfout van de samengestelde trapeziumregel elimineren.

T_2h = I + 4c₁h² + 16c₂h⁴ + 64c₃h⁶ + · · · T_h = I + c₁h² + c₂h⁴ + c₃h⁶ + · · · aftrekken geeft

T_2h− T_h = 3c₁h² + 15c₂h⁴ + 63c₃h⁶ + · · · Eliminatie (h²-extrapolatie) van de term van orde h² geeft dus

T_h¹ := T_h + 1

3 (T_h − T_2h) = I − 4c₂h⁴ − 20c₃h⁶ + · · · = I + c¹₂h⁴ + c¹₃h⁶ + · · · De afbreekfout in T_h¹ is weer een asymptotische reeks in machten van h² beginnend met een term van orde h⁴. Deze term kunnen we op analoge wijze elimineren (h⁴-extrapolatie) uit T_2h¹ en T_h¹,

T_2h¹ − T_h¹ = 15c₂¹h⁴ + 63c¹₃h⁶ + · · · , zodat

T_h² := T_h¹ + 1

15 (T_h¹ − T_2h¹ ) = I − 16

5 h⁶ + · · · = I + c²₃h⁶ + c²₄h⁸ + · · · In het algemeen kunnen we zo de dominante term van orde h^2k elimineren met

T_h^k := T_h^k−1 + 1

4^k− 1 (T_h^k−1 − T_2h^k−1) . (3.22)

Door extrapolaties te berekenen van steeds hogere orde kunnen we het volgende driehoekige tableau van benaderingen opbouwen,

T_4h T_2h T_2h¹ T_h T_h¹ T_h²

T_h/2 T_h/2¹ T_h/2² T_h/2³

· · · ·

(3.23)

De afbreekfout in de k-de kolom is van orde h^2k, mits de ge¨ıntegreerde functie voldoend glad is.

In de praktijk kunnen we de elementen van dit tableau rij voor rij met weinig moeite berekenen uit het resutaat van de samengestelde trapeziumregel voor de stapgrootten h, h/2, h/4, · · · . Als we T_h berekend hebben, kunnen we ons bij de berekening van T_h/2 werk besparen door op te merken, dat T_h/2 berekend kan worden uit T_h en de waarden van f op ”nieuwe” steunpunten:

T_h/2 = h 2



¹₂ f (a) + ¹₂ f (b) +

2n−1

X

j=1

f (a + jh/2)





= h

2



¹₂ f (a) + ¹₂ f (b) +

n−1

X

j=1

f (a + jh) +

n

X

j=1

f (a + jh − h/2)





= ¹₂ T_h + h 2

n−1

X

j=1

f (a + jh − h/2) .

De steunpunten van T_h zijn ook steunpunten van T_h/2. In deze steunpunten behoeven we f niet opnieuw te evalueren (noch de functiewaarden in het geheugen te bewaren) voor de berekening van T_h/2. Door deze eigenschap is een extrapolatieschema gebaseerd op intervalhalvering voor de trapeziumregel zeer aantrekkelijk. Andere strategie¨en voor de keuze van de stapgrootte worden ook wel gebruikt, b.v. de rij h, h/2, h/3, h/4, h/6, h/8, h/12, h/16, · · · .

Bij het rekenen op een rekenmachine in de gebruikelijke “floating point arithmetiek” wordt de

“diepte” van de extrapolatie natuurlijk beperkt door (afrond- en afbreek)fouten bij de evaluatie van de te integreren functie en door afrondfouten in de berekende trapeziumsommen. Als de correctie T_h^k−1 − T_2h^k−1 kleiner wordt, dan de afrondfout in T_h^k−1 , dan heeft verder extrapoleren natuurlijk weinig zin, daar formule (3.22) dan haar praktische geldigheid heeft verloren. In het tableau (3.23) kunnen we dit i.h.a. goed zien aan het feit, dat de geschatte fout in de k-de kolom vanaf zekere k niet meer met een factor 2^2k afneemt. Tenslotte merken we op, dat we de extrapolaties in (3.23) beter met behulp van formule (3.22) kunnen berekenen dan met de (equivalente) formule

T_h^k = 4^k T_h^k−1 − T_2h^k−1 4^k − 1 , omdat (3.22) i.h.a. een kleinere afrondfout induceert in T_h^k.

3.e Voorbeeld van het Rombergschema

Een voorbeeld van de effectiviteit van het schema van Romberg voor het benaderen van de integraal van een (gladde) functie vinden we bij de berekening van een benadering van de integraal voor π,

Z 1 0

4

1 + x² dx = π ,

met Rombergintegratie. In de eerste tabel zijn de benaderingen afgedrukt, verkregen met de trape- ziumregel op 2^k deelintervallen, en achtereenvolgens de h²-, h⁴-, h⁶-, h⁸- en h¹⁰-extrapolaties hier- van.

k t[0, k] t[1, k] t[2, k] t[3, k] t[4, k] t[5, k]

1 3.000000000000000

2 3.100000000000009 3.133333333333340

3 3.131176470588244 3.141568627450994 3.142117647058839

4 3.138988494491102 3.141592502458721 3.141594094125907 3.141585783761897

5 3.140941612041402 3.141592651224840 3.141592661142582 3.141592638396816 3.141592665277742

6 3.141429893174987 3.141592653552848 3.141592653708045 3.141592653590038 3.141592653649624 3.141592653638256 7 3.141551963485654 3.141592653589214 3.141592653591644 3.141592653589797 3.141592653589797 3.141592653589740

In de tweede tabel zijn de verschillen afgedrukt tussen de bovenstaande berekende waarden en de exakte waarde π en de ²log ervan. Het is duidelijk, dat de fouten in de eerste kolom met een faktor 4 afnemen bij verdubbeling van de stapgrootte h. Dit bevestigt het feit, dat de overwegende term in de fout gelijk is aan een constante maal h². In de tweede kolom verwachten we een afname van de fout met een faktor 16 (van de ²log vermindert dan met 4), maar we vinden een afname met 2⁶, hetgeen erop wijst, dat de co¨effici¨ent van h⁴ in de asymptotische reeks voor de fout nul is. Inderdaad geldt f^′′′(1) − f^′′′(0) = 0 voor f (x) = 4/(1 + x²). Desalniettemin wordt in het (standaard) Romberg schema toch h⁴-extrapolatie toegepast, met als resultaat, dat de fout niet daalt. De h⁶-afhankelijkheid blijft echter wel bestaan, zodat de h⁶-extrapolatie in de volgende kolom de fout inderdaad verder reduceert. In de vierde kolom zien we een h⁸-afhankelijkheid en in de vijfde een (vage) h¹⁰-afhankelijkheid. Extrapolatie van nog hogere orde heeft weinig zin, daar de verschillen in grootte vergelijkbaar zijn geworden met de afrondfouten.

j 0 1 2 3 4 5

k fout ²log fout ²log fout ²log fout ²log fout ²log fout ²log 1 -1.42e-001 -2.8

2 -4.16e-002 -4.6 -8.26e-003 -6.9

3 -1.04e-002 -6.6 -2.40e-005 -15.3 5.25e-004 -10.9

4 -2.60e-003 -8.6 -1.51e-007 -22.7 1.44e-006 -19.4 -6.87e-006 -17.2

5 -6.51e-004 -10.6 -2.36e-009 -28.7 7.55e-009 -27.0 -1.52e-008 -26.0 1.17e-008 -26.4

6 -1.63e-004 -12.6 -3.69e-011 -34.7 1.18e-010 -33.0 2.42e-013 -41.9 5.98e-011 -34.0 4.85e-011 -34.3 7 -4.07e-005 -14.6 -5.83e-013 -40.6 1.85e-012 -39.0 0.00e+000 -332.2 0.00e+000 -332.2 -5.68e-014 -44.0

In de praktijk beschikken we niet over de exakte waarde, maar we kunnen wel de fout schatten met de formule (t[j, i] − t[j, i − 1])/(4^j+1− 1). In de derde tabel worden deze schatting en de²log ervan afgedrukt. In deze tabel vertoont de afname van de fout hetzelfde patroon als in de vorige.

j 0 1 2 3 4 5

i schatting ²log schatting ²log schatting ²log schatting ²log schatting ²log schatting ²log 2 3.33e-002 -4.9

3 1.04e-002 -6.6 5.49e-004 -10.8

4 2.60e-003 -8.6 1.59e-006 -19.3 -8.31e-006 -16.9

5 6.51e-004 -10.6 9.92e-009 -26.6 -2.27e-008 -25.4 2.69e-008 -25.1

6 1.63e-004 -12.6 1.55e-010 -32.6 -1.18e-010 -33.0 5.96e-011 -34.0 -1.14e-011 -36.4

7 4.07e-005 -14.6 2.42e-012 -38.6 -1.85e-012 -39.0 -9.47e-016 -49.9 -5.85e-014 -44.0 -1.18e-014 -46.3

In de vierde tabel worden het resultaat van de trapeziumregel op n deelintervallen, de fout en de

2log van de fout nog eens bijeen gezet.

n trap-regel fout ²log

1 3.000000000000000 -1.42e-001 -2.8 2 3.100000000000009 -4.16e-002 -4.6 4 3.131176470588244 -1.04e-002 -6.6 8 3.138988494491102 -2.60e-003 -8.6 16 3.140941612041402 -6.51e-004 -10.6 32 3.141429893174987 -1.63e-004 -12.6 64 3.141551963485654 -4.07e-005 -14.6

In de vijfde tabel worden het resultaat van de trapeziumregel op n deelintervallen, de h-kwadraat extrapolatie ervan, de fout in deze extrapolatie en de ²log van deze fout nog eens bijeen gezet.

Zoals reeds opgemerkt, is deze fout toevallig van orde h⁶.

n trap-regel extrapolatie fout ²log 1 3.000000000000000

2 3.100000000000009 3.133333333333340 -8.26e-003 -6.9 4 3.131176470588244 3.141568627450994 -2.40e-005 -15.3 8 3.138988494491102 3.141592502458721 -1.51e-007 -22.7 16 3.140941612041402 3.141592651224840 -2.36e-009 -28.7 32 3.141429893174987 3.141592653552848 -3.69e-011 -34.7 64 3.141551963485654 3.141592653589214 -5.83e-013 -40.6

3.f Integratie met veranderlijke stapgrootte

Stel, dat we op het interval [0, 1] de integraal willen bepalen van de functie sin (_.01+x¹ ) met een afbreekfout van ten hoogste ε := 10⁻⁶. Als we hiervoor de samengestelde middelpuntsregel willen gebruiken met stapgrootte h (op n = 1/h deelintervallen), dan moeten we volgens (3.14) in ieder deelinterval eisen:

| h³ 24

n

X

j=1

f^′′(ξ_j) | ≤ 10⁻⁶ , ξ_j ∈ (jh − h, jh).

Aan deze voorwaarde is zeker voldaan, als h²

24 max

0 ≤ x ≤ 1 | f^′′(x) | ≤ 10⁻⁶ . Met

f^′′(x) = − 2(.01 + x)⁻³ cos 1

.01 + x − (.01 + x)⁻⁴sin 1 .01 + x betekent dit

k f^′′ k_{∞, [0, 1]} ≤ 1.02₁₀8 , en dus h ≤ 5₁₀− 7 .

Met zo’n equidistant integratieschema zouden we dus twee millioen deelintervallen moeten ge- bruiken en evenveel evaluaties van sin (_.01+x¹ ) ! Op het deelinterval [¹₂, 1] is de afbreekfout echter veel kleiner en kunnen we volstaan met een veel grotere deelintervallengte h. Omdat

k f^′′ k

∞, [1 2 , 1]

≤ 32 ,

vinden we, dat h := .001 voldoend klein is om een fout van 510− 7 (de helft van de toegestane fout over het gehele interval) te garanderen. Evenzo is op [¹₄, ¹₂] een stapgrootte h := .0001 voldoende om een fout van 2.5₁₀− 7 over dat deelinterval te garanderen.

Opgave3.6: Bepaal de stapgrootte h, waarmee de integraal,^R₀² | cos x | dx berekend moet worden om een afbreekfout van ten hoogste ₁₀− 8 te garanderen, als we een samengestelde tweepunts Gauss-formule gebruiken op equidistante deelintervallen. Aanwijzing: neem de foutschatting voor het deelinterval dat π/2 bevat apart.

Beantwoord dezelfde vraag voor het geval dat het interval eerst wordt opgebroken in de stukken [0, π/2] en [π/2, 2].

In feite is het gebruik van een uniforme stapgrootte dus geen slimme strategie voor integralen met een sterk wisselende gladheid. Veel beter lijkt het om de stapgrootte aan te passen aan de plaatselijke gladheid van de integrand: verdeel het integratie-interval in stukken van lengte h_i die zo gekozen worden, dat de fout op het i-de deelinterval kleiner dan εh_i is; de totale fout is dan hoogstens gelijk aan ε maal de intervallengte. In de praktijk willen we natuurlijk niet, dat we voor iedere integrand opnieuw een opdeling van het integratie-interval moeten opgeven. Een automatische strategie voor de keuze van de stapgrootte kunnen we baseren op het idee, dat we de afbreekfout kunnen schatten door twee verschillende integratieformules te gebruiken (of twee verschillende stapgrootten), zoals we reeds zagen in 3c. Als we dan een integraal over het interval [a, b] moeten benaderen met tolerantie van ten hoogste ε(b − a), dan kunnen we op een deelinterval [y, y + h] een fout toestaan van ten hoogste εh. Als de geschatte fout inderdaad kleiner dan εh is

op dat interval, dan accepteren we de gevonden benadering, tellen deze op bij het reeds berekende stuk en gaan verder met het volgende deelinterval. Als de gevraagde tolerantie niet gehaald wordt, verdelen we het interval in stukken, b.v. in [y, y + h/2] en [y + h/2, y + h], en passen de procedure toe op ieder van de stukken.

Een voorbeeld van zo’n formulepaar zou trapezium- en middelpuntsregel kunnen zijn. Hiervan weten we immers, cf. (3.2-3),

Z _y+h

y

f (x) dx = h f (y + h/2) + h³

24 f^′′(ξ) = h

2 (f (y) + f (y + h)) − h³

12 f^′′(η) , zodat ^h₆ (f (y) − 2f (y + h/2) + f (y + h)) een goede schatting is van de afbreekfout in de middelpuntsregel als benadering van de integraal over dit interval. We accepteren h f (y + h/2) dus als voldoend goede benadering voor de integraal over dit interval als deze (lokale) foutschatting kleiner is dan de (lokale) tolerantie εh en tellen dit op by het reeds berekende stuk van de integraal.

Als de lokale fouttolerantie niet gehaald wordt, delen we het interval in de stukken [y, y + h/2]

en [y + h/2, y + h] en passen we de formules toe op beide deelintervallen apart. Merk op dat we in een goede implementatie van deze methode bij intervalhalvering de functie alleen opnieuw hoeven te evalueren in de nieuwe middelpunten, omdat het oude midden, waar we de functie reeds kenden, randpunt van de beide deelintervallen geworden is. Overigens let niets ons om i.p.v. het resultaat van de middelpuntsregel, h f (y + h/2), het (nauwkeuriger) resultaat van Simpson’s regel,

h

6 (f (y) + 4f (y +h/2) + f (y +h)), toe te voegen aan het reeds berekende stuk van de integraal, als de lokale fouttolerantie gehaald wordt (de benodigde functiewaarden hebben we toch al berekend).

De implementatie van een dergelijke adaptieve integratiemethode is niet geheel triviaal. In de eerste plaats willen we het aantal functie¨evaluaties minimaliseren. Dus moeten we een formule (of formulepaar) kiezen, waarvan de steunpunten zo geplaatst zijn, dat ze bij verfijning steunpunt blijven en de reeds berekende functiewaarde opnieuw gebruikt kunnen worden. Een adaptieve routine kan slechts aan ´e´en deelinterval tegelijk werken. Als de routine dus besluit een interval op te delen in stukken, dan moet de informatie betreffende de rest opgeborgen kunnen worden op een stapel. Ook moeten er grenzen worden ingebouwd, zodat er niet tot in het oneindige verfijnd wordt, als de lokale fouttest telkens negatief blijkt uit te vallen.

De keuze voor een integratieroutine met zelfzoekende stapgroote of een samengestelde formule met vaste stapgrootte zal afhangen van de integrand, maar ook van het doel waarvoor ge¨ıntegreerd wordt. Als er ´e´en integraal berekend moet worden, dan ligt een adaptieve methode voor de hand, waarbij alleen maar de maximaal toelaatbare fout behoeft te worden gespecificeerd. Als de in- tegrand echter nog van andere parameters afhangt en als deze integraal voor een groot aantal waarden van deze parameters berekend moet worden, dan zou integratie met vaste stapgroote wel eens de voorkeur kunnen verdienen, omdat de afbreekfout dan een veel gladdere functie is van deze parameters. Bij een adaptieve methode immers kan de beslissing verfijnen of niet door een kleine parameterwijziging net anders uitvallen, hetgeen een sprong (discontinu¨ıteit) veroorzaakt in de afbreekfout.

Opgave 3.6: Schrijf een procedure voor de integratie van een gegeven functie op een gegeven interval met een voorgeschreven (geschatte) tolerantie met behulp van een adaptieve methode, gebaseerd op intervalhalvering en de regel van Simpson. (Schrijf eerst een recursieve versie; daarvoor is geen (expliciete) stapel nodig, omdat deze reeds aanwezig is in het recursiemechanisme.)

3.g Numerieke stabiliteit

Als we een functie f in ieder punt x van het integratie-interval kunnen berekenen met een absolute fout δf (x), waarvoor geldt

| δf (x) | ≤ d ∀ x ∈ [a, b] , (3.24)

dan vinden we voor de integraal over dat interval een (onvermijdelijke) absolute fout δI waarvoor we de volgende bovengrens vinden,

| δI | = | Z _b

a

(f (x) + δf (x)) dx − Z _b

a

f (x) dx | = | Z _b

a

δf (x) dx | ≤ d | b − a | . (3.25) Als we de integraal benaderen met een integratieformule (eventueel een samengestelde formule) met positieve gewichten, die exact is voor constanten,

Z b

a f (x) dx =

n

X

i=1

w_i f (t_i) + rest met w_i > 0 ∀ i en

n

X

i=1

w_i = b − a ,

dan geldt voor het verschil tussen berekende en exacte waarde (afgezien van afrondfouten bij het sommeren)

|

n

X

i=1

w_i δf (t_i) | ≤ d | b − a | . (3.26) Deze fout is van gelijke orde als de onvermijdelijke fout (3.25). Integratie met een formule met positieve gewichten is dus een numeriek stabiel proces.

Opgave3.7: In het integratieschema van Romberg bouwen we in feite een serie integratieformules T_h, T_h¹, T_h², · · · van steeds hogere orde op de steunpunten {a + jh | j = 0, 1, · · · , n} met h :=

(b − a)/n. Laat zien dat T_h¹ overeenkomt met de samengestelde Simpsonregel en dus positieve gewichten heeft. Laat vervolgens zien dat alle extrapolaties in dit schema formules zijn met positieve gewichten:

T_h^j =

n

X

i=1

w_i^j f (t_i) met w^j_i > 0 ,

voor alle i en j . (Gebruik volledige induktie naar j.) We zien hieruit dat het Romberg schema numeriek stabiel is.

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

De Bernoulli-polynomen B_i (i = 0, 1, · · · ) worden gedefinieerd door

B0 = 1 , (3.27a)

B_k^′ = k B_k−1, k ≥ 1 , (3.27b)

Z 1 0

B_k(x) dx = 0 , k ≥ 1 . (3.27c)

Ze voldoen aan de volgende eigenschappen:

B_k(1) − B_k(0) = Z ₁

0

k B_k−1(x) dx = 0 , k ≥ 2 , (3.28a) Bk(¹₂ − x) = (−1)^k Bk(¹₂+ x) , k ≥ 0 , (3.28b)

B_2k+1(1) = 0 = B_2k+1(0) , k ≥ 1 . (3.28c)

Deze laatste eigenschap (3.28c) kan natuurlijk niet gelden voor B₁, omdat dit polynoom van graad 1 dan minstens twee nulpunten zou hebben. Voor even k zijn de Bk in 0 en 1 gelijk en niet nul.

Deze getallen worden gedefinieerd als de Bernoulli-getallen B_k,

B_k := (−1)^k+1B_2k(0) , k = 1, 2, · · · . (3.29)

Een bewijs voor (3.28a) vinden we door (3.27b) te integreren en (3.27c) toe te passen. We bewijzen de eigenschappen (3.28b-c) met volledige induktie naar de graad van B_k. Uit de definitie volgt onmiddellijk

B₁(x) = x −¹₂, (3.30)

hetgeen voldoet aan (3.28b) (induktiebasis). Stel nu, dat (3.28b) waar is voor B_2k−1 (induktiehy- pothese), dan kunnen we dit polynoom integreren,

B_2k(x) = Z x

1 2

2k B_2k−1(x) dx + c_2k ,

met een (vooralsnog) onbekende integratieconstante c_2k. Omdat B_2k−1 antisymmetrisch is t.o.v.

het punt x = ¹₂, is B_2k symmetrisch is t.o.v. dit punt ongeacht de waarde van c_2k. Voor B_2k+1 geldt dan

B_2k+1(x) = Z _x

1 2

(2k + 1) B_2k(t) dt + c_2k+1= Z _x

1 2

Z _t

1 2

2k (2k + 1) B_2k−1(s) ds dt + c_2k(x −¹₂) + c_2k+1 Als we vervolgens eis (3.27c) toepassen en deze expressie integreren over het interval [0, 1], vinden we dat de integraal van de eerste twee termen nul is omdat de integranden antisymmetrisch zijn t.o.v. het punt x = ¹₂. De integratieconstante c_2k+1 moet dus nul zijn om de integraal over het geheel nul te maken. Dit bewijst (3.28b) voor B_2k en B_2k+1. Tenslotte is (3.28c) een gevolg van de antisymmetrie van B_2k+1en de gelijkheid van de functiewaarden in x = 0 en x = 1. Hiermee is (3.28) bewezen.

Met behulp van de polynomen en de getallen van Bernoulli kunnen we nu een asymptotische ontwikkeling geven voor de afbreekfout in de samengestelde trapeziumregel:

Stelling (Sommatieformule van Euler-McLaurin). Zij f ∈ C^2m+1([a, b]) en h := (b − a)/n met m en n geheel, dan is er een ξ ∈ [a, b] zodat

Z _b

a f (x) dx = h

2(f (a) + f (b)) + h

n−1

X

j=1

f (a + jh) +

m

X

j−1

(−1)^j

(2j)! B_jh^2j(f^(2j−1)(b) − f^(2j−1)(a)) + O(h^2m+1f^(2m+1)(ξ)).

(3.31)

Bewijs.Bekijk de integraal ^R f dx op een deelinterval van lengte h, zeg op [0, h], en integreer partieel met gebruik van (3.27):

Z _h

0

f (x) dx = Z _h

0

B₀(x

h) f (x) dx =

·

h B₁(x h) f (x)

¸h 0

− Z _h

0

h B₁(x

h) f^′(x) dx

= h

2(f (0) + f (h)) −

"

h² 2 B₂(x

h) f^′(x)

#h

0

+ Z h

0

h² 2 B₂(x

h) f^′′(x) dx . Als we doorgaan met het parti¨eel integreren van de integraal vinden we als algemene formule:

h^2k (2k)!

Z h 0

B_2k(x

h) f^(2k)(x) dx = h^2k+1 (2k + 1)!

·

B_2k+1(x

h) f^(2k)(x)

¸h 0

− h^2k+1 (2k + 1)!

Z h

0 B_2k+1(x

h) f^(2k+1)(x) dx , waarin de stokterm nul is wegens (3.28c). Nogmaals parti¨eel integreren geeft

= − h^2k+2 (2k + 2)!

·

B_2k+2(x

h) f^(2k+1)(x)

¸h 0

+ h^2k+2 (2k + 2)!

Z _h

0 B_2k+2(x

h) f^(2k+2)(x) dx

en door invullen van het Bernoulligetal wordt dit

= −B_k+1h^2k+2

(2k + 2)! (f^(2k+1)(h) − f^(2k+1)(0)) + h^2k+2 (2k + 2)!

Z h

0 B_2k+2(x

h) f^(2k+2)(x) dx . In het totaal vinden we dus de asymptotische reeks

Z _h

0

f (x) dx = h

2(f (0) + f (h)) +

m

X

k=1

(−1)^kB_kh^2k

(2k)! (f^(2k+1)(h) − f^(2k+1)(0)) + h^2m+1

(2m + 1)!

Z _h

0

B_2m+1(x

h) f^(2m+1)(x)dx .

Als we deze expressie uitschrijven voor ieder deelinterval (a + jh , a + jh + h), dan zien we dat alle afgeleiden op de interne deelpunten tegen elkaar wegvallen, zodat de gevraagde formule (3.31) overblijft.

Opmerking: Formule (3.31) maakt de samengestelde trapeziumregel bijzonder aantrekkelijk voor periodieke functies, immers de afbreekfout gaat dan harder naar nul dan iedere macht van h. Deze eigenschap hebben we in feite al gebruikt bij de discretisatie van de Fouriertransformatie, zie (2.4).

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.