Interpolacija i diferenciranje

1. zadatak

Neka je funkcija $f$ zadata tablično M-fajlom $tablica.m$ koji generiše dva niza $X=[x_{1},x_{2}...,x_{n}]$ i $F=[f_{1},f_{2}...,f_{n}]$ (od kojih je prvi strogo rastući) za tu tablično zadatu funkciju. Tablica ne mora biti ekvidistantna.

$tablica.m$

​x
% Komandni fajl sa vektorima X i F
​
% X=[0 0.2 0.4 0.6 0.8 1];
% Umesto da navodimo sve ove elemente, možemo pisati
% X je vektor vrednosti od 0 : sa korakom 0.2 : do 1
X = 0:0.2:1;
​
% F=[1 2 3 4 5 6];
F = 1:6;
​
% primer sa vežbi:
% X=[100 121 144];
% F=[10 11 12];

• Napisati M-fajl $novatablica.m$ u kom se prethodna tablica proširuje do nove dodavanjem čvorova $\frac{x_{i}+x_{i+1}}{2}, i=1,...n-1$ i računanjem vrednosti funkcije $f$ u njima korišćenjem formule: $f(\frac{x_{i}+x_{i+1}}{2})=\frac{f(x_{i})+f(x_{i+1})}{2}, i=1,...n-1.$

   

  xxxxxxxxxx
  % Komandani fajl koji od stare tablice u komandnom fajlu tablica.m (X i F)
  % formira novu tablicu (X1 i F1)
  ​
  tablica; % Učitavamo podatke iz tablica.m (X i F)
  n = length(X); % Dužina niza X
  ​
  % Poželjno je formirati nula nizove dužine 2*n-1 (dužina proširenog niza) kako ne bismo u %svakom prolasku for petlje menjali dimenziju nizova 
  X1 = zeros(1, 2*n-1);
  F1 = zeros(1, 2*n-1);
  ​
  for i = 1:2:2*n-1
      % indeks mora biti ceo broj
      X1(i) = X((i + 1)/2);
      F1(i) = F((i + 1)/2);
      % ili    X1(i)=X(round(i/2));
      %        F1(i)=F(round(i/2));
  end
  ​
  % Računamo nove elemente
  for i = 2:2:2*n-1
      X1(i)=(X1(i+1)+X1(i-1))/2;
      F1(i)=(F1(i+1)+F1(i-1))/2;
      % ili    X1(i) = (X(i/2) + X(i/2+1)) / 2;
      %        F1(i) = (F(i/2) + F(i/2+1)) / 2;
  end
  ​
  % Funkcije za zaokruživanje: round(), ceil() i floor()
  % round(1.8)=2
  % round(1.3)=1
  % round(1.5)=2
  % ceil(1.3)=2
  % floor(1.8)=1
  

• Napisati M-fajl $Lagr1.m$ sa funkcijom $L = Lagr1(x)$ koja za uneti argument $x$ vraća približnu vrednost funkcije $f$ u toj tački izračunatu pomoću Lagranžovog interpolacionog polinoma $L$ korišćenjem svih vrednosti iz nove tablice.

  xxxxxxxxxx
  function L = Lagr1(x)
  ​
  % Računamo vrednost Lagranžovog interpolacionog polinoma u tački x
  % Argumenti funkcije:
  %   x - (broj) tačka u kojoj računamo vrednost Lagranžovog int. polinoma
  % Funkcija vraća:
  %   L - vrednost Lagranžovog int. polinoma u tački x 
  % NAPOMENA: ne konstruišemo polinom, ne nalazimo koeficijente polinoma!
  ​
  novatablica;
  L = 0;
  n = length(X1);
  for i = 1 : n
      p = 1;
      for j = 1 : n
          if i ~= j
              p = p * (x - X1(j)) / (X1(i) - X1(j));
          end
      end
      L = L + p * F1(i);
  end
  

2. zadatak

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

$funkcija.m$

% Primer funkcije kod koje se sa povećanjem broja čvorova povećava i greška (Rungeov fenomen)
f = @(x) 1 ./ (1 + 25 * x.^2);
​
% Drugi način definisanja funkcije 
% f=inline('1./(1+25*x.^2)');
​
% neki drugi primeri funkcija (anonimne funkcije)
%f=@(x) sin(x)/3+cos(x).^2;
%f=@(x) 45./(10.*x.^2+1);

• Napisati M-fajl $tablica.m$ sa funkcijom $[X, Y ] = tablica(a, b, n)$ koja tabelira zadatu funkciju $f$ na intervalu $[a, b]$ sa $n$ čvorova.

  xxxxxxxxxx
  function [X, Y] = tablica(a, b, n)
  ​
  % Agrumenti funkcije:
  %   a i b - granice intervala
  %   n - broj tačaka na koji delimo interval [a,b]
  % Funkcija vraća:
  %   X - čvorove interpolacije
  %   Y - vrednosti funkcije u tim čvorovima
  ​
  % Poziv funkcije: [X Y] = tablica(-2,2,5)
  ​
  funkcija;
  ​
  % I način
  % n ravnomerno raspoređenih tačaka od a do b => broj podintervala  je n-1 
  % dužina svakog od njih (korak h) iznosi    
  h = (b - a) / (n - 1);
  ​
  for i = 1 : n
      X(i) = a + (i - 1) * h;
      Y(i) = f(X(i));
      %Y(i)=feval(f,X(i));
  end
  ​
  % II način - ugrađena funkcija linspace()
  % X = linspace(a, b, n);
  % Y = f(X);
  ​
  

   

