Kvadratur förklarad: numeriska metoder för att approximera integraler

Det finns en kategori av bestämda integraler som helt enkelt inte låter sig lösas med papper-och-penna-kalkyl. Inte för att du gör något fel. Vissa funktioner, som e^(−x²) eller sin(x)/x, har över huvud taget inga primitiva funktioner i sluten form. När det händer behöver du ett annat angreppssätt: utvärdera funktionen i ett antal noggrant valda punkter och kombinera dessa värden för att uppskatta arean under kurvan. Det är kvadratur.

Den här artikeln går igenom vad kvadratur är, hur den allmänna ramen fungerar och de viktigaste metoderna, från den intuitiva mittpunktsregeln till den förvånansvärt skarpa Gaussiska metoden. Vi tar också upp sammansatta metoder, felkällor, adaptiva strategier och hur du väljer rätt metod för problemet framför dig.

Vad är kvadratur?

Ordet "kvadratur" är gammalt. Det kommer från det klassiska geometriska problemet att konstruera en kvadrat med samma area som en given figur, vilket bara är ett annat sätt att säga "hitta arean". I dag betyder termen den numeriska approximationen av en bestämd integral.

Kärnidén är en substitution: i stället för att hitta en exakt primitiv funktion, samplar du funktionen f(x) i en uppsättning punkter som kallas noder, multiplicerar varje värde med en tillhörande vikt och summerar resultaten:

Termerna kvadratur och numerisk integrering betyder i huvudsak samma sak. Kvadratur dyker oftare upp i mer matematiska sammanhang; numerisk integrering är vardagstermer inom vetenskaplig beräkning.

Hur numerisk kvadratur fungerar

Varje kvadraturregel verkar på ett intervall [a, b]. Inom det intervallet väljer den en uppsättning noder xᵢ och tilldelar var och en en vikt wᵢ. Approximationen är den viktade summan ovan. Det som skiljer metoderna åt är var noderna placeras och hur vikterna bestäms. Det är allt. Trapetsregeln, Simpsons regel och Gaussisk kvadratur är alla varianter på samma mall.

Approximationsfelet är gapet mellan den viktade summan och den sanna integralen. Varje metod introducerar ett sådant fel, och hur snabbt det minskar när du lägger till fler noder (eller gör intervallen mindre) är det som skiljer en bra metod från en utmärkt. Ha denna utgångspunkt i åtanke: att välja en kvadraturmetod är i grunden att välja hur du använder dina funktionsutvärderingar.

Vanliga kvadraturmetoder

Metoderna nedan går från mest intuitiva till mest sofistikerade.

Visuell jämförelse av mittpunkts-, trapets- och Simpsons integreringsregler. Bild av författaren.

Mittpunktsregeln

Ta mittpunkten av [a, b], utvärdera f där och multiplicera med bredden (b − a). Antagandet är att funktionsvärdet i mitten representerar medelvärdet över hela intervallet.

Det låter nästan för enkelt, och för ett enskilt intervall är det det. Men här kommer den första överraskningen: mittpunktsregeln är exakt för konstanta och linjära funktioner, och dess fel skalar som O(h³) per intervall där h = b − a. Det är faktiskt bättre än trapetsregeln för släta funktioner, vilket brukar överraska folk.

Trapetsregeln

Använd båda ändpunkterna i stället för bara mitten. Förbind dem med en rät linje och beräkna arean av den resulterande trapetsen: (b − a) · [f(a) + f(b)] / 2. Exakt för linjära funktioner, med O(h³) fel per intervall, samma ordning som mittpunktsregeln men med en större konstant.

Trapetsregeln är förmodligen den mest undervisade kvadraturmetoden, vilket har gett den ett oproportionerligt rykte. För välbeteende funktioner är den hyfsad. Men sammansatt Simpsons regel gör nästan alltid ett bättre jobb för samma antal utvärderingar.

Simpsons regel

Anpassa ett andragradspolynom genom tre punkter: de två ändpunkterna och mittpunkten. Integrera det polynomet exakt:

Den extra utvärderingen i mitten kostar ett ytterligare funktionsanrop, men utdelningen är verklig: Simpsons regel är exakt för polynom upp till grad tre, och felet per intervall skalar som O(h⁵), två ordningar bättre än trapetsregeln. För släta funktioner är detta oftast den första metoden att ta till.

Gaussisk kvadratur

Här blir det på riktigt intressant. Mittpunkts-, trapets- och Simpsons regler använder jämnt fördelade noder. Gaussisk kvadratur ställer en annan fråga: vad händer om vi väljer både nodernas positioner och vikterna för att maximera noggrannheten?

