Simpsons regel

Från Wikipedia
Hoppa till: navigering, sök

Simpsons regel, efter Thomas Simpson, används för att approximera en integral.

 \int_{a}^{b} f(x) \, dx \approx \frac{b-a}{6}\left[f(a) + 4f\left(\frac{a+b}{2}\right)+f(b)\right].


Förklaring[redigera | redigera wikitext]

En funktion f(x) och dess interpolationspolynom P(x) som har samma funktionsvärde som f(x) i a, b och m

Simpsons regel är ett sätt att uppskatta en integral till en funktion f(x) genom att ersätta detta med ett andragradspolynom som antar samma värde som f(x) i ändpunkterna a, b och mittpunkten m. För att kunna simpsons regel så måste alltså värdet på f(x) vara känt i dessa tre punkter.

Det andragradspolynom som uppfyller dessa krav kallar vi P(x) och detta kan vi få fram med bland annat Newtons eller Lagranges interpolationspolynom där den senare ser ut som följer:

P(x)= f(a)\frac{(x-m)(x-b)}{(a-m)(a-b)} +f(m)\frac{(x-a)(x-b)}{(m-a)(m-b)} +f(b)\frac{(x-a)(x-m)}{(b-a)(b-m)}

Där alltså m är

m = \frac {a + b}{2}

och

 \int_{a}^{b} f(x) \, dx \approx \int_{a}^{b} P(x) \, dx

Ur detta ses att om funtionen f(x) är ett andragradspolynom så kommer P(x) vara exakt lika med f(x).

Härledning[redigera | redigera wikitext]

Härledning med hjälp av Lagranges interpolationspolynom[redigera | redigera wikitext]

För att härleda att regeln ser ut som den gör behöver man bara integrera det uttrycket som står ovan. Om man nu sätter att

h = \frac {b - a}{2},

dvs. den längd det mellan a och m eller m och b kan man skriva om uttrycket på det något trevligare utseendet:

P(x)= f(a)\frac{(x-(a + h))(x-(a + 2h))}{2h^{2}} +f(a + h)\frac{(x-a)(x-(a + h))}{-h^{2}} +f(a + 2h)\frac{(x-a)(x-(a + h))}{2h^{2}}

Om man nu integrerar P(x) med avseende på x så kommer man få en enkel integralberäkning som dock är väldigt lång och därför kommer inte hela presenteras. Om man bara tittar på kvoten efter f(a) och integrerar denna så kan man uppfattning om att det ändå stämmer, så att:

P_{1}(x) = f(a)\frac{(x-(a + h))(x-(a + 2h))}{2h^{2}}

Och integralen av P_{1}(x) blir då alltså

\int_{a}^{a+2h} P_{1}(x) \, dx = \int_{a}^{a+2h} f(a)\frac{(x-(a + h))(x-(a + 2h))}{2h^{2}} \, dx 
=
  \frac {f(a)}{2h^{2}} \left[ x (a^2 + 3ah + 2h^2)- x^2 \left( \frac {2a + 3h}{2} \right) + \frac {x^3}{3} \right]_{a}^{a+2h} =
 f(a) \frac {h}{3} = f(a) \frac {b-a}6

Vilket stämmer med ovan. De andra två räknas ut på likartat sätt och resultatet ovan kommer att erhållas. Nämnaren 3 gör att uttrycket ibland kallas Simpsons 1/3 regel.

Restterm[redigera | redigera wikitext]

Integralen av interpolationspolynomet kan även skrivas med en restterm R_{n} och får då utseendet

\int_{a}^{a+2h} f(x) \, dx = \int_{a}^{a+2h} P(x) \, dx + R_{n}

Där

R_{n} = \frac {h^5}{90} \ f^{(4)}(\xi),\ \xi \in \left[a, a+2h\right] [1]

Det största möjliga felet som kan inträffa gör alltså det när absolutbeloppet av fjärdederivatan av f är så stort som möjligt. Det största felet kan skrivas som

\left| f(\xi)^{(4)} \right| \le M, \ \xi \in \left[a, a+2h\right]
\left|R_{n}\right| \le \frac {h^5}{90} M

Simpsons 3/8 regel[redigera | redigera wikitext]

Om man istället vet värdet hos en funktion f(x) i fyra skilda punkter kan man göra som ovan men ersätta funktionen med ett tredjegradspolynom. Härledningen ser väldigt likartad ut som den ovan. Lagranges interpolationspolynom blir då

 P(x) = f(a)\frac {(x-(a+h))(x-(a+2h))(x-(a+3h))}{-6h^3} + f(a+h)\frac {(x-a)(x-(a+2h))(x-(a+3h))}{2h^3}
 + f(a+2h)\frac {(x-a)(x-(a+h))(x-(a+3h))}{-2h^3} + f(a+3h)\frac {(x-a)(x-(a+h))(x-(a+2h))}{6h^3}

Där

h = \frac{b-a}{3}

Simpsons 3/8 regel säger då att

 \int_{a}^{b} f(x) \, dx \approx \int_{a}^{b} P(x) \, dx

Och att detta i sin tur är

 \int_{a}^{a+2h} f(x) \, dx \approx \frac {3h}{8} \left[ f(a) + 3f(a+h) + 3f(a+2h) + f(a+3h) \right]

MATLAB-implementation[redigera | redigera wikitext]

Simpsons regel kan implementeras i MATLAB enligt följande:

function [P] = simpson(f,a,b,n)
%f=funktionens namn, a=startvärde, b=slutvärde, 
%n=antal iterationer
h=(b-a)/n;
S=f(a);
i=1:2:n-1;
x=a+h.*i;
y=f(x);
S=S+4*sum(y);
i=2:2:n-2;
x=a+h.*i;
y=f(x);
S=S+2*sum(y);
S=S+f(b);
P=h*S/3;
end

Källor[redigera | redigera wikitext]

  1. ^ http://planetmath.org/?op=getobj&from=objects&id=7979 Planet Math: bound on error of Simpson's rule