• Napisati M-fajl $Lagr1b.m$ sa funkcijom $[L, y] = Lagr1b(x, a, b, n)$ koji formira i vraća koeficijente Lagranžovog interpolacionog polinoma $L$ formiranog koristeći sve vrednosti iz tablice, kao i vrednost formiranog polinoma u tački $x$.

  xxxxxxxxxx
  function [L, y] = Lagr1b(x, a, b, n)
  ​
  % Poziv funkcije, npr: [L,y]=Lagr1b(0.5,-2,2,5)
  % Argumenti funkcije:
  %   x - tačka u kojoj se računa vrednost Lagranžovog int. polinoma
  %   a i b - granice intervala [a,b]
  %   n - broj čvorova
  % Funkcija vraća:
  %   L - koeficijente odgovarajućeg Lagranzovog int. polinoma 
  %   y - vrednost polinoma u tački x
  %   NAPOMENA: sada konstruišemo polinom!
  ​
  [X, Y] = tablica(a, b, n);
  n = length(X);
  L = zeros(1, n); % L je sada polinom
  ​
  for i = 1 : n
      p = 1; % p je trenutno polinom nultog stepena, tj. p=1*x^0
      for j = 1 : n
          if i ~= j
              % Množimo 2 polinoma
              % [1 -X(j)] je polinom 1*x^1 - X(j)*x^0
              % NAPOMENA: paziti na razmak
              % [1 - X(j)] nije isto što i [1 -X(j)], prvi vektor je polinom (1-X(j))*x^0
              p = conv(p, [1, -X(j)]/(X(i) - X(j)));
          end
      end
      L = L + p * Y(i);
  end
  ​
  % Računamo vrednost polinoma L u tački x
  y = polyval(L, x);
  

• Uporediti grafike funkcije $f$ i formiranog interpolacionog polinoma.

  xxxxxxxxxx
  function grafik(a, b, n)
  ​
  % Poziv funkcije:
  % grafik(-1,1,5) interpolacija sa 5 tačaka
  % grafik(-1,1,10) interpolacija sa 10 tačaka
  % grafik(-1,1,11) interpolacija sa 11 tačaka
  ​
  [L, y] = Lagr1b(1, a, b, n); %x je nebitan za grafik
  funkcija;
  ​
  % Delimo interval [a,b] na 100 tačaka
  x = linspace(a, b);
  ​
  % Želimo više grafika na jednoj slici (moramo uključiti hold on i kasnije isključiti sa
  % hold off). Svaki poziv plot() izmedju ove dve naredbe crtaće novi
  % grafik i pritom se grafici prethodnog pozivanja funkcije plot() neće obrisati
  ​
  hold on
  ​
  % Crtamo prvi grafik funkcije
  plot(x, f(x));
  % Crtamo drugi grafik polinoma
  plot(x, polyval(L, x), 'r');
  ​
  hold off
  ​
  legend('funkcija f', 'interpolacioni polinom')
  % Izjednačavamo ose
  axis equal
  

3. zadatak

Neka je funkcija $f$ zadata tablično M-fajlom $tablica.m$ koji generiše dva niza $X=[x_{1},x_{2}...,x_{n}]$ i $Y=[y_{1},y_{2}...,y_{n}]$ za tu tablično zadatu funkciju.

$tablica.m$

xxxxxxxxxx
X = 15:5:55;
Y = [0.2588, 0.3420, 0.4226, 0.5000, 0.5736, 0.6428, 0.7071, 0.7660, 0.8192];

• Napisati M-fajl $tablicaCheck.m$ sa funkcijom $t = tablicaCheck()$ koja vrši proveru da li je tablica u komandnom fajlu $tablica.m$ ekvidistantna i da li je niz $X$ zadat u strogo rastućem poretku. Ukoliko su oba uslova ispunjena funkcija vraća vrednost $1$, u suprotnom vraća vrednost $0$ i u oba slučaja ispisuje odgovarajuću poruku.

  xxxxxxxxxx
  function t = tablicaCheck()
  ​
  % Funkcija vraća:
  %       t=1 - ukoliko je sve u redu
  %       t=0 - ukoliko nije ispunjen jedan od ova dva uslova
  ​
  tablica;
  r = all(diff(X) > 0);
  h = X(2) - X(1);
  % ekv=all(diff(X)==h); Radi samo za celobrojne vrednosti koraka h!
  % zbog netačne reprezentacije relane brojeve nikada ne treba porediti na jednakost
  % već smatramo da su jednaki ako je apsolutna vrednost njihove razlike manja od nekog %dovoljno malog broja
  ekv = all(abs(diff(X)-h) <= 1e-10);
  t = r && ekv;
  if t == 1
      disp('Tablica je ekvidistantna i vektor X je zadat u strogo rastucem poretku.')
  else
      disp('Nisu ispunjeni uslovi zadatka!')
  end
  

