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)
Kvadraturne formule Njutn-Kotesovog tipa
Pretpostavimo da imamo ravnomernu podelu intervala sa čvorovima i korakom .
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:
- formula levih pravougaonika
- formula desnih pravougaonika
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.
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)
Formula trapeza
Kvadraturna formula glasi:
Ocena greške:
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)
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):
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)
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);