Integracija

10. zadatak

Neka je funkcija $f$ zadata eksplicitno komadnim M-fajlom $funkcija.m$.

$funkcija.m$

f = inline('x.^2+1');
%f=@(x) x.^2+1;

• Napisati M-fajl $trapez.m$ sa funkcijom $I = trapez(a, b)$ koja približno računa i vraća vrednost određenog integrala funkcije $f$ (granice integracije su $a$ i $b$) korišćenjem uopštene Trapezne kvadraturne formule sa $n = 9$ čvorova.

  ​x
  function I = trapez(a, b)
  ​
  % Korišćenjem Trapezne kvadraturne formule računamo približne vrednosti integrala
  % od a do b funkcije f
  % Poziv funkcije: I = trapez(1, 5)
  ​
  funkcija;
  n = 9;
  % čvorovi kvadraturne formule
  X = linspace(a, b, n);
  % korak
  h = (b - a) / (n - 1); % ili h = X(2) - X(1)
  ​
  Y = f(X);
  %vrednost integrala
  I = (h / 2) * (Y(1) + Y(n) + 2 * sum(Y(2:n-1)));
  

• Napisati M-fajl $simps.m$ sa funkcijom $I = simps(a, b)$ koja približno računa i vraća vrednost određenog integrala funkcije $f$ (granice integracije su $a$ i $b$) korišćenjem uopštene Simpsonove kvadraturne formule sa $n = 9$ čvorova.

  xxxxxxxxxx
  function I = simps(a, b)
  ​
  % Korišćenjem Simpsonove kvadraturne formule računamo vrednosti integrala
  % od a do b funkcije f
  % Poziv funkcije: I = simps(1, 5)
  ​
  funkcija;
  n = 9;
  ​
  % čvorovi kvadraturne formule
  X = linspace(a, b, n);
  ​
  % korak
  h = (b - a) / (n - 1); % ili h = X(2) - X(1)
  Y = f(X);
  ​
  % vrednost integrala
  I = (h / 3) * (Y(1) + Y(n) + 2 * sum(Y(3:2:n-1)) + 4 * sum(Y(2:2:n-1)));
  % U Simpsonovoj formuli čvorove sa neparnim indeksom množimo sa 4 (formula sa vežbi) 
  % Pošto indeksiranje vektora u MATLAB-u kreće od 1 (a ne od 0), ovde ćemo sa 4 množiti
  % elemente vektora sa parnim indeksima.
  

• Napisati M-fajl $vredfunk.m$ sa funkcijom $[X, Y] = vredfunk(k, p)$ koja približno izračunava vrednost funkcije $I(x) = \int_{1}^{x}f(t)dt,$ kada se $x$ kreće od $2$ do $k \in \mathbb{N}, k\geq2$ sa korakom $1$. Ukoliko je $p = 1$ integrale računati koristeći uopštenu Simpsonovu kvadraturnu formulu (sa $n = 9$ čvorova), a u slučaju kada je $p = 2$ uopštenu Trapeznu kvadraturnu formulu (sa $n = 9$ čvorova). Funkcija vraća dva vektora: $X$ sa vrednostima $x_i$ i $Y$ sa izračunatim vrednostima funkcije $I(x)$ u tačkama $x_i$.

  xxxxxxxxxx
  function [X, I] = vredfunk(k, p)
  ​
  % Iz postavke zadatka, vektor X uzima vrednosti od 2 do k, sa korakom 1, tj.
  % X = 2:1:k (a možemo ga formirati i u okviru for petlje kako je urađeno u nastavku)
  % Vektor I sadrži vrednosti integrala od 1 do X(i) funkcije f
  % Pozivi funkcija: [X, I]=vredfunk(5,1), [X, I]=vredfunk(5,2)
  ​
  for i = 1 : k-1
      % X=[2 3 4 5 .... k]
      X(i) = i + 1;
      if p == 1
          % Svaki element vektora I predstavlja približnu vrednost integrala dobijenu
          % Simpsonovom kvadraturnom formulom
          I(i) = simps(1, X(i));
      else
          % Svaki element vektora I predstavlja približnu vrednost integrala dobijenu
          % Trapeznom kvadraturnom formulom
          I(i) = trapez(1, X(i));
      end
  end
  

11. zadatak