• Napisati M-fajl $polozaj.m$ sa funkcijom $polozaj(x)$ koja za uneti argument $x$ vraća vrednost $1$ ukoliko je $x < x_{2}$, $2$ ukoliko je $x > x_{n-1}$ i $0$ inače.

  xxxxxxxxxx
  % Određuje položaj tačke x u odnosu na čvorove interpolacije
  function pol = polozaj(x)
  ​
  tablica;
  n = length(X);
  ​
  if x < X(2)
      pol = 1;
  else if x > X(n-1)
          pol = 2;
      else
          pol = 0;
      end
  end
  

• Napisati M-fajl $Njutn.m$ sa funkcijom $Njutn(x)$ koja ukoliko su svi uslovi ispunjeni, vraća približnu vrednost funkcije $f$ u tački $x$ izračunatu korišćenjem I (II) Njutnovog interpolacionog polinoma, ako je vrednost funkcije $polozaj$ u tački $x$ jednaka $1(2)$, odnosno izdaje odgovarajuću poruku ukoliko je $polozaj(x) = 0.$

  xxxxxxxxxx
  function nj = Njutn(x)
  ​
  if tablicaCheck() == 1 % da li su ispunjeni uslovi zadatka
      if polozaj(x) == 1 % da li je x na početku tablice
          nj = Njutn1(x);
      else
          if polozaj(x) == 2 % da li je x na kraju tablice
              nj = Njutn2(x);
          else
              error('Tacka se nalazi u sredini');
          end
      end
  else
      error('Nisu ispunjeni uslovi zadatka');
  end
  ​
    % Ugrađena funkcija error() prekida izvršavanje progrma i vraća poruku (string koji dobija 
    % kao argument) kojom nas obaveštava zašto je prekinutno izvršavanje.
  ​
  

  xxxxxxxxxx
  ​
  function nj1 = Njutn1(x)
  ​
  % Pomoćna funkcija koju koristi Njutn.m za računanje vrednosti I Njutnovog
  % interpolacionog polinoma u tački x
  ​
  tablica;
  n = length(X);
  ​
  % KONAČNE RAZLIKE:
  k_razlike = zeros(n, n-1);
  ​
  % For petlja koja popunjava prvu kolonu tabele konačnih razlika, jer se
  % one računaju oduzimanjem odgovarajućih elemenata vektora Y
  for i = 1 : n-1
      k_razlike(i, 1) = Y(i+1) - Y(i);
  end
  ​
  % For petlja koja računa konačne razlike od drugog do n-1 reda 
  % Svako izračunavanje se vrši oduzimanjem odgovarajućih vrednosti
  % iz prethodne kolone tabele konačnih razlika
  for j = 2 : n-1 % po kolonama (red konačne razlike)
      for i = 1 : n-j % po vrstama (poslednji čvor za koje se može naći k.r. reda j je n-j)
          k_razlike(i, j) = k_razlike(i+1, j-1) - k_razlike(i, j-1);
      end
  end
  ​
  disp(k_razlike);
  % Želimo da matricu k_razlike nadovežemo na čvorove i vrednosti funk. u čvorovima i da 
  % dobijemo tablicu kao na vežbama
  % disp([X' Y' k_razlike])
  ​
  % I NJUTNOV INTERPOLACIONI POLINOM:
  y = Y(1);
  h = X(2) - X(1);
  q = (x - X(1)) / h;
  % čuvamo početnu vrednost za q jer nam treba za formiranje svakog sledećeg sabirka 
  Q = q;
  ​
  for j = 1 : n-1
      y = y + q * k_razlike(1, j) / factorial(j);
      q = q * (Q - j);
  end
  ​
  nj1 = y;
  % Samo računamo vrednost polinoma u tački, ne konstruišemo polinom!
  

  xxxxxxxxxx
  function nj2 = Njutn2(x)
  ​
  % Pomoćna funkcija koju koristi Njutn.m za racunanje II Njutnovog
  % interpolacionog polinoma u tački x
  ​
  tablica;
  n = length(X);
  ​
  % KONAČNE RAZLIKE:
  k_razlike = zeros(n, n-1);
  ​
  for i = 1 : n-1
      k_razlike(i, 1) = Y(i+1) - Y(i);
  end
  ​
  for j = 2 : n-1
      for i = 1 : n-j
          k_razlike(i, j) = k_razlike(i+1, j-1) - k_razlike(i, j-1);
      end
  end
  ​
  disp(k_razlike);
  ​
  % II NJUTNOV INTERPOLACIONI POLINOM:
  y = Y(end);
  h = X(2) - X(1);
  q = (x - X(end)) / h;
  Q = q;
  ​
  for j = 1 : n-1
      y = y + q * k_razlike(n-j, j) / factorial(j);
      q = q * (Q + j);
  end
  ​
  nj2 = y;
  