Med n noder kan Gaussisk kvadratur (specifikt Gauss–Legendre-varianten) integrera vilket polynom som helst av grad upp till 2n − 1 exakt. Två noder ger exakthet för kubiker; fem noder hanterar polynom av grad 9. Noderna är rötterna till Legendres polynom, och de klustrar mot intervallets ändpunkter i stället för att ligga jämnt fördelade. Den extra noggrannheten kommer av att låta noderna flytta sig, inte av var de råkar hamna. Att de klustrar nära ändpunkterna är en biprodukt av matematiken, och vi ska se varför nedan.

Så här jämförs alla fyra metoderna:

Sammansatta kvadraturmetoder

Att tillämpa en regel en gång över hela intervallet [a, b] räcker vanligtvis inte. Standardlösningen: dela upp intervallet i n lika stora delintervall med bredd h = (b − a)/n och tillämpa grundregeln på vart och ett.

Sammansatta mittpunkts- och trapetsregler

För den sammansatta trapetsregeln utvärderar du f i n+1 jämnt fördelade punkter och summerar n trapetser. Totalt fel sjunker till O(h²). Halvera h, och felet minskar med en faktor fyra. Den sammansatta mittpunktsregeln följer samma logik och har samma globala felordning.

Sammansatt Simpsons regel

Den sammansatta varianten kräver ett jämnt antal delintervall. Utvärdera f i n+1 punkter och få O(h⁴) globalt fel. Att halvera h reducerar felet med en faktor sexton. Det är därför sammansatt Simpsons regel ofta är standardvalet när du behöver måttlig noggrannhet för en slät funktion.

Pythonkoden nedan visar hur approximationen skärps när du lägger till delintervall för ∫[0 till π] sin(x) dx (exakt värde: 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}")

Med bara 32 intervall sjunker felet under 10⁻⁶. Det är O(h⁴)-konvergens.

Gaussisk kvadratur

Idén bakom Gaussisk kvadratur förtjänar mer än en tabellrad.

Varför jämnt fördelade noder inte alltid är optimala

När noderna ligger jämnt fördelade lägger du utvärderingsbudgeten på ett sätt som kanske inte matchar integrandens svårighetsgrad. Vissa områden behöver högre upplösning; andra slösar utvärderingar på flacka partier. Gaussisk kvadratur undviker detta genom att välja noder som är bevisat optimala för den polynomgrad du kan integrera exakt.

För n noder väljer Gauss–Legendre-regeln de n rötterna till det n:te Legendres polynom Pₙ(x) som noder och beräknar sedan vikter så att varje polynom av grad ≤ 2n − 1 integreras utan fel. I praktiken returnerar numpy.polynomial.legendre.leggauss(n) båda. Noderna klustrar mot ändpunkterna av [−1, 1], standardintervallet, där interpolationsfelet för polynom tenderar att vara störst och där extra sampeltäthet ger mest utdelning.

När du ska använda den

Gaussisk kvadratur är som bäst när du har en slät, analytisk funktion och behöver hög noggrannhet med få utvärderingar. Fem Gauss–Legendre-punkter kan överträffa sammansatt Simpsons regel med dussintals intervall på välbeteende integrander. Men den passar dåligt för funktioner med knyckar, diskontinuiteter eller data som bara finns i fasta, jämnt fördelade punkter. I sådana fall motverkar den noggrant valda nodplaceringen ditt syfte.

Kvadraturfel och noggrannhet

Ingen numerisk metod slipper undan fel. Två separata källor påverkar varje kvadraturresultat, och att förstå båda är avgörande för att veta när du kan lita på ditt svar.

Felkällor inkluderar:

• Diskretiseringsfel (även kallat trunkeringsfel) är gapet mellan den exakta integralen och vad din kvadraturregel ger i perfekt aritmetik. Det är detta som O(hᵏ)-formlerna fångar. Det krymper när du lägger till fler punkter eller minskar intervallbredden.
• Flyttalsfel är ett annat problem. Att summera många små kvantiteter i ändlig precision introducerar avrundningar som inte krymper med h och som faktiskt kan växa. Med sammansatta regler som använder tusentals delintervall kan du nå en nivå där flyttalsbruset dominerar och fler intervall inte hjälper. För de flesta vetenskapliga beräkningar når du din önskade tolerans långt innan den nivån. Men den finns.

Funktionens släthet, intervallbredd och metodval påverkar alla diskretiseringsfelet. En grov funktion med skarpa drag kräver fler intervall eller ett adaptivt angreppssätt.

Adaptiv kvadratur

Ett enhetligt rutnät slösar utvärderingar när integranden är slät på de flesta ställen men komplicerad på några få.

Så fungerar det

Använd en kvadraturregel på ett intervall, använd sedan en högre ordningens regel (eller samma regel på två halvor) och jämför. Skillnaden uppskattar det lokala felet. Om det är under toleransen, acceptera intervallet. Annars, dela det och rekursivt behandla vardera halvan. Upprepa tills varje delintervall uppfyller toleransen.

