Kwadratury Gaussa - Legendre są to kwadratury postaci:
\(\displaystyle{ \int_{-1}^{1} g(z)dz = \sum_{i=1}^{n} w_{i}g(z_{i}) \ \ (1)}\)
gdzie wagi \(\displaystyle{ w_{j} }\) i węzły \(\displaystyle{ z_{j} }\) tak dobieramy, aby błąd przybliżenia
\(\displaystyle{ E_{n}(f) = \int_{-1}^{1} g(z)dz - \sum_{i=1}^{n} w_{i}g(z_{i}) }\) był równy zero.
Do obliczenia wartości węzłów i wag kwadratury korzystamy z równania:
\(\displaystyle{ E_{n}(a_{0}+a_{1}z +...+a_{m}z^{m}) = a_{0}E_{n}(1) + a_{1}E_{n}(z) +...+ a_{m}E_{n}(z^{m}). }\)
Z równania tego wynika, że \(\displaystyle{ E_{n}(f) = 0 }\) dla każdego wielomianu stopnia \(\displaystyle{ \leq m }\), wtedy i tylko wtedy, gdy
\(\displaystyle{ E_{n}(z^{i})=0 }\) dla \(\displaystyle{ i = 0,1,2, ..., m. }\)
Przypadek 1
\(\displaystyle{ n = 1 }\)
Otrzymujemy dwa parametry: \(\displaystyle{ w_{1}, z_{1}.}\)
Wymagając, aby \(\displaystyle{ E_{1} = 0, \ \ E_{n}(g)= 0 }\) otrzymujemy
\(\displaystyle{ \int_{-1}^{1} 1dz -w_{1}=0, \ \ \int_{-1}^{1}z dz - w_{1}z_{1}=0 }\)
Z tych równań wynika, że \(\displaystyle{ w_{1}=2 }\) i \(\displaystyle{ z_{1} =0 }\)
Na podstawie \(\displaystyle{ (1) \ \ \int_{-1}^{1} g(z)dz \approx 2g(0) }\) (reguła punktu środkowego)
Przypadek 2
\(\displaystyle{ n=2 }\)
Mamy cztery parametry: \(\displaystyle{ w_{1}, x_{1}, z_{1}, z_{2} }\) i cztery ograniczenia:
\(\displaystyle{ E_{n}(z^{i}) = \int_{-1}^{1} z^{i}dx - [w_{1}\cdot z^{i}_{1} + w_{2}\cdot z^{i}_{2}] = 0, \ \ i = 0,1,2,3. }\)
Otrzymujemy układ równań:
\(\displaystyle{ w_{1}+ w_{2}= 2,}\)
\(\displaystyle{ w_{1}\cdot z_{1} +w_{2}\cdot z_{2} = 0, }\)
\(\displaystyle{ w_{1}\cdot z^2_{1} + w_{2}\cdot z^2_{2} = \frac{2}{3} }\)
\(\displaystyle{ w_{1}\cdot z^3_{1} + w_{2}\cdot z^3_{2} = 0,}\)
z którego \(\displaystyle{ w_{1}= 1, \ \ w_{2} = 1, \ \ z_{1}=-\frac{\sqrt{3}}{3},\ \ z_{2}= \frac{\sqrt{3}}{3}. }\)
Podstawiając te wartości do \(\displaystyle{ (1) }\) otrzymujemy kwadraturę Gaussa rzędu trzeciego:
\(\displaystyle{ \int_{-1}^{1} g(x) dx \approx g\left(-\frac{\sqrt{3}}{3}\right) + g\left(\frac{\sqrt{3}}{3}\right). }\)
Przypadek 3
ogólnie dla \(\displaystyle{ n }\) mamy \(\displaystyle{ 2n }\) parametrów \(\displaystyle{ (z_{i}), \ \ (w_{i}) }\)
i ograniczeń:
\(\displaystyle{ E_{n}(z^{i}) = 0, \ \ i=0, 1,2, ..., 2n-1 }\)
\(\displaystyle{ \sum_{j=1}^{n} w_{j} z^{i}_{j} = \begin{cases} 0 \ \ \mbox{dla} \ \ i=1,3,...,2n-1 \\ \frac{2}{i+1} \ \ \mbox{dla} \ \ i= 0,2,,...,2n-2 \end{cases} }\)
Otrzymujemy kwadraturę Gaussa-Legendre rzędu \(\displaystyle{ 2n-1. }\)
Można udowodnić (*) że błąd przybliżenia kwadraturą Gaussa spełnia równość
\(\displaystyle{ E_{n} = \frac{g^{(2n)}(\zeta)}{2n!} \int_{a}^{b} p^2_{n}, }\) dla \(\displaystyle{ \zeta \in [a, b] }\)
(*) DAVID KINCAID WARD CHENEY analiza numeryczna. WNT Warszawa 2006.
Program w MATLAB
Kod: Zaznacz cały
function gauss_legendre(f,a,b,n)
% Obliczenie kwadratury Gaussa-Legendre rzędu piątego
fprintf('\n')
disp('Trzypunktowa kwadratura Gaussa-Legendre')
disp('_______________________________________________________')
disp(' i zi wi g(zi) wi*g(zi) ')
disp('_______________________________________________________ ')
if (n==2)
z(1)=-sqrt(1/3);
z(2)=-z(1);
w(1)=1; w(2)=1;
end
if (n==3)
z(1)=-sqrt(3/5);
z(2)=0;
z(3)=-z(1);
w(1)=5/9; w(2)= 8/9; w(3)=w(1);
end
if (n==4)
z(1)=-sqrt(1/7*(3-4*sqrt(0.3)));
z(2)=-sqrt(1/7*(3+4*sqrt(0.3)));
z(3)=-z(1);
z(4)=-z(2);
w(1)=1/2+1/12*sqrt(10/3);
w(2)=1/2-1/12*sqrt(10/3);
w(3)=w(1); w(4)=w(2);
end
if (n==5)
z(1)=-sqrt(1/9*(5-2*sqrt(10/7)));
z(2)=-sqrt(1/9*(5+2*sqrt(10/7)));
z(3)=0;
z(4)=-z(1);
z(5)=-z(2);
w(1)=0.3*((-0.7+5*sqrt(0.7))/(-2+5*sqrt(0.7)));
w(2)=0.3*((0.7+5*sqrt(0.7))/(2+5*sqrt(0.7)));
w(3)=128/225;
w(4)=w(1); w(5)=w(2);
end;
S=0;
for i=1:n
x=((b-a)*z(i)+a+b)/2;
g=feval(f,x);
S=S+w(i)*g;
fprintf('%2.0f %12.4f %12.4f% 12.4f %12.4f\n',i-1,z(i),w(i),g,g*w(i));
end;
INT=S*(b-a)/2;
fprintf('\n Wartość całki funkcji f(x) jest =%16.8f\n',INT);
Kod: Zaznacz cały
function f=f1(x)
f=7+14*x.^6;
>> gauss_legendre('f1',0,1,3)
Trzypunktowa kwadratura Gaussa-Legendre
________________________________________________________
i zi wi g(zi) wi*g(zi)
________________________________________________________
0 -0.7746 0.5556 7.0000 3.8889
1 0.0000 0.8889 7.2188 6.4167
2 0.7746 0.5556 13.8320 7.6844
Wartość całki funkcji f(x) jest = 8.99500000
Kod: Zaznacz cały
Czteropunktowa kwadratura Gaussa-Legendre
________________________________________________________
i zi wi g(zi) wi*g(zi)
________________________________________________________
0 -0.3400 0.6521 7.0181 4.5768
1 -0.8611 0.3479 7.0000 2.4350
2 0.3400 0.6521 8.2663 5.3908
3 0.8611 0.3479 16.0911 5.5974
Wartość całki funkcji f(x) jest = 9.00000000
Programy MATLAB i OCTAVE mają dwie wbudowane funkcje "quad" i "quadl" do obliczeń numerycznych kwadratur, aproksymujących całkę z określoną tolerancją:
Kod: Zaznacz cały
>>quad(fun,a,b,tol,)
Kod: Zaznacz cały
function f=f1(x)
f=7+14*x.^6;
>> quad('f1',0,1,10^(-3))
ans =
9.0000