Kwadratury Gaussa-Legendre

Sigma-ciała i zbiory borelowskie. Miary, miary zewnętrze i miara Lebesgue'a. Funkcje mierzalne. Całka Lebesgue'a. Inne zagadnienia analizy rzeczywistej.
janusz47
Użytkownik
Użytkownik
Posty: 7910
Rejestracja: 18 mar 2009, o 16:24
Płeć: Mężczyzna
Podziękował: 30 razy
Pomógł: 1670 razy

Kwadratury Gaussa-Legendre

Post autor: janusz47 »

Kwadratury Gaussa-Legendre

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); 
Przykład 1

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
Przykład 2

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
 
Dodano po 12 godzinach 30 minutach 24 sekundach:
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,) 
 
Przykład

Kod: Zaznacz cały

function f=f1(x)
f=7+14*x.^6;

>> quad('f1',0,1,10^(-3))

ans =

    9.0000
 
ODPOWIEDZ