Neka je funkcija $f$ zadata eksplicitno komadnim M-fajlom $funkcija.m$.

$funkcija.m$

xxxxxxxxxx
f = @(x) x.^2 + 1 + sin(x);
%f=inline('x.^2+1+sin(x)');

• Napisati M-fajl $integralt.m$ sa funkcijom $[I, briter] = integralt(a, b, tol)$ koja sa tačnošću $tol$ računa i vraća približnu vrednost određenog integrala funkcije $f$ (granice integracije su $a$ i $b$) korišćenjem uopštene Trapezne kvadraturne formule. Funkcija vraća vrednost integrala i broj iteracija.

x
function [I, briter] = integralt(a, b, tol)
​
% Poziv funkcije: [I,briter] = integralt(0,2,1e-4)
​
n1 = 3; % Krećemo od 3 cvora, 2 podintervala
I1 = trapez(a, b, n1);
​
% Potrebna nam je vrednost kvadraturne formule sa duplo manjim korakom, a to ćemo     
% dobiti dodavanjem središnje tačke unutar svakog podintervala. Broj podintervala  
% je n1-1 što je ujedno i broj novododatih tačaka, pa je ukupan broj tačaka 2*n1-1
n2 = 2 * n1 - 1; 
I2 = trapez(a, b, n2);
​
runge = abs(I2-I1) / 3; % Rungeova ocena greške Trapezne kv. formule 
briter = 2; % Već su izračunate 2 približne vr. integrala I1 i I2
​
% Dokle god je Rungeova ocena greške veća od dozvoljene tačnosti, polovimo
% korak i računamo novu vrednost integrala
while runge > tol
    I1 = I2;
    n2 = 2 * n2 - 1; 
    I2 = trapez(a, b, n2);
    runge = abs(I2-I1) / 3;
    briter = briter + 1;
end
I = I2;

x
function I = trapez(a, b, n)
​
% Pomoćni fajl za integralt.m 
% Uopštenom Trapeznom kv. formulom računa približnu vrednost integrala 
% funkcije f na segmetu [a,b] koristeći n čvorova. 
% Kako ćemo u svakoj iteraciji poloviti korak, tj. povećavati
% broj čvorova, potrebno je broj čvorova n prosleđivati kao argument
​
funkcija;
X = linspace(a, b, n);
h = (b - a) / (n - 1); % korak (n - cvorova => n-1 intervala izmedju njih)
Y = f(X);
​
% Trapezna kvadraturna formula
I = (h / 2) * (Y(1) + Y(end) + 2 * sum(Y(2:end-1)));

• Napisati M-fajl $integrals.m$ sa funkcijom $[I, briter] = integrals(a, b, tol)$ koja sa tačnošću $tol$ računa i vraća približnu vrednost određenog integrala funkcije $f$ (granice integracije su $a$ i $b$) korišćenjem uopštene Simpsonove kvadraturne formule. Funkcija vraća vrednost integrala i broj iteracija.

x
  function [I, briter] = integrals(a, b, tol)
  
% Poziv funkcije: [I,briter] = integrals(0,2,1e-4)
  
% NAPOMENA: Ovde je neophodno krenuti sa neparnim brojem čvorova!!! 
  n1 = 3; % Krećemo od 3 cvora, 2 intervala
  I1 = simps(a, b, n1);
  n2 = 2 * n1 - 1;
  I2 = simps(a, b, n2);
​
  runge = abs(I2-I1) / 15; % Rungeova ocena greške Simpsonove kv. formule
  briter = 2;
​
  while runge > tol
      I1 = I2;
      n2 = 2 * n2 - 1;
      I2 = simps(a, b, n2);
      runge = abs(I2-I1) / 15;
      briter = briter + 1;
  end
  I = I2;

xxxxxxxxxx
  function I = simps(a, b, n)
​
  % Pomoćni fajl za integrals.m 
  % Uopštenom Simpsonovom kv. formulom računa približnu vrednost integrala 
  % funkcije f na segmetu [a,b] koristeći n čvorova. 
  % Kako ćemo u svakoj iteraciji poloviti korak, tj. povećavati
  % broj čvorova, potrebno je broj čvorova n prosleđivati kao argument
​
  funkcija;
  X = linspace(a, b, n); 
  h = (b - a) / (n - 1); 
  Y = f(X);
