Kwadratuur uitgelegd: numerieke methoden om integralen te benaderen

Er is een categorie bepaalde integralen die zich simpelweg niet laat oplossen met pen-en-papierrekenwerk. Niet omdat je iets verkeerd doet. Sommige functies, zoals e^(−x²) of sin(x)/x, hebben helemaal geen gesloten primitieve. In dat geval heb je een andere aanpak nodig: evalueer de functie op een handvol zorgvuldig gekozen punten en combineer die waarden om de onderliggende oppervlakte te schatten. Dat is kwadratuur.

Dit artikel legt uit wat kwadratuur is, hoe het algemene raamwerk werkt en de belangrijkste methoden, van de intuïtieve middenpuntsregel tot de verrassend scherpe Gaussiaanse aanpak. We behandelen ook samengestelde methoden, foutbronnen, adaptieve strategieën en hoe je de juiste methode kiest voor het probleem dat voor je ligt.

Wat is kwadratuur?

Het woord "kwadratuur" is oud. Het komt van het klassieke meetkundige probleem om een vierkant te construeren met dezelfde oppervlakte als een gegeven vorm, wat gewoon een andere manier is om te zeggen: "vind de oppervlakte." Tegenwoordig betekent de term de numerieke benadering van een bepaalde integraal.

Het kernidee is een vervanging: in plaats van een exacte primitieve te vinden, neem je steekproeven van de functie f(x) op een set punten, de zogeheten knooppunten (nodes), vermenigvuldig je elke waarde met een bijbehorende weging en tel je de resultaten op:

De termen kwadratuur en numerieke integratie betekenen in wezen hetzelfde. Kwadratuur duikt vaker op in wiskundige contexten; numerieke integratie is de alledaagse term in wetenschappelijke computing.

Hoe numerieke kwadratuur werkt

Elke kwadratuurregel werkt op een interval [a, b]. Binnen dat interval kiest hij een set knooppunten xᵢ en kent aan elk een gewicht wᵢ toe. De benadering is de gewogen som hierboven. Wat tussen methoden verschilt, is waar de knooppunten komen en hoe de gewichten worden toegewezen. Dat is alles. De trapeziumregel, de regel van Simpson en Gaussiaanse kwadratuur zijn allemaal variaties op hetzelfde sjabloon.

De benaderingsfout is het verschil tussen de gewogen som en de echte integraal. Elke methode introduceert er wat van, en hoe snel die fout krimpt naarmate je meer knooppunten toevoegt (of intervallen kleiner maakt), is wat een goede methode van een geweldige onderscheidt. Houd dit kader in gedachten: een kwadratuurmethode kiezen is eigenlijk kiezen hoe je je functiewaarderingen uitgeeft.

Veelgebruikte kwadratuurmethoden

De onderstaande methoden lopen van het meest intuïtief naar het meest verfijnd.

Visuele vergelijking van de middenpunts-, trapezium- en Simpson-integratieregels. Afbeelding door de auteur.

Middenpuntsregel

Neem het midden van [a, b], evalueer f daar en vermenigvuldig met de breedte (b − a). De aanname is dat de functiewaarde in het midden het gemiddelde over het hele interval vertegenwoordigt.

Het klinkt bijna té simpel, en voor één interval is dat het ook. Maar hier is de eerste verrassing: de middenpuntsregel is exact voor constante en lineaire functies, en de fout schaalt als O(h³) per interval, waarbij h = b − a. Dat is eigenlijk beter dan de trapeziumregel voor gladde functies, wat vaak tegen de intuïtie ingaat.

Trapeziumregel

Gebruik beide eindpunten in plaats van alleen het midden. Verbind ze met een rechte lijn en bereken de oppervlakte van het resulterende trapezium: (b − a) · [f(a) + f(b)] / 2. Exact voor lineaire functies, met O(h³)-fout per interval, dezelfde orde als de middenpuntsregel, zij het met een grotere constante.

De trapeziumregel is waarschijnlijk de meest onderwezen kwadratuurmethode, wat haar een onevenredige reputatie heeft gegeven. Voor goedgedragende functies is ze prima. Maar samengestelde Simpson doet bijna altijd beter voor hetzelfde aantal evaluaties.

Regel van Simpson

