Polinomska interpolacija

Aproksimacija

Često se u praksi javlja potreba da se umesto neke funkcije koristi neka njoj bliska funkcija. Na primer, može se desiti da je funkcija previše komplikovana za neku upotrebu ili nema neka zgodna matematička svojstva, poput diferencijabilnosti ili konveksnosti. Ono što je češći problem je da su vrednosti funkcije poznate samo u nekim tačkama odnosno sama funkcija koja predstavlja vezu između nekih podataka uopšte nije poznata.

Razmotrimo primer broja stanovnika Beograda kroz godine:

Želimo da procenimo broj stanovnika 1994. godine. Posmatraćemo funkciju broja stanovnika u zavisnosti od vremena. Pretpostavićemo da je funkcija neprekidna pa u slučaju da vrednost funkcije u nekom trenutku nije ceo broj, samo ćemo uzeti ceo deo te vrednosti. Ovo je upravo situacija koju smo malopre spomenuli, gde funkcija ima poznate vrednosti samo u konačno mnogo tačaka (zadata je tabelarno).

Kriterijumi kvaliteta aproksimacije mogu biti različiti, i zavise od konkretnog slučaja. Razmatraćemo probleme nalaženje aproksimacije funkcije na nekom intervalu definisanosti te funkcije tako da se optimalnost aproksimacije posmatra na celom intervalu definisanosti, ili na nekom njegovom podskupu. Prirodna je ideja da se za aproksimaciju funkcije koristi neki polinom (npr. na ranijim kursevima smo se već sretali sa Tejlorovim polinomom).

Traženu funkciju broja stanovnika ćemo aproksimirati polinomom. Smatraćemo da je naša aproksimacija kvalitetna ako se vrednosti tog polinoma poklapaju sa vrednostima funkcije u datim tačkama (godinama). Naredna teorema je ključna za celu temu kojom ćemo se u nastavku baviti.

Teorema: Ako su nam poznate vrednosti (), , tada postoji tačno jedan polinom stepena ne većeg od , , za koji važi .

Vrednosti nazivamo čvorovima interpolacije, a polinom interpolacionim polinomom funkcije . Za grešku interpolacije, , važi narena ocena:

gde je .

Tačka u kojoj određujemo vrednost interpolacionog polinoma treba da bude između prvog i poslednjeg čvora interpolacije - inače govorimo o problemu ekstrapolacije.

Iako je sam polinom sa ovim svojstvom jedinstven, postoji više načina na koji se on može zapisati, te različiti načini donose interpolacione polinome različitih imena.

Lagranžov interpolacioni polinom

Interpolacioni polinom dat u obliku

nosi ime Lagranžov interpolacioni polinom. Ako bismo uvrstili vrednosti iz tabele stanovnika Beograda, interpolacijom bismo došli do zaključka da je 1994. godine u Beogradu bilo približno 1126212 stanovnika.

U nastavku je dat kod kojim smo dobili prethodni rezultat.

% function P = LagrangePoly(X,Y,x)

% n = length(X);

% P = zeros(1,n);

% for i = 1:n

% L = 1;

% prod = 1;

% for j = 1:n

% if j~=i

% prod = prod*(X(i)-X(j));

% L = conv(L,[1 -X(j)]);

% end

% P = P + L*Y(i)/prod;

% end

% val = polyval(P,x);

% end

X = [1953,1961,1971,1981,1991,2002,2011];

Y = [477982,657362,899094,1087915,1133146,1119642,1233796];

god = 1994;

brStanovnika = LagrangeVal(X,Y,god);

fprintf('Broj stanovika %d. godine je priblizno %d.', god, floor(brStanovnika));

Broj stanovika 1994. godine je priblizno 1126212.

Malom modifikacijom prethodnog programa, možemo odrediti i sam interpolacioni polinom. Polinomi se u MATLAB-u zadaju kao vektori koeficijenata, te se tako polinon zadaje kao Množenje dva polinoma se obavlja funkcijom Rezultati dobijeni na ovaj način su nešto manje precizni nego oni dobijeni prethodnom implementacijom zbog nagomilavanja računske greške, ali ovako možemo odrediti grafički prikaz interpolacionog polinoma. Kodovi su u nastavku:

% function P = LagrangePoly(X,Y,x)

% n = length(X);

% P = zeros(1,n);

% for i = 1:n

% l = 1;

% prod = 1;

% for j = 1:n

% if j~=i

% prod = prod*(X(i)-X(j));

% l = conv(l,[1 -X(j)]);

% end

% P = P + l*Y(i)/prod;

% end

% val = polyval(P,x);

% end

%

% function LagrangePrikaz(X,Y,P,x)

% plot(X,Y,'bo');

% hold on

% grid on

% Xrange = X(1):0.1:X(end);

% plot(Xrange,polyval(P,Xrange),'-r');

% plot(x,polyval(P,x),'r*');

% hold off

% end

LagPoly = LagrangePoly(X,Y,god);

LagrangePrikaz(X,Y,LagPoly,god);

Njutnov interpolacioni polinom

Pre nego što navedemo Njutnov oblik interpolacionog polinoma, moramo uvesti pojam podeljenih razlika. Definišemo redom:

- podeljene razlike prvog reda za čvor ;
- podeljene razliike drugog reda za čvor ;
- -tog reda za čvor .

Njutnov interpolacioni polinom sada definišemo na sledeći način:

Podeljene razlike često zapisujemo tabelarno (obojene su one vrednosti koje se pojavljuju z zapisu interpolacionog polinoma):

U nastavku je data implementacija metode kojom određujemo tabelu podeljenih razlika i Njutnov interpolacioni polinom, kao i vrednost polinoma u datoj tački.

% function [P,val] = Newton(X,Y,x)

% n = length(X);

% pRaz = zeros(n,n+1);

% pRaz(:,1) = X;

% pRaz(:,2) = Y;

% for i = 3:n+1

% for j = 1:n-i+2

% pRaz(j,i) = (pRaz(j+1,i-1)-pRaz(j,i-1))/(X(j+i-2)-X(j));

% end

% disp(pRaz);

% P = pRaz(1,2);

% L = 1;

% for i = 1:n-1

% L = conv(L,[1 -X(i)]);

% P = [0 P] + L*pRaz(1,i+2);

% end

% val = polyval(P,x);

% fprintf('Trazena vrednost je %f',val);

% end

X = [1 2 4 5 8];

Y = [1.5811 1.8708 2.3452 2.5495 3.0822];

[P,val] = Newton(X,Y,3.6);

  Columns 1 through 5

   1.0000e+00   1.5811e+00   2.8970e-01  -1.7500e-02   1.6333e-03
   2.0000e+00   1.8708e+00   2.3720e-01  -1.0967e-02   7.1389e-04
   4.0000e+00   2.3452e+00   2.0430e-01  -6.6833e-03            0
   5.0000e+00   2.5495e+00   1.7757e-01            0            0
   8.0000e+00   3.0822e+00            0            0            0

  Column 6

  -1.3135e-04
            0
            0
            0
            0

Trazena vrednost je 2.258496

Njutnov interpolacioni polinom za ekvidistantne čvorove

U slučaju kada su čvorovi interpolacije ekvidistantni, odnosno kada su susedni čvorovi uvek na konstantnoj udaljenosti , (imamo da je ) tada se umesto podeljenih najčešće koriste konačne razlike koje ćemo sada redom definisati:

- konačne razlike prvog reda za čvor .
- konačne razlike drugog reda za čvor .
- konačne razlike -tog reda za čvor .

Nije teško utvrditi vezu između podeljenih i konačnih razlika: .

Prvi Njutnov interpolacioni polinom

Ovaj oblik se najčešće koristi kada je tačka u kojoj želimo da odredimo vrednost interpolacionog polinoma blizu prvog čvora interpolacije. Ako je , tada se interpolacioni polinom može zapisati kao

Za grešku inteprolacije važi ocena

međutim, nekad se koristi i aproksimativna ocena greške (nekad ne koristimo konačne razlike svakog reda):

Drugi Njutnov interpolacioni polinom

Ovaj oblik se najčešće koristi kada je tačka u kojoj želimo da odredimo vrednost interpolacionog polinoma blizu poslednjeg čvora interpolacije. Ako je , tada se interpolacioni polinom može zapisati kao

Ocena greške:

Aproksimativna ocena greške:

Kao i kod podeljenih, obično se formira tabela konačnih razlika:

Prvi Njutnov polinom

1. Njutnov.png

Drugi Njutnov polinom

2. Njutnov.png

U nastavku je data implementacija metode za određivanje tabele konačnih razlika i vrednosti prvog Njutnovog interpolacionog polinoma u datoj tački (za zapis rezultata je korišćen format short e).

% function p = NewtonFirst(X,Y,x)

% n = length(X);

% kRaz = zeros(n,n+1);

% kRaz(:,1) = X;

% kRaz(:,2) = Y;

% for i = 3:n+1

% for j = 1:n-i+2

% kRaz(j,i) = kRaz(j+1,i-1)-kRaz(j,i-1);

% end

% disp(kRaz);

% h = X(2)-X(1);

% s = (x-X(1))/h;

% p = 0;

% l = 1;

% for i = 0:n-1

% p = p + l*kRaz(1,i+2)/factorial(i);

% l = l*(s-i);

% end

% fprintf('Trazena vrednost je %f',p);

% end

X = [150, 160, 170, 180, 190, 200];

Y = [2.17609, 2.20412, 2.23045, 2.25527, 2.27875, 2.30103];

p = NewtonFirst(X,Y,154);

  Columns 1 through 5

   1.5000e+02   2.1761e+00   2.8030e-02  -1.7000e-03   1.9000e-04
   1.6000e+02   2.2041e+00   2.6330e-02  -1.5100e-03   1.7000e-04
   1.7000e+02   2.2304e+00   2.4820e-02  -1.3400e-03   1.4000e-04
   1.8000e+02   2.2553e+00   2.3480e-02  -1.2000e-03            0
   1.9000e+02   2.2788e+00   2.2280e-02            0            0
   2.0000e+02   2.3010e+00            0            0            0

  Columns 6 through 7

  -2.0000e-05  -1.0000e-05
  -3.0000e-05            0
            0            0
            0            0
            0            0
            0            0

Trazena vrednost je 2.187519