​
  % Simpsonova kvadraturna formula
  I = (h / 3) * (Y(1) + Y(end) + 2 * sum(Y(3:2:end-1)) + 4 * sum(Y(2:2:end-1)));
  

• Napisati M-fajl $grafik.m$ sa funkcijom $grafik(a, b)$ koja prikazuje grafik zavisnosti brzine konvergencije Simpsonove kvadraturne formule (plavo) i Trapezne kvadraturne formule (crveno), za različite tolerancije $(tol = 10^{-1},\ldots, 10^{-6}).$

  xxxxxxxxxx
  function grafik(a, b)
  ​
  % Poziv funkcije: grafik(0,10)
  ​
  tol = [1e-1 1e-2 1e-3 1e-4 1e-5 1e-6];
  %tol=1./10.^(1:6);
  ​
  briter_t = zeros(1, 6);
  briter_s = zeros(1, 6);
  ​
  ​
  for i = 1 : 6
      [~, briter] = integralt(a, b, tol(i)); % Vrednost I nam ne treba, zato ~
      briter_t(i) = briter;
      [~, briter] = integrals(a, b, tol(i));
      briter_s(i) = briter;
  end
  ​
  plot(tol, briter_t, 'r', tol, briter_s, 'b');
  legend('Trapezna', 'Simpsonova')
  xlabel('tolerancije')
  ylabel('broj iteracija')
  

12. zadatak

Neka je funkcija $f$ (koja ne mora biti samo pozitivna) zadata eksplicitno komandnim M-fajlom $funkcija.m$.

xxxxxxxxxx
% Komandni fajl
%f = inline('1./(x+cos(x))');
f = inline('sin(x)');
​
% Funkciju možemo zadati i u okviru funkcijskog fajla
%function y = funkcija(x) 
% y = sin(x);

• Napisati M-fajl $Runge.m$ sa funkcijom $r=Runge(S1, S2)$ koja vraća vrednost Rungeove ocene greške uopštene Simpsonove kvadraturne formule, ako su $S1$ i $S2$ njene vrednosti od kojih je jedna izračunata sa dvostruko manjim korakom u odnosu na drugu.

  xxxxxxxxxx
  function r = Runge(S1, S2)
  ​
  % Računamo Rungeovu ocenu greške Simpsonove kvadraturne formule
  % S1 - vrednost izračunata Simpsonovom kvadraturnom formulom sa korakom H
  % S2 - vrednost izračunata Simpsonovom kvadraturnom formulom sa korakom h=H/2
  ​
  k = 4;
  r = abs(S1-S2) / (2^k - 1);
  

• Napisati M-fajl $zapremina.m$ sa funkcijom $V=zapremina(a, b, tol)$ koja koristeći uopštenu Simp- sonovu kvadraturnu formulu vraća zapreminu tela nastalog obrtanjem figure ograničene pravama $y = 0, x = a, x = b$ i funkcijom $f$ oko ose $Ox$ izračunatu sa tačnošću $tol$. (Za ocenu tačnosti koristiti funkciju $Runge$.)

  xxxxxxxxxx
  function V = zapremina(a, b, tol)
  ​
  % Poziv funkcije: V = zapremina(1, 3, 1e-4)
  % Provera za f(x)=sin(x): V=pi*integral(g,1,3) gde je
  % g = @(x) f(x).^2 (ugrađenoj funkciji integral kao prvi argument
  % treba proslediti anonimnu funkciju)
  ​
  n = 3;
  h = (b - a) / (n - 1);
  ​
  S1 = Simpson_zap(a, b, h);
  H = h / 2; %polovimo korak
  S2 = Simpson_zap(a, b, H);
  ​
  % Dokle god je Rungeova ocena greške veća od dozvoljene tačnosti, polovimo
  % korak i računamo novu vrednost integrala
  while Runge(S1, S2) > tol
      S1 = S2;
      H = H / 2;
      S2 = Simpson_zap(a, b, H);
  end
  ​
  V = pi * S2;
  

  xxxxxxxxxx
  function I = Simpson_zap(a, b, h)
  ​
  % Pomoćni fajl za zapremina.m
  % Računamo integral od a do b funkcije f(x)^2
  ​
  funkcija;
  % Potrebni su nam čvorovi da bismo dobili vektor vrednosti 
  % kvadrata funkcije f u zadatim čvorovima. Nađimo najpre broj podintervala
  % dužine h na segmentu [a,b]
  % n - broj podintervala.
  n = (b - a) / h;
  ​
  % Za n podintervala, imaćemo n+1 čvor (n+1 mora biti neparan)
  X = linspace(a, b, n+1);
  ​
  Y = f(X).^2;
  I = (h / 3) * (Y(1) + 2 * sum(Y(3:2:end-1)) + 4 * sum(Y(2:2:end-1)) + Y(end));
  

