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 p0 (x) een polynoom als functie van x, en het spreekt voor zich dat p0 (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 c1 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 qn (x):











Bernoulli

De functiewaarden van de polynomen qn (x) voor x = 0 zijn de Bernoulli-getallen:











Merk op dat de functiewaarden van de polynomen qn (x) voor x = 1 hieraan gelijk zijn (behalve voor n = 1, dat scheelt een minteken). Dat qn (0) en qn (1) gelijk zijn aan elkaar komt uiteraard door de integraal-eis van vergelijking (20). Dit zijn de qn (1) waarden:










Het Bernoulli-getal B1 valt om twee redenen uit de toon. Ten eerste omdat q1 (0) en q1 (1) een minteken respectievelijk een plusteken hebben, maar ook omdat B1 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: 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 B1):
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:

Euler

Maclaurin

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/x2 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/x4 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/x2:
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!