Effektivitetsvinsten jämfört med enhetlig förfining kan vara stor när integranden har isolerade särdrag: en smal topp, en knyck, snabba svängningar omgivna av släta partier. Adaptiva metoder koncentrerar utvärderingar exakt där funktionen kräver dem.

scipy.integrate.quad()

I Python implementerar scipy.integrate.quad() adaptiv Gaussisk kvadratur (specifikt QUADPACK:s DQAGSE-algoritm). Skicka in din funktion och gränser; få tillbaka uppskattningen och en felgräns. För de flesta endimensionella integraler är detta rätt verktyg, utan manuell trimning.

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() hanterar denna smala topp utan problem. En enhetlig Simpsons regel skulle behöva hundratals delintervall för att matcha den.

Hur du väljer en kvadraturmetod

Börja med scipy.integrate.quad() för endimensionella integraler i Python om du inte har en anledning att låta bli. Den är adaptiv, noggrann och vältestad. Allt nedan är relevant när du implementerar något själv, arbetar i ett språk utan ett bra adaptivt bibliotek eller optimerar för ett specifikt fall.

• Släta, analytiska funktioner med känd formel: Gaussisk kvadratur ger störst noggrannhet per utvärdering. Om du behöver måttlig noggrannhet och enkelhet är sammansatt Simpsons regel ett stabilt standardval.
• Funktioner med knyckar eller diskontinuiteter: Jämnt fördelade metoder (trapets, Simpson) tenderar att överträffa Gaussisk kvadratur här, eftersom fasta nodpositioner inte störs av oregelbundenheter. Dela upp integrationsintervallet vid känd diskontinuitet och hantera varje del separat.
• Data i fasta, jämnt fördelade punkter: Gaussisk kvadratur är utesluten. Du kan inte utvärdera funktionen i godtyckliga noder, så använd i stället sammansatt trapets- eller sammansatt Simpsons regel.
• Högre-dimensionella integraler: Kvadratur drabbas snabbt av dimensionsförbannelsen. Med n punkter per dimension och d dimensioner behöver du nᵈ totala utvärderingar. Tio dimensioner med tio punkter vardera är 10¹⁰ utvärderingar, vilket inte är realistiskt. Monte Carlo skalar mycket bättre över ungefär fem dimensioner, så där sker bytet.

Kvadratur vs. relaterade numeriska metoder

Kvadratur existerar inte i ett vakuum. Det hjälper att se var det placerar sig i förhållande till andra verktyg du sannolikt överväger.

Kvadratur vs. analytisk integrering

Analytisk integrering hittar en exakt primitiv; kvadratur approximerar ett tal. Om en exakt primitiv finns och du kan hitta den, använd den. Kvadratur förtjänar sin plats när exakt integrering är omöjlig eller opraktisk.

Kvadratur vs. Monte Carlo-integrering

Båda approximerar bestämda integraler, men på olika sätt. Kvadratur är deterministiskt: det utvärderar f i fasta, förvalda noder, med konvergens kopplad till funktionens släthet och antalet noder. Monte Carlo drar slumpmässiga stickprov och medelvärdesbildar dem, med fel som skalar som O(1/√N) oavsett dimension. Det är långsammare konvergens i låga dimensioner, men dimensionsoberoende, vilket gör det till det praktiska valet över ungefär fem dimensioner.

Kvadratur vs. interpolation

Vissa kvadraturregler, inklusive Simpsons regel, arbetar genom att anpassa ett polynom genom utvärderingspunkter och integrera det exakt. Men målet skiljer sig från interpolation. Interpolation återskapar funktionen; kvadratur vill bara ha ett tal. Polynomet är ett verktyg, inte resultatet.

Slutsats

Det jag tycker är tillfredsställande med kvadratur är hur mycket du får ut av en enda organiserande idé. Varje metod här, mittpunktsregel, trapets, Simpson, Gaussisk, är egentligen bara olika svar på samma fråga: givet en budget av funktionsutvärderingar, var placerar du dem och hur viktar du dem?

Att anpassa bättre polynom genom utvärderingspunkter driver snabbt ner felet. Att flytta till optimalt valda noder pressar ut ännu mer noggrannhet ur samma budget. Adaptiva metoder knyter sedan ihop allt och dirigerar utvärderingar dit funktionen faktiskt behöver dem.

I Python är scipy.integrate.quad() rätt standardval för 1D-problem. Allt annat i den här artikeln är till för när du behöver förstå vad som händer under huven eller när du implementerar något själv. Om du vill fortsätta bygga numeriska färdigheter med SciPy, gå vår kurs Introduction to Optimization in Python.