13. zadatak - domaći

• Formirati M-fajl $sistem.m$ sa funkcijom $[B, b]=sistem(d, t, n)$ koja formira sistem linearnih jednačina koji se dobija prilikom nalaženja koeficijenata kvadraturne formule oblika $\int_{0}^{d}t(x)f(x)dx\approx \sum_{i=0}^{n} A_if(\frac{i*d}{n})$ koja treba da je tačna za polinome što je moguće većeg stepena. Funkcija treba da vraća matricu sistema i vektor desne strane.

• Formirati M-fajl $koeficijenti.m$ sa funkcijom $koeficijenti(d, t, n)$ koja određuje koeficijente $A_i$ gore napisane kvadraturne formule. Dozvoljeno je korišćenje operatora \ za rešavanje sistema. (* Nakon sistema linearnih jednačina, zadatak se može rešavati i nekom od metoda za sisteme linearnih jednačina: LU, iterativna,...).

  test primer:

  xxxxxxxxxx
  >> t=@(x) sin(x);
  ​
  >> [B,b]=sistem(1,t,5)
  B =
   1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
   0 0.2000 0.4000 0.6000 0.8000 1.0000
   0 0.0400 0.1600 0.3600 0.6400 1.0000
   0 0.0080 0.0640 0.2160 0.5120 1.0000
   0 0.0016 0.0256 0.1296 0.4096 1.0000
   0 0.0003 0.0102 0.0778 0.3277 1.0000
  b =
   0.4597
   0.3012
   0.2232
   0.1771
   0.1467
   0.1251
  >> koeficijenti(1,t,5)
  ans =
   0.0051
   0.0270
   0.1151
   0.0531
   0.2077
   0.0517
  

14. zadatak

• Napisati M-fajl $legendre\_poly.m$ sa funkcijom $L = legendre\_poly(n)$ koja formira i vraća niz SVIH Ležandrovih polinoma do stepena $n$ na intervalu [-1, 1]. Nacrtati grafik svih formiranih Ležandrovih polinoma.

  xxxxxxxxxx
  function L = legendre_poly(n)
  % Funkcija formira i vraća niz svih Ležandrovih polinoma do stepena n
  ​
  % Funkcija vraća cell array L, tako da je L{i}=Li-1
  % tj. na i-toj poziciji nalazi se Ležandrov polinom stepena i-1
  ​
  L=cell(1,n+1);
  L{1} = 1; % L0 = 1
  L{2} = [1, 0]; % L1 = x
  ​
  % Ostali polinomi se računaju na osnovu rekurentne formule
  % Li(x)=((2i-1)*x*Li-1(x) - (i-1)*Li-2(x))/i
  ​
  % Množenje polinoma sa x ekvivalentno je dodavanju nule na kraj
  % Dodavanje nule na početak ne menja vektor i vrši se kako bi se omogućilo 
  % sabiranje vektora, tj. kako bismo izjednačili njihove dimenzije 
  for i = 2:n
      L{i+1} = ((2 * i - 1) * [L{i}, 0] - (i - 1)* [0, 0, L{i-1}])/i;
  end
  celldisp(L);
  ​
  % Crtamo grafik sa svim polinomima u različitim bojama
  X = linspace(-1, 1);
  colors = 'bgrcmyk'; % 7 različitih boja
  axis equal
  ​
  hold on
  for i = 1:length(L)
      plot(X, polyval(L{i}, X), colors(mod(i, 7)+1));
  end
  hold off
  ​
  