Pas een kwadratisch polynoom door drie punten: de twee eindpunten en het midden. Integreer dat polynoom exact:

De extra evaluatie in het midden kost één extra functiewaarde, maar de opbrengst is reëel: de regel van Simpson is exact voor polynomen tot en met graad drie, en de fout per interval schaalt als O(h⁵), twee ordes beter dan de trapeziumregel. Voor gladde functies is dit meestal de eerste methode om naar te grijpen.

Gaussiaanse kwadratuur

Hier wordt het echt interessant. De middenpunt-, trapezium- en Simpson-regels gebruiken allemaal gelijkmatig verdeelde knooppunten. Gaussiaanse kwadratuur stelt een andere vraag: wat als we zowel de locaties van de knooppunten als de gewichten zo kiezen dat de nauwkeurigheid maximaal is?

Met n knooppunten kan Gaussiaanse kwadratuur (specifiek de Gauss-Legendre-variant) elk polynoom tot en met graad 2n − 1 exact integreren. Twee knooppunten geven je exactheid voor derdegraads polynomen; vijf knooppunten dekken graad 9. De knooppunten zijn de nulpunten van Legendre-polynomen en clusteren naar de eindpunten van het interval in plaats van gelijkmatig verdeeld te staan. De extra nauwkeurigheid komt doordat de knooppunten mogen bewegen, niet door waar ze toevallig belanden. Hun clustering nabij de eindpunten is een bijproduct van de wiskunde, en hieronder zie je waarom.

Zo vergeleken de vier methoden met elkaar:

Samengestelde kwadratuurmethoden

Een regel één keer toepassen over het volledige interval [a, b] is meestal niet genoeg. De standaardoplossing: splits het interval in n gelijke subintervallen met breedte h = (b − a)/n en pas de basisregel op elk daarvan toe.

Samengestelde middenpunt- en trapeziumregels

Voor de samengestelde trapeziumregel evalueer je f op n+1 gelijkmatig verdeelde punten en som je n trapezia. De totale fout daalt naar O(h²). Halveer h, en de fout daalt met een factor vier. De samengestelde middenpuntsregel volgt dezelfde logica en heeft dezelfde globale foutorde.

Samengestelde regel van Simpson

De samengestelde versie vereist een even aantal subintervallen. Evalueer f op n+1 punten en krijg O(h⁴) globale fout. Halvering van h reduceert de fout met een factor zestien. Daarom is de samengestelde regel van Simpson vaak de standaard wanneer je middelhoge nauwkeurigheid op een gladde functie nodig hebt.

De onderstaande Python-code laat zien hoe de benadering strakker wordt naarmate je subintervallen toevoegt voor ∫[0 tot π] sin(x) dx (exacte waarde: 2.0):

import numpy as np

def composite_simpson(func, lower, upper, num_intervals):
    """Approximate integral using composite Simpson's rule."""
    if num_intervals % 2 != 0:
        raise ValueError("num_intervals must be even")
    nodes = np.linspace(lower, upper, num_intervals + 1)
    values = func(nodes)
    h = (upper - lower) / num_intervals
    weights = np.ones(num_intervals + 1)
    weights[1:-1:2] = 4   # odd-indexed interior points
    weights[2:-2:2] = 2   # even-indexed interior points
    return (h / 3) * np.dot(weights, values)

exact = 2.0
for n in [4, 8, 16, 32]:
    approx = composite_simpson(np.sin, 0, np.pi, n)
    print(f"n={n:3d}: approx={approx:.10f}, error={abs(approx - exact):.2e}")

Met slechts 32 intervallen daalt de fout onder 10⁻⁶. Dat is O(h⁴)-convergentie.

Gaussiaanse kwadratuur

Het idee achter Gaussiaanse kwadratuur verdient meer dan een tabelregel.

Waarom gelijkmatig verdeelde knooppunten niet altijd optimaal zijn

Wanneer knooppunten op gelijkmatig verdeelde punten staan, besteed je je evaluatiebudget op een manier die mogelijk niet past bij de moeilijkheid van de integrand. Sommige regio's hebben meer resolutie nodig; andere verspillen evaluaties aan vlakke stukken. Gaussiaanse kwadratuur omzeilt dit door knooppunten te kiezen die aantoonbaar optimaal zijn voor de polynomiale graad die je exact kunt integreren.

