Integratiemethode: 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):





















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
21 0 0.00000
22 854513/138 6192.12319
23 0 0.00000
24 −236364091/2730 −86580.25311
25 0 0.00000
26 8553103/6 1425517.16667
27 0 0.00000
28 −23749461029/870 −27298231.06782
29 0 0.00000
30 8615841276005/14322 601580873.90064
31 0 0.00000
32 −7709321041217/510 −15116315767.09216
33 0 0.00000
34 2577687858367/6 429614643061.16667
35 0 0.00000
36 −26315271553053477373/1919190 −13711655205088.33277
37 0 0.00000
38 2929993913841559/6 488332318973593.16667
39 0 0.00000
40 −261082718496449122051/13530 −19296579341940068.14863
41 0.00000 ∙ 100
42 8.41693 ∙ 1017
43 0.00000 ∙ 100
44 −4.03381 ∙ 1019
45 0.00000 ∙ 100
46 2.11507 ∙ 1021
47 0.00000 ∙ 100
48 −1.20866 ∙ 1023
49 0.00000 ∙ 100
50 7.50087 ∙ 1024
51 0.00000 ∙ 100
52 −5.03878 ∙ 1026
53 0.00000 ∙ 100
54 3.65288 ∙ 1028
55 0.00000 ∙ 100
56 −2.84988 ∙ 1030
57 0.00000 ∙ 100
58 2.38654 ∙ 1032
59 0.00000 ∙ 100
60 −2.13999 ∙ 1034
61 0.00000 ∙ 100
62 2.05010 ∙ 1036
63 0.00000 ∙ 100
64 −2.09380 ∙ 1038
65 0.00000 ∙ 100
66 2.27527 ∙ 1040
67 0.00000 ∙ 100
68 −2.62577 ∙ 1042
69 0.00000 ∙ 100
70 3.21251 ∙ 1044
71 0.00000 ∙ 100
72 −4.15983 ∙ 1046
73 0.00000 ∙ 100
74 5.69207 ∙ 1048
75 0.00000 ∙ 100
76 −8.21836 ∙ 1050
77 0.00000 ∙ 100
78 1.25029 ∙ 1053
79 0.00000 ∙ 100
80 −2.00156 ∙ 1055
81 0.00000 ∙ 100
82 3.36750 ∙ 1057
83 0.00000 ∙ 100
84 −5.94710 ∙ 1059
85 0.00000 ∙ 100
86 1.10119 ∙ 1062
87 0.00000 ∙ 100
88 −2.13553 ∙ 1064
89 0.00000 ∙ 100
90 4.33289 ∙ 1066
91 0.00000 ∙ 100
92 −9.18855 ∙ 1068
93 0.00000 ∙ 100
94 2.03469 ∙ 1071
95 0.00000 ∙ 100
96 −4.70038 ∙ 1073
97 0.00000 ∙ 100
98 1.13180 ∙ 1076
99 0.00000 ∙ 100
100 −2.83822 ∙ 1078
... ...
Tabel 1: Bernoulli-getallen
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:
  1. ten eerste moeten de grenzen a en b gehele getallen (integers) zijn,
  2. 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 B1):



En dan hebben we hier de trapezium-methode 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:



Er van 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 trapezium-methode. De rechterterm met de afgeleiden compenseert alle stukjes die we met de trapezium-methode 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 trapezium-methode 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 trapezium-methode 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 trapezium-methode, 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 trapezium-methode levert dat het volgende op (alles afgerond op veertien decimalen):
1/2 f (1.0) 0.50000000000000
f (1.1) 0.90909090909091
f (1.2) 0.83333333333333
f (1.3) 0.76923076923077
f (1.4) 0.71428571428571
f (1.5) 0.66666666666667
f (1.6) 0.62500000000000
f (1.7) 0.58823529411765
f (1.8) 0.55555555555556
f (1.9) 0.52631578947368
1/2 f (2.0) 0.25000000000000
Totaal 6.93771403175428
Tabel 2: trapezium-methode f (x) = 1/x
Wanneer we dit totaal met de stapgrootte 0.1 vermenigvuldigen hebben we 0.69377140317543. Het exacte antwoord van de integraal is:



Het trapeziumdeel zit hier dus 0.09 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.69314640317543, een afwijking van −0.0001 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.69314718442543, een afwijking van nog maar 0.0000006 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 meer dan tien 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.00000000000000 0
10 1.54976773116654 1
100 1.63498390018489 3
1000 1.64393456668156 4
10000 1.64483407184806 4
100000 1.64492406689824 5
1000000 1.64493306684877 6
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.6447677311665406903502141597xx
(1/6) ∙ 10−3 = 0.0001666666666666666666666667xx
-------------------------------------------------- +
1.6449343978332073670168808264xx
(1/30) ∙ 10−5 = 0.0000003333333333333333333333xx
-------------------------------------------------- −
1.6449340644998740336835474931xx
(1/42) ∙ 10−7 = 0.0000000023809523809523809524xx
-------------------------------------------------- +
1.6449340668808264146359284455xx
(1/30) ∙ 10−9 = 0.0000000000333333333333333333xx
-------------------------------------------------- −
1.6449340668474930813025951122xx
(5/66) ∙ 10−11 = 0.0000000000007575757575757576xx
-------------------------------------------------- +
1.6449340668482506570601708698xx
(691/2730) ∙ 10−13 = 0.0000000000000253113553113553xx
-------------------------------------------------- −
1.6449340668482253457048595145xx
(7/6) ∙ 10−15 = 0.0000000000000011666666666667xx
-------------------------------------------------- +
1.6449340668482265123715261812xx
(3617/510) ∙ 10−17 = 0.0000000000000000709215686275xx
-------------------------------------------------- −
1.6449340668482264414499575537xx
Tabel 4: termen met afgeleiden van 1/x2
Hiermee hebben we het antwoord, namelijk π2/6 = 1.644934066848226436472415166646, tot op zestien cijfers achter de komma benaderd (dus zeventien 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!