De Euler-Maclaurin-formule
Stel je hebt de
integraal
van een willekeurige functie f (x):
Deze
integraal kan ik ook als volgt schrijven:
Hierin is p
0 (x) een polynoom als functie van x, en het spreekt voor zich dat p
0 (x) in
bovenstaande vergelijking gelijk is aan één.
Vervolgens ga ik gebruik maken van
partieel integreren:
Voor onze
integraal, zie vergelijking (2),
ziet dat er zo uit:
Er geldt dus:
Hierin is c
1 een
integratieconstante.
Het uitvoeren van het
partieel integreren levert aldus het volgende op:
We gaan gewoon nog een keer
partieel integreren:
Nu geldt er dus:
Het uitvoeren van de tweede keer
partieel integreren levert het volgende op:
Je voelde het waarschijnlijk al aankomen, we gaan nog een keer
partieel integreren:
Ditmaal geldt het volgende:
De derde keer
partieel integreren brengt ons dit resultaat:
Ik breng even een kleine verandering aan in de notatie waarbij ik de functie zelf als de nulde
afgeleide aangeef:
Het
partieel integreren
kunnen we uiteraard eindeloos herhalen:
En dit kan ik natuurlijk ook compact opschrijven als volgt:
Waarbij voor de polynomen geldt:
Ik stel:
Hiermee verandert vergelijking (15) in:
Door de vergelijkingen (16) en (17) te combineren ontstaat:
Tot zo ver heb ik alleen nog maar een heleboel
afgeleiden, polynomen en een zooi onbekende
integratieconstanten gegenereerd.
Laten we maar eens op weg gaan om hier iets zinvols van te maken.
Om te beginnen eis ik:
Dat levert het volgende op:
Daarmee worden de polynomen q
n (x):
De functiewaarden van de polynomen qn (x) voor x = 0 zijn de Bernoulli-getallen:
Merk op dat de functiewaarden van de polynomen q
n (x) voor x = 1 hieraan gelijk zijn (behalve voor n = 1,
dat scheelt een minteken).
Dat q
n (0) en q
n (1) gelijk zijn aan elkaar komt uiteraard door de
integraal-eis van vergelijking (20).
Dit zijn de q
n (1) waarden:
Het Bernoulli-getal B
1 valt om twee redenen uit de toon.
Ten eerste omdat q
1 (0) en q
1 (1) een minteken respectievelijk een plusteken hebben,
maar ook omdat B
1 het enige oneven Bernoulli-getal is ongelijk aan nul.
Alle overige oneven Bernoulli-getallen zijn nul.
Uit het voorgaande, met name vergelijking (19l), is duidelijk geworden:
Met behulp van vergelijking (21k) kan ik ook schrijven:
Hiermee kan ik rechtstreeks de Bernoulli-getallen bepalen:
En zo kunnen we eindeloos doorgaan:
n |
Bn (als breuk) |
Bn (decimaal) |
0 |
1 |
1.00000 |
1 |
−1/2 |
−0.50000 |
2 |
1/6 |
0.16667 |
3 |
0 |
0.00000 |
4 |
−1/30 |
−0.03333 |
5 |
0 |
0.00000 |
6 |
1/42 |
0.02381 |
7 |
0 |
0.00000 |
8 |
−1/30 |
−0.03333 |
9 |
0 |
0.00000 |
10 |
5/66 |
0.07576 |
11 |
0 |
0.00000 |
12 |
−691/2730 |
−0.25311 |
13 |
0 |
0.00000 |
14 |
7/6 |
1.16667 |
15 |
0 |
0.00000 |
16 |
−3617/510 |
−7.09216 |
17 |
0 |
0.00000 |
18 |
43867/798 |
54.97118 |
19 |
0 |
0.00000 |
20 |
−174611/330 |
−529.12424 |
... |
... |
Tabel 1: Bernoulli-getallen (voor meer Bernoulli-getallen zie deze pagina) |
We begonnen dit verhaal bij:
En dat bracht ons bij:
Oftewel, door de vergelijkingen (1) en (18) te combineren:
Ik ga nu de
integraal
grenzen meegeven:
Ik ga nog twee eisen toevoegen:
- ten eerste moeten de grenzen
a en b gehele getallen (integers) zijn,
- ten tweede moeten de polynomen qn (x) periodiek zijn met periode één.
Met andere woorden, de polynomen zijn nog wel een functie van x maar om precies te zijn zijn ze nu een
functie van (x − int (x)) (de int-functie geeft het hoogste gehele getal, de grootste integer,
die kleiner of gelijk aan x is).
Het zal zo duidelijk worden waarom, maar ik neem de term k = 0 even apart:
Ik ga nu bij de rechterterm de
grenzen invullen:
De tweede term aan de rechterkant neem ik apart onder de loep en ik ga wat met de
grenzen
spelen (en alle halfjes die er bij in komen zijn de Bernoulli-getallen B
1):
En dan hebben we hier de
trapeziummethode
staan voor het numeriek benaderen van een
integraal.
Ik ga (32) nog iets anders opschrijven door er een half f (b) van af te trekken en er ook bij op te tellen:
Vergelijking (33) stop ik in vergelijking (31):
Ik ga de termen even wat reorganiseren en wat verduidelijking erbij schrijven:
Ervan uitgaande dat de restterm nul wordt indien n naar oneindig gaat (let op: dit is absoluut niet
vanzelfsprekend!) krijgen we:
En omdat de oneven Bernoulli-getallen nul zijn brengt ons dit uiteindelijk bij de Euler-Maclaurin-formule
voor het benaderen van integralen die niet
rechtstreeks op te lossen zijn:
Dit is in feite een verbeterde versie van de
trapeziummethode.
De rechterterm met de
afgeleiden
compenseert alle stukjes die we met de
trapeziummethode
te veel of te weinig in rekening brengen.
Echter, we hebben onderweg een beperking aangebracht door te stellen dat a en b gehele getallen zijn en daar
zouden we eigenlijk wel weer vanaf willen (want we houden niet van beperkingen).
De
trapeziummethode is niet gebonden aan gehele getallen
en dus is er geen enkele reden om nu wel vast te houden aan gehele getallen.
Ik noem de stapgrootte s en ik schrijf het deel
trapeziummethode
weer op volgens vergelijking (32):
Maar nu zit ik nog met het deel met de
afgeleiden.
Dat los ik als volgt op, ik stel:
Hierin is x een geheel getal, s de stapgrootte en t een variabele die niet meer een geheel getal behoeft te zijn.
Dit leidt middels de
kettingregel tot het volgende:
Dus voor iedere keer
differentiëren
komt er een s bij en er komt, net als bij de
trapeziummethode,
nog een s voor.
Hiermee verandert vergelijking (38) in:
Hoe pakt dit nou uit in de praktijk?
Laten we daarvoor de functie f (x) = 1/x eens gaan
integreren van x = 1 tot x = 2 met stapjes van 0.1.
Voor het deel
trapeziummethode
levert dat het volgende op:
1/2 f (1.0) |
0.50000000000000000000 |
f (1.1) |
0.90909090909090909091 |
f (1.2) |
0.83333333333333333333 |
f (1.3) |
0.76923076923076923077 |
f (1.4) |
0.71428571428571428571 |
f (1.5) |
0.66666666666666666667 |
f (1.6) |
0.62500000000000000000 |
f (1.7) |
0.58823529411764705882 |
f (1.8) |
0.55555555555555555556 |
f (1.9) |
0.52631578947368421053 |
1/2 f (2.0) |
0.25000000000000000000 |
Totaal |
6.93771403175427943230 |
Tabel 2: trapeziummethode f (x) = 1/x |
Wanneer we dit totaal met de stapgrootte 0.1 vermenigvuldigen hebben we 0.69377140317542794323.
Het exacte antwoord van de
integraal is:
Het trapeziumdeel zit hier dus 0.090 procent naast (boven).
Nu gaan we de schoonheid van de wiskunde in actie zien.
De
afgeleide
van 1/x is −1/x
2 en daarmee wordt de term k = 1 van het deel met
afgeleiden:
Dit trekken we van het trapeziumdeel af en dan krijgen we 0.69314640317542794323, een afwijking van −0.00011 procent.
De derde
afgeleide
van 1/x is −6/x
4 en daarmee wordt de term k = 2 van het deel met
afgeleiden:
Dit trekken we af van het vorige tussenresultaat en dan krijgen we 0.69314718442542794323, een afwijking van nog maar
0.00000056 procent.
Door slechts twee termen met
afgeleiden
erin te betrekken is het antwoord met meer dan een factor honderdduizend nauwkeuriger geworden en hebben we
de
natuurlijke logaritme
van twee verkregen tot op acht decimalen nauwkeurig!
Wat ook spectaculair is, is dat de Euler-Maclaurin-formule een verband legt tussen
integralen en somreeksen.
Je kunt ermee
integralen
oplossen middels somreeksen, maar ook somreeksen oplossen middels
integralen!
Daarom schrijf ik vergelijking (37) even iets anders op:
Volgens de overlevering gebruikte Euler bovenstaande vergelijking om het
Bazel-probleem op te lossen, een probleem
dat de wiskundigen destijds al heel lang bezig hield:
Om dit uit te rekenen met behulp van alleen pen en papier is een onmogelijke opgave omdat deze reeks heel langzaam
convergeert (men wist toen zelfs nog niet of de reeks überhaupt convergeert) zoals onderstaande tabel laat zien:
Aantal termen |
Resultaat |
Aantal cijfers goed (afgerond) |
1 |
1.00000000000000000000 |
0 |
10 |
1.54976773116654069035 |
1 |
100 |
1.63498390018489286508 |
3 |
1000 |
1.64393456668155980314 |
4 |
10000 |
1.64483407184805976981 |
4 |
100000 |
1.64492406689822626981 |
5 |
1000000 |
1.64493306684872643631 |
6 |
∞ |
1.64493406684822630823 |
∞ |
Tabel 3: somreeks 1/x2 |
Zoals je ziet gaat de nauwkeurigheid van de uitkomst (het aantal goede cijfers) ongeveer met de
logaritme van het aantal termen dus dat schiet niet op.
Met behulp van de Euler-Maclaurin-formule is dit probleem uiteraard een fluitje van een cent, maar toch zit hier nog
een addertje onder het gras.
Ik zal eerst de
afgeleiden
opschrijven van 1/x
2:
Dit stop ik in vergelijking (45) en ik vul gelijk de
grenzen in:
Dit ziet er heel simpel uit en dat is het ook, maar zoals tabel 1 laat zien worden de Bernoulli-getallen steeds groter
en gaan uiteindelijk naar +∞ of −∞.
Met andere woorden: deze reeks divergeert!
De oplossing vereist enige creativiteit.
We gaan de somreeks opdelen in twee delen:
De overblijvende somreeks, hierboven de rechterterm, heeft nu als ondergrens tien in plaats van één en daarmee wordt
(48) wél convergerend:
Dat gaan we maar eens uitwerken:
1.644767731166540690350214159738xx |
(1/6) ∙ 10−3 = |
0.000166666666666666666666666667xx |
-------------------------------------------------- + |
1.644934397833207357016880826405xx |
(1/30) ∙ 10−5 = |
0.000000333333333333333333333333xx |
-------------------------------------------------- − |
1.644934064499874023683547493071xx |
(1/42) ∙ 10−7 = |
0.000000002380952380952380952381xx |
-------------------------------------------------- + |
1.644934066880826404635928445452xx |
(1/30) ∙ 10−9 = |
0.000000000033333333333333333333xx |
-------------------------------------------------- − |
1.644934066847493071302595112119xx |
(5/66) ∙ 10−11 = |
0.000000000000757575757575757576xx |
-------------------------------------------------- + |
1.644934066848250647060170869695xx |
(691/2730) ∙ 10−13 = |
0.000000000000025311355311355311xx |
-------------------------------------------------- − |
1.644934066848225335704859514383xx |
(7/6) ∙ 10−15 = |
0.000000000000001166666666666667xx |
-------------------------------------------------- + |
1.644934066848226502371526181050xx |
(3617/510) ∙ 10−17 = |
0.000000000000000070921568627451xx |
-------------------------------------------------- − |
1.644934066848226431449957553599xx |
Tabel 4: termen met
afgeleiden van 1/x2 |
Hiermee hebben we het antwoord, namelijk π
2/6 = 1.644934066848226308227702251860, tot op vijftien
cijfers achter de komma benaderd (dus zestien cijfers goed)!
Wanneer we meer termen gaan toevoegen dan raken we toch weer in de problemen omdat de Bernoulli-getallen sterk gaan toenemen.
Daarom is dit ook geen echt bewijs dat de uitkomst gelijk is aan π
2/6, maar het stuurde Euler wel
in de goede richting (hij kwam een aantal jaren later alsnog met een waterdicht bewijs).
Het belangrijkste hier is dat het benaderen van deze reeks tot op meer dan vijftien decimalen nauwkeurig, door alle
termen uit te rekenen, handmatig een onmogelijke opgave is en dat zelfs een computer daar wel enige tijd mee bezig is
(er zouden vele
biljarden termen opgeteld moeten worden).
Echter, via bovenstaande methode kost het slechts een paar uurtjes werk met enkel en alleen pen en papier!