Numerička integracija

Kvadraturne formule

Problem kojim se bavimo je približno izračunavanje određenog integrala . Tom problemu možemo pristupiti tako što ćemo odrediti zgodnu funkciju koja aproksimira funkciju tako da važi . Prirodna je ideja da se za uzme interpolacioni polinom funkcije na . Ako je , tada je

gde su čvorovi interpolacije i .

Formule oblika

zovemo kvadraturnim formulama. Greška kvadraturne formule za funkciju je tada

Ukoliko je , tada kažemo da je kvadraturna formula tačna za funkciju .

Ukoliko je najveći broj takav da je kvadraturna formula tačna za svako , , , tada za broj kažemo da je algebarski stepen tačnosti te kvadraturne formule. U tom slučaju, kvadraturna formula je tačna za svaki polinom stepena .

U nastavku je data skripta za određivanje koeficijenata kvadraturne formule za zadate čvorove, kao i testiranje dobijene formule za zadatu funkciju .

% function fQuad = QuadFormula(x,a,b,f)

% fIntegral = integral(@(x)f(x),a,b); %ovo cemo smatrati tacnom vrednoscu integrala

% n = length(x);

% L = zeros(n,n);

% d = zeros(n,1);

% for i = 1:n

% L(i,:) = x.^(i-1);

% d(i) = integral(@(x)x.^(i-1),a,b);

% end

% A = L\d;

% fprintf('Koeficijenti kvadraturne formule:');

% disp(A);