Voor n knooppunten kiest de Gauss-Legendre-regel de n nulpunten van het n-de Legendre-polynoom Pₙ(x) als knooppunten en berekent vervolgens gewichten zodat elk polynoom van graad ≤ 2n − 1 zonder fout wordt geïntegreerd. In de praktijk geeft numpy.polynomial.legendre.leggauss(n) beide terug. De knooppunten clusteren richting de eindpunten van [−1, 1], het standaardreferentie-interval, waar de interpolatiefout van polynomen meestal het grootst is en waar extra steekproefdichtheid het meest loont.

Wanneer gebruiken

Gaussiaanse kwadratuur is op zijn best wanneer je een gladde, analytische functie hebt en hoge nauwkeurigheid met weinig evaluaties nodig is. Vijf Gauss-Legendre-punten kunnen samengestelde Simpson met tientallen intervallen overtreffen op goedgedragende integranden. Maar de methode is minder geschikt voor functies met knikken, discontinuïteiten of data die alleen op vaste, gelijkmatig verdeelde punten beschikbaar zijn. In die gevallen werkt de zorgvuldig gekozen knooppuntplaatsing juist tegen je.

Kwadratuurfout en nauwkeurigheid

Geen enkele numerieke methode krijgt vrijstelling van fouten. Twee verschillende bronnen beïnvloeden elk kwadratuurresultaat, en begrip van beide is essentieel om te weten wanneer je je antwoord kunt vertrouwen.

Bronnen van fouten zijn onder meer:

• Discretisatiefout (ook truncatiefout genoemd) is het verschil tussen de exacte integraal en wat je kwadratuurregel oplevert bij perfecte rekenkunde. Dit is wat de O(hᵏ)-formules vastleggen. Het krimpt naarmate je meer punten toevoegt of de intervalbreedte verkleint.
• Floating-pointfout is een ander probleem. Het optellen van veel kleine hoeveelheden in eindige-precisierekenkunde introduceert afrondingen die niet met h krimpen en zelfs kunnen groeien. Bij samengestelde regels met duizenden subintervallen kun je een bodem bereiken waar floating-pointruis domineert en toevoegen van intervallen niet meer helpt. Voor de meeste wetenschappelijke toepassingen bereik je je gewenste tolerantie ruim vóór die bodem. Maar hij bestaat.

Functiesmoothness, intervalbreedte en methodeselectie bepalen samen de discretisatiefout. Een ruwe functie met scherpe kenmerken heeft meer intervallen of een adaptieve aanpak nodig.

Adaptieve kwadratuur

Een uniform rooster verspilt evaluaties wanneer de integrand op de meeste plaatsen glad is maar op een paar ingewikkeld.

Hoe het werkt

Pas een kwadratuurregel toe op een interval, pas vervolgens een hoger-orde regel toe (of dezelfde regel op twee helften) en vergelijk. Het verschil schat de lokale fout. Is die onder de tolerantie, accepteer dan het interval. Zo niet, splits het en herhaal recursief op elke helft. Ga door totdat elk subinterval aan de tolerantie voldoet.

De efficiëntiewinst ten opzichte van uniforme verfijning kan groot zijn wanneer de integrand geïsoleerde kenmerken heeft: een smalle piek, een knik, snelle oscillatie omringd door glad gebied. Adaptieve methoden concentreren evaluaties precies waar de functie ze nodig heeft.

scipy.integrate.quad()