4. zadatak

Neka je funkcija $f$ zadata tablično M-fajlom $tablica.m$ koji generiše dva niza $X=[x_{1},x_{2}...,x_{n}]$ i

$Y=[y_{1},y_{2}...,y_{n}]$ za tu tablično zadatu funkciju. Tablica ne mora biti ekvidistantna.

$tablica.m$

xxxxxxxxxx
X = [2 2.5 3.5 4];
Y = sin(X);
​
% ili npr.
% Y=X+1
%Y=[3 3.5 4.5 5];

• Napisati M-fajl $tablicaCheck.m$ sa funkcijom $t = tablicaCheck()$ koja vrši proveru da li je niz $X$ zadat u strogo rastućem poretku i da li je niz $Y$ monoton. Ukoliko su oba uslova ispunjena funkcija vraća vrednost $1$, u suprotnom vraća vrednost $0$. Ukoliko neki od uslova nije ispunjen, funkcija ispisuje odgovarajuću poruku.

  xxxxxxxxxx
  function t = tablicaCheck()
  ​
  % Funkcija vraća:
  %       t=1 - ukoliko je sve u redu
  %       t=0 - ukoliko nije ispunjen jedan od ova dva uslova
  ​
  tablica;
  n = length(X);
  ​
  provera_x = 1;
  provera_y = 1;
  ​
  % Da li je X strogo rastući?
  for i = 1 : n-1
      if X(i) >= X(i+1)
          provera_x = 0;
          disp('Niz X nije rastuci');
          break
      end
  end
  ​
  % Da li je Y monoton?
  i = 1;
  if Y(2) > Y(1) % Pretenduje da bude strogo rastući
      % Operator && ne računa logičku vrednost drugog operanda, ako prvi ima vr. 0
      while (i < n && Y(i+1) > Y(i)) % Bitan redosled, kako se ne bi pristupalo Y(n+1)
          i = i + 1;
      end
  else if Y(2) < Y(1) % Pretenduje da bude strogo opadajući
          while (i < n && Y(i+1) < Y(i))
              i = i + 1;
          end
      end
  end
  ​
  if i < n % prekinuto je izvršavanje jedne od while petlji jer nije ispunjen uslov za Y
      provera_y = 0;
      disp('Niz Y nije monoton');
  end
  ​
  if provera_x == 1 && provera_y == 1
      % if (provera_x && provera_y)
      t = 1;
  else
      t = 0;
  end
  % Ili umesto if-a:  t=proverax && proveray
  ​
  % II Nacin:
  ​
    % provera_x=all(diff(X)>0);
    % provera_y=all(diff(Y)>0) || all(diff(Y)<0);
    % t=provera_x && provera_y;
  ​
  

• Napisati M-fajl $vredfunk.m$ sa funkcijom $y = vredfunk(x)$ koja za uneti argument $x$ vraća

približnu vrednost funkcije $f$ u toj tački izračunatu pomoću Njutnovog interpolacionog polinoma sa podeljenim razlikama korišćenjem svih vrednosti iz tablice.

xxxxxxxxxx
​
  function y = vredfunk(x)
​
  tablica;
  n = length(X);
​
  % Da li su ispunjeni uslovi zadatka?
  if tablicaCheck() == 0
      error('Nisu ispunjeni uslovi zadatka');
  end
​
  % Podeljene razlike
  p_razlike = zeros(n, n-1);
​
  % For petlja koja popunjava prvu kolonu tabele podeljenih razlika, jer se
  % one racunaju oduzimanjem odgovarajucih elemenata vektora Y
  for i = 1 : n-1
      p_razlike(i, 1) = (Y(i+1) - Y(i)) / (X(i+1) - X(i));
  end
​
  % For petlja koja računa podeljene razlike od drugog do n-1 reda
  % Svako izračunavanje se vrši oduzimanjem odgovarajućih vrednosti
  % iz prethodne kolone tabele podeljenih razlika 
  for j = 2 : n-1 % po kolonama (red podeljene razlike)
      for i = 1 : n-j % po vrstama
          p_razlike(i, j) = (p_razlike(i+1, j-1) - p_razlike(i, j-1)) / (X(i+j) - X(i));
      end
  end
​
  disp('Stampamo tablicu po kolonama:');
  % "lepimo" vektore X i Y na početak tablice podeljenih razlika
  disp([X', Y', p_razlike]);
​
  % Njutnov interpolacioni polinom sa podeljenim razlikama
  y = Y(1);
  p = 1;
  for i = 1 : n-1
      p = p * (x - X(i));
      y = y + p * p_razlike(1, i);
  end
​
  disp(y);