% fQuad = sum((A').*f(x));

% fprintf('Greska kvadraturne formule:');

% disp(abs(fIntegral-fQuad));

% end

x = [0 0.5 1];

a = 0;

b = 1;

f = @(x)exp(-sin(x));

F = QuadFormula(x,a,b,f)

Koeficijenti kvadraturne formule:
   0.166666666666667
   0.666666666666666
   0.166666666666667
Greska kvadraturne formule:
     5.294665310706659e-05
F =
   0.651271965839420

Kvadraturne formule Njutn-Kotesovog tipa

Pretpostavimo da imamo ravnomernu podelu intervala sa čvorovima i korakom .

Podela.png

Ukoliko koristimo te čvorove da aproksimiramo funkciju na segmentima podele interpolacionim polinomom nultog, prvog i drugog stepena, dobijamo redom kvadraturne formula pravougaonika, trapeza i Simpsona.

Formula pravougaonika

Kao što smo već naveli, ovo je slučaj kada za aproksimaciju koristimo polinom nultog stepena, odnosno konstantu. Imamo dva slučaja:

Ocena greške formule pravougaonika je:

Geometrijski, određeni integral odgovara, do na znak, površini između grafika funkcije i ose. Metoda prevougaonika aproksimira tu površinu površinom odgovarajućeg pravougaonika. Slična će biti situacija i sa metodom trapeza.

FormulaPravougaonika.png

Navedimo kratku implementaciju metode pravougaonika. Za približno računanje integrala funkcije na intervalu se može koristiti ugrađena funkcija integral(f,a,b).

% function I = RectFormula(f,a,b,n)

% figure()

% fplot(f,[a,b],'-r');

% hold on

% grid on

% h = (b-a)/n;

% x = linspace(a,b,n+1);

% s = 0;

% for i = 1:n

% s = s+f(x(i));

% line([x(i) x(i)],[0 f(x(i))]);

% line([x(i) x(i+1)], [f(x(i)) f(x(i))]);

% end

% I = h*s;

% fIntegral = integral(f,a,b);

% fprintf('Greska formule pravougaonika za datu funkciju:');

% disp(abs(I-fIntegral));

% hold off

% end

a = 0;

b = 1;

n = 10;

f = @(x)exp(x.^2);

I = RectFormula(f,a,b,n)

Greska formule pravougaonika za datu funkciju:
   0.081391144591335
I =
   1.381260601315846

Formula trapeza

Kvadraturna formula glasi:

Ocena greške:

FormulaTrapeza.png

Implementacija:

% function I = TrapFormula(f,a,b,n)

% figure()

% fplot(f,[a,b],'-r');

% hold on

% grid on

% h = (b-a)/n;

% x = linspace(a,b,n+1);

% s = (f(x(1))+f(x(end)))/2;

% line([x(1) x(1)],[0 f(x(1))]);

% line([x(1) x(2)], [f(x(1)) f(x(2))]);

% for i = 2:n

% s = s+f(x(i));

% line([x(i) x(i)],[0 f(x(i))]);

% line([x(i) x(i+1)], [f(x(i)) f(x(i+1))]);

% end

% I = h*s;

% fIntegral = integral(f,a,b);

% fprintf('Greska formule trapeza za datu funkciju:');

% disp(abs(I-fIntegral));

% hold off

% end

a = 0;

b = 1;

n = 10;

f = @(x)exp(x.^2);

I = TrapFormula(f,a,b,n)

Greska formule trapeza za datu funkciju:
   0.004522946831617
I =
   1.467174692738799

Simpsonova formula

Kod Simpsonove formule se koriste po 3 susedna čvora podele, pa broj čvorova mora biti neparan, odnosno broj segmenata paran . Tada kvadraturna formula glasi:

Ocena greške:

Geometrijska interpretacija (3 tačke nam u opštem slučaju određuju parabolu):

FormulaSimpsona.png

Implementacija Simpsonove metode:

% function I = SimpFormula(f,a,b,n)

% h = (b-a)/n;

% x = linspace(a,b,n+1);

% F = f(x);

% s = F(1)+F(end)+4*sum(F(2:2:end-1))+2*sum(F(3:2:end-1));

% I = (h/3)*s;

% fIntegral = integral(f,a,b);

% fprintf('Greska formule Simpsona za datu funkciju:');

% disp(abs(I-fIntegral));

% end

a = 0;

b = 1;

n = 8;

f = @(x)exp(x.^2);

I = SimpFormula(f,a,b,n)

Greska formule Simpsona za datu funkciju:
     7.166876608710737e-05
I =
   1.462723414673269

Rungeova ocena greške

Za ocene greške prethodnih formula je neophodno računanje izvoda višeg reda i njihovo ocenjivanje, što je često komplikovano. Rungeova ocena greška je aproksimacionog karaktera i korisna je jer izbegava rad sa izvodima:

Ako je kvadraturna formula za funkciju , sa korakom podele intervala , a formula kad je korak , tada važi procena greške

gde je za formulu pravougaonika, za formulu trapeza i za Simpsonovu formulu. Da bismo sa korakom imali ekvidistantnu podelu na intervalu integracije, broj podeoka mora biti paran.

U nastavku je data implementacija metoda kojom određujemo Rungeovu ocenu greške. Primetimo da u slučaju Simpsonove metode broj podeoka mora biti deljiv sa 4 da bismo imali paran broj podeoka kad primenjujemo formulu sa korakom .

% function RungeError(f,a,b,n,p)

% switch p

% case 1

% E = RectFormula(f,a,b,n)-Rect(f,a,b,n/2);

% fprintf('Rungeova ocena greske za formulu pravougaonika: E = %f', E);

% case 2

% E = (TrapFormula(f,a,b,n)-TrapFormula(f,a,b,n/2))/3;

% fprintf('Rungeova ocena greske za formulu trapeza: E = %f', E);

% case 4

% E = (SimpFormula(f,a,b,n)-SimpFormula(f,a,b,4))/15;

% fprintf('Rungeova ocena greske za formulu Simpsona: E = %f', E);

% otherwise

% error('p mora biti iz skupa {1,2,4}!');

% end

% end

a = 0;

b = 1;

n = 8;

p = 4;

f = @(x)exp(x.^2);

RungeError(f,a,b,n,p);

Greska formule Simpsona za datu funkciju:
     7.166876608710737e-05
Greska formule Simpsona za datu funkciju:
   0.001059014538415
Rungeova ocena greske za formulu Simpsona: E = -0.000066