In Python implementeert scipy.integrate.quad() adaptieve Gaussiaanse kwadratuur (specifiek QUADPACK's DQAGSE-algoritme). Geef je functie en grenzen door; je krijgt de schatting en een foutgrens terug. Voor de meeste eendimensionale integralen is dit het juiste gereedschap, zonder handmatig afstellen.

from scipy import integrate
import numpy as np

# Integrate a function with a peak near x=0.5
def peaked_func(x):
    return np.exp(-100 * (x - 0.5)**2)

result, error_bound = integrate.quad(peaked_func, 0, 1)
print(f"Estimate: {result:.10f}")
print(f"Error bound: {error_bound:.2e}")

quad() verwerkt deze smalle piek moeiteloos. Een uniforme Simpson-regel zou honderden subintervallen nodig hebben om dit te evenaren.

Hoe kies je een kwadratuurmethode

Begin met scipy.integrate.quad() voor eendimensionale integralen in Python, tenzij je een reden hebt om dat niet te doen. Het is adaptief, nauwkeurig en goed getest. Alles hieronder is relevant wanneer je zelf iets implementeert, werkt in een taal zonder goede adaptieve bibliotheek of optimaliseert voor een specifieke casus.

• Gladde, analytische functies met een bekende formule: Gaussiaanse kwadratuur geeft de meeste nauwkeurigheid per evaluatie. Heb je middelhoge nauwkeurigheid en eenvoud nodig, dan is de samengestelde regel van Simpson een solide standaard.
• Functies met knikken of discontinuïteiten: Gelijkmatig verdeelde methoden (trapezium, Simpson) presteren hier doorgaans beter dan Gaussiaanse kwadratuur, omdat vaste knooppuntplaatsing niet wordt verstoord door onregelmatigheden. Splits het integratie-interval bij de bekende discontinuïteit en behandel elk stuk apart.
• Data op vaste, gelijkmatig verdeelde punten: Gaussiaanse kwadratuur valt af. Je kunt de functie niet op willekeurige knooppunten evalueren, dus gebruik in plaats daarvan de samengestelde trapezium- of samengestelde Simpson-regel.
• Hoger-dimensionale integralen: Kwadratuur loopt al snel tegen de vloek van dimensionaliteit aan. Met n punten per dimensie en d dimensies heb je in totaal nᵈ evaluaties nodig. Tien dimensies met tien punten elk is 10¹⁰ evaluaties, wat niet gaat gebeuren. Monte Carlo schaalt veel beter boven ruwweg vijf dimensies, dus daar maak je de overstap.

Kwadratuur vs. verwante numerieke methoden

Kwadratuur staat niet op zichzelf. Het helpt om te zien waar het staat ten opzichte van andere tools die je waarschijnlijk overweegt.

Kwadratuur vs. analytische integratie

Analytische integratie vindt een exacte primitieve; kwadratuur benadert een getal. Als er een exacte primitieve bestaat en je kunt die vinden, gebruik die dan. Kwadratuur verdient zijn plek wanneer exacte integratie onmogelijk of onpraktisch is.

Kwadratuur vs. Monte Carlo-integratie

Beide benaderen bepaalde integralen, maar op verschillende manieren. Kwadratuur is deterministisch: het evalueert f op vaste, vooraf gekozen knooppunten, met convergentie gekoppeld aan functiesmoothness en het aantal knooppunten. Monte Carlo trekt willekeurige monsters en neemt het gemiddelde, met fout die schaalt als O(1/√N) ongeacht de dimensie. Dat is langzamere convergentie in lage dimensies, maar dimensie-onafhankelijk, wat het de praktische keuze maakt boven ruwweg vijf dimensies.

Kwadratuur vs. interpolatie

Sommige kwadratuurregels, waaronder die van Simpson, werken door een polynoom door evaluatiepunten te passen en dat exact te integreren. Maar het doel verschilt van interpolatie. Interpolatie reconstrueert de functie; kwadratuur wil alleen een getal. Het polynoom is een hulpmiddel, niet de uitkomst.

Conclusie

Wat ik bevredigend vind aan kwadratuur is hoeveel je haalt uit één organiserend idee. Elke methode hier, middenpuntsregel, trapezium, Simpson, Gaussiaans, is eigenlijk gewoon een ander antwoord op dezelfde vraag: gegeven een budget aan functiewaarderingen, waar plaats je ze en hoe weeg je ze?

Betere polynomen passen door evaluatiepunten drijft de fout snel omlaag. Verschuiven naar optimaal gekozen knooppunten perst nog meer nauwkeurigheid uit hetzelfde budget. Adaptieve methoden brengen vervolgens alles samen en sturen evaluaties naar precies daar waar de functie ze daadwerkelijk nodig heeft.

In Python is scipy.integrate.quad() de juiste standaard voor 1D-problemen. Al het andere in dit artikel is voor wanneer je wilt begrijpen wat er onder de motorkap gebeurt, of wanneer je zelf iets implementeert. Wil je je numerieke vaardigheden met SciPy verder opbouwen, volg dan onze cursus Introduction to Optimization in Python.