Tak jest na przykład , gdy chcemy obliczyć całkę \(\displaystyle{ \int_{0}^{2} e^{x^2} dx.}\)
Musimy odwoływać się do metod przybliżonych. Jeśli mamy znaleźć przybliżoną wartość całki oznaczonej
\(\displaystyle{ \int_{a}^{b} f(x)dx, }\)
to sposób postępowania polega na zastąpieniu funkcji \(\displaystyle{ f }\) inną bliską funkcją \(\displaystyle{ g(x),}\) dla której całkę możemy obliczyć.
Stosujemy więc przybliżenie:
\(\displaystyle{ \int_{a}^{b} f(x)dx \approx \int_{a}^{b} g(x) dx. }\)
Dobrym przybliżeniem (bo łatwo całkowalnym) jest na przykład suma częściowa szeregu Taylora.
Przykład:
\(\displaystyle{ \int_{0}^{1}e^{x^2}dx \approx \int_{0}^{1}\sum_{i=1}^{n}\frac{1}{i!}x^{2i}dx = \sum_{i=0}^{n}\frac{1}{(2i+1)\cdot i!}.}\)
Kwadraturę Simpsona konstruujemy, dzieląc przedział całkowania \(\displaystyle{ [a, b] }\) na parzystą liczbę \(\displaystyle{ n }\) podprzedziałów tej samej długości \(\displaystyle{ h.}\)
Pole zawarte między wykresem funkcji \(\displaystyle{ f(x) }\) i osią \(\displaystyle{ Ox }\) przybliżamy polem, pomiędzy wykresem wielomianu
kwadratowego (paraboli) i osią \(\displaystyle{ Ox }\) w punktach (węzłach): \(\displaystyle{ x_{i-1}, x_{i}, x_{i+1} }\) przedziału \(\displaystyle{ [ a,b]. }\)
Z podstawowego twierdzenia rachunku całkowego wynika, że
\(\displaystyle{ I_{i} = \int_{x_{i-1}}^{x_{i}} f(x) dx = F(x_{i+1}) - F(x_{i}) \ \ (1) }\)
Rozwijając funkcję \(\displaystyle{ F(x) }\) w szereg Taylora mamy
\(\displaystyle{ F(x_{i+1}) = F(x_{i}+h)= F(x_{i}) +hF^{'}(x{i}) + \frac{h^2}{2!}F^{''}(x_{i}) + \frac{h^3}{3!}F^{'''}(x_{i}) + \frac{h^4}{4!}F^{(4)}(x_{i})+ \frac{h^5}{5!}f^{(5)}(x_{i})+ \ \ ...}\)
\(\displaystyle{ = F(x_{i}) + hf(x_{i}) + \frac{h^2}{2!}f'(x_{i}) + \frac{h^3}{3!}f^{''}(x_{i})+ \frac{h^4}{4!}f^{'''}(x_{i}) +\frac{h^5}{5!}f^{(4)}(x_{i})+ \ \ ... \ \ (2) }\)
Podobnie
\(\displaystyle{ F(x_{i-1}) = F(x_{i}) - hf(x_{i}) + \frac{h^2}{2!}f'(x_{i}) - \frac{h^3}{3!}f^{''}(x_{i})+ \frac{h^4}{4!}f^{'''}(x_{i}) - \ \ ... \ \ (3) }\)
Podstawiając \(\displaystyle{ (2), (3) }\) do \(\displaystyle{ (1) }\) otrzymujemy
\(\displaystyle{ I_{i} = 2h f(x_{i}) +\frac{h^3}{3!} f^{''}(x_{i}) +\frac{h^5}{90}f(x_{i}) + \ \ ... (3) }\)
Przybliżenie różnicowe drugiej pochodnej centralnej
\(\displaystyle{ f^{''}(x_{i}) = \frac{f(x_{i+1}) -2f(x_{i})+f(x_{i-1})}{h^2} - \frac{h^2}{12}f^{(4)}(x_{i})+ \ \ ... }\)
w \(\displaystyle{ (3) }\), daje ostatecznie wzór na kwadraturę Simpsona zwaną też kwadraturą parabol:
\(\displaystyle{ I_{i} = \frac{h}{3} \left[ f(x_{i+1}) + 4 f(x_{i}) +f(x_{i-1}) \right] - \frac{h^{5}}{90}f^{(4)}(x_{i}) + \ \ ... }\)
Możemy przenieść ten wzór na przedział \(\displaystyle{ [x_{0}, x_{2}] }\), otrzymując:
\(\displaystyle{ \int_{x_{0}}^{x_{2}} f(x)dx = \frac{h}{3}\left[ f(x_{0}) + 4f(x_{1}) +f(x_{2}) \right ] - \frac{h^5}{90}f^{(4)}(\xi) \ \ (4) }\)
gdzie
\(\displaystyle{ x_{0}< \xi < x_{2}.}\)
Błąd tej kwadratury \(\displaystyle{ E_{i} = -\frac{h^5}{90} f^{(4)}(x_{i}),}\) jeśli \(\displaystyle{ |f^{(4)}(x_{i})|\leq M, }\) możemy zapisać w postaci:
\(\displaystyle{ |E_{i}|\leq \frac{h^5}{90}\cdot M \ \ (5) }\)
Przybliżenie całki \(\displaystyle{ (4) }\) jest niedokładne, zwłaszcza gdy mamy doczynienia z dużymi przedziałami całkowania.
Dlatego przedział \(\displaystyle{ [a, b] }\) dzielimy na \(\displaystyle{ n = 2m, \ \ m\in \ZZ }\) podprzedziałów i sumując
\(\displaystyle{ I = \sum_{i=1}^{m} \int_{x_{2i-2}}^{x_{2i}} f(x) dx = \frac{h}{3}\sum_{i=1}^{m} \left[ f(x_{2i-2}) +4f(x_{2i-1})+f(x_{2i})\right] - \frac{h^5}{90}\sum_{i=1}^{m}f^{(4)}(\xi_{i}),}\)
otrzymujemy równanie złożonej kwadratury Simpsona, które redukuje się do równania:
\(\displaystyle{ I = \frac{h}{3}\left[ f(x_{0})+2\sum_{i=1}^{m-1}f(x_{2i})+ 4 \sum_{i=1}^{m} f(x_{2i-1}) +f(x_{2m}) \right] - \frac{h^5}{90}\sum_{i=1}^{m}f^{(4)}(\xi).}\)
Jeśli błąd kwadratury \(\displaystyle{ (5) }\) zastosujemy do \(\displaystyle{ n }\) podprzedziałów przedziału \(\displaystyle{ [a,b], }\) to błąd przybliżenia
złożoną kwadraturą Simpsona możemy napisać w postaci nierówności:
\(\displaystyle{ |E_{T}| \leq \frac{n}{2}\cdot \frac{h^2}{90}\cdot M = (b-a)\cdot \frac{h^4}{180}\cdot M.}\)
Program w OCTAVE złożonej kwadratury Simpsona
Kod: Zaznacz cały
function simpson(f,a,b,n)
h = (b-a)/n;
disp('___________________________________________________')
disp([' i xi f(xi) h=',num2str(h)])
disp('___________________________________________________')
S=feval(f,a);
fprintf('%2.0f%12.4f%14.6f\n',0,a,S);
for i=1:n/2
m=2*i-1;
x=a+h*m;
g=feval(f,x);
S=S+4*g;
fprintf('%2.0f%12.4f % 14.6f\n',m,x,g);
m=2*i;
x=a+h*m;
g=feval(f,x);
if(i==n/2)
S=S+g;
else
S=S+2*g;
end;
fprintf('%2.0f%12.4f%14.6f\n',m,x,g);
end
INT=h*S/3;
fprintf('\n Całka funkcji f(x) =% 16.8f\n',INT);
function f = f1(x)
f= -2*x*exp(0.7*x);
Wywołanie programu
>> simpson('f1',2,7,4)
________________________________
i xi f(xi) h = 1.25
________________________________
0 2.0000 -16.220800
1 3.2500 -63.231474
2 4.5000 -210.024581
3 5.7500 -643.773551
4 7.0000 -1880.056916
Całka funkcji f(x) = -2143.47790675