​
   
  % II nacin za formiranje matrice podeljenih razlika - funkcija diff()
​
  % prazlike(1:n-1,1)=diff(Y)./diff(X);
  % for j=2:n-1
  %     prazlike(1:n-j,j)=diff(prazlike(1:n-j+1,j-1))./(X(1+j:n)-X(1:n-j))';
  % end

​

5. zadatak

Neka su u komandnom fajlu $podaci.m$ dati funkcija $f$ i vektor $X$ koji sadrži samo celobrojne vrednosti.

$podaci.m$

f = @(x) (x + 1) / 3;
X = [-5 -4 -2 -1 2 3 5 6 7 8 9 10]
f = @(x) (x + 1) / 3;

• Napisati M-fajl $tablica.m$ sa funkcijom $[X1, Y 1] = tablica()$ koja formira tablicu gde se vektor $X1$ sastoji samo od parnih vrednosti vektora $X$, a vektor $Y 1$ su vrednosti eksplicitno zadate funkcije $f$ u elementima vektora $X1$ zaokruženi na $3$ decimale.

  xxxxxxxxxx
  ​
  function [X1, Y1] = tablica()
  ​
  podaci;
  ​
  % Ugradjena funkcija find(uslov) vraća vektor pozicija onih elemenata u vektoru (koji je sadržan u uslovu) koji ispunjavaju zadati uslov
  % find(mod(X,2)==0) vratiće vektor pozicija (indeksa) parnih elemenata vektora X
  % X(find(mod(X,2)==0)) vraća parne elemente vektora X
  ​
  X1 = X(find(mod(X, 2) == 0));
  Y1 = f(X1);
  Y1 = round(Y1.*1000) / 1000;
  % Ugrađenoj funkciji round() kao drugi argument možete proslediti na koliko decimala %želite da izvršite zaokruživanje
  % Y1=round(Y1,3);
  