• Napisati M-fajl $Cebisev\_poly.m$ sa funkcijom $C= Cebisev\_poly(n)$ koja formira i vraća niz SVIH Čebiševljevih polinoma do stepena $n$ na intervalu [-1, 1]. Nacrtati grafik svih formiranih Čebiševljevih polinoma.

  xxxxxxxxxx
  ​
  function C = Cebisev_poly(n)
  % Funkcija formira i vraća niz svih Čebiševljevih polinoma do stepena n
  ​
  % Funkcija vraća cell array C, tako da je C{i}=Ci-1
  % tj. na i-toj poziciji nalazi se Čebiševljev polinom stepena i-1
  ​
  C=cell(1,n+1);
  C{1} = 1; %C0=1
  C{2} = [1, 0]; %C1=x
  ​
  % Ostali polinomi se računaju na osnovu rekurentne formule
  %  Ci(x) = 2*x*Ci-1(x) - Ci-2(x)
  for i = 2:n
      C{i+1} = 2 * [C{i}, 0] - [0, 0, C{i-1}];
  end
  celldisp(C);
  ​
  %Crtmo grafik
  X = linspace(-1, 1);
  colors = 'bgrcmyk';
  ​
  hold on
  for i = 1:length(C)
      plot(X, polyval(C{i}, X), colors(mod(i, 7)+1));
  end
  hold off
  ​
  

• Napisati M-fajl $integrali.m$ sa funkcijom $integrali(f)$ koja korišćenjem ugrađene MATLAB funkcije integral() računa i štampa vrednosti sledećih integrala: $\int_{-1}^{1}f(x)dx, \int_{-1}^{1}f(x)\cdot sin(x)dx, \int_{-1}^{1}f(x)\cdot L_5(x)dx, \int_{-1}^{1}L_3(x)\cdot L_5(x)dx$ gde je $L_i(x)$ Ležandrov polinom stepena $i$. Prosleđena funkicija $f$ može biti složena funkcija.

  xxxxxxxxxx
  function integrali(f)
  ​
  % Ako je funkcija prosleđena kao string, tada moramo kreirati inline
  % objekat (ne možemo direktno anonimnu funkciju)
  ​
  % Ako je funkcija prosleđena kao anonimna, onda je direktno možemo proslediti 
  % ugrađenoj funkciji integral()
  f = inline(f);
  g = @(x) sin(x);
  L = legendre_poly(5);
  ​
  disp('Vrednost integrala funkcije f(x) je:');
  integral(@(x) f(x),-1,1)
  disp('Vrednost integrala funkcije f(x)*sin(x) je:');
  integral(@(x) f(x).*g(x),-1,1)
  disp('Vrednost integrala funkcije f(x)*L_5 je:');
  integral(@(x) f(x).*polyval(L{6},x),-1,1)
  disp('Vrednost integrala funkcije L_3*L_5 je:');
  integral(@(x) polyval(L{4},x).*polyval(L{6},x),-1,1)
  %quad(@(x) polyval(L{4},x).*polyval(L{6},x),-1,1)
  % U prethodnim verzijama MATLAB-a quad() se koristio umesto integral()
  ​
  

15. zadatak - domaći

• Napisati M-fajl $legendre.m$ sa funkcijom $L = legendre(n)$ koja kao rezultat vraća Ležandrov polinom $L$ stepena $n$ na intervalu $[-1, 1]$.

• Napisati M-fajl $polinom.m$ sa funkcijom $P = polinom(n,m)$ koja kao rezultat vraća polinom $P$ dobijen preko formule $P(x)=(1-x^2)\frac{d^m}{dx^m}L_n(x)$, gde je $L_n(x)$ Ležandrov polinom stepena $n$ na intervalu $[-1, 1]$.

• Napisati M-fajl $integral.m$ sa funkcijom $I = integral(n, m, tol)$ koja sa tačnošću $tol$ približno određuje i kao rezultat vraća vrednost integrala $\int_{-1}^{1}P(x)e^xdx$. Integral računati korišćenjem uopštene Simpsonove formule. Polinom $P(x)$ je polinom dobijen pod (2).

  test primer

  xxxxxxxxxx
  >> L=legendre(5)
  ​
  L =
  ​
      7.8750         0   -8.7500         0    1.8750         0
  ​
  >> P=polinom(5,3)
  ​
  P =
  ​
   -472.5000         0  525.0000         0  -52.5000
  ​
  >> I = integral(5,3,1e-4)
  ​
  I =
  ​
     77.0221
  ​
  Provera:
  >> integral(@(x) polyval(P,x).*exp(x), -1,1)
  ​
  ans =
  ​
     77.0222