Ben's Hobbyhoekje -- Numerieke wiskunde XVIII: Integreren met de Simpsonregel

Numerieke wiskunde XVIII:
Integreren met de Simpsonregel

Stel: je wilt de oppervlakte berekenen onder de grafiek van een functie. Zie de figuur hiernaast. De oppervlakte die we zoeken is gearceerd.
Het domein (dat is de waarde van x) loopt van x₀ tot x₁.
Bij een wat meer ingewikkeld functievoorschrift is analytisch integreren meestal lastig (lees: moeilijk). Numeriek Integreren is dan een mogelijke oplossing.

Op deze pagina gaat het over de Regel van Simpson. Met deze techniek kan in een beperkt aantal integratiestappen een goede nauwkeurigheid worden behaald.

Het domein wordt verdeeld in een eindig aantal stapjes. Uit de functiewaarden aan het begin en eind van elke stap, en de stapgrootte h, wordt het oppervlak van een vierhoek berekend. Één zijde van die vierhoek volgt de grafiek via een driepunts interpolatiepolynoom (een parabool). De som van alle 'plakjes oppervlak' is een benadering van de totale oppervlakte onder de grafiek.
Hierbij geldt: Hoe kleiner de stap h, hoe nauwkeuriger de oplossing.
De oppervlakte van één plakje is gelijk aan:
a = (y_i + 4·y_i+h/2 + y_i+1) · h / 6.
Het totale oppervlak is dan: A = ∑ (a) = ∑ (y_i + 4·y_i+h/2 + y_i+1) · h / 6, waarbij i loopt van 0 tot en met n-1 en h = (x₁ - x₀) / n.
De stapgrootte h hoeft tijdens het rekenproces niet constant te zijn. Het is in principe mogelijk om h tijdens het rekenproces te vergroten of te verkleinen, afhankelijk van het verloop van de functie.
Een eenvoudiger aanpak is door de berekening meerdere maken uit te voeren, met kleinere, constante, stapgrootte en de opeenvolgende resultaten te vergelijken. De resultaten zullen steeds minder verschillen. Als de verschillen klein genoeg zijn, heb je een behoorlijke benadering van de oppervlakte. Maar het is nog steeds een benadering!

Een uitgewerkt voorbeeld
Vraag: Wat is het oppervlak onder de grafiek van y = sin(x) tussen x₀ = 0° en x₁ = 180°?

Uitwerking: In de meeste (alle?) computertalen rekent de ingebouwde sinus-functie in radialen. Het is daarom het handigste om de integraal uit te rekenen van x₀ = 0 tot x₁ = π

Als eerste tekenen we de grafiek van y = sin(x) tussen en π, zie de figuur hiernaast.
Om het voorbeeld wat overzichtelijk te houden lossen we dit op 10 stappen, dus n = 10. y = sin(x) wordt dus berekend voor x = 0; π/10; 2π/10; 3π/10; .... π;/, in totaal dus 11 keer. Zie onderstaande tabel.
De grafiek van sin(x) verloopt glad over het hele domein. De figuur is gebaseerd op de tabel hieronder, bij de hele waarden van i, waarbij de elf punten zijn verbonden door rechte lijnen.

i	x	y	a
0	`0.00000`	`0.00000`	`0.04894365`
0+h/2	`0.15708`	`0.15643`
1	`0.31416`	`0.30902`	`0.142040004`
1+h/2	`0.47124`	`0.45399`
2	`0.62832`	`0.58779`	`0.221232493`
2+h/2	`0.78540`	`0.70711`
3	`0.94248`	`0.80902`	`0.278769204`
3+h/2	`1.09956`	`0.89101`
4	`1.25664`	`0.95106`	`0.309018043`
4+h/2	`1.41372`	`0.98769`
5	`1.57080`	`1.00000`	`0.309018043`
5+h/2	`1.72788`	`0.98769`
6	`1.88496`	`0.95106`	`0.278769204`
6+h/2	`2.04204`	`0.89101`
7	`2.19911`	`0.80902`	`0.221232493`
7+h/2	`2.35619`	`0.70711`
8	`2.51327`	`0.58779`	`0.142040004`
8+h/2	`2.67035`	`0.45339`
9	`2.82743`	`0.30902`	`0.04894365`
9+h/2	`2.98451`	`0.15643`
10	`3.14159`	`0.00000`
	`A = ∑ (a) =`		`2.000006784`

De exacte oplossing van ∫ sin(x)dx tussen 0 en π is gelijk aan 2. De afwijking is een gevolg van de afwijking in het bepalen van de driepunts-polynomen, die in dit geval klein is.
Natuurlijk is er ook een applicatie om zelf met deze techniek aan de slag te gaan. Klik op de knop hiernaast.
De applicatie berekent de stapgrootte uit de integratiegrenzen en het aantal stappen. Dat is constant gedurende de berekening.
De applicatie toont alle beschikbare decimalen in de uitkomst. Er wordt dus niet afgerond!
De code van de applicatie kun je downloaden om zelf aan door te ontwikkelen.
Als je verder wilt werken aan de applicatie, download je de .zip-file en pak je hem uit. Je hebt dan meteen een werkend voorbeeld.

Downloaden:

Druk op de knop: File: voorb668.zip, 3217 bytes.