• Napisati M-fajl $inverz.m$ sa funkcijom $x=inverz(y)$ koja za zadatu vrednost $y$ inverznom interpolacijom približno određuje $x$ za koje je $f(x) = y$. (*Tablica neće biti ekvidistantna, pa koristimo Lagranžov interpolacioni polinom)

  x
  ​
  function x = inverz(y)
  ​
  [X, Y] = tablica();
  n = length(X);
  ​
  % Pošto iz postavke zadatka tablica ne mora biti ekvidistantna, za inverznu
  % interpolaciju najpre ćemo invertovati tablicu, a zatim odrediti vrednost
  % Lagranžovog interpolacionog polinoma u tački y
  ​
  inv_Tablica = zeros(2, n);
  % U prvu vrstu upisujemo vektor Y
  inv_Tablica(1, :) = Y;
  % U drugu vrstu upisujemo vektor X
  inv_Tablica(2, :) = X;
  ​
  % Tablicu sortiramo po Y u rastućem poretku
  % Ugradjena funkcija sortrows(M,br_kol) sortiraće matricu M po koloni br_kol
  % Želimo sortiranje po prvoj vrsti matrice inv_Tablica, tj. po prvoj koloni
  % transponovane matrice
  inv_Tablica = (sortrows(inv_Tablica', 1))';
  ​
  % Prikazujemo kako izgleda naša tablica nakon invertovanja i sortiranja
  disp('Invertovana tablica:');
  disp(inv_Tablica);
  ​
  % Izvlačimo prvu vrstu invertovane tablice (čvorovi za interpolaciju inverza funk. f) % u vektor Y1 i drugu vrstu (vrednosti inverza funkcije f) u vektor X1
  Y1 = inv_Tablica(1, :);
  X1 = inv_Tablica(2, :);
  ​
  % Računamo vrednost Lagranžovog interpolacionog polinoma  u tački y (ne konstruišemo polinom!)
  L = 0;
  for i = 1 : n
      p = 1;
      for j = 1 : n
          if i ~= j
              p = p * (y - Y1(j)) / (Y1(i) - Y1(j)); % u vektoru Y1 se nalaze čvorovi
          end
      end
      L = L + p * X1(i); % u vektoru X1 se nalaze vrednosti inverza funkcije f
  end
  

6. zadatak

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

$funkcija.m$

xxxxxxxxxx
f = inline('x.^2+x/2-exp(x)/4');

• Napisati M-fajl $tablica.m$ sa funkcijom $[X, Y ] = tablica(a, b, n)$ koja formira ekvidistantnu tabelu funkcije $f$ na segmentu $[a, b]$ sa $n$ čvorova.

  x
  funkcija;
  function [X, Y] = tablica(a, b, n)
  ​
  % Poziv funkcije: [X Y] = tablica(1,4.5,10)
  ​
  funkcija; % učitavanje komandnog fajla kako bismo imali na raspolaganju funk. f 
  ​
  % Korak
  h = (b - a) / (n - 1);
  ​
  for i = 1 : n
      X(i) = a + (i - 1) * h;
      Y(i) = f(X(i));
      % računanje vrednosti funk. korišćenjem ugrađene funkcije feval()
      % Y(i) = feval(f, X(i)); 
  ​
  end
  ​
  % Ili u dva reda:
  ​
  % X=linspace(a,b,n) 
  % Y=f(X);
  ​
  % NAPOMENA: Funkciji f prosleđujemo vektor i želimo da kao rezultat dobijemo vektor %vrednosti funkcije f u svim elementima vektora X, stoga je neophodno je uvesti '.' %notaciju tj. '.^'  pri definisanju funkcije f (svaki element iz X treba kvadrirati). 
  

• Napisati M-fajl $promenaZnaka.m$ sa funkcijom $[c, d] = promenaZnaka(a, b, n)$ koja na osnovu nizova $X$ i $Y$ dobijenih pozivanjem funkcije $tablica(a, b, n)$ pronalazi i kao rezultat vraća prvi interval $[x_{i}, x_{i+1}]$ u kome funkcija menja znak $(c = x_{i}, d = x_{i+1})$. Pretpostavlja se da takav interval postoji.

  x
  ​
  function [c, d] = promenaZnaka(a, b, n)
  ​
  % Poziv funkcije: [c,d] = promenaZnaka(1,5,10)
  ​
  [X, Y] = tablica(a, b, n);
  n = length(X);
  for i = 1 : n-1
      if Y(i) * Y(i+1) < 0
          break
      end
  end
  ​
  % II način - korišćenjem ugrađene funk. find()
  % koordinatno množimo vektore [y1 y2 ... yn-1] i [y2 y3 ... yn] i u vekoru rezultata
  % tražimo poziciju prvog negativnog elementa. Kada funk. find() kao drugi
  % argument prosledite 1, ona će vratiti poziciju prvog elementa koji ispunjava uslov
  ​
  %i=find(Y(1:n-1).*Y(2:n)<0,1);
  ​
  c = X(i);
  d = X(i+1);
  

• Napisati M-fajl $nula.m$ sa funkcijom $nula(a, b, n)$ koja metodom inverzne interpolacije približno određuje nulu funkcije $f$ na intervalu $[c, d],$ koristeći II Njutnov interpolacioni polinom zaključno sa konačnim razlikama reda $3$. Kriterijum zaustavljanja iterativnog niza: $|q_{i}-q_{i-1}|\leq 10^{-4}, i = 2...$

  x
  % Nakon tabeliranja na segmentu [c,d] nula ce biti na kraju tablice
  function x = nula(a, b, n)
  ​
  % Poziv funkcije: x=nula(1,4.5,10)
  % Nakon tabeliranja na segmentu [c,d] nula će se nalaziti na kraju tablice
  ​
  [c, d] = promenaZnaka(a, b, n);
  % tabeliramo funk. na podintervalu [c,d] u dodatnih n tačaka
  [X1, Y1] = tablica(c, d, n);
  ​
  h = X1(2) - X1(1);
  % x - nula funkcije
  ​
  % TABLICA KONACNIH RAZLIKA
  % X(1)   Y(1)   kR(1,1)     kR(1,2)     kR(1,3)
  % X(2)   Y(2)   kR(2,1)     kR(2,2)     kR(2,3)
  % X(3)   Y(3)   kR(3,1)     kR(3,2)     kR(3,3)
  %  .       .       .           .           .
  %  .       .       .           .           .
  %  .       .       .           .           .
  % X(n-3) Y(n-3) kR(n-3,1)   kR(n-3,2)   kR(n-3,3)
  % X(n-2) Y(n-2) kR(n-2,1)   kR(n-2,2)
  % X(n-1) Y(n-1) kR(n-1,1)
  % X(n)   Y(n)
  ​
  k_razlike = zeros(n, 3); % potrebne su nam konačne razlike do 3. reda
  ​
  for i = 1 : n-1
      k_razlike(i, 1) = Y1(i+1) - Y1(i);
  end
  ​
  for j = 2 : 3
      for i = 1 : n-j
          k_razlike(i, j) = k_razlike(i+1, j-1) - k_razlike(i, j-1);
      end
  end
  ​
  % Štampamo tabelu konačnih razlika
  disp([X1', Y1', k_razlike]);
  ​
  % II NJUTNOV INTERPOLACIONI POLINOM
  % y = Y(n)+q*k_razlike(n-1,1)+q*(q+1)*k_razlike(n-2,2)/2+...
  %                             q*(q+1)*(q+2)*k_razlike(n-3,3)/6
  ​
  % ODGOVARAJUAĆA ITERATIVNA FORMULA
  % q =(y-Y(n)-q*(q+1)*k_razlike(n-2,2)/2-
  %           -q*(q+1)*(q+2)*k_razlike(n-3,3)/6)/k_razlike(n-1,1)
  ​
  y = 0; % tražimo nulu funkcije
  % q0=0, q1=slobodan član prethodno napisanog polinoma po q
  % u vektor q smeštamo elemente iterativnog niza
  q = [0, (y-Y1(end)) / k_razlike(n-1, 1)];
  ​
  while (abs(q(end)-q(end-1)) > 0.0001)
      q1 = q(end);
      q = [q, (y-Y1(end)-k_razlike(n-2, 2)*q1*(q1 + 1)/2- ...  % za prelazak u novi red
          k_razlike(n-3, 3)*q1*(q1 + 1)*(q1 + 2)/6) / k_razlike(n-1, 1)];
  end
  disp('Nula funkcije je: ');
  x = X1(end) + q(end) * h;
  disp('iterativni niz q je: ');
  q'
  

7. zadatak - domaći

Neka je funkcija $f$ zadata tablično M-fajlom $tablica.m$ koji generiše dva niza $X=[x_{1},x_{2}...,x_{n}]$ i $Y=[y_{1},y_{2}...,y_{n}]$ (od kojih je prvi strogo rastući) za tu tablično zadatu funkciju. Tablica ne mora biti ekvidistantna.

$tablica.m$

x
X = [-2, -1.5, 0, 0.7, 0.9, 1.1];
Y = 2 * X - exp(X);

• Napisati M-fajl $izvod.m$ sa funkcijom $[X_{i}, Y_{i}] = izvod()$ u kom se na osnovu prethodne tablice formira tablica prvog izvoda funkcije $f$ u tačkama $x_{2},...,x_{n-1}$ korišćenjem sledeće formule: $f'(x_{i})=\frac{f(x_{i+1})-f(x_{i-1})}{x_{i+1}-x_{i-1}}$, gde je $Y_{i}=[f'(x_{2}),..., f'(x_{n-1})]$, a $X_{i}=[x_{2},..., x_{n-1}]$.

• Napisati M-fajl $vredizvod.m$ sa funkcijom $y=vredizvod(x)$ koja za uneti argument $x$ vraća približnu vrednost prvog izvoda funkcije $f$ izračunatu korišćenjem Njutnovog interpolacionog polinoma sa podeljenim razlikama konstruisanog na osnovu svih vrednosti iz tablice formirane korišćenjem funkcije $izvod()$.

• Napisati M-fajl $nula.m$ sa funkcijom $nul=nula()$ koja metodom inverzne interpolacije približno određuje i vraća jednu nulu prvog izvoda funkcije $f$ korišćenjem Njutnovog interpolacionog polinoma sa podeljenim razlikama (pretpostavka je da je prvi izvod monotona funkcija).

  test primer

  x
      >> [Xi,Yi] = izvod()
        
      X =
        
         -1.5000         0    0.7000    0.9000   
  ​
        
      Yi =
        
          1.5677    1.1861    0.3782   -0.4760
        
      >> y = vredizvod(0.4)
      >>   Tablica podeljenih razlika:
         -1.5000    1.5677   -0.2544   -0.4090   -1.2727
               0    1.1861   -1.1541   -3.4635         0
          0.7000    0.3782   -4.2713         0         0
          0.9000   -0.4760         0         0         0
      
      y =
        
          1.0637
        
      >> nul = nula()
      >>   Tablica podeljenih razlika:
         -0.4760    0.9000   -0.2341   -0.3805   -1.0745
          0.3782    0.7000   -0.8665   -2.5764         0
          1.1861         0   -3.9310         0         0
          1.5677   -1.5000         0         0         0
      
      nul =
        
          0.6276
  

8. zadatak - domaći

Neka je funkcija $f$ zadata tablično M-fajlom $tablica.m$ koji generiše dva niza $X=[x_{1},x_{2}...,x_{n}]$ i $Y=[y_{1},y_{2}...,y_{n}]$ (od kojih je prvi strogo rastući) za tu tablično zadatu funkciju. Tablica mora biti ekvidistantna ( sa korakom h ).

$tablica.m$

xxxxxxxxxx
X = linspace(1,2,6);
f = @(x) exp(x) - 15 * x;
Y = f(X);

• Napisati M-fajl $drugiizvod.m$ sa funkcijom $[X2_{i}, Y2_{i}] = drugiizvod()$ u kom se na osnovu prethodne tablice formira tablica drugog izvoda funkcije $f$ u tačkama $x_{2},...,x_{n-1}$ korišćenjem sledeće formule: $f''(x_{i})=\frac{f(x_{i-1})-2f(x_{i})+f(x_{i+1})}{h^2}$, gde je $Y2_{i}=[f''(x_{2}),..., f''(x_{n-1})],$ a $X2_{i}=[x_{2},..., x_{n-1}].$

• Napisati M-fajl $vred2izvod.m$ sa funkcijom $y=vred2izvod(x)$ koja za uneti argument $x$ vraća približnu vrednost drugog izvoda funkcije $f$ izračunatu korišćenjem I Njutnovog interpolacionog polinoma konstruisanog na osnovu svih vrednosti iz tablice formirane korišćenjem funkcije $drugiizvod()$.

• Napisati M-fajl $nula.m$ sa funkcijom $y=nula()$ koja metodom inverzne interpolacije približno određuje nulu drugog izvoda funkcije $f$ (pretpostavka je da je drugi izvod monotona funkcija) koristeći I Njutnov interpolacioni polinom zaključno sa konačnim razlikama reda $3$. Kriterijum zaustavljanja iterativnog niza: $|q_{i}-q_{i-1}|\leq 10^{-4}, i = 2...$

  test primer

  xxxxxxxxxx
      -12.2817  -14.6799  -16.9448  -19.0470  -20.9504  -22.6109
    >> [X Y Y2] = drugiizvod()
  ​
    X =
  ​
        1.2000    1.4000    1.6000    1.8000   
      
    Y2 =
      
        3.3312    4.0687    4.9696    6.0698
      
    >> y = vred2izvod(1.66)
    >>   Tablica konacnih razlika:
        1.2000    3.3312    0.7375    0.1633    0.0362
        1.4000    4.0687    0.9008    0.1994         0
        1.6000    4.9696    1.1003         0         0
        1.8000    6.0698         0         0         0
  ​
    y =
      
        5.2771
      
    >> y = nula()
    >>   Tablica konacnih razlika:
        1.2000    3.3312    0.7375    0.1633    0.0362
        1.4000    4.0687    0.9008    0.1994         0
        1.6000    4.9696    1.1003         0         0
        1.8000    6.0698         0         0         0
  ​
    y =
      
       -0.1068
  

9. zadatak - domaći

Neka je funkcija $f$ zadata tablično M-fajlom $tablica.m$ koji generiše dva niza $X=[x_{1},x_{2}...,x_{n}]$ i $Y=[y_{1},y_{2}...,y_{n}]$ (od kojih je prvi strogo rastući) za tu tablično zadatu funkciju. Tablica mora biti ekvidistantna ( sa korakom h ).

$tablica.m$

xxxxxxxxxx
X = linspace(2, 4, 10);
Y = exp(X);

• Napisati M-fajl $izvod1.m$ sa funkcijom $izvod1(x)$ koja računa vrednost prvog izvoda tabelirane funkcije u tački $x$ koristeći diferenciranje I Njutnovog interpolacionog polinoma zaključno sa konačnim razlikama reda $4$.

  xxxxxxxxxx
      % TABLICA KONACNIH RAZLIKA
      function izvod1(x)
      % Poziv funkcije: izvod1(2.1)
    
      tablica;
    
      n = length(X);
      % TABLICA KONAČNIH RAZLIKA
    
      k_razlike = zeros(n, 4);
      for i = 1 : n-1
          k_razlike(i, 1) = Y(i+1) - Y(i);
      end
      for j = 2 : 4
          for i = 1 : n-j
              k_razlike(i, j) = k_razlike(i+1, j-1) - k_razlike(i, j-1);
          end
      end
    
      disp([X', Y', k_razlike]);
    
      % Izvod I Njutnovog interpolacionog polinoma
      % y' = (k_razlike(1,1) + (q-0.5)*k_razlike(1,2) +
      % (3*q^2-6*q+2)*k_razlike(1,3)/6+(4*q^3-18*q^2+22*q-6)/24*krazlike(1,4))/h
    
      h = X(2) - X(1);
      Q = (x - X(1)) / h;
      yi = (k_razlike(1, 1) + (Q - 0.5) * k_razlike(1, 2) + (3 * Q^2 - 6 * Q + 2) * k_razlike(1, 3) / 6 + (4 * Q^3 - 18 * Q^2 + 22 * Q - 6) / 24 * k_razlike(1, 4)) / h;
    
      disp('Vrednost prvog izvoda funkcije f u tacki ');
      x
      disp('je: ');
      yi
  

• Napisati M-fajl $izvod2.m$ sa funkcijom $izvod2(x)$ koja računa vrednost drugog izvoda tabelirane funkcije u tački $x$ koristeći diferenciranje I Njutnovog interpolacionog polinoma zaključno sa konačnim razlikama reda $4$.

  test primer

  xxxxxxxxxx
    >>izvod2(2.3)
  ​
    Tablica konacnih razlika:
        2.0000    7.3891    1.8388    0.4576    0.1139    0.0283
        2.2222    9.2278    2.2963    0.5714    0.1422    0.0354
        2.4444   11.5241    2.8678    0.7136    0.1776    0.0442
        2.6667   14.3919    3.5814    0.8912    0.2218    0.0552
        2.8889   17.9733    4.4726    1.1130    0.2770    0.0689
        3.1111   22.4460    5.5857    1.3900    0.3459    0.0861
        3.3333   28.0316    6.9756    1.7359    0.4320         0
        3.5556   35.0073    8.7115    2.1679         0         0
        3.7778   43.7188   10.8794         0         0         0
        4.0000   54.5982         0         0         0         0
  ​
    Vrednost drugog izvoda funkcije f u tacki
  ​
    x =
  ​
        2.3000
  ​
    je:
  ​
    yi2 =
  ​
       9.9598